1! Tests DMDAGetVecGetArray() 2#include <petsc/finclude/petscdmda.h> 3program main 4 use petscdmda 5 implicit none 6 7 type(tVec) g 8 type(tDM) ada 9 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 16 17 PetscInt nen, nel 18 PetscInt, pointer :: elements(:) 19 20 PetscInt nfields 21 character(80), pointer :: namefields(:) 22 IS, pointer :: isfields(:) 23 DM, pointer :: dmfields(:) 24 PetscInt zero, one 25 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)) 49 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)) 65 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)) 72 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)) 91 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)) 110 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)) 119 120 PetscCallA(DMDestroy(ada, ierr)) 121 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)) 140 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)) 162 163 PetscCallA(PetscFinalize(ierr)) 164end program 165 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*/ 176