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