1! Tests DMDAGetVecGetArray() 2 3 program 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 enddo 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 enddo 64 enddo 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) .ge. 0,PETSC_COMM_SELF,PETSC_ERR_PLIB,'Error getting DMDA elements') 72 enddo 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 enddo 88 enddo 89 enddo 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 enddo 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 enddo 138 enddo 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 enddo 159 enddo 160 enddo 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)) 167 END 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