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)      {u[0] = coords[0]*coords[1]; u[1] = coords[1]*coords[2]; u[2] = coords[2]*coords[0];}
 69:   else if (dim > 1) {u[0] = coords[0]*coords[0]; u[1] = coords[0]*coords[1];}
 70:   else if (dim > 0) {u[0] = coords[0]*coords[0];}
 71:   return 0;
 72: }
 73: PetscErrorCode quadraticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 74: {
 75:   if (dim > 2)      {u[0] = coords[1]*n[0] + coords[0]*n[1]; u[1] = coords[2]*n[1] + coords[1]*n[2]; u[2] = coords[2]*n[0] + coords[0]*n[2];}
 76:   else if (dim > 1) {u[0] = 2.0*coords[0]*n[0]; u[1] = coords[1]*n[0] + coords[0]*n[1];}
 77:   else if (dim > 0) {u[0] = 2.0*coords[0]*n[0];}
 78:   return 0;
 79: }

 81: /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
 82: PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 83: {
 84:   if (dim > 2)      {u[0] = coords[0]*coords[0]*coords[1]; u[1] = coords[1]*coords[1]*coords[2]; u[2] = coords[2]*coords[2]*coords[0];}
 85:   else if (dim > 1) {u[0] = coords[0]*coords[0]*coords[0]; u[1] = coords[0]*coords[0]*coords[1];}
 86:   else if (dim > 0) {u[0] = coords[0]*coords[0]*coords[0];}
 87:   return 0;
 88: }
 89: PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 90: {
 91:   if (dim > 2)      {u[0] = 2.0*coords[0]*coords[1]*n[0] + coords[0]*coords[0]*n[1]; u[1] = 2.0*coords[1]*coords[2]*n[1] + coords[1]*coords[1]*n[2]; u[2] = 2.0*coords[2]*coords[0]*n[2] + coords[2]*coords[2]*n[0];}
 92:   else if (dim > 1) {u[0] = 3.0*coords[0]*coords[0]*n[0]; u[1] = 2.0*coords[0]*coords[1]*n[0] + coords[0]*coords[0]*n[1];}
 93:   else if (dim > 0) {u[0] = 3.0*coords[0]*coords[0]*n[0];}
 94:   return 0;
 95: }

 97: /* u = tanh(x) */
 98: PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 99: {
100:   PetscInt d;
101:   for (d = 0; d < dim; ++d) u[d] = PetscTanhReal(coords[d] - 0.5);
102:   return 0;
103: }
104: PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
105: {
106:   PetscInt d;
107:   for (d = 0; d < dim; ++d) u[d] = 1.0/PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d];
108:   return 0;
109: }

111: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
112: {
113:   PetscInt       n = 3;

117:   options->useDA           = PETSC_FALSE;
118:   options->shearCoords     = PETSC_FALSE;
119:   options->nonaffineCoords = PETSC_FALSE;
120:   options->qorder          = 0;
121:   options->numComponents   = PETSC_DEFAULT;
122:   options->porder          = 0;
123:   options->convergence     = PETSC_FALSE;
124:   options->convRefine      = PETSC_TRUE;
125:   options->constraints     = PETSC_FALSE;
126:   options->tree            = PETSC_FALSE;
127:   options->treeCell        = 0;
128:   options->testFEjacobian  = PETSC_FALSE;
129:   options->testFVgrad      = PETSC_FALSE;
130:   options->testInjector    = PETSC_FALSE;
131:   options->constants[0]    = 1.0;
132:   options->constants[1]    = 2.0;
133:   options->constants[2]    = 3.0;

135:   PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
136:   PetscOptionsBool("-use_da", "Flag for DMDA mesh", "ex3.c", options->useDA, &options->useDA, NULL);
137:   PetscOptionsBool("-shear_coords", "Transform coordinates with a shear", "ex3.c", options->shearCoords, &options->shearCoords, NULL);
138:   PetscOptionsBool("-non_affine_coords", "Transform coordinates with a non-affine transform", "ex3.c", options->nonaffineCoords, &options->nonaffineCoords, NULL);
139:   PetscOptionsBoundedInt("-qorder", "The quadrature order", "ex3.c", options->qorder, &options->qorder, NULL,0);
140:   PetscOptionsBoundedInt("-num_comp", "The number of field components", "ex3.c", options->numComponents, &options->numComponents, NULL,PETSC_DEFAULT);
141:   PetscOptionsBoundedInt("-porder", "The order of polynomials to test", "ex3.c", options->porder, &options->porder, NULL,0);
142:   PetscOptionsBool("-convergence", "Check the convergence rate", "ex3.c", options->convergence, &options->convergence, NULL);
143:   PetscOptionsBool("-conv_refine", "Use refinement for the convergence rate", "ex3.c", options->convRefine, &options->convRefine, NULL);
144:   PetscOptionsBool("-constraints", "Test local constraints (serial only)", "ex3.c", options->constraints, &options->constraints, NULL);
145:   PetscOptionsBool("-tree", "Test tree routines", "ex3.c", options->tree, &options->tree, NULL);
146:   PetscOptionsBoundedInt("-tree_cell", "cell to refine in tree test", "ex3.c", options->treeCell, &options->treeCell, NULL,0);
147:   PetscOptionsBool("-test_fe_jacobian", "Test finite element Jacobian assembly", "ex3.c", options->testFEjacobian, &options->testFEjacobian, NULL);
148:   PetscOptionsBool("-test_fv_grad", "Test finite volume gradient reconstruction", "ex3.c", options->testFVgrad, &options->testFVgrad, NULL);
149:   PetscOptionsBool("-test_injector","Test finite element injection", "ex3.c", options->testInjector, &options->testInjector,NULL);
150:   PetscOptionsRealArray("-constants","Set the constant values", "ex3.c", options->constants, &n,NULL);
151:   PetscOptionsEnd();

153:   return 0;
154: }

