Actual source code: ex3.c

  1: static char help[] = "Check that a DM can accurately represent and interpolate functions of a given polynomial order\n\n";

  3: #include <petscdmplex.h>
  4: #include <petscdm.h>
  5: #include <petscdmda.h>
  6: #include <petscfe.h>
  7: #include <petscds.h>
  8: #include <petscksp.h>
  9: #include <petscsnes.h>

 11: typedef struct {
 12:   /* Domain and mesh definition */
 13:   PetscBool useDA;           /* Flag DMDA tensor product mesh */
 14:   PetscBool shearCoords;     /* Flag for shear transform */
 15:   PetscBool nonaffineCoords; /* Flag for non-affine transform */
 16:   /* Element definition */
 17:   PetscInt qorder;        /* Order of the quadrature */
 18:   PetscInt numComponents; /* Number of field components */
 19:   PetscFE  fe;            /* The finite element */
 20:   /* Testing space */
 21:   PetscInt  porder;         /* Order of polynomials to test */
 22:   PetscBool convergence;    /* Test for order of convergence */
 23:   PetscBool convRefine;     /* Test for convergence using refinement, otherwise use coarsening */
 24:   PetscBool constraints;    /* Test local constraints */
 25:   PetscBool tree;           /* Test tree routines */
 26:   PetscBool testFEjacobian; /* Test finite element Jacobian assembly */
 27:   PetscBool testFVgrad;     /* Test finite difference gradient routine */
 28:   PetscBool testInjector;   /* Test finite element injection routines */
 29:   PetscInt  treeCell;       /* Cell to refine in tree test */
 30:   PetscReal constants[3];   /* Constant values for each dimension */
 31: } AppCtx;

 33: /* u = 1 */
 34: PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 35: {
 36:   AppCtx  *user = (AppCtx *)ctx;
 37:   PetscInt d;
 38:   for (d = 0; d < dim; ++d) u[d] = user->constants[d];
 39:   return PETSC_SUCCESS;
 40: }
 41: PetscErrorCode constantDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 42: {
 43:   PetscInt d;
 44:   for (d = 0; d < dim; ++d) u[d] = 0.0;
 45:   return PETSC_SUCCESS;
 46: }

 48: /* u = x */
 49: PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 50: {
 51:   PetscInt d;
 52:   for (d = 0; d < dim; ++d) u[d] = coords[d];
 53:   return PETSC_SUCCESS;
 54: }
 55: PetscErrorCode linearDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 56: {
 57:   PetscInt d, e;
 58:   for (d = 0; d < dim; ++d) {
 59:     u[d] = 0.0;
 60:     for (e = 0; e < dim; ++e) u[d] += (d == e ? 1.0 : 0.0) * n[e];
 61:   }
 62:   return PETSC_SUCCESS;
 63: }

 65: /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */
 66: PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 67: {
 68:   if (dim > 2) {
 69:     u[0] = coords[0] * coords[1];
 70:     u[1] = coords[1] * coords[2];
 71:     u[2] = coords[2] * coords[0];
 72:   } else if (dim > 1) {
 73:     u[0] = coords[0] * coords[0];
 74:     u[1] = coords[0] * coords[1];
 75:   } else if (dim > 0) {
 76:     u[0] = coords[0] * coords[0];
 77:   }
 78:   return PETSC_SUCCESS;
 79: }
 80: PetscErrorCode quadraticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 81: {
 82:   if (dim > 2) {
 83:     u[0] = coords[1] * n[0] + coords[0] * n[1];
 84:     u[1] = coords[2] * n[1] + coords[1] * n[2];
 85:     u[2] = coords[2] * n[0] + coords[0] * n[2];
 86:   } else if (dim > 1) {
 87:     u[0] = 2.0 * coords[0] * n[0];
 88:     u[1] = coords[1] * n[0] + coords[0] * n[1];
 89:   } else if (dim > 0) {
 90:     u[0] = 2.0 * coords[0] * n[0];
 91:   }
 92:   return PETSC_SUCCESS;
 93: }

 95: /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
 96: PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 97: {
 98:   if (dim > 2) {
 99:     u[0] = coords[0] * coords[0] * coords[1];
100:     u[1] = coords[1] * coords[1] * coords[2];
101:     u[2] = coords[2] * coords[2] * coords[0];
102:   } else if (dim > 1) {
103:     u[0] = coords[0] * coords[0] * coords[0];
104:     u[1] = coords[0] * coords[0] * coords[1];
105:   } else if (dim > 0) {
106:     u[0] = coords[0] * coords[0] * coords[0];
107:   }
108:   return PETSC_SUCCESS;
109: }
110: PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
111: {
112:   if (dim > 2) {
113:     u[0] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1];
114:     u[1] = 2.0 * coords[1] * coords[2] * n[1] + coords[1] * coords[1] * n[2];
115:     u[2] = 2.0 * coords[2] * coords[0] * n[2] + coords[2] * coords[2] * n[0];
116:   } else if (dim > 1) {
117:     u[0] = 3.0 * coords[0] * coords[0] * n[0];
118:     u[1] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1];
119:   } else if (dim > 0) {
120:     u[0] = 3.0 * coords[0] * coords[0] * n[0];
121:   }
122:   return PETSC_SUCCESS;
123: }

125: /* u = tanh(x) */
126: PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
127: {
128:   PetscInt d;
129:   for (d = 0; d < dim; ++d) u[d] = PetscTanhReal(coords[d] - 0.5);
130:   return PETSC_SUCCESS;
131: }
132: PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
133: {
134:   PetscInt d;
135:   for (d = 0; d < dim; ++d) u[d] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d];
136:   return PETSC_SUCCESS;
137: }

