1 static char help[] = "Test VTK structured grid (.vts) viewer support\n\n"; 2 3 #include <petscdm.h> 4 #include <petscdmda.h> 5 6 /* 7 Write 3D DMDA vector with coordinates in VTK .vts format 8 9 */ 10 PetscErrorCode test_3d(const char filename[]) 11 { 12 MPI_Comm comm = MPI_COMM_WORLD; 13 const PetscInt M=10,N=15,P=30,dof=1,sw=1; 14 const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0; 15 DM da; 16 Vec v; 17 PetscViewer view; 18 DMDALocalInfo info; 19 PetscScalar ***va; 20 PetscInt i,j,k; 21 22 PetscCall(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)); 23 PetscCall(DMSetFromOptions(da)); 24 PetscCall(DMSetUp(da)); 25 26 PetscCall(DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz)); 27 PetscCall(DMDAGetLocalInfo(da,&info)); 28 PetscCall(DMCreateGlobalVector(da,&v)); 29 PetscCall(DMDAVecGetArray(da,v,&va)); 30 for (k=info.zs; k<info.zs+info.zm; k++) { 31 for (j=info.ys; j<info.ys+info.ym; j++) { 32 for (i=info.xs; i<info.xs+info.xm; i++) { 33 PetscScalar x = (Lx*i)/M; 34 PetscScalar y = (Ly*j)/N; 35 PetscScalar z = (Lz*k)/P; 36 va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2); 37 } 38 } 39 } 40 PetscCall(DMDAVecRestoreArray(da,v,&va)); 41 PetscCall(PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view)); 42 PetscCall(VecView(v,view)); 43 PetscCall(PetscViewerDestroy(&view)); 44 PetscCall(VecDestroy(&v)); 45 PetscCall(DMDestroy(&da)); 46 return 0; 47 } 48 49 /* 50 Write 2D DMDA vector with coordinates in VTK .vts format 51 52 */ 53 PetscErrorCode test_2d(const char filename[]) 54 { 55 MPI_Comm comm = MPI_COMM_WORLD; 56 const PetscInt M=10,N=20,dof=1,sw=1; 57 const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0; 58 DM da; 59 Vec v; 60 PetscViewer view; 61 DMDALocalInfo info; 62 PetscScalar **va; 63 PetscInt i,j; 64 65 PetscCall(DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR, M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da)); 66 PetscCall(DMSetFromOptions(da)); 67 PetscCall(DMSetUp(da)); 68 PetscCall(DMDASetUniformCoordinates(da,0.0,Lx,0.0,Ly,0.0,Lz)); 69 PetscCall(DMDAGetLocalInfo(da,&info)); 70 PetscCall(DMCreateGlobalVector(da,&v)); 71 PetscCall(DMDAVecGetArray(da,v,&va)); 72 for (j=info.ys; j<info.ys+info.ym; j++) { 73 for (i=info.xs; i<info.xs+info.xm; i++) { 74 PetscScalar x = (Lx*i)/M; 75 PetscScalar y = (Ly*j)/N; 76 va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2); 77 } 78 } 79 PetscCall(DMDAVecRestoreArray(da,v,&va)); 80 PetscCall(PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view)); 81 PetscCall(VecView(v,view)); 82 PetscCall(PetscViewerDestroy(&view)); 83 PetscCall(VecDestroy(&v)); 84 PetscCall(DMDestroy(&da)); 85 return 0; 86 } 87 88 /* 89 Write 2D DMDA vector without coordinates in VTK .vts format 90 91 */ 92 PetscErrorCode test_2d_nocoord(const char filename[]) 93 { 94 MPI_Comm comm = MPI_COMM_WORLD; 95 const PetscInt M=10,N=20,dof=1,sw=1; 96 const PetscScalar Lx=1.0,Ly=1.0; 97 DM da; 98 Vec v; 99 PetscViewer view; 100 DMDALocalInfo info; 101 PetscScalar **va; 102 PetscInt i,j; 103 104 PetscCall(DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,dof,sw,NULL,NULL,&da)); 105 PetscCall(DMSetFromOptions(da)); 106 PetscCall(DMSetUp(da)); 107 PetscCall(DMDAGetLocalInfo(da,&info)); 108 PetscCall(DMCreateGlobalVector(da,&v)); 109 PetscCall(DMDAVecGetArray(da,v,&va)); 110 for (j=info.ys; j<info.ys+info.ym; j++) { 111 for (i=info.xs; i<info.xs+info.xm; i++) { 112 PetscScalar x = (Lx*i)/M; 113 PetscScalar y = (Ly*j)/N; 114 va[j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2); 115 } 116 } 117 PetscCall(DMDAVecRestoreArray(da,v,&va)); 118 PetscCall(PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view)); 119 PetscCall(VecView(v,view)); 120 PetscCall(PetscViewerDestroy(&view)); 121 PetscCall(VecDestroy(&v)); 122 PetscCall(DMDestroy(&da)); 123 return 0; 124 } 125 126 /* 127 Write 3D DMDA vector without coordinates in VTK .vts format 128 129 */ 130 PetscErrorCode test_3d_nocoord(const char filename[]) 131 { 132 MPI_Comm comm = MPI_COMM_WORLD; 133 const PetscInt M=10,N=20,P=30,dof=1,sw=1; 134 const PetscScalar Lx=1.0,Ly=1.0,Lz=1.0; 135 DM da; 136 Vec v; 137 PetscViewer view; 138 DMDALocalInfo info; 139 PetscScalar ***va; 140 PetscInt i,j,k; 141 142 PetscCall(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)); 143 PetscCall(DMSetFromOptions(da)); 144 PetscCall(DMSetUp(da)); 145 146 PetscCall(DMDAGetLocalInfo(da,&info)); 147 PetscCall(DMCreateGlobalVector(da,&v)); 148 PetscCall(DMDAVecGetArray(da,v,&va)); 149 for (k=info.zs; k<info.zs+info.zm; k++) { 150 for (j=info.ys; j<info.ys+info.ym; j++) { 151 for (i=info.xs; i<info.xs+info.xm; i++) { 152 PetscScalar x = (Lx*i)/M; 153 PetscScalar y = (Ly*j)/N; 154 PetscScalar z = (Lz*k)/P; 155 va[k][j][i] = PetscPowScalarInt(x-0.5*Lx,2)+PetscPowScalarInt(y-0.5*Ly,2)+PetscPowScalarInt(z-0.5*Lz,2); 156 } 157 } 158 } 159 PetscCall(DMDAVecRestoreArray(da,v,&va)); 160 PetscCall(PetscViewerVTKOpen(comm,filename,FILE_MODE_WRITE,&view)); 161 PetscCall(VecView(v,view)); 162 PetscCall(PetscViewerDestroy(&view)); 163 PetscCall(VecDestroy(&v)); 164 PetscCall(DMDestroy(&da)); 165 return 0; 166 } 167 168 int main(int argc, char *argv[]) 169 { 170 171 PetscCall(PetscInitialize(&argc,&argv,0,help)); 172 PetscCall(test_3d("3d.vts")); 173 PetscCall(test_2d("2d.vts")); 174 PetscCall(test_2d_nocoord("2d_nocoord.vts")); 175 PetscCall(test_3d_nocoord("3d_nocoord.vts")); 176 PetscCall(PetscFinalize()); 177 return 0; 178 } 179 180 /*TEST 181 182 build: 183 requires: !complex 184 185 test: 186 nsize: 2 187 188 TEST*/ 189