Actual source code: ex3.c

petsc-3.10.3 2018-12-18
Report Typos and Errors
  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:   PetscInt  debug;             /* The debugging level */
 13:   /* Domain and mesh definition */
 14:   PetscInt  dim;               /* The topological mesh dimension */
 15:   PetscBool simplex;           /* Flag for simplex or tensor product mesh */
 16:   PetscBool useDA;             /* Flag DMDA tensor product mesh */
 17:   PetscBool interpolate;       /* Generate intermediate mesh elements */
 18:   PetscReal refinementLimit;   /* The largest allowable cell volume */
 19:   PetscBool shearCoords;       /* Flag for shear transform */
 20:   PetscBool nonaffineCoords;   /* Flag for non-affine transform */
 21:   /* Element definition */
 22:   PetscInt  qorder;            /* Order of the quadrature */
 23:   PetscInt  numComponents;     /* Number of field components */
 24:   PetscFE   fe;                /* The finite element */
 25:   /* Testing space */
 26:   PetscInt  porder;            /* Order of polynomials to test */
 27:   PetscBool convergence;       /* Test for order of convergence */
 28:   PetscBool convRefine;        /* Test for convergence using refinement, otherwise use coarsening */
 29:   PetscBool constraints;       /* Test local constraints */
 30:   PetscBool tree;              /* Test tree routines */
 31:   PetscBool testFEjacobian;    /* Test finite element Jacobian assembly */
 32:   PetscBool testFVgrad;        /* Test finite difference gradient routine */
 33:   PetscBool testInjector;      /* Test finite element injection routines */
 34:   PetscInt  treeCell;          /* Cell to refine in tree test */
 35:   PetscReal constants[3];      /* Constant values for each dimension */
 36: } AppCtx;

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

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

 71: /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */
 72: PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 73: {
 74:   AppCtx *user = (AppCtx *) ctx;
 75:   if (user->dim > 2)      {u[0] = coords[0]*coords[1]; u[1] = coords[1]*coords[2]; u[2] = coords[2]*coords[0];}
 76:   else if (user->dim > 1) {u[0] = coords[0]*coords[0]; u[1] = coords[0]*coords[1];}
 77:   else if (user->dim > 0) {u[0] = coords[0]*coords[0];}
 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:   AppCtx *user = (AppCtx *) ctx;
 83:   if (user->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];}
 84:   else if (user->dim > 1) {u[0] = 2.0*coords[0]*n[0]; u[1] = coords[1]*n[0] + coords[0]*n[1];}
 85:   else if (user->dim > 0) {u[0] = 2.0*coords[0]*n[0];}
 86:   return 0;
 87: }

 89: /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
 90: PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 91: {
 92:   AppCtx *user = (AppCtx *) ctx;
 93:   if (user->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];}
 94:   else if (user->dim > 1) {u[0] = coords[0]*coords[0]*coords[0]; u[1] = coords[0]*coords[0]*coords[1];}
 95:   else if (user->dim > 0) {u[0] = coords[0]*coords[0]*coords[0];}
 96:   return 0;
 97: }
 98: PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 99: {
100:   AppCtx *user = (AppCtx *) ctx;
101:   if (user->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];}
102:   else if (user->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];}
103:   else if (user->dim > 0) {u[0] = 3.0*coords[0]*coords[0]*n[0];}
104:   return 0;
105: }

107: /* u = tanh(x) */
108: PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
109: {
110:   AppCtx   *user = (AppCtx *) ctx;
111:   PetscInt d;
112:   for (d = 0; d < user->dim; ++d) u[d] = PetscTanhReal(coords[d] - 0.5);
113:   return 0;
114: }
115: PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
116: {
117:   AppCtx   *user = (AppCtx *) ctx;
118:   PetscInt d;
119:   for (d = 0; d < user->dim; ++d) u[d] = 1.0/PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d];
120:   return 0;
121: }

123: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
124: {

128:   options->debug           = 0;
129:   options->dim             = 2;
130:   options->simplex         = PETSC_TRUE;
131:   options->useDA           = PETSC_TRUE;
132:   options->interpolate     = PETSC_TRUE;
133:   options->refinementLimit = 0.0;
134:   options->shearCoords     = PETSC_FALSE;
135:   options->nonaffineCoords = PETSC_FALSE;
136:   options->qorder          = 0;
137:   options->numComponents   = PETSC_DEFAULT;
138:   options->porder          = 0;
139:   options->convergence     = PETSC_FALSE;
140:   options->convRefine      = PETSC_TRUE;
141:   options->constraints     = PETSC_FALSE;
142:   options->tree            = PETSC_FALSE;
143:   options->treeCell        = 0;
144:   options->testFEjacobian  = PETSC_FALSE;
145:   options->testFVgrad      = PETSC_FALSE;
146:   options->testInjector    = PETSC_FALSE;

148:   PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
149:   PetscOptionsInt("-debug", "The debugging level", "ex3.c", options->debug, &options->debug, NULL);
150:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex3.c", options->dim, &options->dim, NULL);
151:   PetscOptionsBool("-simplex", "Flag for simplices or hexhedra", "ex3.c", options->simplex, &options->simplex, NULL);
152:   PetscOptionsBool("-use_da", "Flag for DMDA mesh", "ex3.c", options->useDA, &options->useDA, NULL);
153:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex3.c", options->interpolate, &options->interpolate, NULL);
154:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex3.c", options->refinementLimit, &options->refinementLimit, NULL);
155:   PetscOptionsBool("-shear_coords", "Transform coordinates with a shear", "ex3.c", options->shearCoords, &options->shearCoords, NULL);
156:   PetscOptionsBool("-non_affine_coords", "Transform coordinates with a non-affine transform", "ex3.c", options->nonaffineCoords, &options->nonaffineCoords, NULL);
157:   PetscOptionsInt("-qorder", "The quadrature order", "ex3.c", options->qorder, &options->qorder, NULL);
158:   PetscOptionsInt("-num_comp", "The number of field components", "ex3.c", options->numComponents, &options->numComponents, NULL);
159:   PetscOptionsInt("-porder", "The order of polynomials to test", "ex3.c", options->porder, &options->porder, NULL);
160:   PetscOptionsBool("-convergence", "Check the convergence rate", "ex3.c", options->convergence, &options->convergence, NULL);
161:   PetscOptionsBool("-conv_refine", "Use refinement for the convergence rate", "ex3.c", options->convRefine, &options->convRefine, NULL);
162:   PetscOptionsBool("-constraints", "Test local constraints (serial only)", "ex3.c", options->constraints, &options->constraints, NULL);
163:   PetscOptionsBool("-tree", "Test tree routines", "ex3.c", options->tree, &options->tree, NULL);
164:   PetscOptionsInt("-tree_cell", "cell to refine in tree test", "ex3.c", options->treeCell, &options->treeCell, NULL);
165:   PetscOptionsBool("-test_fe_jacobian", "Test finite element Jacobian assembly", "ex3.c", options->testFEjacobian, &options->testFEjacobian, NULL);
166:   PetscOptionsBool("-test_fv_grad", "Test finite volume gradient reconstruction", "ex3.c", options->testFVgrad, &options->testFVgrad, NULL);
167:   PetscOptionsBool("-test_injector","Test finite element injection", "ex3.c", options->testInjector, &options->testInjector,NULL);
168:   PetscOptionsEnd();

170:   options->numComponents = options->numComponents < 0 ? options->dim : options->numComponents;

172:   return(0);
173: }

