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 m, n, p, dof, s, i, j, k, xs, xl
14: PetscInt ys, yl
15: PetscInt zs, zl, sw
17: PetscInt nen, nel
18: PetscInt, pointer :: elements(:)
20: PetscInt nfields
21: character(80), pointer :: namefields(:)
22: IS, pointer :: isfields(:)
23: DM, pointer :: dmfields(:)
24: PetscInt zero, one
26: m = 5
27: n = 6
28: p = 4
29: s = 1
30: dof = 1
31: sw = 1
32: zero = 0
33: one = 1
34: PetscCallA(PetscInitialize(ierr))
35: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, m, dof, sw, PETSC_NULL_INTEGER_ARRAY, ada, ierr))
36: PetscCallA(DMSetUp(ada, ierr))
37: PetscCallA(DMGetGlobalVector(ada, g, ierr))
38: PetscCallA(DMDAGetCorners(ada, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xl, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
39: PetscCallA(DMDAVecGetArray(ada, g, x1, ierr))
40: do i = xs, xs + xl - 1
41: ! CHKMEMQ
42: x1(i) = i
43: ! CHKMEMQ
44: end do
45: PetscCallA(DMDAVecRestoreArray(ada, g, x1, ierr))
46: PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
47: PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
48: PetscCallA(DMDestroy(ada, ierr))
50: 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))
51: PetscCallA(DMSetUp(ada, ierr))
52: PetscCallA(DMGetGlobalVector(ada, g, ierr))
53: PetscCallA(DMDAGetCorners(ada, xs, ys, PETSC_NULL_INTEGER, xl, yl, PETSC_NULL_INTEGER, ierr))
54: PetscCallA(DMDAVecGetArray(ada, g, x2, ierr))
55: do i = xs, xs + xl - 1
56: do j = ys, ys + yl - 1
57: ! CHKMEMQ
58: x2(i, j) = i + j
59: ! CHKMEMQ
60: end do
61: end do
62: PetscCallA(DMDAVecRestoreArray(ada, g, x2, ierr))
63: PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
64: PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
66: PetscCallA(DMDAGetElements(ada, nen, nel, elements, ierr))
67: do i = 1, nen*nel
68: PetscCheckA(elements(i) >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Error getting DMDA elements')
69: end do
70: PetscCallA(DMDARestoreElements(ada, nen, nel, elements, ierr))
71: PetscCallA(DMDestroy(ada, ierr))
73: 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))
74: PetscCallA(DMSetUp(ada, ierr))
75: PetscCallA(DMGetGlobalVector(ada, g, ierr))
76: PetscCallA(DMDAGetCorners(ada, xs, ys, zs, xl, yl, zl, ierr))
77: PetscCallA(DMDAVecGetArray(ada, g, x3, ierr))
78: do i = xs, xs + xl - 1
79: do j = ys, ys + yl - 1
80: do k = zs, zs + zl - 1
81: ! CHKMEMQ
82: x3(i, j, k) = i + j + k
83: ! CHKMEMQ
84: end do
85: end do
86: end do
87: PetscCallA(DMDAVecRestoreArray(ada, g, x3, ierr))
88: PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
89: PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
90: PetscCallA(DMDestroy(ada, ierr))
92: !
93: ! Same tests but now with DOF > 1, so dimensions of array are one higher
94: !
95: dof = 2
96: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, m, dof, sw, PETSC_NULL_INTEGER_ARRAY, ada, ierr))
97: PetscCallA(DMSetUp(ada, ierr))
98: PetscCallA(DMGetGlobalVector(ada, g, ierr))
99: PetscCallA(DMDAGetCorners(ada, xs, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, xl, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
100: PetscCallA(DMDAVecGetArray(ada, g, x2, ierr))
101: do i = xs, xs + xl - 1
102: ! CHKMEMQ
103: x2(0, i) = i
104: x2(1, i) = -i
105: ! CHKMEMQ
106: end do
107: PetscCallA(DMDAVecRestoreArray(ada, g, x1, ierr))
108: PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
109: PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
111: ! some testing unrelated to the example
112: PetscCallA(DMDASetFieldName(ada, zero, 'Field 0', ierr))
113: PetscCallA(DMDASetFieldName(ada, one, 'Field 1', ierr))
114: PetscCallA(DMCreateFieldDecomposition(ada, nfields, namefields, PETSC_NULL_IS_POINTER, PETSC_NULL_DM_POINTER, ierr))
115: ! print*,nfields,trim(namefields(1)),trim(namefields(2))
116: PetscCallA(DMDestroyFieldDecomposition(ada, nfields, namefields, PETSC_NULL_IS_POINTER, PETSC_NULL_DM_POINTER, ierr))
117: PetscCallA(DMCreateFieldDecomposition(ada, nfields, namefields, isfields, dmfields, ierr))
118: PetscCallA(DMDestroyFieldDecomposition(ada, nfields, namefields, isfields, dmfields, ierr))
120: PetscCallA(DMDestroy(ada, ierr))
122: dof = 2
123: 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))
124: PetscCallA(DMSetUp(ada, ierr))
125: PetscCallA(DMGetGlobalVector(ada, g, ierr))
126: PetscCallA(DMDAGetCorners(ada, xs, ys, PETSC_NULL_INTEGER, xl, yl, PETSC_NULL_INTEGER, ierr))
127: PetscCallA(DMDAVecGetArray(ada, g, x3, ierr))
128: do i = xs, xs + xl - 1
129: do j = ys, ys + yl - 1
130: ! CHKMEMQ
131: x3(0, i, j) = i + j
132: x3(1, i, j) = -(i + j)
133: ! CHKMEMQ
134: end do
135: end do
136: PetscCallA(DMDAVecRestoreArray(ada, g, x3, ierr))
137: PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
138: PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
139: PetscCallA(DMDestroy(ada, ierr))
141: dof = 3
142: 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))
143: PetscCallA(DMSetUp(ada, ierr))
144: PetscCallA(DMGetGlobalVector(ada, g, ierr))
145: PetscCallA(DMDAGetCorners(ada, xs, ys, zs, xl, yl, zl, ierr))
146: PetscCallA(DMDAVecGetArray(ada, g, x4, ierr))
147: do i = xs, xs + xl - 1
148: do j = ys, ys + yl - 1
149: do k = zs, zs + zl - 1
150: ! CHKMEMQ
151: x4(0, i, j, k) = i + j + k
152: x4(1, i, j, k) = -(i + j + k)
153: x4(2, i, j, k) = i + j + k
154: ! CHKMEMQ
155: end do
156: end do
157: end do
158: PetscCallA(DMDAVecRestoreArray(ada, g, x4, ierr))
159: PetscCallA(VecView(g, PETSC_VIEWER_STDOUT_WORLD, ierr))
160: PetscCallA(DMRestoreGlobalVector(ada, g, ierr))
161: PetscCallA(DMDestroy(ada, ierr))
163: PetscCallA(PetscFinalize(ierr))
164: end program
166: !
167: !/*TEST
168: !
169: ! build:
170: ! requires: !complex
171: !
172: ! test:
173: ! filter: Error: grep -v "Vec Object" | grep -v "Warning: ieee_inexact is signaling"
174: !
175: !TEST*/