1 static char help[] = "Test GLVis high-order support with DMDAs\n\n"; 2 3 #include <petscdm.h> 4 #include <petscdmda.h> 5 #include <petscdmplex.h> 6 #include <petscdt.h> 7 8 static PetscErrorCode MapPoint(PetscScalar xyz[], PetscScalar mxyz[]) 9 { 10 PetscScalar x, y, z, x2, y2, z2; 11 12 x = xyz[0]; 13 y = xyz[1]; 14 z = xyz[2]; 15 x2 = x * x; 16 y2 = y * y; 17 z2 = z * z; 18 mxyz[0] = x * PetscSqrtScalar(1.0 - y2 / 2.0 - z2 / 2.0 + y2 * z2 / 3.0); 19 mxyz[1] = y * PetscSqrtScalar(1.0 - z2 / 2.0 - x2 / 2.0 + z2 * x2 / 3.0); 20 mxyz[2] = z * PetscSqrtScalar(1.0 - x2 / 2.0 - y2 / 2.0 + x2 * y2 / 3.0); 21 return PETSC_SUCCESS; 22 } 23 24 static PetscErrorCode test_3d(PetscInt cells[], PetscBool plex, PetscBool ho) 25 { 26 DM dm; 27 Vec v; 28 PetscScalar *c; 29 PetscInt nl, i; 30 PetscReal u[3] = {1.0, 1.0, 1.0}, l[3] = {-1.0, -1.0, -1.0}; 31 32 PetscFunctionBeginUser; 33 if (ho) { 34 u[0] = 2. * cells[0]; 35 u[1] = 2. * cells[1]; 36 u[2] = 2. * cells[2]; 37 l[0] = 0.; 38 l[1] = 0.; 39 l[2] = 0.; 40 } 41 if (plex) { 42 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 43 PetscCall(DMSetType(dm, DMPLEX)); 44 PetscCall(DMSetFromOptions(dm)); 45 } else { 46 PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, cells[0] + 1, cells[1] + 1, cells[2] + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, NULL, &dm)); 47 } 48 PetscCall(DMSetUp(dm)); 49 if (!plex) PetscCall(DMDASetUniformCoordinates(dm, l[0], u[0], l[1], u[1], l[2], u[2])); 50 if (ho) { /* each element mapped to a sphere */ 51 DM cdm; 52 Vec cv; 53 PetscSection cSec; 54 DMDACoor3d ***_coords; 55 PetscScalar shift[3], *cptr; 56 PetscInt nel, dof = 3, nex, ney, nez, gx = 0, gy = 0, gz = 0; 57 PetscInt i, j, k, pi, pj, pk; 58 PetscReal *nodes, *weights; 59 char name[256]; 60 61 PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &dof, NULL)); 62 dof += 1; 63 64 PetscCall(PetscMalloc1(dof, &nodes)); 65 PetscCall(PetscMalloc1(dof, &weights)); 66 PetscCall(PetscDTGaussLobattoLegendreQuadrature(dof, PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, nodes, weights)); 67 PetscCall(DMGetCoordinatesLocal(dm, &cv)); 68 PetscCall(DMGetCoordinateDM(dm, &cdm)); 69 if (plex) { 70 PetscInt cEnd, cStart; 71 72 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 73 PetscCall(DMGetCoordinateSection(dm, &cSec)); 74 nel = cEnd - cStart; 75 nex = nel; 76 ney = 1; 77 nez = 1; 78 } else { 79 PetscCall(DMDAVecGetArray(cdm, cv, &_coords)); 80 PetscCall(DMDAGetElementsCorners(dm, &gx, &gy, &gz)); 81 PetscCall(DMDAGetElementsSizes(dm, &nex, &ney, &nez)); 82 nel = nex * ney * nez; 83 } 84 PetscCall(VecCreate(PETSC_COMM_WORLD, &v)); 85 PetscCall(VecSetSizes(v, 3 * dof * dof * dof * nel, PETSC_DECIDE)); 86 PetscCall(VecSetType(v, VECSTANDARD)); 87 PetscCall(VecGetArray(v, &c)); 88 cptr = c; 89 for (k = gz; k < gz + nez; k++) { 90 for (j = gy; j < gy + ney; j++) { 91 for (i = gx; i < gx + nex; i++) { 92 if (plex) { 93 PetscScalar *t = NULL; 94 95 PetscCall(DMPlexVecGetClosure(dm, cSec, cv, i, NULL, &t)); 96 shift[0] = t[0]; 97 shift[1] = t[1]; 98 shift[2] = t[2]; 99 PetscCall(DMPlexVecRestoreClosure(dm, cSec, cv, i, NULL, &t)); 100 } else { 101 shift[0] = _coords[k][j][i].x; 102 shift[1] = _coords[k][j][i].y; 103 shift[2] = _coords[k][j][i].z; 104 } 105 for (pk = 0; pk < dof; pk++) { 106 PetscScalar xyz[3]; 107 108 xyz[2] = nodes[pk]; 109 for (pj = 0; pj < dof; pj++) { 110 xyz[1] = nodes[pj]; 111 for (pi = 0; pi < dof; pi++) { 112 xyz[0] = nodes[pi]; 113 PetscCall(MapPoint(xyz, cptr)); 114 cptr[0] += shift[0]; 115 cptr[1] += shift[1]; 116 cptr[2] += shift[2]; 117 cptr += 3; 118 } 119 } 120 } 121 } 122 } 123 } 124 if (!plex) PetscCall(DMDAVecRestoreArray(cdm, cv, &_coords)); 125 PetscCall(VecRestoreArray(v, &c)); 126 PetscCall(PetscSNPrintf(name, sizeof(name), "FiniteElementCollection: L2_T1_3D_P%" PetscInt_FMT, dof - 1)); 127 PetscCall(PetscObjectSetName((PetscObject)v, name)); 128 PetscCall(PetscObjectCompose((PetscObject)dm, "_glvis_mesh_coords", (PetscObject)v)); 129 PetscCall(VecDestroy(&v)); 130 PetscCall(PetscFree(nodes)); 131 PetscCall(PetscFree(weights)); 132 PetscCall(DMViewFromOptions(dm, NULL, "-view")); 133 } else { /* map the whole domain to a sphere */ 134 PetscCall(DMGetCoordinates(dm, &v)); 135 PetscCall(VecGetLocalSize(v, &nl)); 136 PetscCall(VecGetArray(v, &c)); 137 for (i = 0; i < nl / 3; i++) PetscCall(MapPoint(c + 3 * i, c + 3 * i)); 138 PetscCall(VecRestoreArray(v, &c)); 139 PetscCall(DMSetCoordinates(dm, v)); 140 PetscCall(DMViewFromOptions(dm, NULL, "-view")); 141 } 142 PetscCall(DMDestroy(&dm)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 int main(int argc, char *argv[]) 147 { 148 PetscBool ho = PETSC_TRUE; 149 PetscBool plex = PETSC_FALSE; 150 PetscInt cells[3] = {2, 2, 2}; 151 152 PetscFunctionBeginUser; 153 PetscCall(PetscInitialize(&argc, &argv, 0, help)); 154 PetscCall(PetscOptionsGetBool(NULL, NULL, "-ho", &ho, NULL)); 155 PetscCall(PetscOptionsGetBool(NULL, NULL, "-plex", &plex, NULL)); 156 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nex", &cells[0], NULL)); 157 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ney", &cells[1], NULL)); 158 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nez", &cells[2], NULL)); 159 PetscCall(test_3d(cells, plex, ho)); 160 PetscCall(PetscFinalize()); 161 return 0; 162 } 163 164 /*TEST 165 166 testset: 167 nsize: 1 168 args: -view glvis: 169 170 test: 171 suffix: dmda_glvis_ho 172 args: -nex 1 -ney 1 -nez 1 173 174 test: 175 suffix: dmda_glvis_lo 176 args: -ho 0 177 178 test: 179 suffix: dmplex_glvis_ho 180 args: -nex 1 -ney 1 -nez 1 181 182 test: 183 suffix: dmplex_glvis_lo 184 args: -ho 0 185 186 TEST*/ 187