1! Tests DMDAGetVecGetArray() 2 3program main 4#include <petsc/finclude/petscdm.h> 5#include <petsc/finclude/petscdmda.h> 6 use petscdmda 7 use petsc 8 implicit none 9 10 Type(tVec) g 11 Type(tDM) ada 12 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 19 20 PetscInt nen, nel 21 PetscInt, pointer :: elements(:) 22 23 PetscInt nfields 24 character(80), pointer :: namefields(:) 25 IS, pointer :: isfields(:) 26 DM, pointer :: dmfields(:) 27 PetscInt zero, one 28 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)) 52 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)) 68 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)) 75 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)) 94 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)) 113 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)) 122 123 PetscCallA(DMDestroy(ada, ierr)) 124 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)) 143 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)) 165 166 PetscCallA(PetscFinalize(ierr)) 167END PROGRAM 168 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*/ 179