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>
 10: #include <petscsf.h>

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

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

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

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

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

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

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

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

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

576:     PetscCall(DMCreateLabel(dmRedist, "Face Sets"));
577:     PetscCall(DMGetLabel(dmRedist, "Face Sets", &label));
578:     PetscCall(DMGetPointSF(dmRedist, &pointSF));
579:     PetscCall(PetscObjectReference((PetscObject)pointSF));
580:     PetscCall(DMSetPointSF(dmRedist, NULL));
581:     PetscCall(DMPlexMarkBoundaryFaces(dmRedist, 1, label));
582:     PetscCall(DMSetPointSF(dmRedist, pointSF));
583:     PetscCall(PetscSFDestroy(&pointSF));
584:   }
585:   PetscCall(DMPlexConstructGhostCells(dmRedist, NULL, NULL, &dmfv));
586:   PetscCall(DMDestroy(&dmRedist));
587:   PetscCall(DMSetNumFields(dmfv, 1));
588:   PetscCall(DMSetField(dmfv, 0, NULL, (PetscObject)fv));
589:   PetscCall(DMCreateDS(dmfv));
590:   PetscCall(DMPlexGetReferenceTree(dm, &refTree));
591:   if (refTree) PetscCall(DMCopyDisc(dmfv, refTree));
592:   PetscCall(DMPlexGetGradientDM(dmfv, fv, &dmgrad));
593:   PetscCall(DMPlexGetHeightStratum(dmfv, 0, &cStart, &cEnd));
594:   nvecs = dim * (dim + 1) / 2;
595:   PetscCall(DMPlexGetGeometryFVM(dmfv, NULL, &cellgeom, NULL));
596:   PetscCall(VecGetDM(cellgeom, &dmCell));
597:   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
598:   PetscCall(DMGetGlobalVector(dmgrad, &grad));
599:   PetscCall(DMGetLocalVector(dmgrad, &locGrad));
600:   PetscCall(DMPlexGetCellTypeStratum(dmgrad, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
601:   cEndInterior = (cEndInterior < 0) ? cEnd : cEndInterior;
602:   for (v = 0; v < nvecs; v++) {
603:     Vec                locX;
604:     PetscInt           c;
605:     PetscScalar        trueGrad[3][3] = {{0.}};
606:     const PetscScalar *gradArray;
607:     PetscReal          maxDiff, maxDiffGlob;

609:     PetscCall(DMGetLocalVector(dmfv, &locX));
610:     /* get the local projection of the rigid body mode */
611:     for (c = cStart; c < cEnd; c++) {
612:       PetscFVCellGeom *cg;
613:       PetscScalar      cx[3] = {0., 0., 0.};

615:       PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
616:       if (v < dim) {
617:         cx[v] = 1.;
618:       } else {
619:         PetscInt w = v - dim;

621:         cx[(w + 1) % dim] = cg->centroid[(w + 2) % dim];
622:         cx[(w + 2) % dim] = -cg->centroid[(w + 1) % dim];
623:       }
624:       PetscCall(DMPlexVecSetClosure(dmfv, NULL, locX, c, cx, INSERT_ALL_VALUES));
625:     }
626:     /* TODO: this isn't in any header */
627:     PetscCall(DMPlexReconstructGradientsFVM(dmfv, locX, grad));
628:     PetscCall(DMGlobalToLocalBegin(dmgrad, grad, INSERT_VALUES, locGrad));
629:     PetscCall(DMGlobalToLocalEnd(dmgrad, grad, INSERT_VALUES, locGrad));
630:     PetscCall(VecGetArrayRead(locGrad, &gradArray));
631:     /* compare computed gradient to exact gradient */
632:     if (v >= dim) {
633:       PetscInt w = v - dim;

635:       trueGrad[(w + 1) % dim][(w + 2) % dim] = 1.;
636:       trueGrad[(w + 2) % dim][(w + 1) % dim] = -1.;
637:     }
638:     maxDiff = 0.;
639:     for (c = cStart; c < cEndInterior; c++) {
640:       PetscScalar *compGrad;
641:       PetscInt     i, j, k;
642:       PetscReal    FrobDiff = 0.;

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

646:       for (i = 0, k = 0; i < dim; i++) {
647:         for (j = 0; j < dim; j++, k++) {
648:           PetscScalar diff = compGrad[k] - trueGrad[i][j];
649:           FrobDiff += PetscRealPart(diff * PetscConj(diff));
650:         }
651:       }
652:       FrobDiff = PetscSqrtReal(FrobDiff);
653:       maxDiff  = PetscMax(maxDiff, FrobDiff);
654:     }
655:     PetscCall(MPIU_Allreduce(&maxDiff, &maxDiffGlob, 1, MPIU_REAL, MPIU_MAX, comm));
656:     allVecMaxDiff = PetscMax(allVecMaxDiff, maxDiffGlob);
657:     PetscCall(VecRestoreArrayRead(locGrad, &gradArray));
658:     PetscCall(DMRestoreLocalVector(dmfv, &locX));
659:   }
660:   if (allVecMaxDiff < fvTol) {
661:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Finite volume gradient reconstruction: PASS\n"));
662:   } else {
663:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Finite volume gradient reconstruction: FAIL at tolerance %g with max difference %g\n", (double)fvTol, (double)allVecMaxDiff));
664:   }
665:   PetscCall(DMRestoreLocalVector(dmgrad, &locGrad));
666:   PetscCall(DMRestoreGlobalVector(dmgrad, &grad));
667:   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
668:   PetscCall(DMDestroy(&dmfv));
669:   PetscCall(PetscFVDestroy(&fv));
670:   PetscFunctionReturn(PETSC_SUCCESS);
671: }

673: 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)
674: {
675:   Vec       u;
676:   PetscReal n[3] = {1.0, 1.0, 1.0};

678:   PetscFunctionBeginUser;
679:   PetscCall(DMGetGlobalVector(dm, &u));
680:   /* Project function into FE function space */
681:   PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u));
682:   PetscCall(VecViewFromOptions(u, NULL, "-projection_view"));
683:   /* Compare approximation to exact in L_2 */
684:   PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error));
685:   PetscCall(DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer));
686:   PetscCall(DMRestoreGlobalVector(dm, &u));
687:   PetscFunctionReturn(PETSC_SUCCESS);
688: }