139: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
140: {
141:   PetscInt n = 3;

143:   PetscFunctionBeginUser;
144:   options->useDA           = PETSC_FALSE;
145:   options->shearCoords     = PETSC_FALSE;
146:   options->nonaffineCoords = PETSC_FALSE;
147:   options->qorder          = 0;
148:   options->numComponents   = PETSC_DEFAULT;
149:   options->porder          = 0;
150:   options->convergence     = PETSC_FALSE;
151:   options->convRefine      = PETSC_TRUE;
152:   options->constraints     = PETSC_FALSE;
153:   options->tree            = PETSC_FALSE;
154:   options->treeCell        = 0;
155:   options->testFEjacobian  = PETSC_FALSE;
156:   options->testFVgrad      = PETSC_FALSE;
157:   options->testInjector    = PETSC_FALSE;
158:   options->constants[0]    = 1.0;
159:   options->constants[1]    = 2.0;
160:   options->constants[2]    = 3.0;

162:   PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
163:   PetscCall(PetscOptionsBool("-use_da", "Flag for DMDA mesh", "ex3.c", options->useDA, &options->useDA, NULL));
164:   PetscCall(PetscOptionsBool("-shear_coords", "Transform coordinates with a shear", "ex3.c", options->shearCoords, &options->shearCoords, NULL));
165:   PetscCall(PetscOptionsBool("-non_affine_coords", "Transform coordinates with a non-affine transform", "ex3.c", options->nonaffineCoords, &options->nonaffineCoords, NULL));
166:   PetscCall(PetscOptionsBoundedInt("-qorder", "The quadrature order", "ex3.c", options->qorder, &options->qorder, NULL, 0));
167:   PetscCall(PetscOptionsBoundedInt("-num_comp", "The number of field components", "ex3.c", options->numComponents, &options->numComponents, NULL, PETSC_DEFAULT));
168:   PetscCall(PetscOptionsBoundedInt("-porder", "The order of polynomials to test", "ex3.c", options->porder, &options->porder, NULL, 0));
169:   PetscCall(PetscOptionsBool("-convergence", "Check the convergence rate", "ex3.c", options->convergence, &options->convergence, NULL));
170:   PetscCall(PetscOptionsBool("-conv_refine", "Use refinement for the convergence rate", "ex3.c", options->convRefine, &options->convRefine, NULL));
171:   PetscCall(PetscOptionsBool("-constraints", "Test local constraints (serial only)", "ex3.c", options->constraints, &options->constraints, NULL));
172:   PetscCall(PetscOptionsBool("-tree", "Test tree routines", "ex3.c", options->tree, &options->tree, NULL));
173:   PetscCall(PetscOptionsBoundedInt("-tree_cell", "cell to refine in tree test", "ex3.c", options->treeCell, &options->treeCell, NULL, 0));
174:   PetscCall(PetscOptionsBool("-test_fe_jacobian", "Test finite element Jacobian assembly", "ex3.c", options->testFEjacobian, &options->testFEjacobian, NULL));
175:   PetscCall(PetscOptionsBool("-test_fv_grad", "Test finite volume gradient reconstruction", "ex3.c", options->testFVgrad, &options->testFVgrad, NULL));
176:   PetscCall(PetscOptionsBool("-test_injector", "Test finite element injection", "ex3.c", options->testInjector, &options->testInjector, NULL));
177:   PetscCall(PetscOptionsRealArray("-constants", "Set the constant values", "ex3.c", options->constants, &n, NULL));
178:   PetscOptionsEnd();

180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: static PetscErrorCode TransformCoordinates(DM dm, AppCtx *user)
184: {
185:   PetscSection coordSection;
186:   Vec          coordinates;
187:   PetscScalar *coords;
188:   PetscInt     vStart, vEnd, v;

190:   PetscFunctionBeginUser;
191:   if (user->nonaffineCoords) {
192:     /* x' = r^(1/p) (x/r), y' = r^(1/p) (y/r), z' = z */
193:     PetscCall(DMGetCoordinateSection(dm, &coordSection));
194:     PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
195:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
196:     PetscCall(VecGetArray(coordinates, &coords));
197:     for (v = vStart; v < vEnd; ++v) {
198:       PetscInt  dof, off;
199:       PetscReal p = 4.0, r;

201:       PetscCall(PetscSectionGetDof(coordSection, v, &dof));
202:       PetscCall(PetscSectionGetOffset(coordSection, v, &off));
203:       switch (dof) {
204:       case 2:
205:         r               = PetscSqr(PetscRealPart(coords[off + 0])) + PetscSqr(PetscRealPart(coords[off + 1]));
206:         coords[off + 0] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 0];
207:         coords[off + 1] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 1];
208:         break;
209:       case 3:
210:         r               = PetscSqr(PetscRealPart(coords[off + 0])) + PetscSqr(PetscRealPart(coords[off + 1]));
211:         coords[off + 0] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 0];
212:         coords[off + 1] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 1];
213:         coords[off + 2] = coords[off + 2];
214:         break;
215:       }
216:     }
217:     PetscCall(VecRestoreArray(coordinates, &coords));
218:   }
219:   if (user->shearCoords) {
220:     /* x' = x + m y + m z, y' = y + m z,  z' = z */
221:     PetscCall(DMGetCoordinateSection(dm, &coordSection));
222:     PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
223:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
224:     PetscCall(VecGetArray(coordinates, &coords));
225:     for (v = vStart; v < vEnd; ++v) {
226:       PetscInt  dof, off;
227:       PetscReal m = 1.0;

229:       PetscCall(PetscSectionGetDof(coordSection, v, &dof));
230:       PetscCall(PetscSectionGetOffset(coordSection, v, &off));
231:       switch (dof) {
232:       case 2:
233:         coords[off + 0] = coords[off + 0] + m * coords[off + 1];
234:         coords[off + 1] = coords[off + 1];
235:         break;
236:       case 3:
237:         coords[off + 0] = coords[off + 0] + m * coords[off + 1] + m * coords[off + 2];
238:         coords[off + 1] = coords[off + 1] + m * coords[off + 2];
239:         coords[off + 2] = coords[off + 2];
240:         break;
241:       }
242:     }
243:     PetscCall(VecRestoreArray(coordinates, &coords));
244:   }
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

248: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
249: {
250:   PetscInt  dim = 2;
251:   PetscBool simplex;

253:   PetscFunctionBeginUser;
254:   if (user->useDA) {
255:     switch (dim) {
256:     case 2:
257:       PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 2, 2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, dm));
258:       PetscCall(DMSetFromOptions(*dm));
259:       PetscCall(DMSetUp(*dm));
260:       PetscCall(DMDASetVertexCoordinates(*dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
261:       break;
262:     default:
263:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot create structured mesh of dimension %" PetscInt_FMT, dim);
264:     }
265:     PetscCall(PetscObjectSetName((PetscObject)*dm, "Hexahedral Mesh"));
266:   } else {
267:     PetscCall(DMCreate(comm, dm));
268:     PetscCall(DMSetType(*dm, DMPLEX));
269:     PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
270:     PetscCall(DMSetFromOptions(*dm));

272:     PetscCall(DMGetDimension(*dm, &dim));
273:     PetscCall(DMPlexIsSimplex(*dm, &simplex));
274:     PetscCallMPI(MPI_Bcast(&simplex, 1, MPIU_BOOL, 0, comm));
275:     if (user->tree) {
276:       DM refTree, ncdm = NULL;

278:       PetscCall(DMPlexCreateDefaultReferenceTree(comm, dim, simplex, &refTree));
279:       PetscCall(DMViewFromOptions(refTree, NULL, "-reftree_dm_view"));
280:       PetscCall(DMPlexSetReferenceTree(*dm, refTree));
281:       PetscCall(DMDestroy(&refTree));
282:       PetscCall(DMPlexTreeRefineCell(*dm, user->treeCell, &ncdm));
283:       if (ncdm) {
284:         PetscCall(DMDestroy(dm));
285:         *dm = ncdm;
286:         PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_FALSE));
287:       }
288:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "tree_"));
289:       PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
290:       PetscCall(DMSetFromOptions(*dm));
291:       PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
292:     } else {
293:       PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_TRUE));
294:     }
295:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "dist_"));
296:     PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
297:     PetscCall(DMSetFromOptions(*dm));
298:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, NULL));
299:     if (simplex) PetscCall(PetscObjectSetName((PetscObject)*dm, "Simplicial Mesh"));
300:     else PetscCall(PetscObjectSetName((PetscObject)*dm, "Hexahedral Mesh"));
301:   }
302:   PetscCall(DMSetFromOptions(*dm));
303:   PetscCall(TransformCoordinates(*dm, user));
304:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
305:   PetscFunctionReturn(PETSC_SUCCESS);
306: }