175: static PetscErrorCode TransformCoordinates(DM dm, AppCtx *user)
176: {
177:   PetscSection   coordSection;
178:   Vec            coordinates;
179:   PetscScalar   *coords;
180:   PetscInt       vStart, vEnd, v;

184:   if (user->nonaffineCoords) {
185:     /* x' = r^(1/p) (x/r), y' = r^(1/p) (y/r), z' = z */
186:     DMGetCoordinateSection(dm, &coordSection);
187:     DMGetCoordinatesLocal(dm, &coordinates);
188:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
189:     VecGetArray(coordinates, &coords);
190:     for (v = vStart; v < vEnd; ++v) {
191:       PetscInt  dof, off;
192:       PetscReal p = 4.0, r;

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

222:       PetscSectionGetDof(coordSection, v, &dof);
223:       PetscSectionGetOffset(coordSection, v, &off);
224:       switch (dof) {
225:       case 2:
226:         coords[off+0] = coords[off+0] + m*coords[off+1];
227:         coords[off+1] = coords[off+1];
228:         break;
229:       case 3:
230:         coords[off+0] = coords[off+0] + m*coords[off+1] + m*coords[off+2];
231:         coords[off+1] = coords[off+1] + m*coords[off+2];
232:         coords[off+2] = coords[off+2];
233:         break;
234:       }
235:     }
236:     VecRestoreArray(coordinates, &coords);
237:   }
238:   return(0);
239: }

241: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
242: {
243:   PetscInt       dim             = user->dim;
244:   PetscBool      interpolate     = user->interpolate;
245:   PetscReal      refinementLimit = user->refinementLimit;
246:   PetscBool      isPlex;

250:   if (user->simplex) {
251:     DM refinedMesh = NULL;

253:     DMPlexCreateBoxMesh(comm, dim, PETSC_TRUE, NULL, NULL, NULL, NULL, interpolate, dm);
254:     /* Refine mesh using a volume constraint */
255:     DMPlexSetRefinementLimit(*dm, refinementLimit);
256:     DMRefine(*dm, comm, &refinedMesh);
257:     if (refinedMesh) {
258:       DMDestroy(dm);
259:       *dm  = refinedMesh;
260:     }
261:     DMPlexSetRefinementUniform(*dm, PETSC_TRUE);
262:   } else {
263:     if (user->constraints || user->tree || !user->useDA) {
264:       PetscInt cells[3] = {2, 2, 2};

266:       PetscOptionsGetInt(NULL,NULL,"-da_grid_x",&cells[0],NULL);
267:       PetscOptionsGetInt(NULL,NULL,"-da_grid_y",&cells[1],NULL);
268:       PetscOptionsGetInt(NULL,NULL,"-da_grid_z",&cells[2],NULL);
269:       DMPlexCreateBoxMesh(comm, dim, PETSC_FALSE, cells, NULL, NULL, NULL, PETSC_TRUE, dm);
270:     } else {
271:       switch (user->dim) {
272:       case 2:
273:         DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 2, 2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, dm);
274:         DMSetFromOptions(*dm);
275:         DMSetUp(*dm);
276:         DMDASetVertexCoordinates(*dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
277:         break;
278:       default:
279:         SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot create structured mesh of dimension %d", dim);
280:       }
281:       PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");
282:     }
283:   }
284:   PetscObjectTypeCompare((PetscObject)*dm,DMPLEX,&isPlex);
285:   if (isPlex) {
286:     PetscPartitioner part;
287:     DM               distributedMesh = NULL;

289:     if (user->tree) {
290:       DM refTree;
291:       DM ncdm = NULL;

293:       DMPlexCreateDefaultReferenceTree(comm,user->dim,user->simplex,&refTree);
294:       DMPlexSetReferenceTree(*dm,refTree);
295:       DMDestroy(&refTree);
296:       DMPlexTreeRefineCell(*dm,user->treeCell,&ncdm);
297:       if (ncdm) {
298:         DMDestroy(dm);
299:         *dm = ncdm;
300:         DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
301:       }
302:     } else {
303:       DMPlexSetRefinementUniform(*dm, PETSC_TRUE);
304:     }
305:     /* Distribute mesh over processes */
306:     DMPlexGetPartitioner(*dm,&part);
307:     PetscPartitionerSetFromOptions(part);
308:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
309:     if (distributedMesh) {
310:       DMDestroy(dm);
311:       *dm  = distributedMesh;
312:     }
313:     if (user->simplex) {
314:       PetscObjectSetName((PetscObject) *dm, "Simplicial Mesh");
315:     } else {
316:       PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");
317:     }
318:   }
319:   DMSetFromOptions(*dm);
320:   TransformCoordinates(*dm, user);
321:   DMViewFromOptions(*dm,NULL,"-dm_view");
322:   return(0);
323: }

325: static void simple_mass(PetscInt dim, PetscInt Nf, PetscInt NfAux,
326:                         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
327:                         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
328:                         PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
329: {
330:   PetscInt d, e;
331:   for (d = 0, e = 0; d < dim; d++, e+=dim+1) {
332:     g0[e] = 1.;
333:   }
334: }

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

344:   for (compI = 0; compI < dim; ++compI) {
345:     for (compJ = 0; compJ < dim; ++compJ) {
346:       for (d = 0; d < dim; ++d) {
347:         for (e = 0; e < dim; e++) {
348:           if (d == e && d == compI && d == compJ) {
349:             C[((compI*dim+compJ)*dim+d)*dim+e] = 1.0;
350:           } else if ((d == compJ && e == compI) || (d == e && compI == compJ)) {
351:             C[((compI*dim+compJ)*dim+d)*dim+e] = 0.5;
352:           } else {
353:             C[((compI*dim+compJ)*dim+d)*dim+e] = 0.0;
354:           }
355:         }
356:       }
357:     }
358:   }
359: }

