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 RT;             /* Test for Raviart-Thomas elements */
 24:   PetscBool convergence;    /* Test for order of convergence */
 25:   PetscBool convRefine;     /* Test for convergence using refinement, otherwise use coarsening */
 26:   PetscBool constraints;    /* Test local constraints */
 27:   PetscBool tree;           /* Test tree routines */
 28:   PetscBool testFEjacobian; /* Test finite element Jacobian assembly */
 29:   PetscBool testFVgrad;     /* Test finite difference gradient routine */
 30:   PetscBool testInjector;   /* Test finite element injection routines */
 31:   PetscInt  treeCell;       /* Cell to refine in tree test */
 32:   PetscReal constants[3];   /* Constant values for each dimension */
 33: } AppCtx;

 35: /*
 36: Derivatives are set as n_i \partial u_j / \partial x_i
 37: */

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

 54: /* RT_0: u = (1 + x, 1 + y) or (1 + x, 1 + y, 1 + z) */
 55: PetscErrorCode rt0(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 56: {
 57:   PetscInt d;
 58:   for (d = 0; d < dim; ++d) u[d] = 1.0 + coords[d];
 59:   return PETSC_SUCCESS;
 60: }

 62: PetscErrorCode rt0Der(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 63: {
 64:   PetscInt d, e;
 65:   for (d = 0; d < dim; ++d) {
 66:     u[d] = 0.0;
 67:     for (e = 0; e < dim; ++e) u[d] += (d == e ? 1.0 : 0.0) * n[e];
 68:   }
 69:   return PETSC_SUCCESS;
 70: }

 72: /* u = (x + y, y + x) or (x + z, 2y, z + x) */
 73: PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 74: {
 75:   PetscInt d;
 76:   for (d = 0; d < dim; ++d) u[d] = coords[d] + coords[dim - d - 1];
 77:   return PETSC_SUCCESS;
 78: }
 79: PetscErrorCode linearDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 80: {
 81:   PetscInt d, e;
 82:   for (d = 0; d < dim; ++d) {
 83:     u[d] = 0.0;
 84:     for (e = 0; e < dim; ++e) u[d] += ((d == e ? 1. : 0.) + (d == (dim - e - 1) ? 1. : 0.)) * n[e];
 85:   }
 86:   return PETSC_SUCCESS;
 87: }

 89: /* RT_1: u = (1 + x + y + x^2 + xy, 1 + x + y + xy + y^2) or (1 + x + y + z + x^2 + xy + xz, 1 + x + y + z + xy + y^2 + yz, 1 + x + y + z + xz + yz + z^2) */
 90: PetscErrorCode rt1(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 91: {
 92:   if (dim > 2) {
 93:     u[0] = 1.0 + coords[0] + coords[1] + coords[2] + coords[0] * coords[0] + coords[0] * coords[1] + coords[0] * coords[2];
 94:     u[1] = 1.0 + coords[0] + coords[1] + coords[2] + coords[0] * coords[1] + coords[1] * coords[1] + coords[1] * coords[2];
 95:     u[2] = 1.0 + coords[0] + coords[1] + coords[2] + coords[0] * coords[2] + coords[1] * coords[2] + coords[2] * coords[2];
 96:   } else if (dim > 1) {
 97:     u[0] = 1.0 + coords[0] + coords[1] + coords[0] * coords[0] + coords[0] * coords[1];
 98:     u[1] = 1.0 + coords[0] + coords[1] + coords[0] * coords[1] + coords[1] * coords[1];
 99:   }
100:   return PETSC_SUCCESS;
101: }

103: PetscErrorCode rt1Der(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
104: {
105:   if (dim > 2) {
106:     u[0] = (1.0 + 2.0 * coords[0] + coords[1] + coords[2]) * n[0] + (1.0 + coords[0]) * n[1] + (1.0 + coords[0]) * n[2];
107:     u[1] = (1.0 + coords[1]) * n[0] + (1.0 + coords[0] + 2.0 * coords[1] + coords[2]) * n[1] + (1.0 + coords[1]) * n[2];
108:     u[2] = (1.0 + coords[2]) * n[0] + (1.0 + coords[2]) * n[1] + (1.0 + coords[0] + coords[1] + 2.0 * coords[2]) * n[2];
109:   } else if (dim > 1) {
110:     u[0] = (1.0 + 2.0 * coords[0] + coords[1]) * n[0] + (1.0 + coords[0]) * n[1];
111:     u[1] = (1.0 + coords[1]) * n[0] + (1.0 + coords[0] + 2.0 * coords[1]) * n[1];
112:   }
113:   return PETSC_SUCCESS;
114: }

116: /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */
117: PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
118: {
119:   if (dim > 2) {
120:     u[0] = coords[0] * coords[1];
121:     u[1] = coords[1] * coords[2];
122:     u[2] = coords[2] * coords[0];
123:   } else if (dim > 1) {
124:     u[0] = coords[0] * coords[0];
125:     u[1] = coords[0] * coords[1];
126:   } else if (dim > 0) {
127:     u[0] = coords[0] * coords[0];
128:   }
129:   return PETSC_SUCCESS;
130: }
131: PetscErrorCode quadraticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
132: {
133:   if (dim > 2) {
134:     u[0] = coords[1] * n[0] + coords[0] * n[1];
135:     u[1] = coords[2] * n[1] + coords[1] * n[2];
136:     u[2] = coords[2] * n[0] + coords[0] * n[2];
137:   } else if (dim > 1) {
138:     u[0] = 2.0 * coords[0] * n[0];
139:     u[1] = coords[1] * n[0] + coords[0] * n[1];
140:   } else if (dim > 0) {
141:     u[0] = 2.0 * coords[0] * n[0];
142:   }
143:   return PETSC_SUCCESS;
144: }

146: /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
147: PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
148: {
149:   if (dim > 2) {
150:     u[0] = coords[0] * coords[0] * coords[1];
151:     u[1] = coords[1] * coords[1] * coords[2];
152:     u[2] = coords[2] * coords[2] * coords[0];
153:   } else if (dim > 1) {
154:     u[0] = coords[0] * coords[0] * coords[0];
155:     u[1] = coords[0] * coords[0] * coords[1];
156:   } else if (dim > 0) {
157:     u[0] = coords[0] * coords[0] * coords[0];
158:   }
159:   return PETSC_SUCCESS;
160: }
161: PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
162: {
163:   if (dim > 2) {
164:     u[0] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1];
165:     u[1] = 2.0 * coords[1] * coords[2] * n[1] + coords[1] * coords[1] * n[2];
166:     u[2] = 2.0 * coords[2] * coords[0] * n[2] + coords[2] * coords[2] * n[0];
167:   } else if (dim > 1) {
168:     u[0] = 3.0 * coords[0] * coords[0] * n[0];
169:     u[1] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1];
170:   } else if (dim > 0) {
171:     u[0] = 3.0 * coords[0] * coords[0] * n[0];
172:   }
173:   return PETSC_SUCCESS;
174: }

176: /* u = tanh(x) */
177: PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
178: {
179:   PetscInt d;
180:   for (d = 0; d < dim; ++d) u[d] = PetscTanhReal(coords[d] - 0.5);
181:   return PETSC_SUCCESS;
182: }
183: PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
184: {
185:   PetscInt d;
186:   for (d = 0; d < dim; ++d) u[d] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d];
187:   return PETSC_SUCCESS;
188: }

190: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
191: {
192:   PetscInt n = 3;

194:   PetscFunctionBeginUser;
195:   options->useDA           = PETSC_FALSE;
196:   options->shearCoords     = PETSC_FALSE;
197:   options->nonaffineCoords = PETSC_FALSE;
198:   options->qorder          = 0;
199:   options->numComponents   = PETSC_DEFAULT;
200:   options->porder          = 0;
201:   options->RT              = PETSC_FALSE;
202:   options->convergence     = PETSC_FALSE;
203:   options->convRefine      = PETSC_TRUE;
204:   options->constraints     = PETSC_FALSE;
205:   options->tree            = PETSC_FALSE;
206:   options->treeCell        = 0;
207:   options->testFEjacobian  = PETSC_FALSE;
208:   options->testFVgrad      = PETSC_FALSE;
209:   options->testInjector    = PETSC_FALSE;
210:   options->constants[0]    = 1.0;
211:   options->constants[1]    = 2.0;
212:   options->constants[2]    = 3.0;

214:   PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
215:   PetscCall(PetscOptionsBool("-use_da", "Flag for DMDA mesh", "ex3.c", options->useDA, &options->useDA, NULL));
216:   PetscCall(PetscOptionsBool("-shear_coords", "Transform coordinates with a shear", "ex3.c", options->shearCoords, &options->shearCoords, NULL));
217:   PetscCall(PetscOptionsBool("-non_affine_coords", "Transform coordinates with a non-affine transform", "ex3.c", options->nonaffineCoords, &options->nonaffineCoords, NULL));
218:   PetscCall(PetscOptionsBoundedInt("-qorder", "The quadrature order", "ex3.c", options->qorder, &options->qorder, NULL, 0));
219:   PetscCall(PetscOptionsBoundedInt("-num_comp", "The number of field components", "ex3.c", options->numComponents, &options->numComponents, NULL, PETSC_DEFAULT));
220:   PetscCall(PetscOptionsBoundedInt("-porder", "The order of polynomials to test", "ex3.c", options->porder, &options->porder, NULL, 0));
221:   PetscCall(PetscOptionsBool("-RT", "Use the Raviart-Thomas elements", "ex3.c", options->RT, &options->RT, NULL));
222:   PetscCall(PetscOptionsBool("-convergence", "Check the convergence rate", "ex3.c", options->convergence, &options->convergence, NULL));
223:   PetscCall(PetscOptionsBool("-conv_refine", "Use refinement for the convergence rate", "ex3.c", options->convRefine, &options->convRefine, NULL));
224:   PetscCall(PetscOptionsBool("-constraints", "Test local constraints (serial only)", "ex3.c", options->constraints, &options->constraints, NULL));
225:   PetscCall(PetscOptionsBool("-tree", "Test tree routines", "ex3.c", options->tree, &options->tree, NULL));
226:   PetscCall(PetscOptionsBoundedInt("-tree_cell", "cell to refine in tree test", "ex3.c", options->treeCell, &options->treeCell, NULL, 0));
227:   PetscCall(PetscOptionsBool("-test_fe_jacobian", "Test finite element Jacobian assembly", "ex3.c", options->testFEjacobian, &options->testFEjacobian, NULL));
228:   PetscCall(PetscOptionsBool("-test_fv_grad", "Test finite volume gradient reconstruction", "ex3.c", options->testFVgrad, &options->testFVgrad, NULL));
229:   PetscCall(PetscOptionsBool("-test_injector", "Test finite element injection", "ex3.c", options->testInjector, &options->testInjector, NULL));
230:   PetscCall(PetscOptionsRealArray("-constants", "Set the constant values", "ex3.c", options->constants, &n, NULL));
231:   PetscOptionsEnd();
232:   PetscFunctionReturn(PETSC_SUCCESS);
233: }

235: static PetscErrorCode TransformCoordinates(DM dm, AppCtx *user)
236: {
237:   PetscSection coordSection;
238:   Vec          coordinates;
239:   PetscScalar *coords;
240:   PetscInt     vStart, vEnd, v;

242:   PetscFunctionBeginUser;
243:   if (user->nonaffineCoords) {
244:     /* x' = r^(1/p) (x/r), y' = r^(1/p) (y/r), z' = z */
245:     PetscCall(DMGetCoordinateSection(dm, &coordSection));
246:     PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
247:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
248:     PetscCall(VecGetArray(coordinates, &coords));
249:     for (v = vStart; v < vEnd; ++v) {
250:       PetscInt  dof, off;
251:       PetscReal p = 4.0, r;

253:       PetscCall(PetscSectionGetDof(coordSection, v, &dof));
254:       PetscCall(PetscSectionGetOffset(coordSection, v, &off));
255:       switch (dof) {
256:       case 2:
257:         r               = PetscSqr(PetscRealPart(coords[off + 0])) + PetscSqr(PetscRealPart(coords[off + 1]));
258:         coords[off + 0] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 0];
259:         coords[off + 1] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 1];
260:         break;
261:       case 3:
262:         r               = PetscSqr(PetscRealPart(coords[off + 0])) + PetscSqr(PetscRealPart(coords[off + 1]));
263:         coords[off + 0] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 0];
264:         coords[off + 1] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 1];
265:         coords[off + 2] = coords[off + 2];
266:         break;
267:       }
268:     }
269:     PetscCall(VecRestoreArray(coordinates, &coords));
270:   }
271:   if (user->shearCoords) {
272:     /* x' = x + m y + m z, y' = y + m z,  z' = z */
273:     PetscCall(DMGetCoordinateSection(dm, &coordSection));
274:     PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
275:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
276:     PetscCall(VecGetArray(coordinates, &coords));
277:     for (v = vStart; v < vEnd; ++v) {
278:       PetscInt  dof, off;
279:       PetscReal m = 1.0;

281:       PetscCall(PetscSectionGetDof(coordSection, v, &dof));
282:       PetscCall(PetscSectionGetOffset(coordSection, v, &off));
283:       switch (dof) {
284:       case 2:
285:         coords[off + 0] = coords[off + 0] + m * coords[off + 1];
286:         coords[off + 1] = coords[off + 1];
287:         break;
288:       case 3:
289:         coords[off + 0] = coords[off + 0] + m * coords[off + 1] + m * coords[off + 2];
290:         coords[off + 1] = coords[off + 1] + m * coords[off + 2];
291:         coords[off + 2] = coords[off + 2];
292:         break;
293:       }
294:     }
295:     PetscCall(VecRestoreArray(coordinates, &coords));
296:   }
297:   PetscFunctionReturn(PETSC_SUCCESS);
298: }

