Actual source code: ex11f90.F90

  1: !     Tests DMDAGetVecGetArray()

  3:       program main
  4: #include <petsc/finclude/petscdm.h>
  5: #include <petsc/finclude/petscdmda.h>
  6:       use petscdmda
  7:       use petsc
  8:       implicit none

 10:       Type(tVec)  g
 11:       Type(tDM)   ada

 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

 20:       PetscInt nen,nel
 21:       PetscInt, pointer :: elements(:)

 23:       PetscInt nfields
 24:       character(80), pointer :: namefields(:)
 25:       IS, pointer :: isfields(:)
 26:       DM, pointer :: dmfields(:)
 27:       PetscInt zero, one

 29:       m = 5
 30:       n = 6
 31:       p = 4;
 32:       s = 1
 33:       dof = 1
 34:       sw = 1
 35:       zero = 0
 36:       one = 1
 37:       PetscCallA(PetscInitialize(ierr))
 38:       PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m,dof,sw,PETSC_NULL_INTEGER_ARRAY,ada,ierr))
 39:       PetscCallA(DMSetUp(ada,ierr))
 40:       PetscCallA(DMGetGlobalVector(ada,g,ierr))
 41:       PetscCallA(DMDAGetCorners(ada,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xl,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
 42:       PetscCallA(DMDAVecGetArray(ada,g,x1,ierr))
 43:       do i=xs,xs+xl-1
 44: !         CHKMEMQ
 45:          x1(i) = i
 46: !         CHKMEMQ
 47:       enddo
 48:       PetscCallA(DMDAVecRestoreArray(ada,g,x1,ierr))
 49:       PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr))
 50:       PetscCallA(DMRestoreGlobalVector(ada,g,ierr))
 51:       PetscCallA(DMDestroy(ada,ierr))

 53:       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))
 54:       PetscCallA(DMSetUp(ada,ierr))
 55:       PetscCallA(DMGetGlobalVector(ada,g,ierr))
 56:       PetscCallA(DMDAGetCorners(ada,xs,ys,PETSC_NULL_INTEGER,xl,yl,PETSC_NULL_INTEGER,ierr))
 57:       PetscCallA(DMDAVecGetArray(ada,g,x2,ierr))
 58:       do i=xs,xs+xl-1
 59:         do j=ys,ys+yl-1
 60: !           CHKMEMQ
 61:            x2(i,j) = i + j
 62: !           CHKMEMQ
 63:         enddo
 64:       enddo
 65:       PetscCallA(DMDAVecRestoreArray(ada,g,x2,ierr))
 66:       PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr))
 67:       PetscCallA(DMRestoreGlobalVector(ada,g,ierr))

 69:       PetscCallA(DMDAGetElements(ada,nen,nel,elements,ierr))
 70:       do i=1,nen*nel
 71:          PetscCheckA(elements(i) .ge. 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error getting DMDA elements')
 72:       enddo
 73:       PetscCallA(DMDARestoreElements(ada,nen,nel,elements,ierr))
 74:       PetscCallA(DMDestroy(ada,ierr))

 76:       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))
 77:       PetscCallA(DMSetUp(ada,ierr))
 78:       PetscCallA(DMGetGlobalVector(ada,g,ierr))
 79:       PetscCallA(DMDAGetCorners(ada,xs,ys,zs,xl,yl,zl,ierr))
 80:       PetscCallA(DMDAVecGetArray(ada,g,x3,ierr))
 81:       do i=xs,xs+xl-1
 82:         do j=ys,ys+yl-1
 83:           do k=zs,zs+zl-1
 84: !            CHKMEMQ
 85:             x3(i,j,k) = i + j + k
 86: !            CHKMEMQ
 87:           enddo
 88:         enddo
 89:       enddo
 90:       PetscCallA(DMDAVecRestoreArray(ada,g,x3,ierr))
 91:       PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr))
 92:       PetscCallA(DMRestoreGlobalVector(ada,g,ierr))
 93:       PetscCallA(DMDestroy(ada,ierr))

 95: !
 96: !  Same tests but now with DOF > 1, so dimensions of array are one higher
 97: !
 98:       dof = 2
 99:       PetscCallA(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m,dof,sw,PETSC_NULL_INTEGER_ARRAY,ada,ierr))