308: static void simple_mass(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
309: {
310:   PetscInt d, e;
311:   for (d = 0, e = 0; d < dim; d++, e += dim + 1) g0[e] = 1.;
312: }

314: /* < \nabla v, 1/2(\nabla u + {\nabla u}^T) > */
315: static void symmetric_gradient_inner_product(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar C[])
316: {
317:   PetscInt compI, compJ, d, e;

319:   for (compI = 0; compI < dim; ++compI) {
320:     for (compJ = 0; compJ < dim; ++compJ) {
321:       for (d = 0; d < dim; ++d) {
322:         for (e = 0; e < dim; e++) {
323:           if (d == e && d == compI && d == compJ) {
324:             C[((compI * dim + compJ) * dim + d) * dim + e] = 1.0;
325:           } else if ((d == compJ && e == compI) || (d == e && compI == compJ)) {
326:             C[((compI * dim + compJ) * dim + d) * dim + e] = 0.5;
327:           } else {
328:             C[((compI * dim + compJ) * dim + d) * dim + e] = 0.0;
329:           }
330:         }
331:       }
332:     }
333:   }
334: }

336: static PetscErrorCode SetupSection(DM dm, AppCtx *user)
337: {
338:   PetscFunctionBeginUser;
339:   if (user->constraints) {
340:     /* test local constraints */
341:     DM           coordDM;
342:     PetscInt     fStart, fEnd, f, vStart, vEnd, v;
343:     PetscInt     edgesx = 2, vertsx;
344:     PetscInt     edgesy = 2, vertsy;
345:     PetscMPIInt  size;
346:     PetscInt     numConst;
347:     PetscSection aSec;
348:     PetscInt    *anchors;
349:     PetscInt     offset;
350:     IS           aIS;
351:     MPI_Comm     comm = PetscObjectComm((PetscObject)dm);

353:     PetscCallMPI(MPI_Comm_size(comm, &size));
354:     PetscCheck(size <= 1, comm, PETSC_ERR_SUP, "Local constraint test can only be performed in serial");

356:     /* we are going to test constraints by using them to enforce periodicity
357:      * in one direction, and comparing to the existing method of enforcing
358:      * periodicity */

360:     /* first create the coordinate section so that it does not clone the
361:      * constraints */
362:     PetscCall(DMGetCoordinateDM(dm, &coordDM));

364:     /* create the constrained-to-anchor section */
365:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
366:     PetscCall(DMPlexGetDepthStratum(dm, 1, &fStart, &fEnd));
367:     PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &aSec));
368:     PetscCall(PetscSectionSetChart(aSec, PetscMin(fStart, vStart), PetscMax(fEnd, vEnd)));

370:     /* define the constraints */
371:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-da_grid_x", &edgesx, NULL));
372:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-da_grid_y", &edgesy, NULL));
373:     vertsx   = edgesx + 1;
374:     vertsy   = edgesy + 1;
375:     numConst = vertsy + edgesy;
376:     PetscCall(PetscMalloc1(numConst, &anchors));
377:     offset = 0;
378:     for (v = vStart + edgesx; v < vEnd; v += vertsx) {
379:       PetscCall(PetscSectionSetDof(aSec, v, 1));
380:       anchors[offset++] = v - edgesx;
381:     }
382:     for (f = fStart + edgesx * vertsy + edgesx * edgesy; f < fEnd; f++) {
383:       PetscCall(PetscSectionSetDof(aSec, f, 1));
384:       anchors[offset++] = f - edgesx * edgesy;
385:     }
386:     PetscCall(PetscSectionSetUp(aSec));
387:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numConst, anchors, PETSC_OWN_POINTER, &aIS));

389:     PetscCall(DMPlexSetAnchors(dm, aSec, aIS));
390:     PetscCall(PetscSectionDestroy(&aSec));
391:     PetscCall(ISDestroy(&aIS));
392:   }
393:   PetscCall(DMSetNumFields(dm, 1));
394:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)user->fe));
395:   PetscCall(DMCreateDS(dm));
396:   if (user->constraints) {
397:     /* test getting local constraint matrix that matches section */
398:     PetscSection aSec;
399:     IS           aIS;

401:     PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS));
402:     if (aSec) {
403:       PetscDS         ds;
404:       PetscSection    cSec, section;
405:       PetscInt        cStart, cEnd, c, numComp;
406:       Mat             cMat, mass;
407:       Vec             local;
408:       const PetscInt *anchors;

410:       PetscCall(DMGetLocalSection(dm, &section));
411:       /* this creates the matrix and preallocates the matrix structure: we
412:        * just have to fill in the values */
413:       PetscCall(DMGetDefaultConstraints(dm, &cSec, &cMat, NULL));
414:       PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
415:       PetscCall(ISGetIndices(aIS, &anchors));
416:       PetscCall(PetscFEGetNumComponents(user->fe, &numComp));
417:       for (c = cStart; c < cEnd; c++) {
418:         PetscInt cDof;

420:         /* is this point constrained? (does it have an anchor?) */
421:         PetscCall(PetscSectionGetDof(aSec, c, &cDof));
422:         if (cDof) {
423:           PetscInt cOff, a, aDof, aOff, j;
424:           PetscCheck(cDof == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Found %" PetscInt_FMT " anchor points: should be just one", cDof);

426:           /* find the anchor point */
427:           PetscCall(PetscSectionGetOffset(aSec, c, &cOff));
428:           a = anchors[cOff];

430:           /* find the constrained dofs (row in constraint matrix) */
431:           PetscCall(PetscSectionGetDof(cSec, c, &cDof));
432:           PetscCall(PetscSectionGetOffset(cSec, c, &cOff));

434:           /* find the anchor dofs (column in constraint matrix) */
435:           PetscCall(PetscSectionGetDof(section, a, &aDof));
436:           PetscCall(PetscSectionGetOffset(section, a, &aOff));

438:           PetscCheck(cDof == aDof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point and anchor have different number of dofs: %" PetscInt_FMT ", %" PetscInt_FMT, cDof, aDof);
439:           PetscCheck(cDof % numComp == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point dofs not divisible by field components: %" PetscInt_FMT ", %" PetscInt_FMT, cDof, numComp);

441:           /* put in a simple equality constraint */
442:           for (j = 0; j < cDof; j++) PetscCall(MatSetValue(cMat, cOff + j, aOff + j, 1., INSERT_VALUES));
443:         }
444:       }
445:       PetscCall(MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY));
446:       PetscCall(MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY));
447:       PetscCall(ISRestoreIndices(aIS, &anchors));

449:       /* Now that we have constructed the constraint matrix, any FE matrix
450:        * that we construct will apply the constraints during construction */

452:       PetscCall(DMCreateMatrix(dm, &mass));
453:       /* get a dummy local variable to serve as the solution */
454:       PetscCall(DMGetLocalVector(dm, &local));
455:       PetscCall(DMGetDS(dm, &ds));
456:       /* set the jacobian to be the mass matrix */
457:       PetscCall(PetscDSSetJacobian(ds, 0, 0, simple_mass, NULL, NULL, NULL));
458:       /* build the mass matrix */
459:       PetscCall(DMPlexSNESComputeJacobianFEM(dm, local, mass, mass, NULL));
460:       PetscCall(MatView(mass, PETSC_VIEWER_STDOUT_WORLD));
461:       PetscCall(MatDestroy(&mass));
462:       PetscCall(DMRestoreLocalVector(dm, &local));
463:     }
464:   }
465:   PetscFunctionReturn(PETSC_SUCCESS);
466: }

