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 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)); 25 PetscCall(DMSetFromOptions(da)); 26 PetscCall(DMSetUp(da)); 27 28 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz)); 29 PetscCall(DMDAGetLocalInfo(da, &info)); 30 PetscCall(DMCreateGlobalVector(da, &v)); 31 PetscCall(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 PetscCall(DMDAVecRestoreArray(da, v, &va)); 43 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view)); 44 PetscCall(VecView(v, view)); 45 PetscCall(PetscViewerDestroy(&view)); 46 PetscCall(VecDestroy(&v)); 47 PetscCall(DMDestroy(&da)); 48 return PETSC_SUCCESS; 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 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, &da)); 68 PetscCall(DMSetFromOptions(da)); 69 PetscCall(DMSetUp(da)); 70 PetscCall(DMDASetUniformCoordinates(da, 0.0, Lx, 0.0, Ly, 0.0, Lz)); 71 PetscCall(DMDAGetLocalInfo(da, &info)); 72 PetscCall(DMCreateGlobalVector(da, &v)); 73 PetscCall(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 PetscCall(DMDAVecRestoreArray(da, v, &va)); 82 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view)); 83 PetscCall(VecView(v, view)); 84 PetscCall(PetscViewerDestroy(&view)); 85 PetscCall(VecDestroy(&v)); 86 PetscCall(DMDestroy(&da)); 87 return PETSC_SUCCESS; 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 PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, M, N, PETSC_DECIDE, PETSC_DECIDE, dof, sw, NULL, NULL, &da)); 107 PetscCall(DMSetFromOptions(da)); 108 PetscCall(DMSetUp(da)); 109 PetscCall(DMDAGetLocalInfo(da, &info)); 110 PetscCall(DMCreateGlobalVector(da, &v)); 111 PetscCall(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 PetscCall(DMDAVecRestoreArray(da, v, &va)); 120 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view)); 121 PetscCall(VecView(v, view)); 122 PetscCall(PetscViewerDestroy(&view)); 123 PetscCall(VecDestroy(&v)); 124 PetscCall(DMDestroy(&da)); 125 return PETSC_SUCCESS; 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 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)); 145 PetscCall(DMSetFromOptions(da)); 146 PetscCall(DMSetUp(da)); 147 148 PetscCall(DMDAGetLocalInfo(da, &info)); 149 PetscCall(DMCreateGlobalVector(da, &v)); 150 PetscCall(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 PetscCall(DMDAVecRestoreArray(da, v, &va)); 162 PetscCall(PetscViewerVTKOpen(comm, filename, FILE_MODE_WRITE, &view)); 163 PetscCall(VecView(v, view)); 164 PetscCall(PetscViewerDestroy(&view)); 165 PetscCall(VecDestroy(&v)); 166 PetscCall(DMDestroy(&da)); 167 return PETSC_SUCCESS; 168 } 169 170 int main(int argc, char *argv[]) 171 { 172 PetscFunctionBeginUser; 173 PetscCall(PetscInitialize(&argc, &argv, 0, help)); 174 PetscCall(test_3d("3d.vtr")); 175 PetscCall(test_2d("2d.vtr")); 176 PetscCall(test_2d_nocoord("2d_nocoord.vtr")); 177 PetscCall(test_3d_nocoord("3d_nocoord.vtr")); 178 PetscCall(PetscFinalize()); 179 return 0; 180 } 181 182 /*TEST 183 184 build: 185 requires: !complex 186 187 test: 188 nsize: 2 189 output_file: output/empty.out 190 191 TEST*/ 192