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:   end do
 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:     end do
 64:   end do
 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) >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Error getting DMDA elements')
 72:   end do
 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:       end do
 88:     end do
 89:   end do
 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:   end do
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:     end do
138:   end do
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:       end do
159:     end do
160:   end do
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*/