690: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
691: {
692:   PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
693:   PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
694:   void     *exactCtxs[3];
695:   MPI_Comm  comm;
696:   PetscReal error, errorDer, tol = PETSC_SMALL;

698:   PetscFunctionBeginUser;
699:   exactCtxs[0] = user;
700:   exactCtxs[1] = user;
701:   exactCtxs[2] = user;
702:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
703:   /* Setup functions to approximate */
704:   switch (order) {
705:   case 0:
706:     exactFuncs[0]    = constant;
707:     exactFuncDers[0] = constantDer;
708:     break;
709:   case 1:
710:     exactFuncs[0]    = linear;
711:     exactFuncDers[0] = linearDer;
712:     break;
713:   case 2:
714:     exactFuncs[0]    = quadratic;
715:     exactFuncDers[0] = quadraticDer;
716:     break;
717:   case 3:
718:     exactFuncs[0]    = cubic;
719:     exactFuncDers[0] = cubicDer;
720:     break;
721:   default:
722:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for order %" PetscInt_FMT, order);
723:   }
724:   PetscCall(ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
725:   /* Report result */
726:   if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
727:   else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
728:   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));
729:   else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
730:   PetscFunctionReturn(PETSC_SUCCESS);
731: }

733: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
734: {
735:   PetscErrorCode (*exactFuncs[1])(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
736:   PetscErrorCode (*exactFuncDers[1])(PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx);
737:   PetscReal n[3] = {1.0, 1.0, 1.0};
738:   void     *exactCtxs[3];
739:   DM        rdm, idm, fdm;
740:   Mat       Interp;
741:   Vec       iu, fu, scaling;
742:   MPI_Comm  comm;
743:   PetscInt  dim;
744:   PetscReal error, errorDer, tol = PETSC_SMALL;

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

813: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
814: {
815:   DM odm = dm, rdm = NULL, cdm = NULL;
816:   PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)                         = {trig};
817:   PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx) = {trigDer};
818:   void     *exactCtxs[3];
819:   PetscInt  r, c, cStart, cEnd;
820:   PetscReal errorOld, errorDerOld, error, errorDer, rel, len, lenOld;
821:   double    p;

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

882: int main(int argc, char **argv)
883: {
884:   DM        dm;
885:   AppCtx    user; /* user-defined work context */
886:   PetscInt  dim     = 2;
887:   PetscBool simplex = PETSC_FALSE;

889:   PetscFunctionBeginUser;
890:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
891:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
892:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
893:   if (!user.useDA) {
894:     PetscCall(DMGetDimension(dm, &dim));
895:     PetscCall(DMPlexIsSimplex(dm, &simplex));
896:   }
897:   PetscCall(DMPlexMetricSetFromOptions(dm));
898:   user.numComponents = user.numComponents < 0 ? dim : user.numComponents;
899:   PetscCall(PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.numComponents, simplex, NULL, user.qorder, &user.fe));
900:   PetscCall(SetupSection(dm, &user));
901:   if (user.testFEjacobian) PetscCall(TestFEJacobian(dm, &user));
902:   if (user.testFVgrad) PetscCall(TestFVGrad(dm, &user));
903:   if (user.testInjector) PetscCall(TestInjector(dm, &user));
904:   PetscCall(CheckFunctions(dm, user.porder, &user));
905:   {
906:     PetscDualSpace dsp;
907:     PetscInt       k;

909:     PetscCall(PetscFEGetDualSpace(user.fe, &dsp));
910:     PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
911:     if (dim == 2 && simplex == PETSC_TRUE && user.tree == PETSC_FALSE && k == 0) {
912:       PetscCall(CheckInterpolation(dm, PETSC_FALSE, user.porder, &user));
913:       PetscCall(CheckInterpolation(dm, PETSC_TRUE, user.porder, &user));
914:     }
915:   }
916:   PetscCall(CheckConvergence(dm, 3, &user));
917:   PetscCall(PetscFEDestroy(&user.fe));
918:   PetscCall(DMDestroy(&dm));
919:   PetscCall(PetscFinalize());
920:   return 0;
921: }