300: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
301: {
302:   PetscInt  dim = 2;
303:   PetscBool simplex;

305:   PetscFunctionBeginUser;
306:   if (user->useDA) {
307:     switch (dim) {
308:     case 2:
309:       PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 2, 2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, dm));
310:       PetscCall(DMSetFromOptions(*dm));
311:       PetscCall(DMSetUp(*dm));
312:       PetscCall(DMDASetVertexCoordinates(*dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
313:       break;
314:     default:
315:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot create structured mesh of dimension %" PetscInt_FMT, dim);
316:     }
317:     PetscCall(PetscObjectSetName((PetscObject)*dm, "Hexahedral Mesh"));
318:   } else {
319:     PetscCall(DMCreate(comm, dm));
320:     PetscCall(DMSetType(*dm, DMPLEX));
321:     PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
322:     PetscCall(DMSetFromOptions(*dm));

324:     PetscCall(DMGetDimension(*dm, &dim));
325:     PetscCall(DMPlexIsSimplex(*dm, &simplex));
326:     PetscCallMPI(MPI_Bcast(&simplex, 1, MPIU_BOOL, 0, comm));
327:     if (user->tree) {
328:       DM refTree, ncdm = NULL;

330:       PetscCall(DMPlexCreateDefaultReferenceTree(comm, dim, simplex, &refTree));
331:       PetscCall(DMViewFromOptions(refTree, NULL, "-reftree_dm_view"));
332:       PetscCall(DMPlexSetReferenceTree(*dm, refTree));
333:       PetscCall(DMDestroy(&refTree));
334:       PetscCall(DMPlexTreeRefineCell(*dm, user->treeCell, &ncdm));
335:       if (ncdm) {
336:         PetscCall(DMDestroy(dm));
337:         *dm = ncdm;
338:         PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_FALSE));
339:       }
340:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "tree_"));
341:       PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
342:       PetscCall(DMSetFromOptions(*dm));
343:       PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
344:     } else {
345:       PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_TRUE));
346:     }
347:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "dist_"));
348:     PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
349:     PetscCall(DMSetFromOptions(*dm));
350:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, NULL));
351:     if (simplex) PetscCall(PetscObjectSetName((PetscObject)*dm, "Simplicial Mesh"));
352:     else PetscCall(PetscObjectSetName((PetscObject)*dm, "Hexahedral Mesh"));
353:   }
354:   PetscCall(DMSetFromOptions(*dm));
355:   PetscCall(TransformCoordinates(*dm, user));
356:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
357:   PetscFunctionReturn(PETSC_SUCCESS);
358: }