468: static PetscErrorCode TestFEJacobian(DM dm, AppCtx *user)
469: {
470:   PetscFunctionBeginUser;
471:   if (!user->useDA) {
472:     Vec          local;
473:     const Vec   *vecs;
474:     Mat          E;
475:     MatNullSpace sp;
476:     PetscBool    isNullSpace, hasConst;
477:     PetscInt     dim, n, i;
478:     Vec          res = NULL, localX, localRes;
479:     PetscDS      ds;

481:     PetscCall(DMGetDimension(dm, &dim));
482:     PetscCheck(user->numComponents == dim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "The number of components %" PetscInt_FMT " must be equal to the dimension %" PetscInt_FMT " for this test", user->numComponents, dim);
483:     PetscCall(DMGetDS(dm, &ds));
484:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, symmetric_gradient_inner_product));
485:     PetscCall(DMCreateMatrix(dm, &E));
486:     PetscCall(DMGetLocalVector(dm, &local));
487:     PetscCall(DMPlexSNESComputeJacobianFEM(dm, local, E, E, NULL));
488:     PetscCall(DMPlexCreateRigidBody(dm, 0, &sp));
489:     PetscCall(MatNullSpaceGetVecs(sp, &hasConst, &n, &vecs));
490:     if (n) PetscCall(VecDuplicate(vecs[0], &res));
491:     PetscCall(DMCreateLocalVector(dm, &localX));
492:     PetscCall(DMCreateLocalVector(dm, &localRes));
493:     for (i = 0; i < n; i++) { /* also test via matrix-free Jacobian application */
494:       PetscReal resNorm;

496:       PetscCall(VecSet(localRes, 0.));
497:       PetscCall(VecSet(localX, 0.));
498:       PetscCall(VecSet(local, 0.));
499:       PetscCall(VecSet(res, 0.));
500:       PetscCall(DMGlobalToLocalBegin(dm, vecs[i], INSERT_VALUES, localX));
501:       PetscCall(DMGlobalToLocalEnd(dm, vecs[i], INSERT_VALUES, localX));
502:       PetscCall(DMSNESComputeJacobianAction(dm, local, localX, localRes, NULL));
503:       PetscCall(DMLocalToGlobalBegin(dm, localRes, ADD_VALUES, res));
504:       PetscCall(DMLocalToGlobalEnd(dm, localRes, ADD_VALUES, res));
505:       PetscCall(VecNorm(res, NORM_2, &resNorm));
506:       if (resNorm > PETSC_SMALL) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Symmetric gradient action null space vector %" PetscInt_FMT " residual: %E\n", i, (double)resNorm));
507:     }
508:     PetscCall(VecDestroy(&localRes));
509:     PetscCall(VecDestroy(&localX));
510:     PetscCall(VecDestroy(&res));
511:     PetscCall(MatNullSpaceTest(sp, E, &isNullSpace));
512:     if (isNullSpace) {
513:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Symmetric gradient null space: PASS\n"));
514:     } else {
515:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Symmetric gradient null space: FAIL\n"));
516:     }
517:     PetscCall(MatNullSpaceDestroy(&sp));
518:     PetscCall(MatDestroy(&E));
519:     PetscCall(DMRestoreLocalVector(dm, &local));
520:   }
521:   PetscFunctionReturn(PETSC_SUCCESS);
522: }

524: static PetscErrorCode TestInjector(DM dm, AppCtx *user)
525: {
526:   DM          refTree;
527:   PetscMPIInt rank;

529:   PetscFunctionBegin;
530:   PetscCall(DMPlexGetReferenceTree(dm, &refTree));
531:   if (refTree) {
532:     Mat inj;

534:     PetscCall(DMPlexComputeInjectorReferenceTree(refTree, &inj));
535:     PetscCall(PetscObjectSetName((PetscObject)inj, "Reference Tree Injector"));
536:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
537:     if (rank == 0) PetscCall(MatView(inj, PETSC_VIEWER_STDOUT_SELF));
538:     PetscCall(MatDestroy(&inj));
539:   }
540:   PetscFunctionReturn(PETSC_SUCCESS);
541: }

