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