360: 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[])
361: {
362:   PetscInt d, e;
363:   for (d = 0, e = 0; d < dim; d++, e += dim + 1) g0[e] = 1.;
364: }

366: /* < \nabla v, 1/2(\nabla u + {\nabla u}^T) > */
367: 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[])
368: {
369:   PetscInt compI, compJ, d, e;

371:   for (compI = 0; compI < dim; ++compI) {
372:     for (compJ = 0; compJ < dim; ++compJ) {
373:       for (d = 0; d < dim; ++d) {
374:         for (e = 0; e < dim; e++) {
375:           if (d == e && d == compI && d == compJ) {
376:             C[((compI * dim + compJ) * dim + d) * dim + e] = 1.0;
377:           } else if ((d == compJ && e == compI) || (d == e && compI == compJ)) {
378:             C[((compI * dim + compJ) * dim + d) * dim + e] = 0.5;
379:           } else {
380:             C[((compI * dim + compJ) * dim + d) * dim + e] = 0.0;
381:           }
382:         }
383:       }
384:     }
385:   }
386: }

388: static PetscErrorCode SetupSection(DM dm, AppCtx *user)
389: {
390:   PetscFunctionBeginUser;
391:   if (user->constraints) {
392:     /* test local constraints */
393:     DM           coordDM;
394:     PetscInt     fStart, fEnd, f, vStart, vEnd, v;
395:     PetscInt     edgesx = 2, vertsx;
396:     PetscInt     edgesy = 2, vertsy;
397:     PetscMPIInt  size;
398:     PetscInt     numConst;
399:     PetscSection aSec;
400:     PetscInt    *anchors;
401:     PetscInt     offset;
402:     IS           aIS;
403:     MPI_Comm     comm = PetscObjectComm((PetscObject)dm);

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

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

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

416:     /* create the constrained-to-anchor section */
417:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
418:     PetscCall(DMPlexGetDepthStratum(dm, 1, &fStart, &fEnd));
419:     PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &aSec));
420:     PetscCall(PetscSectionSetChart(aSec, PetscMin(fStart, vStart), PetscMax(fEnd, vEnd)));

422:     /* define the constraints */
423:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-da_grid_x", &edgesx, NULL));
424:     PetscCall(PetscOptionsGetInt(NULL, NULL, "-da_grid_y", &edgesy, NULL));
425:     vertsx   = edgesx + 1;
426:     vertsy   = edgesy + 1;
427:     numConst = vertsy + edgesy;
428:     PetscCall(PetscMalloc1(numConst, &anchors));
429:     offset = 0;
430:     for (v = vStart + edgesx; v < vEnd; v += vertsx) {
431:       PetscCall(PetscSectionSetDof(aSec, v, 1));
432:       anchors[offset++] = v - edgesx;
433:     }
434:     for (f = fStart + edgesx * vertsy + edgesx * edgesy; f < fEnd; f++) {
435:       PetscCall(PetscSectionSetDof(aSec, f, 1));
436:       anchors[offset++] = f - edgesx * edgesy;
437:     }
438:     PetscCall(PetscSectionSetUp(aSec));
439:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numConst, anchors, PETSC_OWN_POINTER, &aIS));

441:     PetscCall(DMPlexSetAnchors(dm, aSec, aIS));
442:     PetscCall(PetscSectionDestroy(&aSec));
443:     PetscCall(ISDestroy(&aIS));
444:   }
445:   PetscCall(DMSetNumFields(dm, 1));
446:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)user->fe));
447:   PetscCall(DMCreateDS(dm));
448:   if (user->constraints) {
449:     /* test getting local constraint matrix that matches section */
450:     PetscSection aSec;
451:     IS           aIS;

453:     PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS));
454:     if (aSec) {
455:       PetscDS         ds;
456:       PetscSection    cSec, section;
457:       PetscInt        cStart, cEnd, c, numComp;
458:       Mat             cMat, mass;
459:       Vec             local;
460:       const PetscInt *anchors;

462:       PetscCall(DMGetLocalSection(dm, &section));
463:       /* this creates the matrix and preallocates the matrix nonzero structure: we
464:        * just have to fill in the values */
465:       PetscCall(DMGetDefaultConstraints(dm, &cSec, &cMat, NULL));
466:       PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
467:       PetscCall(ISGetIndices(aIS, &anchors));
468:       PetscCall(PetscFEGetNumComponents(user->fe, &numComp));
469:       for (c = cStart; c < cEnd; c++) {
470:         PetscInt cDof;

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

478:           /* find the anchor point */
479:           PetscCall(PetscSectionGetOffset(aSec, c, &cOff));
480:           a = anchors[cOff];

482:           /* find the constrained dofs (row in constraint matrix) */
483:           PetscCall(PetscSectionGetDof(cSec, c, &cDof));
484:           PetscCall(PetscSectionGetOffset(cSec, c, &cOff));

486:           /* find the anchor dofs (column in constraint matrix) */
487:           PetscCall(PetscSectionGetDof(section, a, &aDof));
488:           PetscCall(PetscSectionGetOffset(section, a, &aOff));

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

493:           /* put in a simple equality constraint */
494:           for (j = 0; j < cDof; j++) PetscCall(MatSetValue(cMat, cOff + j, aOff + j, 1., INSERT_VALUES));
495:         }
496:       }
497:       PetscCall(MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY));
498:       PetscCall(MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY));
499:       PetscCall(ISRestoreIndices(aIS, &anchors));

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

504:       PetscCall(DMCreateMatrix(dm, &mass));
505:       /* get a dummy local variable to serve as the solution */
506:       PetscCall(DMGetLocalVector(dm, &local));
507:       PetscCall(DMGetDS(dm, &ds));
508:       /* set the jacobian to be the mass matrix */
509:       PetscCall(PetscDSSetJacobian(ds, 0, 0, simple_mass, NULL, NULL, NULL));
510:       /* build the mass matrix */
511:       PetscCall(DMPlexSNESComputeJacobianFEM(dm, local, mass, mass, NULL));
512:       PetscCall(MatView(mass, PETSC_VIEWER_STDOUT_WORLD));
513:       PetscCall(MatDestroy(&mass));
514:       PetscCall(DMRestoreLocalVector(dm, &local));
515:     }
516:   }
517:   PetscFunctionReturn(PETSC_SUCCESS);
518: }