100:       PetscCallA(DMSetUp(ada,ierr))
101:       PetscCallA(DMGetGlobalVector(ada,g,ierr))
102:       PetscCallA(DMDAGetCorners(ada,xs,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,xl,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
103:       PetscCallA(DMDAVecGetArray(ada,g,x2,ierr))
104:       do i=xs,xs+xl-1
105: !         CHKMEMQ
106:          x2(0,i) = i
107:          x2(1,i) = -i
108: !         CHKMEMQ
109:       enddo
110:       PetscCallA(DMDAVecRestoreArray(ada,g,x1,ierr))
111:       PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr))
112:       PetscCallA(DMRestoreGlobalVector(ada,g,ierr))

114:       ! some testing unrelated to the example
115:       PetscCallA(DMDASetFieldName(ada,zero,'Field 0',ierr))
116:       PetscCallA(DMDASetFieldName(ada,one,'Field 1',ierr))
117:       PetscCallA(DMCreateFieldDecomposition(ada, nfields, namefields, PETSC_NULL_IS_POINTER, PETSC_NULL_DM_POINTER, ierr))
118:       ! print*,nfields,trim(namefields(1)),trim(namefields(2))
119:       PetscCallA(DMDestroyFieldDecomposition(ada, nfields, namefields, PETSC_NULL_IS_POINTER, PETSC_NULL_DM_POINTER, ierr))
120:       PetscCallA(DMCreateFieldDecomposition(ada, nfields, namefields, isfields, dmfields, ierr))
121:       PetscCallA(DMDestroyFieldDecomposition(ada, nfields, namefields, isfields, dmfields, ierr))

123:       PetscCallA(DMDestroy(ada,ierr))

125:       dof = 2
126:       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))
127:       PetscCallA(DMSetUp(ada,ierr))
128:       PetscCallA(DMGetGlobalVector(ada,g,ierr))
129:       PetscCallA(DMDAGetCorners(ada,xs,ys,PETSC_NULL_INTEGER,xl,yl,PETSC_NULL_INTEGER,ierr))
130:       PetscCallA(DMDAVecGetArray(ada,g,x3,ierr))
131:       do i=xs,xs+xl-1
132:         do j=ys,ys+yl-1
133: !           CHKMEMQ
134:            x3(0,i,j) = i + j
135:            x3(1,i,j) = -(i + j)
136: !           CHKMEMQ
137:         enddo
138:       enddo
139:       PetscCallA(DMDAVecRestoreArray(ada,g,x3,ierr))
140:       PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr))
141:       PetscCallA(DMRestoreGlobalVector(ada,g,ierr))
142:       PetscCallA(DMDestroy(ada,ierr))

144:       dof = 3
145:       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))
146:       PetscCallA(DMSetUp(ada,ierr))
147:       PetscCallA(DMGetGlobalVector(ada,g,ierr))
148:       PetscCallA(DMDAGetCorners(ada,xs,ys,zs,xl,yl,zl,ierr))
149:       PetscCallA(DMDAVecGetArray(ada,g,x4,ierr))
150:       do i=xs,xs+xl-1
151:         do j=ys,ys+yl-1
152:           do k=zs,zs+zl-1
153: !            CHKMEMQ
154:             x4(0,i,j,k) = i + j + k
155:             x4(1,i,j,k) = -(i + j + k)
156:             x4(2,i,j,k) = i + j + k
157: !            CHKMEMQ
158:           enddo
159:         enddo
160:       enddo
161:       PetscCallA(DMDAVecRestoreArray(ada,g,x4,ierr))
162:       PetscCallA(VecView(g,PETSC_VIEWER_STDOUT_WORLD,ierr))
163:       PetscCallA(DMRestoreGlobalVector(ada,g,ierr))
164:       PetscCallA(DMDestroy(ada,ierr))

166:       PetscCallA(PetscFinalize(ierr))
167:       END PROGRAM

169: !
170: !/*TEST
171: !
172: !   build:
173: !     requires: !complex
174: !
175: !   test:
176: !     filter: Error: grep -v "Vec Object" | grep -v "Warning: ieee_inexact is signaling"
177: !
178: !TEST*/