543: static PetscErrorCode TestFVGrad(DM dm, AppCtx *user)
544: {
545:   MPI_Comm           comm;
546:   DM                 dmRedist, dmfv, dmgrad, dmCell, refTree;
547:   PetscFV            fv;
548:   PetscInt           dim, nvecs, v, cStart, cEnd, cEndInterior;
549:   PetscMPIInt        size;
550:   Vec                cellgeom, grad, locGrad;
551:   const PetscScalar *cgeom;
552:   PetscReal          allVecMaxDiff = 0., fvTol = 100. * PETSC_MACHINE_EPSILON;

554:   PetscFunctionBeginUser;
555:   comm = PetscObjectComm((PetscObject)dm);
556:   PetscCall(DMGetDimension(dm, &dim));
557:   /* duplicate DM, give dup. a FV discretization */
558:   PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
559:   PetscCallMPI(MPI_Comm_size(comm, &size));
560:   dmRedist = NULL;
561:   if (size > 1) PetscCall(DMPlexDistributeOverlap(dm, 1, NULL, &dmRedist));
562:   if (!dmRedist) {
563:     dmRedist = dm;
564:     PetscCall(PetscObjectReference((PetscObject)dmRedist));
565:   }
566:   PetscCall(PetscFVCreate(comm, &fv));
567:   PetscCall(PetscFVSetType(fv, PETSCFVLEASTSQUARES));
568:   PetscCall(PetscFVSetNumComponents(fv, user->numComponents));
569:   PetscCall(PetscFVSetSpatialDimension(fv, dim));
570:   PetscCall(PetscFVSetFromOptions(fv));
571:   PetscCall(PetscFVSetUp(fv));
572:   PetscCall(DMPlexConstructGhostCells(dmRedist, NULL, NULL, &dmfv));
573:   PetscCall(DMDestroy(&dmRedist));
574:   PetscCall(DMSetNumFields(dmfv, 1));
575:   PetscCall(DMSetField(dmfv, 0, NULL, (PetscObject)fv));
576:   PetscCall(DMCreateDS(dmfv));
577:   PetscCall(DMPlexGetReferenceTree(dm, &refTree));
578:   if (refTree) PetscCall(DMCopyDisc(dmfv, refTree));
579:   PetscCall(DMPlexGetGradientDM(dmfv, fv, &dmgrad));
580:   PetscCall(DMPlexGetHeightStratum(dmfv, 0, &cStart, &cEnd));
581:   nvecs = dim * (dim + 1) / 2;
582:   PetscCall(DMPlexGetGeometryFVM(dmfv, NULL, &cellgeom, NULL));
583:   PetscCall(VecGetDM(cellgeom, &dmCell));
584:   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
585:   PetscCall(DMGetGlobalVector(dmgrad, &grad));
586:   PetscCall(DMGetLocalVector(dmgrad, &locGrad));
587:   PetscCall(DMPlexGetCellTypeStratum(dmgrad, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
588:   cEndInterior = (cEndInterior < 0) ? cEnd : cEndInterior;
589:   for (v = 0; v < nvecs; v++) {
590:     Vec                locX;
591:     PetscInt           c;
592:     PetscScalar        trueGrad[3][3] = {{0.}};
593:     const PetscScalar *gradArray;
594:     PetscReal          maxDiff, maxDiffGlob;

596:     PetscCall(DMGetLocalVector(dmfv, &locX));
597:     /* get the local projection of the rigid body mode */
598:     for (c = cStart; c < cEnd; c++) {
599:       PetscFVCellGeom *cg;
600:       PetscScalar      cx[3] = {0., 0., 0.};

602:       PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
603:       if (v < dim) {
604:         cx[v] = 1.;
605:       } else {
606:         PetscInt w = v - dim;

608:         cx[(w + 1) % dim] = cg->centroid[(w + 2) % dim];
609:         cx[(w + 2) % dim] = -cg->centroid[(w + 1) % dim];
610:       }
611:       PetscCall(DMPlexVecSetClosure(dmfv, NULL, locX, c, cx, INSERT_ALL_VALUES));
612:     }
613:     /* TODO: this isn't in any header */
614:     PetscCall(DMPlexReconstructGradientsFVM(dmfv, locX, grad));
615:     PetscCall(DMGlobalToLocalBegin(dmgrad, grad, INSERT_VALUES, locGrad));
616:     PetscCall(DMGlobalToLocalEnd(dmgrad, grad, INSERT_VALUES, locGrad));
617:     PetscCall(VecGetArrayRead(locGrad, &gradArray));
618:     /* compare computed gradient to exact gradient */
619:     if (v >= dim) {
620:       PetscInt w = v - dim;

622:       trueGrad[(w + 1) % dim][(w + 2) % dim] = 1.;
623:       trueGrad[(w + 2) % dim][(w + 1) % dim] = -1.;
624:     }
625:     maxDiff = 0.;
626:     for (c = cStart; c < cEndInterior; c++) {
627:       PetscScalar *compGrad;
628:       PetscInt     i, j, k;
629:       PetscReal    FrobDiff = 0.;

631:       PetscCall(DMPlexPointLocalRead(dmgrad, c, gradArray, &compGrad));

633:       for (i = 0, k = 0; i < dim; i++) {
634:         for (j = 0; j < dim; j++, k++) {
635:           PetscScalar diff = compGrad[k] - trueGrad[i][j];
636:           FrobDiff += PetscRealPart(diff * PetscConj(diff));
637:         }
638:       }
639:       FrobDiff = PetscSqrtReal(FrobDiff);
640:       maxDiff  = PetscMax(maxDiff, FrobDiff);
641:     }
642:     PetscCall(MPIU_Allreduce(&maxDiff, &maxDiffGlob, 1, MPIU_REAL, MPIU_MAX, comm));
643:     allVecMaxDiff = PetscMax(allVecMaxDiff, maxDiffGlob);
644:     PetscCall(VecRestoreArrayRead(locGrad, &gradArray));
645:     PetscCall(DMRestoreLocalVector(dmfv, &locX));
646:   }
647:   if (allVecMaxDiff < fvTol) {
648:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Finite volume gradient reconstruction: PASS\n"));
649:   } else {
650:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Finite volume gradient reconstruction: FAIL at tolerance %g with max difference %g\n", (double)fvTol, (double)allVecMaxDiff));
651:   }
652:   PetscCall(DMRestoreLocalVector(dmgrad, &locGrad));
653:   PetscCall(DMRestoreGlobalVector(dmgrad, &grad));
654:   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
655:   PetscCall(DMDestroy(&dmfv));
656:   PetscCall(PetscFVDestroy(&fv));
657:   PetscFunctionReturn(PETSC_SUCCESS);
658: }

660: static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
661: {
662:   Vec       u;
663:   PetscReal n[3] = {1.0, 1.0, 1.0};

665:   PetscFunctionBeginUser;
666:   PetscCall(DMGetGlobalVector(dm, &u));
667:   /* Project function into FE function space */
668:   PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u));
669:   PetscCall(VecViewFromOptions(u, NULL, "-projection_view"));
670:   /* Compare approximation to exact in L_2 */
671:   PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error));
672:   PetscCall(DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer));
673:   PetscCall(DMRestoreGlobalVector(dm, &u));
674:   PetscFunctionReturn(PETSC_SUCCESS);
675: }

677: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
678: {
679:   PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
680:   PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
681:   void     *exactCtxs[3];
682:   MPI_Comm  comm;
683:   PetscReal error, errorDer, tol = PETSC_SMALL;

685:   PetscFunctionBeginUser;
686:   exactCtxs[0] = user;
687:   exactCtxs[1] = user;
688:   exactCtxs[2] = user;
689:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
690:   /* Setup functions to approximate */
691:   switch (order) {
692:   case 0:
693:     exactFuncs[0]    = constant;
694:     exactFuncDers[0] = constantDer;
695:     break;
696:   case 1:
697:     exactFuncs[0]    = linear;
698:     exactFuncDers[0] = linearDer;
699:     break;
700:   case 2:
701:     exactFuncs[0]    = quadratic;
702:     exactFuncDers[0] = quadraticDer;
703:     break;
704:   case 3:
705:     exactFuncs[0]    = cubic;
706:     exactFuncDers[0] = cubicDer;
707:     break;
708:   default:
709:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for order %" PetscInt_FMT, order);
710:   }
711:   PetscCall(ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
712:   /* Report result */
713:   if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
714:   else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
715:   if (errorDer > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer));
716:   else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
717:   PetscFunctionReturn(PETSC_SUCCESS);
718: }