520: static PetscErrorCode TestFEJacobian(DM dm, AppCtx *user)
521: {
522:   PetscFunctionBeginUser;
523:   if (!user->useDA) {
524:     Vec          local;
525:     const Vec   *vecs;
526:     Mat          E;
527:     MatNullSpace sp;
528:     PetscBool    isNullSpace, hasConst;
529:     PetscInt     dim, n, i;
530:     Vec          res = NULL, localX, localRes;
531:     PetscDS      ds;

533:     PetscCall(DMGetDimension(dm, &dim));
534:     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);
535:     PetscCall(DMGetDS(dm, &ds));
536:     PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, symmetric_gradient_inner_product));
537:     PetscCall(DMCreateMatrix(dm, &E));
538:     PetscCall(DMGetLocalVector(dm, &local));
539:     PetscCall(DMPlexSNESComputeJacobianFEM(dm, local, E, E, NULL));
540:     PetscCall(DMPlexCreateRigidBody(dm, 0, &sp));
541:     PetscCall(MatNullSpaceGetVecs(sp, &hasConst, &n, &vecs));
542:     if (n) PetscCall(VecDuplicate(vecs[0], &res));
543:     PetscCall(DMCreateLocalVector(dm, &localX));
544:     PetscCall(DMCreateLocalVector(dm, &localRes));
545:     for (i = 0; i < n; i++) { /* also test via matrix-free Jacobian application */
546:       PetscReal resNorm;

548:       PetscCall(VecSet(localRes, 0.));
549:       PetscCall(VecSet(localX, 0.));
550:       PetscCall(VecSet(local, 0.));
551:       PetscCall(VecSet(res, 0.));
552:       PetscCall(DMGlobalToLocalBegin(dm, vecs[i], INSERT_VALUES, localX));
553:       PetscCall(DMGlobalToLocalEnd(dm, vecs[i], INSERT_VALUES, localX));
554:       PetscCall(DMSNESComputeJacobianAction(dm, local, localX, localRes, NULL));
555:       PetscCall(DMLocalToGlobalBegin(dm, localRes, ADD_VALUES, res));
556:       PetscCall(DMLocalToGlobalEnd(dm, localRes, ADD_VALUES, res));
557:       PetscCall(VecNorm(res, NORM_2, &resNorm));
558:       if (resNorm > PETSC_SMALL) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Symmetric gradient action null space vector %" PetscInt_FMT " residual: %E\n", i, (double)resNorm));
559:     }
560:     PetscCall(VecDestroy(&localRes));
561:     PetscCall(VecDestroy(&localX));
562:     PetscCall(VecDestroy(&res));
563:     PetscCall(MatNullSpaceTest(sp, E, &isNullSpace));
564:     if (isNullSpace) {
565:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Symmetric gradient null space: PASS\n"));
566:     } else {
567:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Symmetric gradient null space: FAIL\n"));
568:     }
569:     PetscCall(MatNullSpaceDestroy(&sp));
570:     PetscCall(MatDestroy(&E));
571:     PetscCall(DMRestoreLocalVector(dm, &local));
572:   }
573:   PetscFunctionReturn(PETSC_SUCCESS);
574: }

576: static PetscErrorCode TestInjector(DM dm, AppCtx *user)
577: {
578:   DM          refTree;
579:   PetscMPIInt rank;

581:   PetscFunctionBegin;
582:   PetscCall(DMPlexGetReferenceTree(dm, &refTree));
583:   if (refTree) {
584:     Mat inj;

586:     PetscCall(DMPlexComputeInjectorReferenceTree(refTree, &inj));
587:     PetscCall(PetscObjectSetName((PetscObject)inj, "Reference Tree Injector"));
588:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
589:     if (rank == 0) PetscCall(MatView(inj, PETSC_VIEWER_STDOUT_SELF));
590:     PetscCall(MatDestroy(&inj));
591:   }
592:   PetscFunctionReturn(PETSC_SUCCESS);
593: }