361: static PetscErrorCode SetupSection(DM dm, AppCtx *user)
362: {
363:   PetscBool      isPlex;

367:   if (!user->simplex && user->constraints) {
368:     /* test local constraints */
369:     DM            coordDM;
370:     PetscInt      fStart, fEnd, f, vStart, vEnd, v;
371:     PetscInt      edgesx = 2, vertsx;
372:     PetscInt      edgesy = 2, vertsy;
373:     PetscMPIInt   size;
374:     PetscInt      numConst;
375:     PetscSection  aSec;
376:     PetscInt     *anchors;
377:     PetscInt      offset;
378:     IS            aIS;
379:     MPI_Comm      comm = PetscObjectComm((PetscObject)dm);

381:     MPI_Comm_size(comm,&size);
382:     if (size > 1) SETERRQ(comm,PETSC_ERR_SUP,"Local constraint test can only be performed in serial");

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

388:     /* first create the coordinate section so that it does not clone the
389:      * constraints */
390:     DMGetCoordinateDM(dm,&coordDM);

392:     /* create the constrained-to-anchor section */
393:     DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
394:     DMPlexGetDepthStratum(dm,1,&fStart,&fEnd);
395:     PetscSectionCreate(PETSC_COMM_SELF,&aSec);
396:     PetscSectionSetChart(aSec,PetscMin(fStart,vStart),PetscMax(fEnd,vEnd));

398:     /* define the constraints */
399:     PetscOptionsGetInt(NULL,NULL,"-da_grid_x",&edgesx,NULL);
400:     PetscOptionsGetInt(NULL,NULL,"-da_grid_y",&edgesy,NULL);
401:     vertsx = edgesx + 1;
402:     vertsy = edgesy + 1;
403:     numConst = vertsy + edgesy;
404:     PetscMalloc1(numConst,&anchors);
405:     offset = 0;
406:     for (v = vStart + edgesx; v < vEnd; v+= vertsx) {
407:       PetscSectionSetDof(aSec,v,1);
408:       anchors[offset++] = v - edgesx;
409:     }
410:     for (f = fStart + edgesx * vertsy + edgesx * edgesy; f < fEnd; f++) {
411:       PetscSectionSetDof(aSec,f,1);
412:       anchors[offset++] = f - edgesx * edgesy;
413:     }
414:     PetscSectionSetUp(aSec);
415:     ISCreateGeneral(PETSC_COMM_SELF,numConst,anchors,PETSC_OWN_POINTER,&aIS);

417:     DMPlexSetAnchors(dm,aSec,aIS);
418:     PetscSectionDestroy(&aSec);
419:     ISDestroy(&aIS);
420:   }
421:   DMSetNumFields(dm, 1);
422:   DMSetField(dm, 0, (PetscObject) user->fe);
423:   PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
424:   if (!isPlex) {
425:     PetscSection    section;
426:     const PetscInt *numDof;
427:     PetscInt        numComp;

429:     PetscFEGetNumComponents(user->fe, &numComp);
430:     PetscFEGetNumDof(user->fe, &numDof);
431:     DMDACreateSection(dm, &numComp, numDof, NULL, &section);
432:     DMSetSection(dm, section);
433:     PetscSectionDestroy(&section);
434:   }
435:   if (!user->simplex && user->constraints) {
436:     /* test getting local constraint matrix that matches section */
437:     PetscSection aSec;
438:     IS           aIS;

440:     DMPlexGetAnchors(dm,&aSec,&aIS);
441:     if (aSec) {
442:       PetscDS         ds;
443:       PetscSection    cSec, section;
444:       PetscInt        cStart, cEnd, c, numComp;
445:       Mat             cMat, mass;
446:       Vec             local;
447:       const PetscInt *anchors;

449:       DMGetSection(dm,&section);
450:       /* this creates the matrix and preallocates the matrix structure: we
451:        * just have to fill in the values */
452:       DMGetDefaultConstraints(dm,&cSec,&cMat);
453:       PetscSectionGetChart(cSec,&cStart,&cEnd);
454:       ISGetIndices(aIS,&anchors);
455:       PetscFEGetNumComponents(user->fe, &numComp);
456:       for (c = cStart; c < cEnd; c++) {
457:         PetscInt cDof;

459:         /* is this point constrained? (does it have an anchor?) */
460:         PetscSectionGetDof(aSec,c,&cDof);
461:         if (cDof) {
462:           PetscInt cOff, a, aDof, aOff, j;
463:           if (cDof != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Found %d anchor points: should be just one",cDof);

465:           /* find the anchor point */
466:           PetscSectionGetOffset(aSec,c,&cOff);
467:           a    = anchors[cOff];

469:           /* find the constrained dofs (row in constraint matrix) */
470:           PetscSectionGetDof(cSec,c,&cDof);
471:           PetscSectionGetOffset(cSec,c,&cOff);

473:           /* find the anchor dofs (column in constraint matrix) */
474:           PetscSectionGetDof(section,a,&aDof);
475:           PetscSectionGetOffset(section,a,&aOff);

477:           if (cDof != aDof) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point and anchor have different number of dofs: %d, %d\n",cDof,aDof);
478:           if (cDof % numComp) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point dofs not divisible by field components: %d, %d\n",cDof,numComp);

480:           /* put in a simple equality constraint */
481:           for (j = 0; j < cDof; j++) {
482:             MatSetValue(cMat,cOff+j,aOff+j,1.,INSERT_VALUES);
483:           }
484:         }
485:       }
486:       MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
487:       MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
488:       ISRestoreIndices(aIS,&anchors);

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

493:       DMCreateMatrix(dm,&mass);
494:       /* get a dummy local variable to serve as the solution */
495:       DMGetLocalVector(dm,&local);
496:       DMGetDS(dm,&ds);
497:       /* set the jacobian to be the mass matrix */
498:       PetscDSSetJacobian(ds, 0, 0, simple_mass, NULL,  NULL, NULL);
499:       /* build the mass matrix */
500:       DMPlexSNESComputeJacobianFEM(dm,local,mass,mass,NULL);
501:       MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
502:       MatDestroy(&mass);
503:       DMRestoreLocalVector(dm,&local);
504: #if 0
505:       {
506:         /* compare this to periodicity with DMDA: this is broken right now
507:          * because DMCreateMatrix() doesn't respect the default section that I
508:          * set */
509:         DM              dmda;
510:         PetscSection    section;
511:         const PetscInt *numDof;
512:         PetscInt        numComp;

514:                                                               /* periodic x */
515:         DMDACreate2d(PetscObjectComm((PetscObject)dm), DM_BOUNDARY_PERIODIC, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, -2, -2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dmda);
516:         DMSetFromOptions(dmda);
517:         DMSetUp(dmda);
518:         DMDASetVertexCoordinates(dmda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);


521:         PetscFEGetNumComponents(user->fe, &numComp);
522:         PetscFEGetNumDof(user->fe, &numDof);
523:         DMDACreateSection(dmda, &numComp, numDof, NULL, &section);
524:         DMSetSection(dmda, section);
525:         PetscSectionDestroy(&section);
526:         DMCreateMatrix(dmda,&mass);
527:         /* there isn't a DMDA equivalent of DMPlexSNESComputeJacobianFEM()
528:          * right now, but we can at least verify the nonzero structure */
529:         MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
530:         MatDestroy(&mass);
531:         DMDestroy(&dmda);
532:       }
533: #endif
534:     }
535:   }
536:   return(0);
537: }

539: static PetscErrorCode TestFEJacobian(DM dm, AppCtx *user)
540: {
541:   PetscBool      isPlex;

545:   PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
546:   if (isPlex) {
547:     Vec          local;
548:     const Vec    *vecs;
549:     Mat          E;
550:     MatNullSpace sp;
551:     PetscBool    isNullSpace, hasConst;
552:     PetscInt     n, i;
553:     Vec          res, localX, localRes;
554:     PetscDS      ds;

556:     if (user->numComponents != user->dim) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "The number of components %d must be equal to the dimension %d for this test", user->numComponents, user->dim);
557:     DMGetDS(dm,&ds);
558:     PetscDSSetJacobian(ds,0,0,NULL,NULL,NULL,symmetric_gradient_inner_product);
559:     DMCreateMatrix(dm,&E);
560:     DMGetLocalVector(dm,&local);
561:     DMPlexSNESComputeJacobianFEM(dm,local,E,E,NULL);
562:     DMPlexCreateRigidBody(dm,&sp);
563:     MatNullSpaceGetVecs(sp,&hasConst,&n,&vecs);
564:     VecDuplicate(vecs[0],&res);
565:     DMCreateLocalVector(dm,&localX);
566:     DMCreateLocalVector(dm,&localRes);
567:     for (i = 0; i < n; i++) { /* also test via matrix-free Jacobian application */
568:       PetscReal resNorm;

570:       VecSet(localRes,0.);
571:       VecSet(localX,0.);
572:       VecSet(local,0.);
573:       VecSet(res,0.);
574:       DMGlobalToLocalBegin(dm,vecs[i],INSERT_VALUES,localX);
575:       DMGlobalToLocalEnd(dm,vecs[i],INSERT_VALUES,localX);
576:       DMPlexComputeJacobianAction(dm,NULL,0,0,local,NULL,localX,localRes,NULL);
577:       DMLocalToGlobalBegin(dm,localRes,ADD_VALUES,res);
578:       DMLocalToGlobalEnd(dm,localRes,ADD_VALUES,res);
579:       VecNorm(res,NORM_2,&resNorm);
580:       if (resNorm > PETSC_SMALL) {
581:         PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient action null space vector %D residual: %E\n",i,resNorm);
582:       }
583:     }
584:     VecDestroy(&localRes);
585:     VecDestroy(&localX);
586:     VecDestroy(&res);
587:     MatNullSpaceTest(sp,E,&isNullSpace);
588:     if (isNullSpace) {
589:       PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient null space: PASS\n");
590:     } else {
591:       PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient null space: FAIL\n");
592:     }
593:     MatNullSpaceDestroy(&sp);
594:     MatDestroy(&E);
595:     DMRestoreLocalVector(dm,&local);
596:   }
597:   return(0);
598: }

