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