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