600: static PetscErrorCode TestInjector(DM dm, AppCtx *user)
601: {
602:   DM             refTree;
603:   PetscMPIInt    rank;

607:   DMPlexGetReferenceTree(dm,&refTree);
608:   if (refTree) {
609:     Mat inj;

611:     DMPlexComputeInjectorReferenceTree(refTree,&inj);
612:     PetscObjectSetName((PetscObject)inj,"Reference Tree Injector");
613:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
614:     if (!rank) {
615:       MatView(inj,PETSC_VIEWER_STDOUT_SELF);
616:     }
617:     MatDestroy(&inj);
618:   }
619:   return(0);
620: }

622: static PetscErrorCode TestFVGrad(DM dm, AppCtx *user)
623: {
624:   MPI_Comm          comm;
625:   DM                dmRedist, dmfv, dmgrad, dmCell, refTree;
626:   PetscFV           fv;
627:   PetscInt          nvecs, v, cStart, cEnd, cEndInterior;
628:   PetscMPIInt       size;
629:   Vec               cellgeom, grad, locGrad;
630:   const PetscScalar *cgeom;
631:   PetscReal         allVecMaxDiff = 0., fvTol = 100. * PETSC_MACHINE_EPSILON;
632:   PetscErrorCode    ierr;

635:   comm = PetscObjectComm((PetscObject)dm);
636:   /* duplicate DM, give dup. a FV discretization */
637:   DMPlexSetAdjacencyUseCone(dm,PETSC_TRUE);
638:   DMPlexSetAdjacencyUseClosure(dm,PETSC_FALSE);
639:   MPI_Comm_size(comm,&size);
640:   dmRedist = NULL;
641:   if (size > 1) {
642:     DMPlexDistributeOverlap(dm,1,NULL,&dmRedist);
643:   }
644:   if (!dmRedist) {
645:     dmRedist = dm;
646:     PetscObjectReference((PetscObject)dmRedist);
647:   }
648:   PetscFVCreate(comm,&fv);
649:   PetscFVSetType(fv,PETSCFVLEASTSQUARES);
650:   PetscFVSetNumComponents(fv,user->numComponents);
651:   PetscFVSetSpatialDimension(fv,user->dim);
652:   PetscFVSetFromOptions(fv);
653:   PetscFVSetUp(fv);
654:   DMPlexConstructGhostCells(dmRedist,NULL,NULL,&dmfv);
655:   DMDestroy(&dmRedist);
656:   DMSetNumFields(dmfv,1);
657:   DMSetField(dmfv, 0, (PetscObject) fv);
658:   DMPlexGetReferenceTree(dm,&refTree);
659:   if (refTree) {
660:     PetscDS ds;
661:     DMGetDS(dmfv,&ds);
662:     DMSetDS(refTree,ds);
663:   }
664:   DMPlexSNESGetGradientDM(dmfv, fv, &dmgrad);
665:   DMPlexGetHeightStratum(dmfv,0,&cStart,&cEnd);
666:   nvecs = user->dim * (user->dim+1) / 2;
667:   DMPlexSNESGetGeometryFVM(dmfv,NULL,&cellgeom,NULL);
668:   VecGetDM(cellgeom,&dmCell);
669:   VecGetArrayRead(cellgeom,&cgeom);
670:   DMGetGlobalVector(dmgrad,&grad);
671:   DMGetLocalVector(dmgrad,&locGrad);
672:   DMPlexGetHybridBounds(dmgrad,&cEndInterior,NULL,NULL,NULL);
673:   cEndInterior = (cEndInterior < 0) ? cEnd: cEndInterior;
674:   for (v = 0; v < nvecs; v++) {
675:     Vec               locX;
676:     PetscInt          c;
677:     PetscScalar       trueGrad[3][3] = {{0.}};
678:     const PetscScalar *gradArray;
679:     PetscReal         maxDiff, maxDiffGlob;

681:     DMGetLocalVector(dmfv,&locX);
682:     /* get the local projection of the rigid body mode */
683:     for (c = cStart; c < cEnd; c++) {
684:       PetscFVCellGeom *cg;
685:       PetscScalar     cx[3] = {0.,0.,0.};

687:       DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
688:       if (v < user->dim) {
689:         cx[v] = 1.;
690:       } else {
691:         PetscInt w = v - user->dim;

693:         cx[(w + 1) % user->dim] =  cg->centroid[(w + 2) % user->dim];
694:         cx[(w + 2) % user->dim] = -cg->centroid[(w + 1) % user->dim];
695:       }
696:       DMPlexVecSetClosure(dmfv,NULL,locX,c,cx,INSERT_ALL_VALUES);
697:     }
698:     /* TODO: this isn't in any header */
699:     DMPlexReconstructGradientsFVM(dmfv,locX,grad);
700:     DMGlobalToLocalBegin(dmgrad,grad,INSERT_VALUES,locGrad);
701:     DMGlobalToLocalEnd(dmgrad,grad,INSERT_VALUES,locGrad);
702:     VecGetArrayRead(locGrad,&gradArray);
703:     /* compare computed gradient to exact gradient */
704:     if (v >= user->dim) {
705:       PetscInt w = v - user->dim;

707:       trueGrad[(w + 1) % user->dim][(w + 2) % user->dim] =  1.;
708:       trueGrad[(w + 2) % user->dim][(w + 1) % user->dim] = -1.;
709:     }
710:     maxDiff = 0.;
711:     for (c = cStart; c < cEndInterior; c++) {
712:       PetscScalar *compGrad;
713:       PetscInt    i, j, k;
714:       PetscReal   FrobDiff = 0.;

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

718:       for (i = 0, k = 0; i < user->dim; i++) {
719:         for (j = 0; j < user->dim; j++, k++) {
720:           PetscScalar diff = compGrad[k] - trueGrad[i][j];
721:           FrobDiff += PetscRealPart(diff * PetscConj(diff));
722:         }
723:       }
724:       FrobDiff = PetscSqrtReal(FrobDiff);
725:       maxDiff  = PetscMax(maxDiff,FrobDiff);
726:     }
727:     MPI_Allreduce(&maxDiff,&maxDiffGlob,1,MPIU_REAL,MPIU_MAX,comm);
728:     allVecMaxDiff = PetscMax(allVecMaxDiff,maxDiffGlob);
729:     VecRestoreArrayRead(locGrad,&gradArray);
730:     DMRestoreLocalVector(dmfv,&locX);
731:   }
732:   if (allVecMaxDiff < fvTol) {
733:     PetscPrintf(PetscObjectComm((PetscObject)dm),"Finite volume gradient reconstruction: PASS\n");
734:   } else {
735:     PetscPrintf(PetscObjectComm((PetscObject)dm),"Finite volume gradient reconstruction: FAIL at tolerance %g with max difference %g\n",fvTol,allVecMaxDiff);
736:   }
737:   DMRestoreLocalVector(dmgrad,&locGrad);
738:   DMRestoreGlobalVector(dmgrad,&grad);
739:   VecRestoreArrayRead(cellgeom,&cgeom);
740:   DMDestroy(&dmfv);
741:   PetscFVDestroy(&fv);
742:   return(0);
743: }

