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