923: /*TEST

925:   test:
926:     suffix: 1
927:     requires: triangle

929:   # 2D P_1 on a triangle
930:   test:
931:     suffix: p1_2d_0
932:     requires: triangle
933:     args: -petscspace_degree 1 -qorder 1 -convergence
934:   test:
935:     suffix: p1_2d_1
936:     requires: triangle
937:     args: -petscspace_degree 1 -qorder 1 -porder 1
938:   test:
939:     suffix: p1_2d_2
940:     requires: triangle
941:     args: -petscspace_degree 1 -qorder 1 -porder 2
942:   test:
943:     suffix: p1_2d_3
944:     requires: triangle mmg
945:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
946:   test:
947:     suffix: p1_2d_4
948:     requires: triangle mmg
949:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
950:   test:
951:     suffix: p1_2d_5
952:     requires: triangle mmg
953:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0

955:   # 3D P_1 on a tetrahedron
956:   test:
957:     suffix: p1_3d_0
958:     requires: ctetgen
959:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -convergence
960:   test:
961:     suffix: p1_3d_1
962:     requires: ctetgen
963:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 1
964:   test:
965:     suffix: p1_3d_2
966:     requires: ctetgen
967:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 2
968:   test:
969:     suffix: p1_3d_3
970:     requires: ctetgen mmg
971:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
972:   test:
973:     suffix: p1_3d_4
974:     requires: ctetgen mmg
975:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
976:   test:
977:     suffix: p1_3d_5
978:     requires: ctetgen mmg
979:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0

981:   # 2D P_2 on a triangle
982:   test:
983:     suffix: p2_2d_0
984:     requires: triangle
985:     args: -petscspace_degree 2 -qorder 2 -convergence
986:   test:
987:     suffix: p2_2d_1
988:     requires: triangle
989:     args: -petscspace_degree 2 -qorder 2 -porder 1
990:   test:
991:     suffix: p2_2d_2
992:     requires: triangle
993:     args: -petscspace_degree 2 -qorder 2 -porder 2
994:   test:
995:     suffix: p2_2d_3
996:     requires: triangle mmg
997:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
998:   test:
999:     suffix: p2_2d_4
1000:     requires: triangle mmg
1001:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1002:   test:
1003:     suffix: p2_2d_5
1004:     requires: triangle mmg
1005:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0

1007:   # 3D P_2 on a tetrahedron
1008:   test:
1009:     suffix: p2_3d_0
1010:     requires: ctetgen
1011:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -convergence
1012:   test:
1013:     suffix: p2_3d_1
1014:     requires: ctetgen
1015:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 1
1016:   test:
1017:     suffix: p2_3d_2
1018:     requires: ctetgen
1019:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 2
1020:   test:
1021:     suffix: p2_3d_3
1022:     requires: ctetgen mmg
1023:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1024:   test:
1025:     suffix: p2_3d_4
1026:     requires: ctetgen mmg
1027:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1028:   test:
1029:     suffix: p2_3d_5
1030:     requires: ctetgen mmg
1031:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0

1033:   # 2D Q_1 on a quadrilaterial DA
1034:   test:
1035:     suffix: q1_2d_da_0
1036:     requires: broken
1037:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -convergence
1038:   test:
1039:     suffix: q1_2d_da_1
1040:     requires: broken
1041:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 1
1042:   test:
1043:     suffix: q1_2d_da_2
1044:     requires: broken
1045:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 2

1047:   # 2D Q_1 on a quadrilaterial Plex
1048:   test:
1049:     suffix: q1_2d_plex_0
1050:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -convergence
1051:   test:
1052:     suffix: q1_2d_plex_1
1053:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1
1054:   test:
1055:     suffix: q1_2d_plex_2
1056:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2
1057:   test:
1058:     suffix: q1_2d_plex_3
1059:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1 -shear_coords
1060:   test:
1061:     suffix: q1_2d_plex_4
1062:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2 -shear_coords
1063:   test:
1064:     suffix: q1_2d_plex_5
1065:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 0 -non_affine_coords -convergence
1066:   test:
1067:     suffix: q1_2d_plex_6
1068:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 1 -non_affine_coords -convergence
1069:   test:
1070:     suffix: q1_2d_plex_7
1071:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 2 -non_affine_coords -convergence

1073:   # 2D Q_2 on a quadrilaterial
1074:   test:
1075:     suffix: q2_2d_plex_0
1076:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1077:   test:
1078:     suffix: q2_2d_plex_1
1079:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1080:   test:
1081:     suffix: q2_2d_plex_2
1082:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1083:   test:
1084:     suffix: q2_2d_plex_3
1085:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1086:   test:
1087:     suffix: q2_2d_plex_4
1088:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1089:   test:
1090:     suffix: q2_2d_plex_5
1091:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 0 -non_affine_coords -convergence
1092:   test:
1093:     suffix: q2_2d_plex_6
1094:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 1 -non_affine_coords -convergence
1095:   test:
1096:     suffix: q2_2d_plex_7
1097:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 2 -non_affine_coords -convergence

1099:   # 2D P_3 on a triangle
1100:   test:
1101:     suffix: p3_2d_0
1102:     requires: triangle !single
1103:     args: -petscspace_degree 3 -qorder 3 -convergence
1104:   test:
1105:     suffix: p3_2d_1
1106:     requires: triangle !single
1107:     args: -petscspace_degree 3 -qorder 3 -porder 1
1108:   test:
1109:     suffix: p3_2d_2
1110:     requires: triangle !single
1111:     args: -petscspace_degree 3 -qorder 3 -porder 2
1112:   test:
1113:     suffix: p3_2d_3
1114:     requires: triangle !single
1115:     args: -petscspace_degree 3 -qorder 3 -porder 3
1116:   test:
1117:     suffix: p3_2d_4
1118:     requires: triangle mmg
1119:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -convergence -conv_refine 0
1120:   test:
1121:     suffix: p3_2d_5
1122:     requires: triangle mmg
1123:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 1 -conv_refine 0
1124:   test:
1125:     suffix: p3_2d_6
1126:     requires: triangle mmg
1127:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 3 -conv_refine 0

1129:   # 2D Q_3 on a quadrilaterial
1130:   test:
1131:     suffix: q3_2d_0
1132:     requires: !single
1133:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -convergence
1134:   test:
1135:     suffix: q3_2d_1
1136:     requires: !single
1137:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 1
1138:   test:
1139:     suffix: q3_2d_2
1140:     requires: !single
1141:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 2
1142:   test:
1143:     suffix: q3_2d_3
1144:     requires: !single
1145:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 3

1147:   # 2D P_1disc on a triangle/quadrilateral
1148:   test:
1149:     suffix: p1d_2d_0
1150:     requires: triangle
1151:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1152:   test:
1153:     suffix: p1d_2d_1
1154:     requires: triangle
1155:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1156:   test:
1157:     suffix: p1d_2d_2
1158:     requires: triangle
1159:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1160:   test:
1161:     suffix: p1d_2d_3
1162:     requires: triangle
1163:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1164:     filter: sed  -e "s/convergence rate at refinement 0: 2/convergence rate at refinement 0: 1.9/g"
1165:   test:
1166:     suffix: p1d_2d_4
1167:     requires: triangle
1168:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1169:   test:
1170:     suffix: p1d_2d_5
1171:     requires: triangle
1172:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2

1174:   # 2D BDM_1 on a triangle
1175:   test:
1176:     suffix: bdm1_2d_0
1177:     requires: triangle
1178:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1179:           -num_comp 2 -qorder 1 -convergence
1180:   test:
1181:     suffix: bdm1_2d_1
1182:     requires: triangle
1183:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1184:           -num_comp 2 -qorder 1 -porder 1
1185:   test:
1186:     suffix: bdm1_2d_2
1187:     requires: triangle
1188:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1189:           -num_comp 2 -qorder 1 -porder 2

1191:   # 2D BDM_1 on a quadrilateral
1192:   test:
1193:     suffix: bdm1q_2d_0
1194:     requires: triangle
1195:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1196:           -petscdualspace_lagrange_tensor 1 \
1197:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -convergence
1198:   test:
1199:     suffix: bdm1q_2d_1
1200:     requires: triangle
1201:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1202:           -petscdualspace_lagrange_tensor 1 \
1203:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 1
1204:   test:
1205:     suffix: bdm1q_2d_2
1206:     requires: triangle
1207:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1208:           -petscdualspace_lagrange_tensor 1 \
1209:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 2

1211:   # Test high order quadrature
1212:   test:
1213:     suffix: p1_quad_2
1214:     requires: triangle
1215:     args: -petscspace_degree 1 -qorder 2 -porder 1
1216:   test:
1217:     suffix: p1_quad_5
1218:     requires: triangle
1219:     args: -petscspace_degree 1 -qorder 5 -porder 1
1220:   test:
1221:     suffix: p2_quad_3
1222:     requires: triangle
1223:     args: -petscspace_degree 2 -qorder 3 -porder 2
1224:   test:
1225:     suffix: p2_quad_5
1226:     requires: triangle
1227:     args: -petscspace_degree 2 -qorder 5 -porder 2
1228:   test:
1229:     suffix: q1_quad_2
1230:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 2 -porder 1
1231:   test:
1232:     suffix: q1_quad_5
1233:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 5 -porder 1
1234:   test:
1235:     suffix: q2_quad_3
1236:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 3 -porder 1
1237:   test:
1238:     suffix: q2_quad_5
1239:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 5 -porder 1

1241:   # Nonconforming tests
1242:   test:
1243:     suffix: constraints
1244:     args: -dm_coord_space 0 -dm_plex_simplex 0 -petscspace_type tensor -petscspace_degree 1 -qorder 0 -constraints
1245:   test:
1246:     suffix: nonconforming_tensor_2
1247:     nsize: 4
1248:     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
1249:   test:
1250:     suffix: nonconforming_tensor_3
1251:     nsize: 4
1252:     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
1253:   test:
1254:     suffix: nonconforming_tensor_2_fv
1255:     nsize: 4
1256:     args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_simplex 0 -num_comp 2
1257:   test:
1258:     suffix: nonconforming_tensor_3_fv
1259:     nsize: 4
1260:     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
1261:   test:
1262:     suffix: nonconforming_tensor_2_hi
1263:     requires: !single
1264:     nsize: 4
1265:     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
1266:   test:
1267:     suffix: nonconforming_tensor_3_hi
1268:     requires: !single skip
1269:     nsize: 4
1270:     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
1271:   test:
1272:     suffix: nonconforming_simplex_2
1273:     requires: triangle
1274:     nsize: 4
1275:     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
1276:   test:
1277:     suffix: nonconforming_simplex_2_hi
1278:     requires: triangle !single
1279:     nsize: 4
1280:     args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_max_projection_height 1 -petscspace_degree 4 -qorder 4
1281:   test:
1282:     suffix: nonconforming_simplex_2_fv
1283:     requires: triangle
1284:     nsize: 4
1285:     args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -num_comp 2
1286:   test:
1287:     suffix: nonconforming_simplex_3
1288:     requires: ctetgen
1289:     nsize: 4
1290:     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
1291:   test:
1292:     suffix: nonconforming_simplex_3_hi
1293:     requires: ctetgen skip
1294:     nsize: 4
1295:     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
1296:   test:
1297:     suffix: nonconforming_simplex_3_fv
1298:     requires: ctetgen
1299:     nsize: 4
1300:     args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_dim 3 -num_comp 3

1302:   # 3D WXY on a triangular prism
1303:   test:
1304:     suffix: wxy_0
1305:     args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism -qorder 2 -porder 0 \
1306:           -petscspace_type sum \
1307:           -petscspace_variables 3 \
1308:           -petscspace_components 3 \
1309:           -petscspace_sum_spaces 2 \
1310:           -petscspace_sum_concatenate false \
1311:           -sumcomp_0_petscspace_variables 3 \
1312:           -sumcomp_0_petscspace_components 3 \
1313:           -sumcomp_0_petscspace_degree 1 \
1314:           -sumcomp_1_petscspace_variables 3 \
1315:           -sumcomp_1_petscspace_components 3 \
1316:           -sumcomp_1_petscspace_type wxy \
1317:           -petscdualspace_refcell triangular_prism \
1318:           -petscdualspace_form_degree 0 \
1319:           -petscdualspace_order 1 \
1320:           -petscdualspace_components 3

1322: TEST*/

1324: /*
1325:    # 2D Q_2 on a quadrilaterial Plex
1326:   test:
1327:     suffix: q2_2d_plex_0
1328:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1329:   test:
1330:     suffix: q2_2d_plex_1
1331:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1332:   test:
1333:     suffix: q2_2d_plex_2
1334:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1335:   test:
1336:     suffix: q2_2d_plex_3
1337:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1338:   test:
1339:     suffix: q2_2d_plex_4
1340:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1341:   test:
1342:     suffix: q2_2d_plex_5
1343:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 0 -non_affine_coords
1344:   test:
1345:     suffix: q2_2d_plex_6
1346:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 1 -non_affine_coords
1347:   test:
1348:     suffix: q2_2d_plex_7
1349:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 2 -non_affine_coords

1351:   test:
1352:     suffix: p1d_2d_6
1353:     requires: mmg
1354:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1355:   test:
1356:     suffix: p1d_2d_7
1357:     requires: mmg
1358:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1359:   test:
1360:     suffix: p1d_2d_8
1361:     requires: mmg
1362:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1363: */