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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 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 PETSC_SUCCESS; 166 } 167 168 int main(int argc, char *argv[]) 169 { 170 PetscFunctionBeginUser; 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