156: static PetscErrorCode TransformCoordinates(DM dm, AppCtx *user)
157: {
158:   PetscSection   coordSection;
159:   Vec            coordinates;
160:   PetscScalar   *coords;
161:   PetscInt       vStart, vEnd, v;

164:   if (user->nonaffineCoords) {
165:     /* x' = r^(1/p) (x/r), y' = r^(1/p) (y/r), z' = z */
166:     DMGetCoordinateSection(dm, &coordSection);
167:     DMGetCoordinatesLocal(dm, &coordinates);
168:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
169:     VecGetArray(coordinates, &coords);
170:     for (v = vStart; v < vEnd; ++v) {
171:       PetscInt  dof, off;
172:       PetscReal p = 4.0, r;

174:       PetscSectionGetDof(coordSection, v, &dof);
175:       PetscSectionGetOffset(coordSection, v, &off);
176:       switch (dof) {
177:       case 2:
178:         r             = PetscSqr(PetscRealPart(coords[off+0])) + PetscSqr(PetscRealPart(coords[off+1]));
179:         coords[off+0] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p)/(2*p))*coords[off+0];
180:         coords[off+1] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p)/(2*p))*coords[off+1];
181:         break;
182:       case 3:
183:         r             = PetscSqr(PetscRealPart(coords[off+0])) + PetscSqr(PetscRealPart(coords[off+1]));
184:         coords[off+0] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p)/(2*p))*coords[off+0];
185:         coords[off+1] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p)/(2*p))*coords[off+1];
186:         coords[off+2] = coords[off+2];
187:         break;
188:       }
189:     }
190:     VecRestoreArray(coordinates, &coords);
191:   }
192:   if (user->shearCoords) {
193:     /* x' = x + m y + m z, y' = y + m z,  z' = z */
194:     DMGetCoordinateSection(dm, &coordSection);
195:     DMGetCoordinatesLocal(dm, &coordinates);
196:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
197:     VecGetArray(coordinates, &coords);
198:     for (v = vStart; v < vEnd; ++v) {
199:       PetscInt  dof, off;
200:       PetscReal m = 1.0;

202:       PetscSectionGetDof(coordSection, v, &dof);
203:       PetscSectionGetOffset(coordSection, v, &off);
204:       switch (dof) {
205:       case 2:
206:         coords[off+0] = coords[off+0] + m*coords[off+1];
207:         coords[off+1] = coords[off+1];
208:         break;
209:       case 3:
210:         coords[off+0] = coords[off+0] + m*coords[off+1] + m*coords[off+2];
211:         coords[off+1] = coords[off+1] + m*coords[off+2];
212:         coords[off+2] = coords[off+2];
213:         break;
214:       }
215:     }
216:     VecRestoreArray(coordinates, &coords);
217:   }
218:   return 0;
219: }

221: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
222: {
223:   PetscInt       dim = 2;
224:   PetscBool      simplex;

227:   if (user->useDA) {
228:     switch (dim) {
229:     case 2:
230:       DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 2, 2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, dm);
231:       DMSetFromOptions(*dm);
232:       DMSetUp(*dm);
233:       DMDASetVertexCoordinates(*dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
234:       break;
235:     default:
236:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot create structured mesh of dimension %d", dim);
237:     }
238:     PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");
239:   } else {
240:     DMCreate(comm, dm);
241:     DMSetType(*dm, DMPLEX);
242:     DMPlexDistributeSetDefault(*dm, PETSC_FALSE);
243:     DMSetFromOptions(*dm);

245:     DMGetDimension(*dm, &dim);
246:     DMPlexIsSimplex(*dm, &simplex);
247:     MPI_Bcast(&simplex, 1, MPIU_BOOL, 0, comm);
248:     if (user->tree) {
249:       DM refTree, ncdm = NULL;

251:       DMPlexCreateDefaultReferenceTree(comm,dim,simplex,&refTree);
252:       DMViewFromOptions(refTree,NULL,"-reftree_dm_view");
253:       DMPlexSetReferenceTree(*dm,refTree);
254:       DMDestroy(&refTree);
255:       DMPlexTreeRefineCell(*dm,user->treeCell,&ncdm);
256:       if (ncdm) {
257:         DMDestroy(dm);
258:         *dm = ncdm;
259:         DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
260:       }
261:       PetscObjectSetOptionsPrefix((PetscObject) *dm, "tree_");
262:       DMPlexDistributeSetDefault(*dm, PETSC_FALSE);
263:       DMSetFromOptions(*dm);
264:       DMViewFromOptions(*dm,NULL,"-dm_view");
265:     } else {
266:       DMPlexSetRefinementUniform(*dm, PETSC_TRUE);
267:     }
268:     PetscObjectSetOptionsPrefix((PetscObject) *dm, "dist_");
269:     DMPlexDistributeSetDefault(*dm, PETSC_FALSE);
270:     DMSetFromOptions(*dm);
271:     PetscObjectSetOptionsPrefix((PetscObject) *dm, NULL);
272:     if (simplex) PetscObjectSetName((PetscObject) *dm, "Simplicial Mesh");
273:     else         PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");
274:   }
275:   DMSetFromOptions(*dm);
276:   TransformCoordinates(*dm, user);
277:   DMViewFromOptions(*dm,NULL,"-dm_view");
278:   return 0;
279: }