745: static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *),
746:                                    PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *),
747:                                    void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
748: {
749:   Vec            u;
750:   PetscReal      n[3] = {1.0, 1.0, 1.0};

754:   DMGetGlobalVector(dm, &u);
755:   /* Project function into FE function space */
756:   DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u);
757:   VecViewFromOptions(u, NULL, "-projection_view");
758:   /* Compare approximation to exact in L_2 */
759:   DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error);
760:   DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer);
761:   DMRestoreGlobalVector(dm, &u);
762:   return(0);
763: }

765: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
766: {
767:   PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
768:   PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
769:   void            *exactCtxs[3];
770:   MPI_Comm         comm;
771:   PetscReal        error, errorDer, tol = PETSC_SMALL;
772:   PetscErrorCode   ierr;

775:   exactCtxs[0]       = user;
776:   exactCtxs[1]       = user;
777:   exactCtxs[2]       = user;
778:   user->constants[0] = 1.0;
779:   user->constants[1] = 2.0;
780:   user->constants[2] = 3.0;
781:   PetscObjectGetComm((PetscObject)dm, &comm);
782:   /* Setup functions to approximate */
783:   switch (order) {
784:   case 0:
785:     exactFuncs[0]    = constant;
786:     exactFuncDers[0] = constantDer;
787:     break;
788:   case 1:
789:     exactFuncs[0]    = linear;
790:     exactFuncDers[0] = linearDer;
791:     break;
792:   case 2:
793:     exactFuncs[0]    = quadratic;
794:     exactFuncDers[0] = quadraticDer;
795:     break;
796:   case 3:
797:     exactFuncs[0]    = cubic;
798:     exactFuncDers[0] = cubicDer;
799:     break;
800:   default:
801:     SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %d order %d", user->dim, order);
802:   }
803:   ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
804:   /* Report result */
805:   if (error > tol)    {PetscPrintf(comm, "Function tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol,(double) error);}
806:   else                {PetscPrintf(comm, "Function tests pass for order %D at tolerance %g\n", order, (double)tol);}
807:   if (errorDer > tol) {PetscPrintf(comm, "Function tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);}
808:   else                {PetscPrintf(comm, "Function tests pass for order %D derivatives at tolerance %g\n", order, (double)tol);}
809:   return(0);
810: }

812: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
813: {
814:   PetscErrorCode (*exactFuncs[1]) (PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
815:   PetscErrorCode (*exactFuncDers[1]) (PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx);
816:   PetscReal       n[3]         = {1.0, 1.0, 1.0};
817:   void           *exactCtxs[3];
818:   DM              rdm, idm, fdm;
819:   Mat             Interp;
820:   Vec             iu, fu, scaling;
821:   MPI_Comm        comm;
822:   PetscInt        dim  = user->dim;
823:   PetscReal       error, errorDer, tol = PETSC_SMALL;
824:   PetscBool       isPlex, isDA;
825:   PetscErrorCode  ierr;

828:   exactCtxs[0]       = user;
829:   exactCtxs[1]       = user;
830:   exactCtxs[2]       = user;
831:   user->constants[0] = 1.0;
832:   user->constants[1] = 2.0;
833:   user->constants[2] = 3.0;
834:   PetscObjectGetComm((PetscObject)dm,&comm);
835:   PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
836:   PetscObjectTypeCompare((PetscObject) dm, DMDA,   &isDA);
837:   DMRefine(dm, comm, &rdm);
838:   DMSetCoarseDM(rdm, dm);
839:   DMPlexSetRegularRefinement(rdm, user->convRefine);
840:   if (user->tree && isPlex) {
841:     DM refTree;
842:     DMPlexGetReferenceTree(dm,&refTree);
843:     DMPlexSetReferenceTree(rdm,refTree);
844:   }
845:   if (!user->simplex && !user->constraints) {DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
846:   SetupSection(rdm, user);
847:   /* Setup functions to approximate */
848:   switch (order) {
849:   case 0:
850:     exactFuncs[0]    = constant;
851:     exactFuncDers[0] = constantDer;
852:     break;
853:   case 1:
854:     exactFuncs[0]    = linear;
855:     exactFuncDers[0] = linearDer;
856:     break;
857:   case 2:
858:     exactFuncs[0]    = quadratic;
859:     exactFuncDers[0] = quadraticDer;
860:     break;
861:   case 3:
862:     exactFuncs[0]    = cubic;
863:     exactFuncDers[0] = cubicDer;
864:     break;
865:   default:
866:     SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %D order %D", dim, order);
867:   }
868:   idm  = checkRestrict ? rdm :  dm;
869:   fdm  = checkRestrict ?  dm : rdm;
870:   DMGetGlobalVector(idm, &iu);
871:   DMGetGlobalVector(fdm, &fu);
872:   DMSetApplicationContext(dm, user);
873:   DMSetApplicationContext(rdm, user);
874:   DMCreateInterpolation(dm, rdm, &Interp, &scaling);
875:   /* Project function into initial FE function space */
876:   DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu);
877:   /* Interpolate function into final FE function space */
878:   if (checkRestrict) {MatRestrict(Interp, iu, fu);VecPointwiseMult(fu, scaling, fu);}
879:   else               {MatInterpolate(Interp, iu, fu);}
880:   /* Compare approximation to exact in L_2 */
881:   DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error);
882:   DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer);
883:   /* Report result */
884:   if (error > tol)    {PetscPrintf(comm, "Interpolation tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol, (double)error);}
885:   else                {PetscPrintf(comm, "Interpolation tests pass for order %D at tolerance %g\n", order, (double)tol);}
886:   if (errorDer > tol) {PetscPrintf(comm, "Interpolation tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);}
887:   else                {PetscPrintf(comm, "Interpolation tests pass for order %D derivatives at tolerance %g\n", order, (double)tol);}
888:   DMRestoreGlobalVector(idm, &iu);
889:   DMRestoreGlobalVector(fdm, &fu);
890:   MatDestroy(&Interp);
891:   VecDestroy(&scaling);
892:   DMDestroy(&rdm);
893:   return(0);
894: }

896: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
897: {
898:   DM               odm = dm, rdm = NULL, cdm = NULL;
899:   PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {trig};
900:   PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx) = {trigDer};
901:   void            *exactCtxs[3];
902:   PetscInt         r, c, cStart, cEnd;
903:   PetscReal        errorOld, errorDerOld, error, errorDer, rel, len, lenOld;
904:   double           p;
905:   PetscErrorCode   ierr;

908:   if (!user->convergence) return(0);
909:   exactCtxs[0] = user;
910:   exactCtxs[1] = user;
911:   exactCtxs[2] = user;
912:   PetscObjectReference((PetscObject) odm);
913:   if (!user->convRefine) {
914:     for (r = 0; r < Nr; ++r) {
915:       DMRefine(odm, PetscObjectComm((PetscObject) dm), &rdm);
916:       DMDestroy(&odm);
917:       odm  = rdm;
918:     }
919:     SetupSection(odm, user);
920:   }
921:   ComputeError(odm, exactFuncs, exactFuncDers, exactCtxs, &errorOld, &errorDerOld, user);
922:   if (user->convRefine) {
923:     for (r = 0; r < Nr; ++r) {
924:       DMRefine(odm, PetscObjectComm((PetscObject) dm), &rdm);
925:       if (!user->simplex && user->useDA) {DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
926:       SetupSection(rdm, user);
927:       ComputeError(rdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
928:       p    = PetscLog2Real(errorOld/error);
929:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Function   convergence rate at refinement %D: %.2g\n", r, (double)p);
930:       p    = PetscLog2Real(errorDerOld/errorDer);
931:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Derivative convergence rate at refinement %D: %.2g\n", r, (double)p);
932:       DMDestroy(&odm);
933:       odm         = rdm;
934:       errorOld    = error;
935:       errorDerOld = errorDer;
936:     }
937:   } else {
938:     /* ComputeLongestEdge(dm, &lenOld); */
939:     DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd);
940:     lenOld = cEnd - cStart;
941:     for (c = 0; c < Nr; ++c) {
942:       DMCoarsen(odm, PetscObjectComm((PetscObject) dm), &cdm);
943:       if (!user->simplex && user->useDA) {DMDASetVertexCoordinates(cdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
944:       SetupSection(cdm, user);
945:       ComputeError(cdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
946:       /* ComputeLongestEdge(cdm, &len); */
947:       DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
948:       len  = cEnd - cStart;
949:       rel  = error/errorOld;
950:       p    = PetscLogReal(rel) / PetscLogReal(lenOld / len);
951:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Function   convergence rate at coarsening %D: %.2g\n", c, (double)p);
952:       rel  = errorDer/errorDerOld;
953:       p    = PetscLogReal(rel) / PetscLogReal(lenOld / len);
954:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Derivative convergence rate at coarsening %D: %.2g\n", c, (double)p);
955:       DMDestroy(&odm);
956:       odm         = cdm;
957:       errorOld    = error;
958:       errorDerOld = errorDer;
959:       lenOld      = len;
960:     }
961:   }
962:   DMDestroy(&odm);
963:   return(0);
964: }

966: int main(int argc, char **argv)
967: {
968:   DM             dm;
969:   AppCtx         user;                 /* user-defined work context */

972:   PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
973:   ProcessOptions(PETSC_COMM_WORLD, &user);
974:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
975:   PetscFECreateDefault(PETSC_COMM_WORLD, user.dim, user.numComponents, user.simplex, NULL, user.qorder, &user.fe);
976:   SetupSection(dm, &user);
977:   if (user.testFEjacobian) {TestFEJacobian(dm, &user);}
978:   if (user.testFVgrad) {TestFVGrad(dm, &user);}
979:   if (user.testInjector) {TestInjector(dm, &user);}
980:   CheckFunctions(dm, user.porder, &user);
981:   if (user.dim == 2 && user.simplex == PETSC_TRUE && user.tree == PETSC_FALSE) {
982:     CheckInterpolation(dm, PETSC_FALSE, user.porder, &user);
983:     CheckInterpolation(dm, PETSC_TRUE,  user.porder, &user);
984:   }
985:   CheckConvergence(dm, 3, &user);
986:   PetscFEDestroy(&user.fe);
987:   DMDestroy(&dm);
988:   PetscFinalize();
989:   return ierr;
990: }

992: /*TEST

994:   test:
995:     suffix: 1
996:     requires: triangle

998:   # 2D P_1 on a triangle
999:   test:
1000:     suffix: p1_2d_0
1001:     requires: triangle
1002:     args: -petscspace_degree 1 -qorder 1 -convergence
1003:   test:
1004:     suffix: p1_2d_1
1005:     requires: triangle
1006:     args: -petscspace_degree 1 -qorder 1 -porder 1
1007:   test:
1008:     suffix: p1_2d_2
1009:     requires: triangle
1010:     args: -petscspace_degree 1 -qorder 1 -porder 2
1011:   test:
1012:     suffix: p1_2d_3
1013:     requires: triangle pragmatic
1014:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1015:   test:
1016:     suffix: p1_2d_4
1017:     requires: triangle pragmatic
1018:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1019:   test:
1020:     suffix: p1_2d_5
1021:     requires: triangle pragmatic
1022:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0

1024:   # 3D P_1 on a tetrahedron
1025:   test:
1026:     suffix: p1_3d_0
1027:     requires: ctetgen
1028:     args: -dim 3 -petscspace_degree 1 -qorder 1 -convergence
1029:   test:
1030:     suffix: p1_3d_1
1031:     requires: ctetgen
1032:     args: -dim 3 -petscspace_degree 1 -qorder 1 -porder 1
1033:   test:
1034:     suffix: p1_3d_2
1035:     requires: ctetgen
1036:     args: -dim 3 -petscspace_degree 1 -qorder 1 -porder 2
1037:   test:
1038:     suffix: p1_3d_3
1039:     requires: ctetgen pragmatic
1040:     args: -dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1041:   test:
1042:     suffix: p1_3d_4
1043:     requires: ctetgen pragmatic
1044:     args: -dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1045:   test:
1046:     suffix: p1_3d_5
1047:     requires: ctetgen pragmatic
1048:     args: -dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0

1050:   # 2D P_2 on a triangle
1051:   test:
1052:     suffix: p2_2d_0
1053:     requires: triangle
1054:     args: -petscspace_degree 2 -qorder 2 -convergence
1055:   test:
1056:     suffix: p2_2d_1
1057:     requires: triangle
1058:     args: -petscspace_degree 2 -qorder 2 -porder 1
1059:   test:
1060:     suffix: p2_2d_2
1061:     requires: triangle
1062:     args: -petscspace_degree 2 -qorder 2 -porder 2
1063:   test:
1064:     suffix: p2_2d_3
1065:     requires: triangle pragmatic
1066:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1067:   test:
1068:     suffix: p2_2d_4
1069:     requires: triangle pragmatic
1070:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1071:   test:
1072:     suffix: p2_2d_5
1073:     requires: triangle pragmatic
1074:     args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0

1076:   # 3D P_2 on a tetrahedron
1077:   test:
1078:     suffix: p2_3d_0
1079:     requires: ctetgen
1080:     args: -dim 3 -petscspace_degree 2 -qorder 2 -convergence
1081:   test:
1082:     suffix: p2_3d_1
1083:     requires: ctetgen
1084:     args: -dim 3 -petscspace_degree 2 -qorder 2 -porder 1
1085:   test:
1086:     suffix: p2_3d_2
1087:     requires: ctetgen
1088:     args: -dim 3 -petscspace_degree 2 -qorder 2 -porder 2
1089:   test:
1090:     suffix: p2_3d_3
1091:     requires: ctetgen pragmatic
1092:     args: -dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1093:   test:
1094:     suffix: p2_3d_4
1095:     requires: ctetgen pragmatic
1096:     args: -dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1097:   test:
1098:     suffix: p2_3d_5
1099:     requires: ctetgen pragmatic
1100:     args: -dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0

1102:   # 2D Q_1 on a quadrilaterial DA
1103:   test:
1104:     suffix: q1_2d_da_0
1105:     requires: mpi_type_get_envelope broken
1106:     args: -simplex 0 -petscspace_degree 1 -qorder 1 -convergence
1107:   test:
1108:     suffix: q1_2d_da_1
1109:     requires: mpi_type_get_envelope broken
1110:     args: -simplex 0 -petscspace_degree 1 -qorder 1 -porder 1
1111:   test:
1112:     suffix: q1_2d_da_2
1113:     requires: mpi_type_get_envelope broken
1114:     args: -simplex 0 -petscspace_degree 1 -qorder 1 -porder 2

1116:   # 2D Q_1 on a quadrilaterial Plex
1117:   test:
1118:     suffix: q1_2d_plex_0
1119:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 1 -convergence
1120:   test:
1121:     suffix: q1_2d_plex_1
1122:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 1 -porder 1
1123:   test:
1124:     suffix: q1_2d_plex_2
1125:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 1 -porder 2
1126:   test:
1127:     suffix: q1_2d_plex_3
1128:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 1 -porder 1 -shear_coords
1129:   test:
1130:     suffix: q1_2d_plex_4
1131:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 1 -porder 2 -shear_coords
1132:   test:
1133:     suffix: q1_2d_plex_5
1134:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 0 -non_affine_coords -convergence
1135:   test:
1136:     suffix: q1_2d_plex_6
1137:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 1 -non_affine_coords -convergence
1138:   test:
1139:     suffix: q1_2d_plex_7
1140:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 2 -non_affine_coords -convergence

1142:   # 2D Q_2 on a quadrilaterial
1143:   test:
1144:     suffix: q2_2d_plex_0
1145:     requires: mpi_type_get_envelope
1146:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1147:   test:
1148:     suffix: q2_2d_plex_1
1149:     requires: mpi_type_get_envelope
1150:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1151:   test:
1152:     suffix: q2_2d_plex_2
1153:     requires: mpi_type_get_envelope
1154:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1155:   test:
1156:     suffix: q2_2d_plex_3
1157:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1158:   test:
1159:     suffix: q2_2d_plex_4
1160:     requires: mpi_type_get_envelope
1161:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1162:   test:
1163:     suffix: q2_2d_plex_5
1164:     requires: mpi_type_get_envelope
1165:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 0 -non_affine_coords -convergence
1166:   test:
1167:     suffix: q2_2d_plex_6
1168:     requires: mpi_type_get_envelope
1169:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 1 -non_affine_coords -convergence
1170:   test:
1171:     suffix: q2_2d_plex_7
1172:     requires: mpi_type_get_envelope
1173:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 2 -non_affine_coords -convergence


1176:   # 2D P_3 on a triangle
1177:   test:
1178:     suffix: p3_2d_0
1179:     requires: triangle !single
1180:     args: -petscspace_degree 3 -qorder 3 -convergence
1181:   test:
1182:     suffix: p3_2d_1
1183:     requires: triangle !single
1184:     args: -petscspace_degree 3 -qorder 3 -porder 1
1185:   test:
1186:     suffix: p3_2d_2
1187:     requires: triangle !single
1188:     args: -petscspace_degree 3 -qorder 3 -porder 2
1189:   test:
1190:     suffix: p3_2d_3
1191:     requires: triangle !single
1192:     args: -petscspace_degree 3 -qorder 3 -porder 3
1193:   test:
1194:     suffix: p3_2d_4
1195:     requires: triangle pragmatic
1196:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -convergence -conv_refine 0
1197:   test:
1198:     suffix: p3_2d_5
1199:     requires: triangle pragmatic
1200:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 1 -conv_refine 0
1201:   test:
1202:     suffix: p3_2d_6
1203:     requires: triangle pragmatic
1204:     args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 3 -conv_refine 0

1206:   # 2D Q_3 on a quadrilaterial
1207:   test:
1208:     suffix: q3_2d_0
1209:     requires: mpi_type_get_envelope !single
1210:     args: -use_da 0 -simplex 0 -petscspace_degree 3 -qorder 3 -convergence
1211:   test:
1212:     suffix: q3_2d_1
1213:     requires: mpi_type_get_envelope !single
1214:     args: -use_da 0 -simplex 0 -petscspace_degree 3 -qorder 3 -porder 1
1215:   test:
1216:     suffix: q3_2d_2
1217:     requires: mpi_type_get_envelope !single
1218:     args: -use_da 0 -simplex 0 -petscspace_degree 3 -qorder 3 -porder 2
1219:   test:
1220:     suffix: q3_2d_3
1221:     requires: mpi_type_get_envelope !single
1222:     args: -use_da 0 -simplex 0 -petscspace_degree 3 -qorder 3 -porder 3

1224:   # 2D P_1disc on a triangle/quadrilateral
1225:   test:
1226:     suffix: p1d_2d_0
1227:     requires: triangle
1228:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1229:   test:
1230:     suffix: p1d_2d_1
1231:     requires: triangle
1232:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1233:   test:
1234:     suffix: p1d_2d_2
1235:     requires: triangle
1236:     args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1237:   test:
1238:     suffix: p1d_2d_3
1239:     requires: triangle
1240:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1241:     filter: sed  -e "s/convergence rate at refinement 0: 2/convergence rate at refinement 0: 1.9/g"
1242:   test:
1243:     suffix: p1d_2d_4
1244:     requires: triangle
1245:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1246:   test:
1247:     suffix: p1d_2d_5
1248:     requires: triangle
1249:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2

1251:   # Test high order quadrature
1252:   test:
1253:     suffix: p1_quad_2
1254:     requires: triangle
1255:     args: -petscspace_degree 1 -qorder 2 -porder 1
1256:   test:
1257:     suffix: p1_quad_5
1258:     requires: triangle
1259:     args: -petscspace_degree 1 -qorder 5 -porder 1
1260:   test:
1261:     suffix: p2_quad_3
1262:     requires: triangle
1263:     args: -petscspace_degree 2 -qorder 3 -porder 2
1264:   test:
1265:     suffix: p2_quad_5
1266:     requires: triangle
1267:     args: -petscspace_degree 2 -qorder 5 -porder 2
1268:   test:
1269:     suffix: q1_quad_2
1270:     requires: mpi_type_get_envelope
1271:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 2 -porder 1
1272:   test:
1273:     suffix: q1_quad_5
1274:     requires: mpi_type_get_envelope
1275:     args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 5 -porder 1
1276:   test:
1277:     suffix: q2_quad_3
1278:     requires: mpi_type_get_envelope
1279:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 3 -porder 1
1280:   test:
1281:     suffix: q2_quad_5
1282:     requires: mpi_type_get_envelope
1283:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 5 -porder 1


1286:   # Nonconforming tests
1287:   test:
1288:     suffix: constraints
1289:     args: -simplex 0 -petscspace_type tensor -petscspace_degree 1 -qorder 0 -constraints
1290:   test:
1291:     suffix: nonconforming_tensor_2
1292:     nsize: 4
1293:     args: -test_fe_jacobian -test_injector -petscpartitioner_type simple -tree -simplex 0 -dim 2 -dm_plex_max_projection_height 1 -petscspace_type tensor -petscspace_degree 2 -qorder 2 -dm_view ascii::ASCII_INFO_DETAIL
1294:   test:
1295:     suffix: nonconforming_tensor_3
1296:     nsize: 4
1297:     args: -test_fe_jacobian -petscpartitioner_type simple -tree -simplex 0 -dim 3 -dm_plex_max_projection_height 2 -petscspace_type tensor -petscspace_degree 1 -qorder 1 -dm_view ascii::ASCII_INFO_DETAIL
1298:   test:
1299:     suffix: nonconforming_tensor_2_fv
1300:     nsize: 4
1301:     args: -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -simplex 0 -dim 2 -num_comp 2
1302:   test:
1303:     suffix: nonconforming_tensor_3_fv
1304:     nsize: 4
1305:     args: -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -simplex 0 -dim 3 -num_comp 3
1306:   test:
1307:     suffix: nonconforming_tensor_2_hi
1308:     requires: !single
1309:     nsize: 4
1310:     args: -test_fe_jacobian -petscpartitioner_type simple -tree -simplex 0 -dim 2 -dm_plex_max_projection_height 1 -petscspace_type tensor -petscspace_degree 4 -qorder 4
1311:   test:
1312:     suffix: nonconforming_tensor_3_hi
1313:     requires: !single skip
1314:     nsize: 4
1315:     args: -test_fe_jacobian -petscpartitioner_type simple -tree -simplex 0 -dim 3 -dm_plex_max_projection_height 2 -petscspace_type tensor -petscspace_degree 4 -qorder 4
1316:   test:
1317:     suffix: nonconforming_simplex_2
1318:     requires: triangle
1319:     nsize: 4
1320:     args: -test_fe_jacobian -test_injector -petscpartitioner_type simple -tree -simplex 1 -dim 2 -dm_plex_max_projection_height 1 -petscspace_degree 2 -qorder 2 -dm_view ascii::ASCII_INFO_DETAIL
1321:   test:
1322:     suffix: nonconforming_simplex_2_hi
1323:     requires: triangle !single
1324:     nsize: 4
1325:     args: -test_fe_jacobian -petscpartitioner_type simple -tree -simplex 1 -dim 2 -dm_plex_max_projection_height 1 -petscspace_degree 4 -qorder 4
1326:   test:
1327:     suffix: nonconforming_simplex_2_fv
1328:     requires: triangle
1329:     nsize: 4
1330:     args: -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -simplex 1 -dim 2 -num_comp 2
1331:   test:
1332:     suffix: nonconforming_simplex_3
1333:     requires: ctetgen
1334:     nsize: 4
1335:     args: -test_fe_jacobian -test_injector -petscpartitioner_type simple -tree -simplex 1 -dim 3 -dm_plex_max_projection_height 2 -petscspace_degree 2 -qorder 2 -dm_view ascii::ASCII_INFO_DETAIL
1336:   test:
1337:     suffix: nonconforming_simplex_3_hi
1338:     requires: ctetgen skip
1339:     nsize: 4
1340:     args: -test_fe_jacobian -petscpartitioner_type simple -tree -simplex 1 -dim 3 -dm_plex_max_projection_height 2 -petscspace_degree 4 -qorder 4
1341:   test:
1342:     suffix: nonconforming_simplex_3_fv
1343:     requires: ctetgen
1344:     nsize: 4
1345:     args: -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -simplex 1 -dim 3 -num_comp 3

1347: TEST*/

1349: /*
1350:    # 2D Q_2 on a quadrilaterial Plex
1351:   test:
1352:     suffix: q2_2d_plex_0
1353:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1354:   test:
1355:     suffix: q2_2d_plex_1
1356:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1357:   test:
1358:     suffix: q2_2d_plex_2
1359:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1360:   test:
1361:     suffix: q2_2d_plex_3
1362:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1363:   test:
1364:     suffix: q2_2d_plex_4
1365:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1366:   test:
1367:     suffix: q2_2d_plex_5
1368:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 0 -non_affine_coords
1369:   test:
1370:     suffix: q2_2d_plex_6
1371:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 1 -non_affine_coords
1372:   test:
1373:     suffix: q2_2d_plex_7
1374:     args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 2 -non_affine_coords

1376:   test:
1377:     suffix: p1d_2d_6
1378:     requires: pragmatic
1379:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1380:   test:
1381:     suffix: p1d_2d_7
1382:     requires: pragmatic
1383:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1384:   test:
1385:     suffix: p1d_2d_8
1386:     requires: pragmatic
1387:     args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1388: */