720: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
721: {
722:   PetscErrorCode (*exactFuncs[1])(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
723:   PetscErrorCode (*exactFuncDers[1])(PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx);
724:   PetscReal n[3] = {1.0, 1.0, 1.0};
725:   void     *exactCtxs[3];
726:   DM        rdm, idm, fdm;
727:   Mat       Interp;
728:   Vec       iu, fu, scaling;
729:   MPI_Comm  comm;
730:   PetscInt  dim;
731:   PetscReal error, errorDer, tol = PETSC_SMALL;

733:   PetscFunctionBeginUser;
734:   exactCtxs[0] = user;
735:   exactCtxs[1] = user;
736:   exactCtxs[2] = user;
737:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
738:   PetscCall(DMGetDimension(dm, &dim));
739:   PetscCall(DMRefine(dm, comm, &rdm));
740:   PetscCall(DMSetCoarseDM(rdm, dm));
741:   PetscCall(DMPlexSetRegularRefinement(rdm, user->convRefine));
742:   if (user->tree) {
743:     DM refTree;
744:     PetscCall(DMPlexGetReferenceTree(dm, &refTree));
745:     PetscCall(DMPlexSetReferenceTree(rdm, refTree));
746:   }
747:   if (user->useDA) PetscCall(DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
748:   PetscCall(SetupSection(rdm, user));
749:   /* Setup functions to approximate */
750:   switch (order) {
751:   case 0:
752:     exactFuncs[0]    = constant;
753:     exactFuncDers[0] = constantDer;
754:     break;
755:   case 1:
756:     exactFuncs[0]    = linear;
757:     exactFuncDers[0] = linearDer;
758:     break;
759:   case 2:
760:     exactFuncs[0]    = quadratic;
761:     exactFuncDers[0] = quadraticDer;
762:     break;
763:   case 3:
764:     exactFuncs[0]    = cubic;
765:     exactFuncDers[0] = cubicDer;
766:     break;
767:   default:
768:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %" PetscInt_FMT " order %" PetscInt_FMT, dim, order);
769:   }
770:   idm = checkRestrict ? rdm : dm;
771:   fdm = checkRestrict ? dm : rdm;
772:   PetscCall(DMGetGlobalVector(idm, &iu));
773:   PetscCall(DMGetGlobalVector(fdm, &fu));
774:   PetscCall(DMSetApplicationContext(dm, user));
775:   PetscCall(DMSetApplicationContext(rdm, user));
776:   PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling));
777:   /* Project function into initial FE function space */
778:   PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu));
779:   /* Interpolate function into final FE function space */
780:   if (checkRestrict) {
781:     PetscCall(MatRestrict(Interp, iu, fu));
782:     PetscCall(VecPointwiseMult(fu, scaling, fu));
783:   } else PetscCall(MatInterpolate(Interp, iu, fu));
784:   /* Compare approximation to exact in L_2 */
785:   PetscCall(DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error));
786:   PetscCall(DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer));
787:   /* Report result */
788:   if (error > tol) PetscCall(PetscPrintf(comm, "Interpolation tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
789:   else PetscCall(PetscPrintf(comm, "Interpolation tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
790:   if (errorDer > tol) PetscCall(PetscPrintf(comm, "Interpolation tests FAIL for order %" PetscInt_FMT " derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer));
791:   else PetscCall(PetscPrintf(comm, "Interpolation tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
792:   PetscCall(DMRestoreGlobalVector(idm, &iu));
793:   PetscCall(DMRestoreGlobalVector(fdm, &fu));
794:   PetscCall(MatDestroy(&Interp));
795:   PetscCall(VecDestroy(&scaling));
796:   PetscCall(DMDestroy(&rdm));
797:   PetscFunctionReturn(PETSC_SUCCESS);
798: }

800: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
801: {
802:   DM odm = dm, rdm = NULL, cdm = NULL;
803:   PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)                         = {trig};
804:   PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx) = {trigDer};
805:   void     *exactCtxs[3];
806:   PetscInt  r, c, cStart, cEnd;
807:   PetscReal errorOld, errorDerOld, error, errorDer, rel, len, lenOld;
808:   double    p;

810:   PetscFunctionBeginUser;
811:   if (!user->convergence) PetscFunctionReturn(PETSC_SUCCESS);
812:   exactCtxs[0] = user;
813:   exactCtxs[1] = user;
814:   exactCtxs[2] = user;
815:   PetscCall(PetscObjectReference((PetscObject)odm));
816:   if (!user->convRefine) {
817:     for (r = 0; r < Nr; ++r) {
818:       PetscCall(DMRefine(odm, PetscObjectComm((PetscObject)dm), &rdm));
819:       PetscCall(DMDestroy(&odm));
820:       odm = rdm;
821:     }
822:     PetscCall(SetupSection(odm, user));
823:   }
824:   PetscCall(ComputeError(odm, exactFuncs, exactFuncDers, exactCtxs, &errorOld, &errorDerOld, user));
825:   if (user->convRefine) {
826:     for (r = 0; r < Nr; ++r) {
827:       PetscCall(DMRefine(odm, PetscObjectComm((PetscObject)dm), &rdm));
828:       if (user->useDA) PetscCall(DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
829:       PetscCall(SetupSection(rdm, user));
830:       PetscCall(ComputeError(rdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
831:       p = PetscLog2Real(errorOld / error);
832:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Function   convergence rate at refinement %" PetscInt_FMT ": %.2f\n", r, (double)p));
833:       p = PetscLog2Real(errorDerOld / errorDer);
834:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Derivative convergence rate at refinement %" PetscInt_FMT ": %.2f\n", r, (double)p));
835:       PetscCall(DMDestroy(&odm));
836:       odm         = rdm;
837:       errorOld    = error;
838:       errorDerOld = errorDer;
839:     }
840:   } else {
841:     /* PetscCall(ComputeLongestEdge(dm, &lenOld)); */
842:     PetscCall(DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd));
843:     lenOld = cEnd - cStart;
844:     for (c = 0; c < Nr; ++c) {
845:       PetscCall(DMCoarsen(odm, PetscObjectComm((PetscObject)dm), &cdm));
846:       if (user->useDA) PetscCall(DMDASetVertexCoordinates(cdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
847:       PetscCall(SetupSection(cdm, user));
848:       PetscCall(ComputeError(cdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
849:       /* PetscCall(ComputeLongestEdge(cdm, &len)); */
850:       PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
851:       len = cEnd - cStart;
852:       rel = error / errorOld;
853:       p   = PetscLogReal(rel) / PetscLogReal(lenOld / len);
854:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Function   convergence rate at coarsening %" PetscInt_FMT ": %.2f\n", c, (double)p));
855:       rel = errorDer / errorDerOld;
856:       p   = PetscLogReal(rel) / PetscLogReal(lenOld / len);
857:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Derivative convergence rate at coarsening %" PetscInt_FMT ": %.2f\n", c, (double)p));
858:       PetscCall(DMDestroy(&odm));
859:       odm         = cdm;
860:       errorOld    = error;
861:       errorDerOld = errorDer;
862:       lenOld      = len;
863:     }
864:   }
865:   PetscCall(DMDestroy(&odm));
866:   PetscFunctionReturn(PETSC_SUCCESS);
867: }

869: int main(int argc, char **argv)
870: {
871:   DM        dm;
872:   AppCtx    user; /* user-defined work context */
873:   PetscInt  dim     = 2;
874:   PetscBool simplex = PETSC_FALSE;

876:   PetscFunctionBeginUser;
877:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
878:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
879:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
880:   if (!user.useDA) {
881:     PetscCall(DMGetDimension(dm, &dim));
882:     PetscCall(DMPlexIsSimplex(dm, &simplex));
883:   }
884:   PetscCall(DMPlexMetricSetFromOptions(dm));
885:   user.numComponents = user.numComponents < 0 ? dim : user.numComponents;
886:   PetscCall(PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.numComponents, simplex, NULL, user.qorder, &user.fe));
887:   PetscCall(SetupSection(dm, &user));
888:   if (user.testFEjacobian) PetscCall(TestFEJacobian(dm, &user));
889:   if (user.testFVgrad) PetscCall(TestFVGrad(dm, &user));
890:   if (user.testInjector) PetscCall(TestInjector(dm, &user));
891:   PetscCall(CheckFunctions(dm, user.porder, &user));
892:   {
893:     PetscDualSpace dsp;
894:     PetscInt       k;

896:     PetscCall(PetscFEGetDualSpace(user.fe, &dsp));
897:     PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
898:     if (dim == 2 && simplex == PETSC_TRUE && user.tree == PETSC_FALSE && k == 0) {
899:       PetscCall(CheckInterpolation(dm, PETSC_FALSE, user.porder, &user));
900:       PetscCall(CheckInterpolation(dm, PETSC_TRUE, user.porder, &user));
901:     }
902:   }
903:   PetscCall(CheckConvergence(dm, 3, &user));
904:   PetscCall(PetscFEDestroy(&user.fe));
905:   PetscCall(DMDestroy(&dm));
906:   PetscCall(PetscFinalize());
907:   return 0;
908: }

910: /*TEST

912:   test:
913:     suffix: 1
914:     requires: triangle

916:   # 2D P_1 on a triangle
917:   test:
918:     suffix: p1_2d_0
919:     requires: triangle
920:     args: -petscspace_degree 1 -qorder 1 -convergence
921:   test:
922:     suffix: p1_2d_1
923:     requires: triangle
924:     args: -petscspace_degree 1 -qorder 1 -porder 1
925:   test:
926:     suffix: p1_2d_2
927:     requires: triangle
928:     args: -petscspace_degree 1 -qorder 1 -porder 2
929:   test:
930:     suffix: p1_2d_3
931:     requires: triangle mmg
932:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
933:   test:
934:     suffix: p1_2d_4
935:     requires: triangle mmg
936:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
937:   test:
938:     suffix: p1_2d_5
939:     requires: triangle mmg
940:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0

942:   # 3D P_1 on a tetrahedron
943:   test:
944:     suffix: p1_3d_0
945:     requires: ctetgen
946:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -convergence
947:   test:
948:     suffix: p1_3d_1
949:     requires: ctetgen
950:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 1
951:   test:
952:     suffix: p1_3d_2
953:     requires: ctetgen
954:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 2
955:   test:
956:     suffix: p1_3d_3
957:     requires: ctetgen mmg
958:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
959:   test:
960:     suffix: p1_3d_4
961:     requires: ctetgen mmg
962:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
963:   test:
964:     suffix: p1_3d_5
965:     requires: ctetgen mmg
966:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0

968:   # 2D P_2 on a triangle
969:   test:
970:     suffix: p2_2d_0
971:     requires: triangle
972:     args: -petscspace_degree 2 -qorder 2 -convergence
973:   test:
974:     suffix: p2_2d_1
975:     requires: triangle
976:     args: -petscspace_degree 2 -qorder 2 -porder 1
977:   test:
978:     suffix: p2_2d_2
979:     requires: triangle
980:     args: -petscspace_degree 2 -qorder 2 -porder 2
981:   test:
982:     suffix: p2_2d_3
983:     requires: triangle mmg
984:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
985:   test:
986:     suffix: p2_2d_4
987:     requires: triangle mmg
988:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
989:   test:
990:     suffix: p2_2d_5
991:     requires: triangle mmg
992:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0

994:   # 3D P_2 on a tetrahedron
995:   test:
996:     suffix: p2_3d_0
997:     requires: ctetgen
998:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -convergence
999:   test:
1000:     suffix: p2_3d_1
1001:     requires: ctetgen
1002:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 1
1003:   test:
1004:     suffix: p2_3d_2
1005:     requires: ctetgen
1006:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 2
1007:   test:
1008:     suffix: p2_3d_3
1009:     requires: ctetgen mmg
1010:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1011:   test:
1012:     suffix: p2_3d_4
1013:     requires: ctetgen mmg
1014:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1015:   test:
1016:     suffix: p2_3d_5
1017:     requires: ctetgen mmg
1018:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0

1020:   # 2D Q_1 on a quadrilaterial DA
1021:   test:
1022:     suffix: q1_2d_da_0
1023:     requires: broken
1024:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -convergence
1025:   test:
1026:     suffix: q1_2d_da_1
1027:     requires: broken
1028:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 1
1029:   test:
1030:     suffix: q1_2d_da_2
1031:     requires: broken
1032:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 2

1034:   # 2D Q_1 on a quadrilaterial Plex
1035:   test:
1036:     suffix: q1_2d_plex_0
1037:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -convergence
1038:   test:
1039:     suffix: q1_2d_plex_1
1040:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1
1041:   test:
1042:     suffix: q1_2d_plex_2
1043:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2
1044:   test:
1045:     suffix: q1_2d_plex_3
1046:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1 -shear_coords
1047:   test:
1048:     suffix: q1_2d_plex_4
1049:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2 -shear_coords
1050:   test:
1051:     suffix: q1_2d_plex_5
1052:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 0 -non_affine_coords -convergence
1053:   test:
1054:     suffix: q1_2d_plex_6
1055:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 1 -non_affine_coords -convergence
1056:   test:
1057:     suffix: q1_2d_plex_7
1058:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 2 -non_affine_coords -convergence

1060:   # 2D Q_2 on a quadrilaterial
1061:   test:
1062:     suffix: q2_2d_plex_0
1063:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1064:   test:
1065:     suffix: q2_2d_plex_1
1066:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1067:   test:
1068:     suffix: q2_2d_plex_2
1069:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1070:   test:
1071:     suffix: q2_2d_plex_3
1072:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1073:   test:
1074:     suffix: q2_2d_plex_4
1075:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1076:   test:
1077:     suffix: q2_2d_plex_5
1078:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 0 -non_affine_coords -convergence
1079:   test:
1080:     suffix: q2_2d_plex_6
1081:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 1 -non_affine_coords -convergence
1082:   test:
1083:     suffix: q2_2d_plex_7
1084:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 2 -non_affine_coords -convergence

1086:   # 2D P_3 on a triangle
1087:   test:
1088:     suffix: p3_2d_0
1089:     requires: triangle !single
1090:     args: -petscspace_degree 3 -qorder 3 -convergence
1091:   test:
1092:     suffix: p3_2d_1
1093:     requires: triangle !single
1094:     args: -petscspace_degree 3 -qorder 3 -porder 1
1095:   test:
1096:     suffix: p3_2d_2
1097:     requires: triangle !single
1098:     args: -petscspace_degree 3 -qorder 3 -porder 2
1099:   test:
1100:     suffix: p3_2d_3
1101:     requires: triangle !single
1102:     args: -petscspace_degree 3 -qorder 3 -porder 3
1103:   test:
1104:     suffix: p3_2d_4
1105:     requires: triangle mmg
1106:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -convergence -conv_refine 0
1107:   test:
1108:     suffix: p3_2d_5
1109:     requires: triangle mmg
1110:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 1 -conv_refine 0
1111:   test:
1112:     suffix: p3_2d_6
1113:     requires: triangle mmg
1114:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 3 -conv_refine 0

1116:   # 2D Q_3 on a quadrilaterial
1117:   test:
1118:     suffix: q3_2d_0
1119:     requires: !single
1120:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -convergence
1121:   test:
1122:     suffix: q3_2d_1
1123:     requires: !single
1124:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 1
1125:   test:
1126:     suffix: q3_2d_2
1127:     requires: !single
1128:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 2
1129:   test:
1130:     suffix: q3_2d_3
1131:     requires: !single
1132:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 3

1134:   # 2D P_1disc on a triangle/quadrilateral
1135:   test:
1136:     suffix: p1d_2d_0
1137:     requires: triangle
1138:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1139:   test:
1140:     suffix: p1d_2d_1
1141:     requires: triangle
1142:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1143:   test:
1144:     suffix: p1d_2d_2
1145:     requires: triangle
1146:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1147:   test:
1148:     suffix: p1d_2d_3
1149:     requires: triangle
1150:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1151:     filter: sed  -e "s/convergence rate at refinement 0: 2/convergence rate at refinement 0: 1.9/g"
1152:   test:
1153:     suffix: p1d_2d_4
1154:     requires: triangle
1155:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1156:   test:
1157:     suffix: p1d_2d_5
1158:     requires: triangle
1159:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2

1161:   # 2D BDM_1 on a triangle
1162:   test:
1163:     suffix: bdm1_2d_0
1164:     requires: triangle
1165:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1166:           -num_comp 2 -qorder 1 -convergence
1167:   test:
1168:     suffix: bdm1_2d_1
1169:     requires: triangle
1170:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1171:           -num_comp 2 -qorder 1 -porder 1
1172:   test:
1173:     suffix: bdm1_2d_2
1174:     requires: triangle
1175:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1176:           -num_comp 2 -qorder 1 -porder 2

1178:   # 2D BDM_1 on a quadrilateral
1179:   test:
1180:     suffix: bdm1q_2d_0
1181:     requires: triangle
1182:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1183:           -petscdualspace_lagrange_tensor 1 \
1184:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -convergence
1185:   test:
1186:     suffix: bdm1q_2d_1
1187:     requires: triangle
1188:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1189:           -petscdualspace_lagrange_tensor 1 \
1190:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 1
1191:   test:
1192:     suffix: bdm1q_2d_2
1193:     requires: triangle
1194:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1195:           -petscdualspace_lagrange_tensor 1 \
1196:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 2

1198:   # Test high order quadrature
1199:   test:
1200:     suffix: p1_quad_2
1201:     requires: triangle
1202:     args: -petscspace_degree 1 -qorder 2 -porder 1
1203:   test:
1204:     suffix: p1_quad_5
1205:     requires: triangle
1206:     args: -petscspace_degree 1 -qorder 5 -porder 1
1207:   test:
1208:     suffix: p2_quad_3
1209:     requires: triangle
1210:     args: -petscspace_degree 2 -qorder 3 -porder 2
1211:   test:
1212:     suffix: p2_quad_5
1213:     requires: triangle
1214:     args: -petscspace_degree 2 -qorder 5 -porder 2
1215:   test:
1216:     suffix: q1_quad_2
1217:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 2 -porder 1
1218:   test:
1219:     suffix: q1_quad_5
1220:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 5 -porder 1
1221:   test:
1222:     suffix: q2_quad_3
1223:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 3 -porder 1
1224:   test:
1225:     suffix: q2_quad_5
1226:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 5 -porder 1

1228:   # Nonconforming tests
1229:   test:
1230:     suffix: constraints
1231:     args: -dm_coord_space 0 -dm_plex_simplex 0 -petscspace_type tensor -petscspace_degree 1 -qorder 0 -constraints
1232:   test:
1233:     suffix: nonconforming_tensor_2
1234:     nsize: 4
1235:     args: -dist_dm_distribute -test_fe_jacobian -test_injector -petscpartitioner_type simple -tree -dm_plex_simplex 0 -dm_plex_max_projection_height 1 -petscspace_type tensor -petscspace_degree 2 -qorder 2 -dm_view ascii::ASCII_INFO_DETAIL
1236:   test:
1237:     suffix: nonconforming_tensor_3
1238:     nsize: 4
1239:     args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_max_projection_height 2 -petscspace_type tensor -petscspace_degree 1 -qorder 1 -dm_view ascii::ASCII_INFO_DETAIL
1240:   test:
1241:     suffix: nonconforming_tensor_2_fv
1242:     nsize: 4
1243:     args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_simplex 0 -num_comp 2
1244:   test:
1245:     suffix: nonconforming_tensor_3_fv
1246:     nsize: 4
1247:     args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -num_comp 3
1248:   test:
1249:     suffix: nonconforming_tensor_2_hi
1250:     requires: !single
1251:     nsize: 4
1252:     args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_simplex 0 -dm_plex_max_projection_height 1 -petscspace_type tensor -petscspace_degree 4 -qorder 4
1253:   test:
1254:     suffix: nonconforming_tensor_3_hi
1255:     requires: !single skip
1256:     nsize: 4
1257:     args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_max_projection_height 2 -petscspace_type tensor -petscspace_degree 4 -qorder 4
1258:   test:
1259:     suffix: nonconforming_simplex_2
1260:     requires: triangle
1261:     nsize: 4
1262:     args: -dist_dm_distribute -test_fe_jacobian -test_injector -petscpartitioner_type simple -tree -dm_plex_max_projection_height 1 -petscspace_degree 2 -qorder 2 -dm_view ascii::ASCII_INFO_DETAIL
1263:   test:
1264:     suffix: nonconforming_simplex_2_hi
1265:     requires: triangle !single
1266:     nsize: 4
1267:     args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_max_projection_height 1 -petscspace_degree 4 -qorder 4
1268:   test:
1269:     suffix: nonconforming_simplex_2_fv
1270:     requires: triangle
1271:     nsize: 4
1272:     args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -num_comp 2
1273:   test:
1274:     suffix: nonconforming_simplex_3
1275:     requires: ctetgen
1276:     nsize: 4
1277:     args: -dist_dm_distribute -test_fe_jacobian -test_injector -petscpartitioner_type simple -tree -dm_plex_dim 3 -dm_plex_max_projection_height 2 -petscspace_degree 2 -qorder 2 -dm_view ascii::ASCII_INFO_DETAIL
1278:   test:
1279:     suffix: nonconforming_simplex_3_hi
1280:     requires: ctetgen skip
1281:     nsize: 4
1282:     args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_dim 3 -dm_plex_max_projection_height 2 -petscspace_degree 4 -qorder 4
1283:   test:
1284:     suffix: nonconforming_simplex_3_fv
1285:     requires: ctetgen
1286:     nsize: 4
1287:     args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_dim 3 -num_comp 3

1289:   # 3D WXY on a triangular prism
1290:   test:
1291:     suffix: wxy_0
1292:     args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism -qorder 2 -porder 0 \
1293:           -petscspace_type sum \
1294:           -petscspace_variables 3 \
1295:           -petscspace_components 3 \
1296:           -petscspace_sum_spaces 2 \
1297:           -petscspace_sum_concatenate false \
1298:           -sumcomp_0_petscspace_variables 3 \
1299:           -sumcomp_0_petscspace_components 3 \
1300:           -sumcomp_0_petscspace_degree 1 \
1301:           -sumcomp_1_petscspace_variables 3 \
1302:           -sumcomp_1_petscspace_components 3 \
1303:           -sumcomp_1_petscspace_type wxy \
1304:           -petscdualspace_refcell triangular_prism \
1305:           -petscdualspace_form_degree 0 \
1306:           -petscdualspace_order 1 \
1307:           -petscdualspace_components 3

1309: TEST*/

1311: /*
1312:    # 2D Q_2 on a quadrilaterial Plex
1313:   test:
1314:     suffix: q2_2d_plex_0
1315:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1316:   test:
1317:     suffix: q2_2d_plex_1
1318:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1319:   test:
1320:     suffix: q2_2d_plex_2
1321:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1322:   test:
1323:     suffix: q2_2d_plex_3
1324:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1325:   test:
1326:     suffix: q2_2d_plex_4
1327:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1328:   test:
1329:     suffix: q2_2d_plex_5
1330:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 0 -non_affine_coords
1331:   test:
1332:     suffix: q2_2d_plex_6
1333:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 1 -non_affine_coords
1334:   test:
1335:     suffix: q2_2d_plex_7
1336:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 2 -non_affine_coords

1338:   test:
1339:     suffix: p1d_2d_6
1340:     requires: mmg
1341:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1342:   test:
1343:     suffix: p1d_2d_7
1344:     requires: mmg
1345:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1346:   test:
1347:     suffix: p1d_2d_8
1348:     requires: mmg
1349:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1350: */