281: static void simple_mass(PetscInt dim, PetscInt Nf, PetscInt NfAux,
282:                         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
283:                         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
284:                         PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
285: {
286:   PetscInt d, e;
287:   for (d = 0, e = 0; d < dim; d++, e+=dim+1) {
288:     g0[e] = 1.;
289:   }
290: }

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

300:   for (compI = 0; compI < dim; ++compI) {
301:     for (compJ = 0; compJ < dim; ++compJ) {
302:       for (d = 0; d < dim; ++d) {
303:         for (e = 0; e < dim; e++) {
304:           if (d == e && d == compI && d == compJ) {
305:             C[((compI*dim+compJ)*dim+d)*dim+e] = 1.0;
306:           } else if ((d == compJ && e == compI) || (d == e && compI == compJ)) {
307:             C[((compI*dim+compJ)*dim+d)*dim+e] = 0.5;
308:           } else {
309:             C[((compI*dim+compJ)*dim+d)*dim+e] = 0.0;
310:           }
311:         }
312:       }
313:     }
314:   }
315: }

317: static PetscErrorCode SetupSection(DM dm, AppCtx *user)
318: {
320:   if (user->constraints) {
321:     /* test local constraints */
322:     DM            coordDM;
323:     PetscInt      fStart, fEnd, f, vStart, vEnd, v;
324:     PetscInt      edgesx = 2, vertsx;
325:     PetscInt      edgesy = 2, vertsy;
326:     PetscMPIInt   size;
327:     PetscInt      numConst;
328:     PetscSection  aSec;
329:     PetscInt     *anchors;
330:     PetscInt      offset;
331:     IS            aIS;
332:     MPI_Comm      comm = PetscObjectComm((PetscObject)dm);

334:     MPI_Comm_size(comm,&size);

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

341:     /* first create the coordinate section so that it does not clone the
342:      * constraints */
343:     DMGetCoordinateDM(dm,&coordDM);

345:     /* create the constrained-to-anchor section */
346:     DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
347:     DMPlexGetDepthStratum(dm,1,&fStart,&fEnd);
348:     PetscSectionCreate(PETSC_COMM_SELF,&aSec);
349:     PetscSectionSetChart(aSec,PetscMin(fStart,vStart),PetscMax(fEnd,vEnd));

351:     /* define the constraints */
352:     PetscOptionsGetInt(NULL,NULL,"-da_grid_x",&edgesx,NULL);
353:     PetscOptionsGetInt(NULL,NULL,"-da_grid_y",&edgesy,NULL);
354:     vertsx = edgesx + 1;
355:     vertsy = edgesy + 1;
356:     numConst = vertsy + edgesy;
357:     PetscMalloc1(numConst,&anchors);
358:     offset = 0;
359:     for (v = vStart + edgesx; v < vEnd; v+= vertsx) {
360:       PetscSectionSetDof(aSec,v,1);
361:       anchors[offset++] = v - edgesx;
362:     }
363:     for (f = fStart + edgesx * vertsy + edgesx * edgesy; f < fEnd; f++) {
364:       PetscSectionSetDof(aSec,f,1);
365:       anchors[offset++] = f - edgesx * edgesy;
366:     }
367:     PetscSectionSetUp(aSec);
368:     ISCreateGeneral(PETSC_COMM_SELF,numConst,anchors,PETSC_OWN_POINTER,&aIS);

370:     DMPlexSetAnchors(dm,aSec,aIS);
371:     PetscSectionDestroy(&aSec);
372:     ISDestroy(&aIS);
373:   }
374:   DMSetNumFields(dm, 1);
375:   DMSetField(dm, 0, NULL, (PetscObject) user->fe);
376:   DMCreateDS(dm);
377:   if (user->constraints) {
378:     /* test getting local constraint matrix that matches section */
379:     PetscSection aSec;
380:     IS           aIS;

382:     DMPlexGetAnchors(dm,&aSec,&aIS);
383:     if (aSec) {
384:       PetscDS         ds;
385:       PetscSection    cSec, section;
386:       PetscInt        cStart, cEnd, c, numComp;
387:       Mat             cMat, mass;
388:       Vec             local;
389:       const PetscInt *anchors;

391:       DMGetLocalSection(dm,&section);
392:       /* this creates the matrix and preallocates the matrix structure: we
393:        * just have to fill in the values */
394:       DMGetDefaultConstraints(dm,&cSec,&cMat,NULL);
395:       PetscSectionGetChart(cSec,&cStart,&cEnd);
396:       ISGetIndices(aIS,&anchors);
397:       PetscFEGetNumComponents(user->fe, &numComp);
398:       for (c = cStart; c < cEnd; c++) {
399:         PetscInt cDof;

401:         /* is this point constrained? (does it have an anchor?) */
402:         PetscSectionGetDof(aSec,c,&cDof);
403:         if (cDof) {
404:           PetscInt cOff, a, aDof, aOff, j;

407:           /* find the anchor point */
408:           PetscSectionGetOffset(aSec,c,&cOff);
409:           a    = anchors[cOff];

411:           /* find the constrained dofs (row in constraint matrix) */
412:           PetscSectionGetDof(cSec,c,&cDof);
413:           PetscSectionGetOffset(cSec,c,&cOff);

415:           /* find the anchor dofs (column in constraint matrix) */
416:           PetscSectionGetDof(section,a,&aDof);
417:           PetscSectionGetOffset(section,a,&aOff);


422:           /* put in a simple equality constraint */
423:           for (j = 0; j < cDof; j++) {
424:             MatSetValue(cMat,cOff+j,aOff+j,1.,INSERT_VALUES);
425:           }
426:         }
427:       }
428:       MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
429:       MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
430:       ISRestoreIndices(aIS,&anchors);

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

435:       DMCreateMatrix(dm,&mass);
436:       /* get a dummy local variable to serve as the solution */
437:       DMGetLocalVector(dm,&local);
438:       DMGetDS(dm,&ds);
439:       /* set the jacobian to be the mass matrix */
440:       PetscDSSetJacobian(ds, 0, 0, simple_mass, NULL,  NULL, NULL);
441:       /* build the mass matrix */
442:       DMPlexSNESComputeJacobianFEM(dm,local,mass,mass,NULL);
443:       MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
444:       MatDestroy(&mass);
445:       DMRestoreLocalVector(dm,&local);
446:     }
447:   }
448:   return 0;
449: }

