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