Actual source code: ex11f90.F90

  1: !     Tests DMDAGetVecGetArray()
  2: #include <petsc/finclude/petscdmda.h>
  3: program main
  4:   use petscdmda
  5:   implicit none

  7:   type(tVec) g
  8:   type(tDM) ada

 10:   PetscScalar, pointer :: x1(:), x2(:, :)
 11:   PetscScalar, pointer :: x3(:, :, :), x4(:, :, :, :)
 12:   PetscErrorCode ierr
 13:   PetscInt, parameter :: m = 5, n = 6, p = 4, s = 1, sw = 1
 14:   PetscInt i, j, k, dof
 15:   PetscInt xs, xl, ys, yl, zs, zl

 17:   PetscInt nen, nel
 18:   PetscInt, pointer :: elements(:)

 20:   PetscInt nfields
 21:   character(80), pointer :: namefields(:)
 22:   IS, pointer :: isfields(:)
 23:   DM, pointer :: dmfields(:)

 25:   PetscCallA(PetscInitialize(ierr))

 27:   dof = 1
 28:   PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, m, dof, sw, PETSC_NULL_INTEGER_ARRAY, ada, ierr))
 29:   PetscCallA(DMSetUp(ada, ierr))
 30:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
 31:   PetscCallA(DMDAGetCorners(ada, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xl, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
 32:   PetscCallA(DMDAVecGetArray(ada, g, x1, ierr))
 33:   do i = xs, xs + xl - 1
 34:     x1(i) = i
 35:   end do
 36:   PetscCallA(DMDAVecRestoreArray(ada, g, x1, ierr))
 37:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
 38:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
 39:   PetscCallA(DMDestroy(ada, ierr))

 41:   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))
 42:   PetscCallA(DMSetUp(ada, ierr))
 43:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
 44:   PetscCallA(DMDAGetCorners(ada, xs, ys, PETSC_NULL_INTEGER, xl, yl, PETSC_NULL_INTEGER, ierr))
 45:   PetscCallA(DMDAVecGetArray(ada, g, x2, ierr))
 46:   do i = xs, xs + xl - 1
 47:     do j = ys, ys + yl - 1
 48:       x2(i, j) = i + j
 49:     end do
 50:   end do
 51:   PetscCallA(DMDAVecRestoreArray(ada, g, x2, ierr))
 52:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
 53:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))

 55:   PetscCallA(DMDAGetElements(ada, nen, nel, elements, ierr))
 56:   do i = 1, nen*nel
 57:     PetscCheckA(elements(i) >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Error getting DMDA elements')
 58:   end do
 59:   PetscCallA(DMDARestoreElements(ada, nen, nel, elements, ierr))
 60:   PetscCallA(DMDestroy(ada, ierr))

 62:   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))
 63:   PetscCallA(DMSetUp(ada, ierr))
 64:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
 65:   PetscCallA(DMDAGetCorners(ada, xs, ys, zs, xl, yl, zl, ierr))
 66:   PetscCallA(DMDAVecGetArray(ada, g, x3, ierr))
 67:   do i = xs, xs + xl - 1
 68:     do j = ys, ys + yl - 1
 69:       do k = zs, zs + zl - 1
 70:         x3(i, j, k) = i + j + k
 71:       end do
 72:     end do
 73:   end do
 74:   PetscCallA(DMDAVecRestoreArray(ada, g, x3, ierr))
 75:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
 76:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
 77:   PetscCallA(DMDestroy(ada, ierr))

 79: !
 80: !  Same tests but now with DOF > 1, so dimensions of array are one higher
 81: !
 82:   dof = 2
 83:   PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, m, dof, sw, PETSC_NULL_INTEGER_ARRAY, ada, ierr))
 84:   PetscCallA(DMSetUp(ada, ierr))
 85:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
 86:   PetscCallA(DMDAGetCorners(ada, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xl, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
 87:   PetscCallA(DMDAVecGetArray(ada, g, x2, ierr))
 88:   do i = xs, xs + xl - 1
 89:     x2(0, i) = i
 90:     x2(1, i) = -i
 91:   end do
 92:   PetscCallA(DMDAVecRestoreArray(ada, g, x1, ierr))
 93:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
 94:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))

 96:   ! some testing unrelated to the example
 97:   PetscCallA(DMDASetFieldName(ada, 0_PETSC_INT_KIND, 'Field 0', ierr))
 98:   PetscCallA(DMDASetFieldName(ada, 1_PETSC_INT_KIND, 'Field 1', ierr))
 99:   PetscCallA(DMCreateFieldDecomposition(ada, nfields, namefields, PETSC_NULL_IS_POINTER, PETSC_NULL_DM_POINTER, ierr))
100:   ! print*,nfields,trim(namefields(1)),trim(namefields(2))
101:   PetscCallA(DMDestroyFieldDecomposition(ada, nfields, namefields, PETSC_NULL_IS_POINTER, PETSC_NULL_DM_POINTER, ierr))
102:   PetscCallA(DMCreateFieldDecomposition(ada, nfields, namefields, isfields, dmfields, ierr))
103:   PetscCallA(DMDestroyFieldDecomposition(ada, nfields, namefields, isfields, dmfields, ierr))

105:   PetscCallA(DMDestroy(ada, ierr))

107:   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))
108:   PetscCallA(DMSetUp(ada, ierr))
109:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
110:   PetscCallA(DMDAGetCorners(ada, xs, ys, PETSC_NULL_INTEGER, xl, yl, PETSC_NULL_INTEGER, ierr))
111:   PetscCallA(DMDAVecGetArray(ada, g, x3, ierr))
112:   do i = xs, xs + xl - 1
113:     do j = ys, ys + yl - 1
114:       x3(0, i, j) = i + j
115:       x3(1, i, j) = -(i + j)
116:     end do
117:   end do
118:   PetscCallA(DMDAVecRestoreArray(ada, g, x3, ierr))
119:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
120:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
121:   PetscCallA(DMDestroy(ada, ierr))

123:   dof = 3
124:   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))
125:   PetscCallA(DMSetUp(ada, ierr))
126:   PetscCallA(DMGetGlobalVector(ada, g, ierr))
127:   PetscCallA(DMDAGetCorners(ada, xs, ys, zs, xl, yl, zl, ierr))
128:   PetscCallA(DMDAVecGetArray(ada, g, x4, ierr))
129:   do i = xs, xs + xl - 1
130:     do j = ys, ys + yl - 1
131:       do k = zs, zs + zl - 1
132:         x4(0, i, j, k) = i + j + k
133:         x4(1, i, j, k) = -(i + j + k)
134:         x4(2, i, j, k) = i + j + k
135:       end do
136:     end do
137:   end do
138:   PetscCallA(DMDAVecRestoreArray(ada, g, x4, ierr))
139:   PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
140:   PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
141:   PetscCallA(DMDestroy(ada, ierr))

143:   PetscCallA(PetscFinalize(ierr))
144: end program

146: !
147: !/*TEST
148: !
149: !   build:
150: !     requires: !complex
151: !
152: !   test:
153: !     filter: Error: grep -v "Vec Object" | grep -v "Warning: ieee_inexact is signaling"
154: !
155: !TEST*/