451: static PetscErrorCode TestFEJacobian(DM dm, AppCtx *user)
452: {
454:   if (!user->useDA) {
455:     Vec          local;
456:     const Vec    *vecs;
457:     Mat          E;
458:     MatNullSpace sp;
459:     PetscBool    isNullSpace, hasConst;
460:     PetscInt     dim, n, i;
461:     Vec          res = NULL, localX, localRes;
462:     PetscDS      ds;

464:     DMGetDimension(dm, &dim);
466:     DMGetDS(dm,&ds);
467:     PetscDSSetJacobian(ds,0,0,NULL,NULL,NULL,symmetric_gradient_inner_product);
468:     DMCreateMatrix(dm,&E);
469:     DMGetLocalVector(dm,&local);
470:     DMPlexSNESComputeJacobianFEM(dm,local,E,E,NULL);
471:     DMPlexCreateRigidBody(dm,0,&sp);
472:     MatNullSpaceGetVecs(sp,&hasConst,&n,&vecs);
473:     if (n) VecDuplicate(vecs[0],&res);
474:     DMCreateLocalVector(dm,&localX);
475:     DMCreateLocalVector(dm,&localRes);
476:     for (i = 0; i < n; i++) { /* also test via matrix-free Jacobian application */
477:       PetscReal resNorm;

479:       VecSet(localRes,0.);
480:       VecSet(localX,0.);
481:       VecSet(local,0.);
482:       VecSet(res,0.);
483:       DMGlobalToLocalBegin(dm,vecs[i],INSERT_VALUES,localX);
484:       DMGlobalToLocalEnd(dm,vecs[i],INSERT_VALUES,localX);
485:       DMSNESComputeJacobianAction(dm,local,localX,localRes,NULL);
486:       DMLocalToGlobalBegin(dm,localRes,ADD_VALUES,res);
487:       DMLocalToGlobalEnd(dm,localRes,ADD_VALUES,res);
488:       VecNorm(res,NORM_2,&resNorm);
489:       if (resNorm > PETSC_SMALL) {
490:         PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient action null space vector %D residual: %E\n",i,resNorm);
491:       }
492:     }
493:     VecDestroy(&localRes);
494:     VecDestroy(&localX);
495:     VecDestroy(&res);
496:     MatNullSpaceTest(sp,E,&isNullSpace);
497:     if (isNullSpace) {
498:       PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient null space: PASS\n");
499:     } else {
500:       PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient null space: FAIL\n");
501:     }
502:     MatNullSpaceDestroy(&sp);
503:     MatDestroy(&E);
504:     DMRestoreLocalVector(dm,&local);
505:   }
506:   return 0;
507: }