595: static PetscErrorCode TestFVGrad(DM dm, AppCtx *user)
596: {
597:   MPI_Comm           comm;
598:   DM                 dmRedist, dmfv, dmgrad, dmCell, refTree;
599:   PetscFV            fv;
600:   PetscInt           dim, nvecs, v, cStart, cEnd, cEndInterior;
601:   PetscMPIInt        size;
602:   Vec                cellgeom, grad, locGrad;
603:   const PetscScalar *cgeom;
604:   PetscReal          allVecMaxDiff = 0., fvTol = 100. * PETSC_MACHINE_EPSILON;

606:   PetscFunctionBeginUser;
607:   comm = PetscObjectComm((PetscObject)dm);
608:   PetscCall(DMGetDimension(dm, &dim));
609:   /* duplicate DM, give dup. a FV discretization */
610:   PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
611:   PetscCallMPI(MPI_Comm_size(comm, &size));
612:   dmRedist = NULL;
613:   if (size > 1) PetscCall(DMPlexDistributeOverlap(dm, 1, NULL, &dmRedist));
614:   if (!dmRedist) {
615:     dmRedist = dm;
616:     PetscCall(PetscObjectReference((PetscObject)dmRedist));
617:   }
618:   PetscCall(PetscFVCreate(comm, &fv));
619:   PetscCall(PetscFVSetType(fv, PETSCFVLEASTSQUARES));
620:   PetscCall(PetscFVSetNumComponents(fv, user->numComponents));
621:   PetscCall(PetscFVSetSpatialDimension(fv, dim));
622:   PetscCall(PetscFVSetFromOptions(fv));
623:   PetscCall(PetscFVSetUp(fv));
624:   {
625:     PetscSF pointSF;
626:     DMLabel label;

628:     PetscCall(DMCreateLabel(dmRedist, "Face Sets"));
629:     PetscCall(DMGetLabel(dmRedist, "Face Sets", &label));
630:     PetscCall(DMGetPointSF(dmRedist, &pointSF));
631:     PetscCall(PetscObjectReference((PetscObject)pointSF));
632:     PetscCall(DMSetPointSF(dmRedist, NULL));
633:     PetscCall(DMPlexMarkBoundaryFaces(dmRedist, 1, label));
634:     PetscCall(DMSetPointSF(dmRedist, pointSF));
635:     PetscCall(PetscSFDestroy(&pointSF));
636:   }
637:   PetscCall(DMPlexConstructGhostCells(dmRedist, NULL, NULL, &dmfv));
638:   PetscCall(DMDestroy(&dmRedist));
639:   PetscCall(DMSetNumFields(dmfv, 1));
640:   PetscCall(DMSetField(dmfv, 0, NULL, (PetscObject)fv));
641:   PetscCall(DMCreateDS(dmfv));
642:   PetscCall(DMPlexGetReferenceTree(dm, &refTree));
643:   if (refTree) PetscCall(DMCopyDisc(dmfv, refTree));
644:   PetscCall(DMPlexGetGradientDM(dmfv, fv, &dmgrad));
645:   PetscCall(DMPlexGetHeightStratum(dmfv, 0, &cStart, &cEnd));
646:   nvecs = dim * (dim + 1) / 2;
647:   PetscCall(DMPlexGetGeometryFVM(dmfv, NULL, &cellgeom, NULL));
648:   PetscCall(VecGetDM(cellgeom, &dmCell));
649:   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
650:   PetscCall(DMGetGlobalVector(dmgrad, &grad));
651:   PetscCall(DMGetLocalVector(dmgrad, &locGrad));
652:   PetscCall(DMPlexGetCellTypeStratum(dmgrad, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
653:   cEndInterior = (cEndInterior < 0) ? cEnd : cEndInterior;
654:   for (v = 0; v < nvecs; v++) {
655:     Vec                locX;
656:     PetscInt           c;
657:     PetscScalar        trueGrad[3][3] = {{0.}};
658:     const PetscScalar *gradArray;
659:     PetscReal          maxDiff, maxDiffGlob;

661:     PetscCall(DMGetLocalVector(dmfv, &locX));
662:     /* get the local projection of the rigid body mode */
663:     for (c = cStart; c < cEnd; c++) {
664:       PetscFVCellGeom *cg;
665:       PetscScalar      cx[3] = {0., 0., 0.};

667:       PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
668:       if (v < dim) {
669:         cx[v] = 1.;
670:       } else {
671:         PetscInt w = v - dim;

673:         cx[(w + 1) % dim] = cg->centroid[(w + 2) % dim];
674:         cx[(w + 2) % dim] = -cg->centroid[(w + 1) % dim];
675:       }
676:       PetscCall(DMPlexVecSetClosure(dmfv, NULL, locX, c, cx, INSERT_ALL_VALUES));
677:     }
678:     /* TODO: this isn't in any header */
679:     PetscCall(DMPlexReconstructGradientsFVM(dmfv, locX, grad));
680:     PetscCall(DMGlobalToLocalBegin(dmgrad, grad, INSERT_VALUES, locGrad));
681:     PetscCall(DMGlobalToLocalEnd(dmgrad, grad, INSERT_VALUES, locGrad));
682:     PetscCall(VecGetArrayRead(locGrad, &gradArray));
683:     /* compare computed gradient to exact gradient */
684:     if (v >= dim) {
685:       PetscInt w = v - dim;

687:       trueGrad[(w + 1) % dim][(w + 2) % dim] = 1.;
688:       trueGrad[(w + 2) % dim][(w + 1) % dim] = -1.;
689:     }
690:     maxDiff = 0.;
691:     for (c = cStart; c < cEndInterior; c++) {
692:       PetscScalar *compGrad;
693:       PetscInt     i, j, k;
694:       PetscReal    FrobDiff = 0.;

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

698:       for (i = 0, k = 0; i < dim; i++) {
699:         for (j = 0; j < dim; j++, k++) {
700:           PetscScalar diff = compGrad[k] - trueGrad[i][j];
701:           FrobDiff += PetscRealPart(diff * PetscConj(diff));
702:         }
703:       }
704:       FrobDiff = PetscSqrtReal(FrobDiff);
705:       maxDiff  = PetscMax(maxDiff, FrobDiff);
706:     }
707:     PetscCallMPI(MPIU_Allreduce(&maxDiff, &maxDiffGlob, 1, MPIU_REAL, MPIU_MAX, comm));
708:     allVecMaxDiff = PetscMax(allVecMaxDiff, maxDiffGlob);
709:     PetscCall(VecRestoreArrayRead(locGrad, &gradArray));
710:     PetscCall(DMRestoreLocalVector(dmfv, &locX));
711:   }
712:   if (allVecMaxDiff < fvTol) {
713:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Finite volume gradient reconstruction: PASS\n"));
714:   } else {
715:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Finite volume gradient reconstruction: FAIL at tolerance %g with max difference %g\n", (double)fvTol, (double)allVecMaxDiff));
716:   }
717:   PetscCall(DMRestoreLocalVector(dmgrad, &locGrad));
718:   PetscCall(DMRestoreGlobalVector(dmgrad, &grad));
719:   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
720:   PetscCall(DMDestroy(&dmfv));
721:   PetscCall(PetscFVDestroy(&fv));
722:   PetscFunctionReturn(PETSC_SUCCESS);
723: }

725: 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)
726: {
727:   Vec       u;
728:   PetscReal n[3] = {1.0, 1.0, 1.0};

730:   PetscFunctionBeginUser;
731:   PetscCall(DMGetGlobalVector(dm, &u));
732:   /* Project function into FE function space */
733:   PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u));
734:   PetscCall(VecViewFromOptions(u, NULL, "-projection_view"));
735:   /* Compare approximation to exact in L_2 */
736:   PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error));
737:   PetscCall(DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer));
738:   PetscCall(DMRestoreGlobalVector(dm, &u));
739:   PetscFunctionReturn(PETSC_SUCCESS);
740: }

742: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
743: {
744:   PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
745:   PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
746:   void     *exactCtxs[3];
747:   MPI_Comm  comm;
748:   PetscReal error, errorDer, tol = PETSC_SMALL;

750:   PetscFunctionBeginUser;
751:   exactCtxs[0] = user;
752:   exactCtxs[1] = user;
753:   exactCtxs[2] = user;
754:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
755:   /* Setup functions to approximate */
756:   switch (order) {
757:   case 0:
758:     if (user->RT) {
759:       exactFuncs[0]    = rt0;
760:       exactFuncDers[0] = rt0Der;
761:     } else {
762:       exactFuncs[0]    = constant;
763:       exactFuncDers[0] = constantDer;
764:     }
765:     break;
766:   case 1:
767:     if (user->RT) {
768:       exactFuncs[0]    = rt1;
769:       exactFuncDers[0] = rt1Der;
770:     } else {
771:       exactFuncs[0]    = linear;
772:       exactFuncDers[0] = linearDer;
773:     }
774:     break;
775:   case 2:
776:     exactFuncs[0]    = quadratic;
777:     exactFuncDers[0] = quadraticDer;
778:     break;
779:   case 3:
780:     exactFuncs[0]    = cubic;
781:     exactFuncDers[0] = cubicDer;
782:     break;
783:   default:
784:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for order %" PetscInt_FMT, order);
785:   }
786:   PetscCall(ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
787:   /* Report result */
788:   if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
789:   else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
790:   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));
791:   else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
792:   PetscFunctionReturn(PETSC_SUCCESS);
793: }

