1 static char help[] = "Test VTK structured (.vts) and rectilinear (.vtr) viewer support with multi-dof DMDAs.\n\ 2 Supply the -namefields flag to test with field names.\n\n"; 3 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 7 /* Helper function to name DMDA fields */ 8 PetscErrorCode NameFields(DM da,PetscInt dof) 9 { 10 PetscInt c; 11 12 PetscFunctionBeginUser; 13 for (c=0; c<dof; ++c) { 14 char fieldname[256]; 15 CHKERRQ(PetscSNPrintf(fieldname,sizeof(fieldname),"field_%D",c)); 16 CHKERRQ(DMDASetFieldName(da,c,fieldname)); 17 } 18 PetscFunctionReturn(0); 19 } 20 21 /* 22 Write 3D DMDA vector with coordinates in VTK format 23 */ 24 PetscErrorCode test_3d(const char filename[],PetscInt dof,PetscBool namefields) 25 { 26 MPI_Comm comm = MPI_COMM_WORLD; 27 const PetscInt M=10,N=15,P=30,sw=1; 28 const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0; 29 DM da; 30 Vec v; 31 PetscViewer view; 32 DMDALocalInfo info; 33 PetscScalar ****va; 34 PetscInt i,j,k,c; 35 36 CHKERRQ(DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,NULL,&da)); 37 CHKERRQ(DMSetFromOptions(da)); 38 CHKERRQ(DMSetUp(da)); 39 if (namefields) CHKERRQ(NameFields(da,dof)); 40 41 CHKERRQ(DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz)); 42 CHKERRQ(DMDAGetLocalInfo(da,&info)); 43 CHKERRQ(DMCreateGlobalVector(da,&v)); 44 CHKERRQ(DMDAVecGetArrayDOF(da,v,&va)); 45 for (k=info.zs; k<info.zs+info.zm; k++) { 46 for (j=info.ys; j<info.ys+info.ym; j++) { 47 for (i=info.xs; i<info.xs+info.xm; i++) { 48 const PetscScalar x = (Lx*i)/M; 49 const PetscScalar y = (Ly*j)/N; 50 const PetscScalar z = (Lz*k)/P; 51 for (c=0; c<dof; ++ c) { 52 va[k][j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2) + 10.0*c; 53 } 54 } 55 } 56 } 57 CHKERRQ(DMDAVecRestoreArrayDOF(da,v,&va)); 58 CHKERRQ(PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view)); 59 CHKERRQ(VecView(v,view)); 60 CHKERRQ(PetscViewerDestroy(&view)); 61 CHKERRQ(VecDestroy(&v)); 62 CHKERRQ(DMDestroy(&da)); 63 return 0; 64 } 65 66 /* 67 Write 2D DMDA vector with coordinates in VTK format 68 */ 69 PetscErrorCode test_2d(const char filename[],PetscInt dof,PetscBool namefields) 70 { 71 MPI_Comm comm = MPI_COMM_WORLD; 72 const PetscInt M=10,N=20,sw=1; 73 const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0; 74 DM da; 75 Vec v; 76 PetscViewer view; 77 DMDALocalInfo info; 78 PetscScalar ***va; 79 PetscInt i,j,c; 80 81 CHKERRQ(DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da)); 82 CHKERRQ(DMSetFromOptions(da)); 83 CHKERRQ(DMSetUp(da)); 84 if (namefields) CHKERRQ(NameFields(da,dof)); 85 CHKERRQ(DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz)); 86 CHKERRQ(DMDAGetLocalInfo(da,&info)); 87 CHKERRQ(DMCreateGlobalVector(da,&v)); 88 CHKERRQ(DMDAVecGetArrayDOF(da,v,&va)); 89 for (j=info.ys; j<info.ys+info.ym; j++) { 90 for (i=info.xs; i<info.xs+info.xm; i++) { 91 const PetscScalar x = (Lx*i)/M; 92 const PetscScalar y = (Ly*j)/N; 93 for (c=0; c<dof; ++c) { 94 va[j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2) + 10.0*c; 95 } 96 } 97 } 98 CHKERRQ(DMDAVecRestoreArrayDOF(da,v,&va)); 99 CHKERRQ(PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view)); 100 CHKERRQ(VecView(v,view)); 101 CHKERRQ(PetscViewerDestroy(&view)); 102 CHKERRQ(VecDestroy(&v)); 103 CHKERRQ(DMDestroy(&da)); 104 return 0; 105 } 106 107 /* 108 Write a scalar and a vector field from two compatible 3d DMDAs 109 */ 110 PetscErrorCode test_3d_compat(const char filename[],PetscInt dof,PetscBool namefields) 111 { 112 MPI_Comm comm = MPI_COMM_WORLD; 113 const PetscInt M=10,N=15,P=30,sw=1; 114 const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0; 115 DM da,daVector; 116 Vec v,vVector; 117 PetscViewer view; 118 DMDALocalInfo info; 119 PetscScalar ***va,****vVectora; 120 PetscInt i,j,k,c; 121 122 CHKERRQ(DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,P,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,/* dof:*/1,sw,NULL,NULL,NULL,&da)); 123 CHKERRQ(DMSetFromOptions(da)); 124 CHKERRQ(DMSetUp(da)); 125 if (namefields) CHKERRQ(NameFields(da,1)); 126 127 CHKERRQ(DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz)); 128 CHKERRQ(DMDAGetLocalInfo(da,&info)); 129 CHKERRQ(DMDACreateCompatibleDMDA(da,dof,&daVector)); 130 if (namefields) CHKERRQ(NameFields(daVector,dof)); 131 CHKERRQ(DMCreateGlobalVector(da,&v)); 132 CHKERRQ(DMCreateGlobalVector(daVector,&vVector)); 133 CHKERRQ(DMDAVecGetArray(da,v,&va)); 134 CHKERRQ(DMDAVecGetArrayDOF(daVector,vVector,&vVectora)); 135 for (k=info.zs; k<info.zs+info.zm; k++) { 136 for (j=info.ys; j<info.ys+info.ym; j++) { 137 for (i=info.xs; i<info.xs+info.xm; i++) { 138 const PetscScalar x = (Lx*i)/M; 139 const PetscScalar y = (Ly*j)/N; 140 const PetscScalar z = (Lz*k)/P; 141 va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2); 142 for (c=0; c<dof; ++c) { 143 vVectora[k][j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2) + 10.0*c; 144 } 145 } 146 } 147 } 148 CHKERRQ(DMDAVecRestoreArray(da,v,&va)); 149 CHKERRQ(DMDAVecRestoreArrayDOF(da,v,&vVectora)); 150 CHKERRQ(PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view)); 151 CHKERRQ(VecView(v,view)); 152 CHKERRQ(VecView(vVector,view)); 153 CHKERRQ(PetscViewerDestroy(&view)); 154 CHKERRQ(VecDestroy(&v)); 155 CHKERRQ(VecDestroy(&vVector)); 156 CHKERRQ(DMDestroy(&da)); 157 CHKERRQ(DMDestroy(&daVector)); 158 return 0; 159 } 160 161 /* 162 Write a scalar and a vector field from two compatible 2d DMDAs 163 */ 164 PetscErrorCode test_2d_compat(const char filename[],PetscInt dof,PetscBool namefields) 165 { 166 MPI_Comm comm = MPI_COMM_WORLD; 167 const PetscInt M=10,N=20,sw=1; 168 const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0; 169 DM da, daVector; 170 Vec v,vVector; 171 PetscViewer view; 172 DMDALocalInfo info; 173 PetscScalar **va,***vVectora; 174 PetscInt i,j,c; 175 176 CHKERRQ(DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,/* dof:*/ 1,sw,NULL,NULL,&da)); 177 CHKERRQ(DMSetFromOptions(da)); 178 CHKERRQ(DMSetUp(da)); 179 if (namefields) CHKERRQ(NameFields(da,1)); 180 CHKERRQ(DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz)); 181 CHKERRQ(DMDACreateCompatibleDMDA(da,dof,&daVector)); 182 if (namefields) CHKERRQ(NameFields(daVector,dof)); 183 CHKERRQ(DMDAGetLocalInfo(da,&info)); 184 CHKERRQ(DMCreateGlobalVector(da,&v)); 185 CHKERRQ(DMCreateGlobalVector(daVector,&vVector)); 186 CHKERRQ(DMDAVecGetArray(da,v,&va)); 187 CHKERRQ(DMDAVecGetArrayDOF(daVector,vVector,&vVectora)); 188 for (j=info.ys; j<info.ys+info.ym; j++) { 189 for (i=info.xs; i<info.xs+info.xm; i++) { 190 const PetscScalar x = (Lx*i)/M; 191 const PetscScalar y = (Ly*j)/N; 192 va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2); 193 for (c=0; c<dof; ++c) { 194 vVectora[j][i][c] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2) + 10.0*c; 195 } 196 } 197 } 198 CHKERRQ(DMDAVecRestoreArray(da,v,&va)); 199 CHKERRQ(DMDAVecRestoreArrayDOF(daVector,vVector,&vVectora)); 200 CHKERRQ(PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view)); 201 CHKERRQ(VecView(v,view)); 202 CHKERRQ(VecView(vVector,view)); 203 CHKERRQ(PetscViewerDestroy(&view)); 204 CHKERRQ(VecDestroy(&v)); 205 CHKERRQ(VecDestroy(&vVector)); 206 CHKERRQ(DMDestroy(&da)); 207 CHKERRQ(DMDestroy(&daVector)); 208 return 0; 209 } 210 211 int main(int argc, char *argv[]) 212 { 213 PetscErrorCode ierr; 214 PetscInt dof; 215 PetscBool namefields; 216 217 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 218 dof = 2; 219 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-dof",&dof,NULL)); 220 namefields = PETSC_FALSE; 221 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-namefields",&namefields,NULL)); 222 CHKERRQ(test_3d("3d.vtr",dof,namefields)); 223 CHKERRQ(test_2d("2d.vtr",dof,namefields)); 224 CHKERRQ(test_3d_compat("3d_compat.vtr",dof,namefields)); 225 CHKERRQ(test_2d_compat("2d_compat.vtr",dof,namefields)); 226 CHKERRQ(test_3d("3d.vts",dof,namefields)); 227 CHKERRQ(test_2d("2d.vts",dof,namefields)); 228 CHKERRQ(test_3d_compat("3d_compat.vts",dof,namefields)); 229 CHKERRQ(test_2d_compat("2d_compat.vts",dof,namefields)); 230 ierr = PetscFinalize(); 231 return ierr; 232 } 233 234 /*TEST 235 236 build: 237 requires: !complex 238 239 test: 240 suffix: 1 241 nsize: 2 242 args: -dof 2 243 244 test: 245 suffix: 2 246 nsize: 2 247 args: -dof 2 248 249 test: 250 suffix: 3 251 nsize: 2 252 args: -dof 3 253 254 test: 255 suffix: 4 256 nsize: 1 257 args: -dof 2 -namefields 258 259 test: 260 suffix: 5 261 nsize: 2 262 args: -dof 4 -namefields 263 264 TEST*/ 265