509: static PetscErrorCode TestInjector(DM dm, AppCtx *user)
510: {
511:   DM             refTree;
512:   PetscMPIInt    rank;

514:   DMPlexGetReferenceTree(dm,&refTree);
515:   if (refTree) {
516:     Mat inj;

518:     DMPlexComputeInjectorReferenceTree(refTree,&inj);
519:     PetscObjectSetName((PetscObject)inj,"Reference Tree Injector");
520:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
521:     if (rank == 0) {
522:       MatView(inj,PETSC_VIEWER_STDOUT_SELF);
523:     }
524:     MatDestroy(&inj);
525:   }
526:   return 0;
527: }

529: static PetscErrorCode TestFVGrad(DM dm, AppCtx *user)
530: {
531:   MPI_Comm          comm;
532:   DM                dmRedist, dmfv, dmgrad, dmCell, refTree;
533:   PetscFV           fv;
534:   PetscInt          dim, nvecs, v, cStart, cEnd, cEndInterior;
535:   PetscMPIInt       size;
536:   Vec               cellgeom, grad, locGrad;
537:   const PetscScalar *cgeom;
538:   PetscReal         allVecMaxDiff = 0., fvTol = 100. * PETSC_MACHINE_EPSILON;

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

584:     DMGetLocalVector(dmfv,&locX);
585:     /* get the local projection of the rigid body mode */
586:     for (c = cStart; c < cEnd; c++) {
587:       PetscFVCellGeom *cg;
588:       PetscScalar     cx[3] = {0.,0.,0.};

590:       DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
591:       if (v < dim) {
592:         cx[v] = 1.;
593:       } else {
594:         PetscInt w = v - dim;

596:         cx[(w + 1) % dim] =  cg->centroid[(w + 2) % dim];
597:         cx[(w + 2) % dim] = -cg->centroid[(w + 1) % dim];
598:       }
599:       DMPlexVecSetClosure(dmfv,NULL,locX,c,cx,INSERT_ALL_VALUES);
600:     }
601:     /* TODO: this isn't in any header */
602:     DMPlexReconstructGradientsFVM(dmfv,locX,grad);
603:     DMGlobalToLocalBegin(dmgrad,grad,INSERT_VALUES,locGrad);
604:     DMGlobalToLocalEnd(dmgrad,grad,INSERT_VALUES,locGrad);
605:     VecGetArrayRead(locGrad,&gradArray);
606:     /* compare computed gradient to exact gradient */
607:     if (v >= dim) {
608:       PetscInt w = v - dim;

610:       trueGrad[(w + 1) % dim][(w + 2) % dim] =  1.;
611:       trueGrad[(w + 2) % dim][(w + 1) % dim] = -1.;
612:     }
613:     maxDiff = 0.;
614:     for (c = cStart; c < cEndInterior; c++) {
615:       PetscScalar *compGrad;
616:       PetscInt    i, j, k;
617:       PetscReal   FrobDiff = 0.;

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

621:       for (i = 0, k = 0; i < dim; i++) {
622:         for (j = 0; j < dim; j++, k++) {
623:           PetscScalar diff = compGrad[k] - trueGrad[i][j];
624:           FrobDiff += PetscRealPart(diff * PetscConj(diff));
625:         }
626:       }
627:       FrobDiff = PetscSqrtReal(FrobDiff);
628:       maxDiff  = PetscMax(maxDiff,FrobDiff);
629:     }
630:     MPI_Allreduce(&maxDiff,&maxDiffGlob,1,MPIU_REAL,MPIU_MAX,comm);
631:     allVecMaxDiff = PetscMax(allVecMaxDiff,maxDiffGlob);
632:     VecRestoreArrayRead(locGrad,&gradArray);
633:     DMRestoreLocalVector(dmfv,&locX);
634:   }
635:   if (allVecMaxDiff < fvTol) {
636:     PetscPrintf(PetscObjectComm((PetscObject)dm),"Finite volume gradient reconstruction: PASS\n");
637:   } else {
638:     PetscPrintf(PetscObjectComm((PetscObject)dm),"Finite volume gradient reconstruction: FAIL at tolerance %g with max difference %g\n",fvTol,allVecMaxDiff);
639:   }
640:   DMRestoreLocalVector(dmgrad,&locGrad);
641:   DMRestoreGlobalVector(dmgrad,&grad);
642:   VecRestoreArrayRead(cellgeom,&cgeom);
643:   DMDestroy(&dmfv);
644:   PetscFVDestroy(&fv);
645:   return 0;
646: }