795: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
796: {
797:   PetscErrorCode (*exactFuncs[1])(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
798:   PetscErrorCode (*exactFuncDers[1])(PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx);
799:   PetscReal n[3] = {1.0, 1.0, 1.0};
800:   void     *exactCtxs[3];
801:   DM        rdm, idm, fdm;
802:   Mat       Interp;
803:   Vec       iu, fu, scaling;
804:   MPI_Comm  comm;
805:   PetscInt  dim;
806:   PetscReal error, errorDer, tol = PETSC_SMALL;

808:   PetscFunctionBeginUser;
809:   exactCtxs[0] = user;
810:   exactCtxs[1] = user;
811:   exactCtxs[2] = user;
812:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
813:   PetscCall(DMGetDimension(dm, &dim));
814:   PetscCall(DMRefine(dm, comm, &rdm));
815:   PetscCall(DMSetCoarseDM(rdm, dm));
816:   PetscCall(DMPlexSetRegularRefinement(rdm, user->convRefine));
817:   if (user->tree) {
818:     DM refTree;
819:     PetscCall(DMPlexGetReferenceTree(dm, &refTree));
820:     PetscCall(DMPlexSetReferenceTree(rdm, refTree));
821:   }
822:   if (user->useDA) PetscCall(DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
823:   PetscCall(SetupSection(rdm, user));
824:   /* Setup functions to approximate */
825:   switch (order) {
826:   case 0:
827:     exactFuncs[0]    = constant;
828:     exactFuncDers[0] = constantDer;
829:     break;
830:   case 1:
831:     exactFuncs[0]    = linear;
832:     exactFuncDers[0] = linearDer;
833:     break;
834:   case 2:
835:     exactFuncs[0]    = quadratic;
836:     exactFuncDers[0] = quadraticDer;
837:     break;
838:   case 3:
839:     exactFuncs[0]    = cubic;
840:     exactFuncDers[0] = cubicDer;
841:     break;
842:   default:
843:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %" PetscInt_FMT " order %" PetscInt_FMT, dim, order);
844:   }
845:   idm = checkRestrict ? rdm : dm;
846:   fdm = checkRestrict ? dm : rdm;
847:   PetscCall(DMGetGlobalVector(idm, &iu));
848:   PetscCall(DMGetGlobalVector(fdm, &fu));
849:   PetscCall(DMSetApplicationContext(dm, user));
850:   PetscCall(DMSetApplicationContext(rdm, user));
851:   PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling));
852:   /* Project function into initial FE function space */
853:   PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu));
854:   /* Interpolate function into final FE function space */
855:   if (checkRestrict) {
856:     PetscCall(MatRestrict(Interp, iu, fu));
857:     PetscCall(VecPointwiseMult(fu, scaling, fu));
858:   } else PetscCall(MatInterpolate(Interp, iu, fu));
859:   /* Compare approximation to exact in L_2 */
860:   PetscCall(DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error));
861:   PetscCall(DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer));
862:   /* Report result */
863:   if (error > tol) PetscCall(PetscPrintf(comm, "Interpolation tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
864:   else PetscCall(PetscPrintf(comm, "Interpolation tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
865:   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));
866:   else PetscCall(PetscPrintf(comm, "Interpolation tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
867:   PetscCall(DMRestoreGlobalVector(idm, &iu));
868:   PetscCall(DMRestoreGlobalVector(fdm, &fu));
869:   PetscCall(MatDestroy(&Interp));
870:   PetscCall(VecDestroy(&scaling));
871:   PetscCall(DMDestroy(&rdm));
872:   PetscFunctionReturn(PETSC_SUCCESS);
873: }

875: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
876: {
877:   DM odm = dm, rdm = NULL, cdm = NULL;
878:   PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)                         = {trig};
879:   PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx) = {trigDer};
880:   void     *exactCtxs[3];
881:   PetscInt  r, c, cStart, cEnd;
882:   PetscReal errorOld, errorDerOld, error, errorDer, rel, len, lenOld;
883:   double    p;

885:   PetscFunctionBeginUser;
886:   if (!user->convergence) PetscFunctionReturn(PETSC_SUCCESS);
887:   exactCtxs[0] = user;
888:   exactCtxs[1] = user;
889:   exactCtxs[2] = user;
890:   PetscCall(PetscObjectReference((PetscObject)odm));
891:   if (!user->convRefine) {
892:     for (r = 0; r < Nr; ++r) {
893:       PetscCall(DMRefine(odm, PetscObjectComm((PetscObject)dm), &rdm));
894:       PetscCall(DMDestroy(&odm));
895:       odm = rdm;
896:     }
897:     PetscCall(SetupSection(odm, user));
898:   }
899:   PetscCall(ComputeError(odm, exactFuncs, exactFuncDers, exactCtxs, &errorOld, &errorDerOld, user));
900:   if (user->convRefine) {
901:     for (r = 0; r < Nr; ++r) {
902:       PetscCall(DMRefine(odm, PetscObjectComm((PetscObject)dm), &rdm));
903:       if (user->useDA) PetscCall(DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
904:       PetscCall(SetupSection(rdm, user));
905:       PetscCall(ComputeError(rdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
906:       p = PetscLog2Real(errorOld / error);
907:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Function   convergence rate at refinement %" PetscInt_FMT ": %.2f\n", r, (double)p));
908:       p = PetscLog2Real(errorDerOld / errorDer);
909:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Derivative convergence rate at refinement %" PetscInt_FMT ": %.2f\n", r, (double)p));
910:       PetscCall(DMDestroy(&odm));
911:       odm         = rdm;
912:       errorOld    = error;
913:       errorDerOld = errorDer;
914:     }
915:   } else {
916:     /* PetscCall(ComputeLongestEdge(dm, &lenOld)); */
917:     PetscCall(DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd));
918:     lenOld = cEnd - cStart;
919:     for (c = 0; c < Nr; ++c) {
920:       PetscCall(DMCoarsen(odm, PetscObjectComm((PetscObject)dm), &cdm));
921:       if (user->useDA) PetscCall(DMDASetVertexCoordinates(cdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
922:       PetscCall(SetupSection(cdm, user));
923:       PetscCall(ComputeError(cdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
924:       /* PetscCall(ComputeLongestEdge(cdm, &len)); */
925:       PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
926:       len = cEnd - cStart;
927:       rel = error / errorOld;
928:       p   = PetscLogReal(rel) / PetscLogReal(lenOld / len);
929:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Function   convergence rate at coarsening %" PetscInt_FMT ": %.2f\n", c, (double)p));
930:       rel = errorDer / errorDerOld;
931:       p   = PetscLogReal(rel) / PetscLogReal(lenOld / len);
932:       PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Derivative convergence rate at coarsening %" PetscInt_FMT ": %.2f\n", c, (double)p));
933:       PetscCall(DMDestroy(&odm));
934:       odm         = cdm;
935:       errorOld    = error;
936:       errorDerOld = errorDer;
937:       lenOld      = len;
938:     }
939:   }
940:   PetscCall(DMDestroy(&odm));
941:   PetscFunctionReturn(PETSC_SUCCESS);
942: }

944: int main(int argc, char **argv)
945: {
946:   DM        dm;
947:   AppCtx    user; /* user-defined work context */
948:   PetscInt  dim     = 2;
949:   PetscBool simplex = PETSC_FALSE;

951:   PetscFunctionBeginUser;
952:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
953:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
954:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
955:   if (!user.useDA) {
956:     PetscCall(DMGetDimension(dm, &dim));
957:     PetscCall(DMPlexIsSimplex(dm, &simplex));
958:   }
959:   PetscCall(DMPlexMetricSetFromOptions(dm));
960:   user.numComponents = user.numComponents < 0 ? dim : user.numComponents;
961:   PetscCall(PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.numComponents, simplex, NULL, user.qorder, &user.fe));
962:   PetscCall(SetupSection(dm, &user));
963:   if (user.testFEjacobian) PetscCall(TestFEJacobian(dm, &user));
964:   if (user.testFVgrad) PetscCall(TestFVGrad(dm, &user));
965:   if (user.testInjector) PetscCall(TestInjector(dm, &user));
966:   PetscCall(CheckFunctions(dm, user.porder, &user));
967:   {
968:     PetscDualSpace dsp;
969:     PetscInt       k;

971:     PetscCall(PetscFEGetDualSpace(user.fe, &dsp));
972:     PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
973:     if (dim == 2 && simplex == PETSC_TRUE && user.tree == PETSC_FALSE && k == 0) {
974:       PetscCall(CheckInterpolation(dm, PETSC_FALSE, user.porder, &user));
975:       PetscCall(CheckInterpolation(dm, PETSC_TRUE, user.porder, &user));
976:     }
977:   }
978:   PetscCall(CheckConvergence(dm, 3, &user));
979:   PetscCall(PetscFEDestroy(&user.fe));
980:   PetscCall(DMDestroy(&dm));
981:   PetscCall(PetscFinalize());
982:   return 0;
983: }

985: /*TEST

987:   test:
988:     suffix: 1
989:     requires: triangle

991:   # 2D P_1 on a triangle
992:   test:
993:     suffix: p1_2d_0
994:     requires: triangle
995:     args: -petscspace_degree 1 -qorder 1 -convergence
996:   test:
997:     suffix: p1_2d_1
998:     requires: triangle
999:     args: -petscspace_degree 1 -qorder 1 -porder 1
1000:   test:
1001:     suffix: p1_2d_2
1002:     requires: triangle
1003:     args: -petscspace_degree 1 -qorder 1 -porder 2
1004:   test:
1005:     suffix: p1_2d_3
1006:     requires: triangle mmg
1007:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1008:   test:
1009:     suffix: p1_2d_4
1010:     requires: triangle mmg
1011:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1012:   test:
1013:     suffix: p1_2d_5
1014:     requires: triangle mmg
1015:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0

1017:   # 3D P_1 on a tetrahedron
1018:   test:
1019:     suffix: p1_3d_0
1020:     requires: ctetgen
1021:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -convergence
1022:   test:
1023:     suffix: p1_3d_1
1024:     requires: ctetgen
1025:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 1
1026:   test:
1027:     suffix: p1_3d_2
1028:     requires: ctetgen
1029:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 2
1030:   test:
1031:     suffix: p1_3d_3
1032:     requires: ctetgen mmg
1033:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1034:   test:
1035:     suffix: p1_3d_4
1036:     requires: ctetgen mmg
1037:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1038:   test:
1039:     suffix: p1_3d_5
1040:     requires: ctetgen mmg
1041:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0

1043:   # 2D P_2 on a triangle
1044:   test:
1045:     suffix: p2_2d_0
1046:     requires: triangle
1047:     args: -petscspace_degree 2 -qorder 2 -convergence
1048:   test:
1049:     suffix: p2_2d_1
1050:     requires: triangle
1051:     args: -petscspace_degree 2 -qorder 2 -porder 1
1052:   test:
1053:     suffix: p2_2d_2
1054:     requires: triangle
1055:     args: -petscspace_degree 2 -qorder 2 -porder 2
1056:   test:
1057:     suffix: p2_2d_3
1058:     requires: triangle mmg
1059:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1060:   test:
1061:     suffix: p2_2d_4
1062:     requires: triangle mmg
1063:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1064:   test:
1065:     suffix: p2_2d_5
1066:     requires: triangle mmg
1067:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0

1069:   # 3D P_2 on a tetrahedron
1070:   test:
1071:     suffix: p2_3d_0
1072:     requires: ctetgen
1073:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -convergence
1074:   test:
1075:     suffix: p2_3d_1
1076:     requires: ctetgen
1077:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 1
1078:   test:
1079:     suffix: p2_3d_2
1080:     requires: ctetgen
1081:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 2
1082:   test:
1083:     suffix: p2_3d_3
1084:     requires: ctetgen mmg
1085:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1086:   test:
1087:     suffix: p2_3d_4
1088:     requires: ctetgen mmg
1089:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1090:   test:
1091:     suffix: p2_3d_5
1092:     requires: ctetgen mmg
1093:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0

1095:   # 2D Q_1 on a quadrilaterial DA
1096:   test:
1097:     suffix: q1_2d_da_0
1098:     TODO: broken
1099:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -convergence
1100:   test:
1101:     suffix: q1_2d_da_1
1102:     TODO: broken
1103:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 1
1104:   test:
1105:     suffix: q1_2d_da_2
1106:     TODO: broken
1107:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 2

1109:   # 2D Q_1 on a quadrilaterial Plex
1110:   test:
1111:     suffix: q1_2d_plex_0
1112:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -convergence
1113:   test:
1114:     suffix: q1_2d_plex_1
1115:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1
1116:   test:
1117:     suffix: q1_2d_plex_2
1118:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2
1119:   test:
1120:     suffix: q1_2d_plex_3
1121:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1 -shear_coords
1122:   test:
1123:     suffix: q1_2d_plex_4
1124:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2 -shear_coords
1125:   test:
1126:     suffix: q1_2d_plex_5
1127:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 0 -non_affine_coords -convergence
1128:   test:
1129:     suffix: q1_2d_plex_6
1130:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 1 -non_affine_coords -convergence
1131:   test:
1132:     suffix: q1_2d_plex_7
1133:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 2 -non_affine_coords -convergence

1135:   # 2D Q_2 on a quadrilaterial
1136:   test:
1137:     suffix: q2_2d_plex_0
1138:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1139:   test:
1140:     suffix: q2_2d_plex_1
1141:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1142:   test:
1143:     suffix: q2_2d_plex_2
1144:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1145:   test:
1146:     suffix: q2_2d_plex_3
1147:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1148:   test:
1149:     suffix: q2_2d_plex_4
1150:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1151:   test:
1152:     suffix: q2_2d_plex_5
1153:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 0 -non_affine_coords -convergence
1154:   test:
1155:     suffix: q2_2d_plex_6
1156:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 1 -non_affine_coords -convergence
1157:   test:
1158:     suffix: q2_2d_plex_7
1159:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 2 -non_affine_coords -convergence

1161:   # 2D P_3 on a triangle
1162:   test:
1163:     suffix: p3_2d_0
1164:     requires: triangle !single
1165:     args: -petscspace_degree 3 -qorder 3 -convergence
1166:   test:
1167:     suffix: p3_2d_1
1168:     requires: triangle !single
1169:     args: -petscspace_degree 3 -qorder 3 -porder 1
1170:   test:
1171:     suffix: p3_2d_2
1172:     requires: triangle !single
1173:     args: -petscspace_degree 3 -qorder 3 -porder 2
1174:   test:
1175:     suffix: p3_2d_3
1176:     requires: triangle !single
1177:     args: -petscspace_degree 3 -qorder 3 -porder 3
1178:   test:
1179:     suffix: p3_2d_4
1180:     requires: triangle mmg
1181:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -convergence -conv_refine 0
1182:   test:
1183:     suffix: p3_2d_5
1184:     requires: triangle mmg
1185:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 1 -conv_refine 0
1186:   test:
1187:     suffix: p3_2d_6
1188:     requires: triangle mmg
1189:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 3 -conv_refine 0

1191:   # 2D Q_3 on a quadrilaterial
1192:   test:
1193:     suffix: q3_2d_0
1194:     requires: !single
1195:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -convergence
1196:   test:
1197:     suffix: q3_2d_1
1198:     requires: !single
1199:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 1
1200:   test:
1201:     suffix: q3_2d_2
1202:     requires: !single
1203:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 2
1204:   test:
1205:     suffix: q3_2d_3
1206:     requires: !single
1207:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 3

1209:   # 2D P_1disc on a triangle/quadrilateral
1210:   test:
1211:     suffix: p1d_2d_0
1212:     requires: triangle
1213:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1214:   test:
1215:     suffix: p1d_2d_1
1216:     requires: triangle
1217:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1218:   test:
1219:     suffix: p1d_2d_2
1220:     requires: triangle
1221:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1222:   test:
1223:     suffix: p1d_2d_3
1224:     requires: triangle
1225:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1226:     filter: sed  -e "s/convergence rate at refinement 0: 2/convergence rate at refinement 0: 1.9/g"
1227:   test:
1228:     suffix: p1d_2d_4
1229:     requires: triangle
1230:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1231:   test:
1232:     suffix: p1d_2d_5
1233:     requires: triangle
1234:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2

1236:   # 2D BDM_1 on a triangle
1237:   test:
1238:     suffix: bdm1_2d_0
1239:     requires: triangle
1240:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1241:           -num_comp 2 -qorder 1 -convergence
1242:   test:
1243:     suffix: bdm1_2d_1
1244:     requires: triangle
1245:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1246:           -num_comp 2 -qorder 1 -porder 1
1247:   test:
1248:     suffix: bdm1_2d_2
1249:     requires: triangle
1250:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1251:           -num_comp 2 -qorder 1 -porder 2

1253:   # 2D BDM_1 on a quadrilateral
1254:   test:
1255:     suffix: bdm1q_2d_0
1256:     requires: triangle
1257:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1258:           -petscdualspace_lagrange_tensor 1 \
1259:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -convergence
1260:   test:
1261:     suffix: bdm1q_2d_1
1262:     requires: triangle
1263:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1264:           -petscdualspace_lagrange_tensor 1 \
1265:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 1
1266:   test:
1267:     suffix: bdm1q_2d_2
1268:     requires: triangle
1269:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1270:           -petscdualspace_lagrange_tensor 1 \
1271:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 2

1273:   # Test high order quadrature
1274:   test:
1275:     suffix: p1_quad_2
1276:     requires: triangle
1277:     args: -petscspace_degree 1 -qorder 2 -porder 1
1278:   test:
1279:     suffix: p1_quad_5
1280:     requires: triangle
1281:     args: -petscspace_degree 1 -qorder 5 -porder 1
1282:   test:
1283:     suffix: p2_quad_3
1284:     requires: triangle
1285:     args: -petscspace_degree 2 -qorder 3 -porder 2
1286:   test:
1287:     suffix: p2_quad_5
1288:     requires: triangle
1289:     args: -petscspace_degree 2 -qorder 5 -porder 2
1290:   test:
1291:     suffix: q1_quad_2
1292:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 2 -porder 1
1293:   test:
1294:     suffix: q1_quad_5
1295:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 5 -porder 1
1296:   test:
1297:     suffix: q2_quad_3
1298:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 3 -porder 1
1299:   test:
1300:     suffix: q2_quad_5
1301:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 5 -porder 1

1303:   # Nonconforming tests
1304:   test:
1305:     suffix: constraints
1306:     args: -dm_coord_space 0 -dm_plex_simplex 0 -petscspace_type tensor -petscspace_degree 1 -qorder 0 -constraints
1307:   test:
1308:     suffix: nonconforming_tensor_2
1309:     nsize: 4
1310:     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
1311:   test:
1312:     suffix: nonconforming_tensor_3
1313:     nsize: 4
1314:     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
1315:   test:
1316:     suffix: nonconforming_tensor_2_fv
1317:     nsize: 4
1318:     args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_simplex 0 -num_comp 2
1319:   test:
1320:     suffix: nonconforming_tensor_3_fv
1321:     nsize: 4
1322:     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
1323:   test:
1324:     suffix: nonconforming_tensor_2_hi
1325:     requires: !single
1326:     nsize: 4
1327:     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
1328:   test:
1329:     suffix: nonconforming_tensor_3_hi
1330:     requires: !single skip
1331:     nsize: 4
1332:     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
1333:   test:
1334:     suffix: nonconforming_simplex_2
1335:     requires: triangle
1336:     nsize: 4
1337:     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
1338:   test:
1339:     suffix: nonconforming_simplex_2_hi
1340:     requires: triangle !single
1341:     nsize: 4
1342:     args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_max_projection_height 1 -petscspace_degree 4 -qorder 4
1343:   test:
1344:     suffix: nonconforming_simplex_2_fv
1345:     requires: triangle
1346:     nsize: 4
1347:     args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -num_comp 2
1348:   test:
1349:     suffix: nonconforming_simplex_3
1350:     requires: ctetgen
1351:     nsize: 4
1352:     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
1353:   test:
1354:     suffix: nonconforming_simplex_3_hi
1355:     requires: ctetgen skip
1356:     nsize: 4
1357:     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
1358:   test:
1359:     suffix: nonconforming_simplex_3_fv
1360:     requires: ctetgen
1361:     nsize: 4
1362:     args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_dim 3 -num_comp 3

1364:   # 3D WXY on a triangular prism
1365:   test:
1366:     suffix: wxy_0
1367:     args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism -qorder 2 -porder 0 \
1368:           -petscspace_type sum \
1369:           -petscspace_variables 3 \
1370:           -petscspace_components 3 \
1371:           -petscspace_sum_spaces 2 \
1372:           -petscspace_sum_concatenate false \
1373:           -sumcomp_0_petscspace_variables 3 \
1374:           -sumcomp_0_petscspace_components 3 \
1375:           -sumcomp_0_petscspace_degree 1 \
1376:           -sumcomp_1_petscspace_variables 3 \
1377:           -sumcomp_1_petscspace_components 3 \
1378:           -sumcomp_1_petscspace_type wxy \
1379:           -petscdualspace_refcell triangular_prism \
1380:           -petscdualspace_form_degree 0 \
1381:           -petscdualspace_order 1 \
1382:           -petscdualspace_components 3

1384:   # 2D RT_0 on a triangle
1385:   test:
1386:     suffix: rt0_2d_tri
1387:     requires: triangle
1388:     args: -qorder 1 -porder 0 -RT \
1389:           -petscspace_type ptrimmed \
1390:           -petscspace_components 2 \
1391:           -petscspace_ptrimmed_form_degree -1 \
1392:           -petscdualspace_order 1 \
1393:           -petscdualspace_form_degree -1 \
1394:           -petscdualspace_lagrange_trimmed true

1396:   # 2D RT_0 on a quadrilateral
1397:   test:
1398:     suffix: rt0_2d_quad
1399:     requires: triangle
1400:     args: -dm_plex_simplex 0 -qorder 1 -porder 0 -RT \
1401:           -petscspace_degree 1 \
1402:           -petscspace_type sum \
1403:           -petscspace_variables 2 \
1404:           -petscspace_components 2 \
1405:           -petscspace_sum_spaces 2 \
1406:           -petscspace_sum_concatenate true \
1407:           -sumcomp_0_petscspace_variables 2 \
1408:           -sumcomp_0_petscspace_type tensor \
1409:           -sumcomp_0_petscspace_tensor_spaces 2 \
1410:           -sumcomp_0_petscspace_tensor_uniform false \
1411:           -sumcomp_0_tensorcomp_0_petscspace_degree 1 \
1412:           -sumcomp_0_tensorcomp_1_petscspace_degree 0 \
1413:           -sumcomp_1_petscspace_variables 2 \
1414:           -sumcomp_1_petscspace_type tensor \
1415:           -sumcomp_1_petscspace_tensor_spaces 2 \
1416:           -sumcomp_1_petscspace_tensor_uniform false \
1417:           -sumcomp_1_tensorcomp_0_petscspace_degree 0 \
1418:           -sumcomp_1_tensorcomp_1_petscspace_degree 1 \
1419:           -petscdualspace_form_degree -1 \
1420:           -petscdualspace_order 1 \
1421:           -petscdualspace_lagrange_trimmed true

1423: TEST*/

1425: /*
1426:    # 2D Q_2 on a quadrilaterial Plex
1427:   test:
1428:     suffix: q2_2d_plex_0
1429:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1430:   test:
1431:     suffix: q2_2d_plex_1
1432:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1433:   test:
1434:     suffix: q2_2d_plex_2
1435:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1436:   test:
1437:     suffix: q2_2d_plex_3
1438:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1439:   test:
1440:     suffix: q2_2d_plex_4
1441:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1442:   test:
1443:     suffix: q2_2d_plex_5
1444:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 0 -non_affine_coords
1445:   test:
1446:     suffix: q2_2d_plex_6
1447:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 1 -non_affine_coords
1448:   test:
1449:     suffix: q2_2d_plex_7
1450:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 2 -non_affine_coords

1452:   test:
1453:     suffix: p1d_2d_6
1454:     requires: mmg
1455:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1456:   test:
1457:     suffix: p1d_2d_7
1458:     requires: mmg
1459:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1460:   test:
1461:     suffix: p1d_2d_8
1462:     requires: mmg
1463:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1464: */