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