648: static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *),
649:                                    PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *),
650:                                    void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
651: {
652:   Vec            u;
653:   PetscReal      n[3] = {1.0, 1.0, 1.0};

656:   DMGetGlobalVector(dm, &u);
657:   /* Project function into FE function space */
658:   DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u);
659:   VecViewFromOptions(u, NULL, "-projection_view");
660:   /* Compare approximation to exact in L_2 */
661:   DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error);
662:   DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer);
663:   DMRestoreGlobalVector(dm, &u);
664:   return 0;
665: }

667: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
668: {
669:   PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
670:   PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
671:   void            *exactCtxs[3];
672:   MPI_Comm         comm;
673:   PetscReal        error, errorDer, tol = PETSC_SMALL;

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

710: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
711: {
712:   PetscErrorCode (*exactFuncs[1]) (PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
713:   PetscErrorCode (*exactFuncDers[1]) (PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx);
714:   PetscReal       n[3]         = {1.0, 1.0, 1.0};
715:   void           *exactCtxs[3];
716:   DM              rdm, idm, fdm;
717:   Mat             Interp;
718:   Vec             iu, fu, scaling;
719:   MPI_Comm        comm;
720:   PetscInt        dim;
721:   PetscReal       error, errorDer, tol = PETSC_SMALL;

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

788: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
789: {
790:   DM               odm = dm, rdm = NULL, cdm = NULL;
791:   PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {trig};
792:   PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx) = {trigDer};
793:   void            *exactCtxs[3];
794:   PetscInt         r, c, cStart, cEnd;
795:   PetscReal        errorOld, errorDerOld, error, errorDer, rel, len, lenOld;
796:   double           p;

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

857: int main(int argc, char **argv)
858: {
859:   DM             dm;
860:   AppCtx         user;                 /* user-defined work context */
861:   PetscInt       dim = 2;
862:   PetscBool      simplex = PETSC_FALSE;

864:   PetscInitialize(&argc, &argv, NULL, help);
865:   ProcessOptions(PETSC_COMM_WORLD, &user);
866:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
867:   if (!user.useDA) {
868:     DMGetDimension(dm, &dim);
869:     DMPlexIsSimplex(dm, &simplex);
870:   }
871:   DMPlexMetricSetFromOptions(dm);
872:   user.numComponents = user.numComponents < 0 ? dim : user.numComponents;
873:   PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.numComponents, simplex, NULL, user.qorder, &user.fe);
874:   SetupSection(dm, &user);
875:   if (user.testFEjacobian) TestFEJacobian(dm, &user);
876:   if (user.testFVgrad) TestFVGrad(dm, &user);
877:   if (user.testInjector) TestInjector(dm, &user);
878:   CheckFunctions(dm, user.porder, &user);
879:   {
880:     PetscDualSpace dsp;
881:     PetscInt       k;

883:     PetscFEGetDualSpace(user.fe, &dsp);
884:     PetscDualSpaceGetDeRahm(dsp, &k);
885:     if (dim == 2 && simplex == PETSC_TRUE && user.tree == PETSC_FALSE && k == 0) {
886:       CheckInterpolation(dm, PETSC_FALSE, user.porder, &user);
887:       CheckInterpolation(dm, PETSC_TRUE,  user.porder, &user);
888:     }
889:   }
890:   CheckConvergence(dm, 3, &user);
891:   PetscFEDestroy(&user.fe);
892:   DMDestroy(&dm);
893:   PetscFinalize();
894:   return 0;
895: }

897: /*TEST

899:   test:
900:     suffix: 1
901:     requires: triangle

903:   # 2D P_1 on a triangle
904:   test:
905:     suffix: p1_2d_0
906:     requires: triangle
907:     args: -petscspace_degree 1 -qorder 1 -convergence
908:   test:
909:     suffix: p1_2d_1
910:     requires: triangle
911:     args: -petscspace_degree 1 -qorder 1 -porder 1
912:   test:
913:     suffix: p1_2d_2
914:     requires: triangle
915:     args: -petscspace_degree 1 -qorder 1 -porder 2
916:   test:
917:     suffix: p1_2d_3
918:     requires: triangle mmg
919:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
920:   test:
921:     suffix: p1_2d_4
922:     requires: triangle mmg
923:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
924:   test:
925:     suffix: p1_2d_5
926:     requires: triangle mmg
927:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0

929:   # 3D P_1 on a tetrahedron
930:   test:
931:     suffix: p1_3d_0
932:     requires: ctetgen
933:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -convergence
934:   test:
935:     suffix: p1_3d_1
936:     requires: ctetgen
937:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 1
938:   test:
939:     suffix: p1_3d_2
940:     requires: ctetgen
941:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 2
942:   test:
943:     suffix: p1_3d_3
944:     requires: ctetgen mmg
945:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
946:   test:
947:     suffix: p1_3d_4
948:     requires: ctetgen mmg
949:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
950:   test:
951:     suffix: p1_3d_5
952:     requires: ctetgen mmg
953:     args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0

955:   # 2D P_2 on a triangle
956:   test:
957:     suffix: p2_2d_0
958:     requires: triangle
959:     args: -petscspace_degree 2 -qorder 2 -convergence
960:   test:
961:     suffix: p2_2d_1
962:     requires: triangle
963:     args: -petscspace_degree 2 -qorder 2 -porder 1
964:   test:
965:     suffix: p2_2d_2
966:     requires: triangle
967:     args: -petscspace_degree 2 -qorder 2 -porder 2
968:   test:
969:     suffix: p2_2d_3
970:     requires: triangle mmg
971:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
972:   test:
973:     suffix: p2_2d_4
974:     requires: triangle mmg
975:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
976:   test:
977:     suffix: p2_2d_5
978:     requires: triangle mmg
979:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0

981:   # 3D P_2 on a tetrahedron
982:   test:
983:     suffix: p2_3d_0
984:     requires: ctetgen
985:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -convergence
986:   test:
987:     suffix: p2_3d_1
988:     requires: ctetgen
989:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 1
990:   test:
991:     suffix: p2_3d_2
992:     requires: ctetgen
993:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 2
994:   test:
995:     suffix: p2_3d_3
996:     requires: ctetgen mmg
997:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
998:   test:
999:     suffix: p2_3d_4
1000:     requires: ctetgen mmg
1001:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1002:   test:
1003:     suffix: p2_3d_5
1004:     requires: ctetgen mmg
1005:     args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0

1007:   # 2D Q_1 on a quadrilaterial DA
1008:   test:
1009:     suffix: q1_2d_da_0
1010:     requires: broken
1011:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -convergence
1012:   test:
1013:     suffix: q1_2d_da_1
1014:     requires: broken
1015:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 1
1016:   test:
1017:     suffix: q1_2d_da_2
1018:     requires: broken
1019:     args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 2

1021:   # 2D Q_1 on a quadrilaterial Plex
1022:   test:
1023:     suffix: q1_2d_plex_0
1024:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -convergence
1025:   test:
1026:     suffix: q1_2d_plex_1
1027:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1
1028:   test:
1029:     suffix: q1_2d_plex_2
1030:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2
1031:   test:
1032:     suffix: q1_2d_plex_3
1033:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1 -shear_coords
1034:   test:
1035:     suffix: q1_2d_plex_4
1036:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2 -shear_coords
1037:   test:
1038:     suffix: q1_2d_plex_5
1039:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 0 -non_affine_coords -convergence
1040:   test:
1041:     suffix: q1_2d_plex_6
1042:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 1 -non_affine_coords -convergence
1043:   test:
1044:     suffix: q1_2d_plex_7
1045:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 2 -non_affine_coords -convergence

1047:   # 2D Q_2 on a quadrilaterial
1048:   test:
1049:     suffix: q2_2d_plex_0
1050:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1051:   test:
1052:     suffix: q2_2d_plex_1
1053:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1054:   test:
1055:     suffix: q2_2d_plex_2
1056:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1057:   test:
1058:     suffix: q2_2d_plex_3
1059:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1060:   test:
1061:     suffix: q2_2d_plex_4
1062:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1063:   test:
1064:     suffix: q2_2d_plex_5
1065:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 0 -non_affine_coords -convergence
1066:   test:
1067:     suffix: q2_2d_plex_6
1068:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 1 -non_affine_coords -convergence
1069:   test:
1070:     suffix: q2_2d_plex_7
1071:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 2 -non_affine_coords -convergence

1073:   # 2D P_3 on a triangle
1074:   test:
1075:     suffix: p3_2d_0
1076:     requires: triangle !single
1077:     args: -petscspace_degree 3 -qorder 3 -convergence
1078:   test:
1079:     suffix: p3_2d_1
1080:     requires: triangle !single
1081:     args: -petscspace_degree 3 -qorder 3 -porder 1
1082:   test:
1083:     suffix: p3_2d_2
1084:     requires: triangle !single
1085:     args: -petscspace_degree 3 -qorder 3 -porder 2
1086:   test:
1087:     suffix: p3_2d_3
1088:     requires: triangle !single
1089:     args: -petscspace_degree 3 -qorder 3 -porder 3
1090:   test:
1091:     suffix: p3_2d_4
1092:     requires: triangle mmg
1093:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -convergence -conv_refine 0
1094:   test:
1095:     suffix: p3_2d_5
1096:     requires: triangle mmg
1097:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 1 -conv_refine 0
1098:   test:
1099:     suffix: p3_2d_6
1100:     requires: triangle mmg
1101:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 3 -conv_refine 0

1103:   # 2D Q_3 on a quadrilaterial
1104:   test:
1105:     suffix: q3_2d_0
1106:     requires: !single
1107:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -convergence
1108:   test:
1109:     suffix: q3_2d_1
1110:     requires: !single
1111:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 1
1112:   test:
1113:     suffix: q3_2d_2
1114:     requires: !single
1115:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 2
1116:   test:
1117:     suffix: q3_2d_3
1118:     requires: !single
1119:     args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 3

1121:   # 2D P_1disc on a triangle/quadrilateral
1122:   test:
1123:     suffix: p1d_2d_0
1124:     requires: triangle
1125:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1126:   test:
1127:     suffix: p1d_2d_1
1128:     requires: triangle
1129:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1130:   test:
1131:     suffix: p1d_2d_2
1132:     requires: triangle
1133:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1134:   test:
1135:     suffix: p1d_2d_3
1136:     requires: triangle
1137:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1138:     filter: sed  -e "s/convergence rate at refinement 0: 2/convergence rate at refinement 0: 1.9/g"
1139:   test:
1140:     suffix: p1d_2d_4
1141:     requires: triangle
1142:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1143:   test:
1144:     suffix: p1d_2d_5
1145:     requires: triangle
1146:     args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2

1148:   # 2D BDM_1 on a triangle
1149:   test:
1150:     suffix: bdm1_2d_0
1151:     requires: triangle
1152:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1153:           -num_comp 2 -qorder 1 -convergence
1154:   test:
1155:     suffix: bdm1_2d_1
1156:     requires: triangle
1157:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1158:           -num_comp 2 -qorder 1 -porder 1
1159:   test:
1160:     suffix: bdm1_2d_2
1161:     requires: triangle
1162:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1163:           -num_comp 2 -qorder 1 -porder 2

1165:   # 2D BDM_1 on a quadrilateral
1166:   test:
1167:     suffix: bdm1q_2d_0
1168:     requires: triangle
1169:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1170:           -petscdualspace_lagrange_tensor 1 \
1171:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -convergence
1172:   test:
1173:     suffix: bdm1q_2d_1
1174:     requires: triangle
1175:     args: -petscspace_degree 1 -petscdualspace_type bdm \
1176:           -petscdualspace_lagrange_tensor 1 \
1177:           -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 1
1178:   test:
1179:     suffix: bdm1q_2d_2
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 -porder 2

1185:   # Test high order quadrature
1186:   test:
1187:     suffix: p1_quad_2
1188:     requires: triangle
1189:     args: -petscspace_degree 1 -qorder 2 -porder 1
1190:   test:
1191:     suffix: p1_quad_5
1192:     requires: triangle
1193:     args: -petscspace_degree 1 -qorder 5 -porder 1
1194:   test:
1195:     suffix: p2_quad_3
1196:     requires: triangle
1197:     args: -petscspace_degree 2 -qorder 3 -porder 2
1198:   test:
1199:     suffix: p2_quad_5
1200:     requires: triangle
1201:     args: -petscspace_degree 2 -qorder 5 -porder 2
1202:   test:
1203:     suffix: q1_quad_2
1204:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 2 -porder 1
1205:   test:
1206:     suffix: q1_quad_5
1207:     args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 5 -porder 1
1208:   test:
1209:     suffix: q2_quad_3
1210:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 3 -porder 1
1211:   test:
1212:     suffix: q2_quad_5
1213:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 5 -porder 1

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

1276:   # 3D WXY on a triangular prism
1277:   test:
1278:     suffix: wxy_0
1279:     args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism -qorder 2 -porder 0 \
1280:           -petscspace_type sum \
1281:           -petscspace_variables 3 \
1282:           -petscspace_components 3 \
1283:           -petscspace_sum_spaces 2 \
1284:           -petscspace_sum_concatenate false \
1285:           -sumcomp_0_petscspace_variables 3 \
1286:           -sumcomp_0_petscspace_components 3 \
1287:           -sumcomp_0_petscspace_degree 1 \
1288:           -sumcomp_1_petscspace_variables 3 \
1289:           -sumcomp_1_petscspace_components 3 \
1290:           -sumcomp_1_petscspace_type wxy \
1291:           -petscdualspace_refcell triangular_prism \
1292:           -petscdualspace_form_degree 0 \
1293:           -petscdualspace_order 1 \
1294:           -petscdualspace_components 3

1296: TEST*/

1298: /*
1299:    # 2D Q_2 on a quadrilaterial Plex
1300:   test:
1301:     suffix: q2_2d_plex_0
1302:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1303:   test:
1304:     suffix: q2_2d_plex_1
1305:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1306:   test:
1307:     suffix: q2_2d_plex_2
1308:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1309:   test:
1310:     suffix: q2_2d_plex_3
1311:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1312:   test:
1313:     suffix: q2_2d_plex_4
1314:     args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1315:   test:
1316:     suffix: q2_2d_plex_5
1317:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 0 -non_affine_coords
1318:   test:
1319:     suffix: q2_2d_plex_6
1320:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 1 -non_affine_coords
1321:   test:
1322:     suffix: q2_2d_plex_7
1323:     args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 2 -non_affine_coords

1325:   test:
1326:     suffix: p1d_2d_6
1327:     requires: mmg
1328:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1329:   test:
1330:     suffix: p1d_2d_7
1331:     requires: mmg
1332:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1333:   test:
1334:     suffix: p1d_2d_8
1335:     requires: mmg
1336:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1337: */