Actual source code: ex3.c

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

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

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

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

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

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

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

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

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

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

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

180:   return 0;
181: }

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

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

201:       PetscSectionGetDof(coordSection, v, &dof);
202:       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:     VecRestoreArray(coordinates, &coords);
218:   }
219:   if (user->shearCoords) {
220:     /* x' = x + m y + m z, y' = y + m z,  z' = z */
221:     DMGetCoordinateSection(dm, &coordSection);
222:     DMGetCoordinatesLocal(dm, &coordinates);
223:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
224:     VecGetArray(coordinates, &coords);
225:     for (v = vStart; v < vEnd; ++v) {
226:       PetscInt  dof, off;
227:       PetscReal m = 1.0;

229:       PetscSectionGetDof(coordSection, v, &dof);
230:       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:     VecRestoreArray(coordinates, &coords);
244:   }
245:   return 0;
246: }

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

254:   if (user->useDA) {
255:     switch (dim) {
256:     case 2:
257:       DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 2, 2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, dm);
258:       DMSetFromOptions(*dm);
259:       DMSetUp(*dm);
260:       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:     PetscObjectSetName((PetscObject)*dm, "Hexahedral Mesh");
266:   } else {
267:     DMCreate(comm, dm);
268:     DMSetType(*dm, DMPLEX);
269:     DMPlexDistributeSetDefault(*dm, PETSC_FALSE);
270:     DMSetFromOptions(*dm);

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

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

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:     DMGetCoordinateDM(dm, &coordDM);

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

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

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

401:     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:       DMGetLocalSection(dm, &section);
411:       /* this creates the matrix and preallocates the matrix structure: we
412:        * just have to fill in the values */
413:       DMGetDefaultConstraints(dm, &cSec, &cMat, NULL);
414:       PetscSectionGetChart(cSec, &cStart, &cEnd);
415:       ISGetIndices(aIS, &anchors);
416:       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:         PetscSectionGetDof(aSec, c, &cDof);
422:         if (cDof) {
423:           PetscInt cOff, a, aDof, aOff, j;

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

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

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


441:           /* put in a simple equality constraint */
442:           for (j = 0; j < cDof; j++) MatSetValue(cMat, cOff + j, aOff + j, 1., INSERT_VALUES);
443:         }
444:       }
445:       MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY);
446:       MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY);
447:       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:       DMCreateMatrix(dm, &mass);
453:       /* get a dummy local variable to serve as the solution */
454:       DMGetLocalVector(dm, &local);
455:       DMGetDS(dm, &ds);
456:       /* set the jacobian to be the mass matrix */
457:       PetscDSSetJacobian(ds, 0, 0, simple_mass, NULL, NULL, NULL);
458:       /* build the mass matrix */
459:       DMPlexSNESComputeJacobianFEM(dm, local, mass, mass, NULL);
460:       MatView(mass, PETSC_VIEWER_STDOUT_WORLD);
461:       MatDestroy(&mass);
462:       DMRestoreLocalVector(dm, &local);
463:     }
464:   }
465:   return 0;
466: }

468: static PetscErrorCode TestFEJacobian(DM dm, AppCtx *user)
469: {
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:     DMGetDimension(dm, &dim);
483:     DMGetDS(dm, &ds);
484:     PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, symmetric_gradient_inner_product);
485:     DMCreateMatrix(dm, &E);
486:     DMGetLocalVector(dm, &local);
487:     DMPlexSNESComputeJacobianFEM(dm, local, E, E, NULL);
488:     DMPlexCreateRigidBody(dm, 0, &sp);
489:     MatNullSpaceGetVecs(sp, &hasConst, &n, &vecs);
490:     if (n) VecDuplicate(vecs[0], &res);
491:     DMCreateLocalVector(dm, &localX);
492:     DMCreateLocalVector(dm, &localRes);
493:     for (i = 0; i < n; i++) { /* also test via matrix-free Jacobian application */
494:       PetscReal resNorm;

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

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

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

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

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

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

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

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

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

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

630:       DMPlexPointLocalRead(dmgrad, c, gradArray, &compGrad);

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

659: 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)
660: {
661:   Vec       u;
662:   PetscReal n[3] = {1.0, 1.0, 1.0};

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

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

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

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

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

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

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

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

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

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

909: /*TEST

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1308: TEST*/

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

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