1! Tests DMDAGetVecGetArray() 2 3 program main 4#include <petsc/finclude/petscdm.h> 5#include <petsc/finclude/petscdmda.h> 6 use petscdmda 7 use petsc 8 implicit none 9 10 Type(tVec) g 11 Type(tDM) ada 12 13 PetscScalar,pointer :: x1(:),x2(:,:) 14 PetscScalar,pointer :: x3(:,:,:),x4(:,:,:,:) 15 PetscErrorCode ierr 16 PetscInt m,n,p,dof,s,i,j,k,xs,xl 17 PetscInt ys,yl 18 PetscInt zs,zl,sw 19 20 PetscInt nen,nel 21 PetscInt, pointer :: elements(:) 22 23 m = 5 24 n = 6 25 p = 4; 26 s = 1 27 dof = 1 28 sw = 1 29 PetscCallA(PetscInitialize(ierr)) 30 PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m,dof,sw,PETSC_NULL_INTEGER_ARRAY,ada,ierr)) 31 PetscCallA(DMSetUp(ada,ierr)) 32 PetscCallA(DMGetGlobalVector(ada,g,ierr)) 33 PetscCallA(DMDAGetCorners(ada,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xl,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)) 34 PetscCallA(DMDAVecGetArray(ada,g,x1,ierr)) 35 do i=xs,xs+xl-1 36! CHKMEMQ 37 x1(i) = i 38! CHKMEMQ 39 enddo 40 PetscCallA(DMDAVecRestoreArray(ada,g,x1,ierr)) 41 PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr)) 42 PetscCallA(DMRestoreGlobalVector(ada,g,ierr)) 43 PetscCallA(DMDestroy(ada,ierr)) 44 45 PetscCallA(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,PETSC_DECIDE,PETSC_DECIDE,dof,s,PETSC_NULL_INTEGER_ARRAY,PETSC_NULL_INTEGER_ARRAY,ada,ierr)) 46 PetscCallA(DMSetUp(ada,ierr)) 47 PetscCallA(DMGetGlobalVector(ada,g,ierr)) 48 PetscCallA(DMDAGetCorners(ada,xs,ys,PETSC_NULL_INTEGER,xl,yl,PETSC_NULL_INTEGER,ierr)) 49 PetscCallA(DMDAVecGetArray(ada,g,x2,ierr)) 50 do i=xs,xs+xl-1 51 do j=ys,ys+yl-1 52! CHKMEMQ 53 x2(i,j) = i + j 54! CHKMEMQ 55 enddo 56 enddo 57 PetscCallA(DMDAVecRestoreArray(ada,g,x2,ierr)) 58 PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr)) 59 PetscCallA(DMRestoreGlobalVector(ada,g,ierr)) 60 61 PetscCallA(DMDAGetElements(ada,nen,nel,elements,ierr)) 62 do i=1,nen*nel 63 PetscCheckA(elements(i) .ge. 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error getting DMDA elements') 64 enddo 65 PetscCallA(DMDARestoreElements(ada,nen,nel,elements,ierr)) 66 PetscCallA(DMDestroy(ada,ierr)) 67 68 PetscCallA(DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX, m,n,p,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,s,PETSC_NULL_INTEGER_ARRAY,PETSC_NULL_INTEGER_ARRAY,PETSC_NULL_INTEGER_ARRAY,ada,ierr)) 69 PetscCallA(DMSetUp(ada,ierr)) 70 PetscCallA(DMGetGlobalVector(ada,g,ierr)) 71 PetscCallA(DMDAGetCorners(ada,xs,ys,zs,xl,yl,zl,ierr)) 72 PetscCallA(DMDAVecGetArray(ada,g,x3,ierr)) 73 do i=xs,xs+xl-1 74 do j=ys,ys+yl-1 75 do k=zs,zs+zl-1 76! CHKMEMQ 77 x3(i,j,k) = i + j + k 78! CHKMEMQ 79 enddo 80 enddo 81 enddo 82 PetscCallA(DMDAVecRestoreArray(ada,g,x3,ierr)) 83 PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr)) 84 PetscCallA(DMRestoreGlobalVector(ada,g,ierr)) 85 PetscCallA(DMDestroy(ada,ierr)) 86 87! 88! Same tests but now with DOF > 1, so dimensions of array are one higher 89! 90 dof = 2 91 PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m,dof,sw,PETSC_NULL_INTEGER_ARRAY,ada,ierr)) 92 PetscCallA(DMSetUp(ada,ierr)) 93 PetscCallA(DMGetGlobalVector(ada,g,ierr)) 94 PetscCallA(DMDAGetCorners(ada,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xl,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)) 95 PetscCallA(DMDAVecGetArray(ada,g,x2,ierr)) 96 do i=xs,xs+xl-1 97! CHKMEMQ 98 x2(0,i) = i 99 x2(1,i) = -i 100! CHKMEMQ 101 enddo 102 PetscCallA(DMDAVecRestoreArray(ada,g,x1,ierr)) 103 PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr)) 104 PetscCallA(DMRestoreGlobalVector(ada,g,ierr)) 105 PetscCallA(DMDestroy(ada,ierr)) 106 107 dof = 2 108 PetscCallA(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,PETSC_DECIDE,PETSC_DECIDE,dof,s,PETSC_NULL_INTEGER_ARRAY,PETSC_NULL_INTEGER_ARRAY,ada,ierr)) 109 PetscCallA(DMSetUp(ada,ierr)) 110 PetscCallA(DMGetGlobalVector(ada,g,ierr)) 111 PetscCallA(DMDAGetCorners(ada,xs,ys,PETSC_NULL_INTEGER,xl,yl,PETSC_NULL_INTEGER,ierr)) 112 PetscCallA(DMDAVecGetArray(ada,g,x3,ierr)) 113 do i=xs,xs+xl-1 114 do j=ys,ys+yl-1 115! CHKMEMQ 116 x3(0,i,j) = i + j 117 x3(1,i,j) = -(i + j) 118! CHKMEMQ 119 enddo 120 enddo 121 PetscCallA(DMDAVecRestoreArray(ada,g,x3,ierr)) 122 PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr)) 123 PetscCallA(DMRestoreGlobalVector(ada,g,ierr)) 124 PetscCallA(DMDestroy(ada,ierr)) 125 126 dof = 3 127 PetscCallA(DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,m,n,p,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,s,PETSC_NULL_INTEGER_ARRAY,PETSC_NULL_INTEGER_ARRAY,PETSC_NULL_INTEGER_ARRAY,ada,ierr)) 128 PetscCallA(DMSetUp(ada,ierr)) 129 PetscCallA(DMGetGlobalVector(ada,g,ierr)) 130 PetscCallA(DMDAGetCorners(ada,xs,ys,zs,xl,yl,zl,ierr)) 131 PetscCallA(DMDAVecGetArray(ada,g,x4,ierr)) 132 do i=xs,xs+xl-1 133 do j=ys,ys+yl-1 134 do k=zs,zs+zl-1 135! CHKMEMQ 136 x4(0,i,j,k) = i + j + k 137 x4(1,i,j,k) = -(i + j + k) 138 x4(2,i,j,k) = i + j + k 139! CHKMEMQ 140 enddo 141 enddo 142 enddo 143 PetscCallA(DMDAVecRestoreArray(ada,g,x4,ierr)) 144 PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr)) 145 PetscCallA(DMRestoreGlobalVector(ada,g,ierr)) 146 PetscCallA(DMDestroy(ada,ierr)) 147 148 PetscCallA(PetscFinalize(ierr)) 149 END PROGRAM 150 151! 152!/*TEST 153! 154! build: 155! requires: !complex 156! 157! test: 158! filter: Error: grep -v "Vec Object" | grep -v "Warning: ieee_inexact is signaling" 159! 160!TEST*/ 161