Actual source code: ex50.c

  1: static char help[] = "Test GLVis high-order support with DMDAs\n\n";

  3: #include <petscdm.h>
  4: #include <petscdmda.h>
  5: #include <petscdmplex.h>
  6: #include <petscdt.h>

  8: static PetscErrorCode MapPoint(PetscScalar xyz[], PetscScalar mxyz[])
  9: {
 10:   PetscScalar x, y, z, x2, y2, z2;

 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: }

 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};

 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];

 61:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-order", &dof, NULL));
 62:     dof += 1;

 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;

 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;

 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];

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: }

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};

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: }

164: /*TEST

166:    testset:
167:      nsize: 1
168:      args: -view glvis:

170:      test:
171:         suffix: dmda_glvis_ho
172:         args: -nex 1 -ney 1 -nez 1

174:      test:
175:         suffix: dmda_glvis_lo
176:         args: -ho 0

178:      test:
179:         suffix: dmplex_glvis_ho
180:         args: -nex 1 -ney 1 -nez 1

182:      test:
183:         suffix: dmplex_glvis_lo
184:         args: -ho 0

186: TEST*/