Actual source code: ex3.c
1: static char help[] = "Check that a DM can accurately represent and interpolate functions of a given polynomial order\n\n";
3: #include <petscdmplex.h>
4: #include <petscdm.h>
5: #include <petscdmda.h>
6: #include <petscfe.h>
7: #include <petscds.h>
8: #include <petscksp.h>
9: #include <petscsnes.h>
10: #include <petscsf.h>
12: typedef struct {
13: /* Domain and mesh definition */
14: PetscBool useDA; /* Flag DMDA tensor product mesh */
15: PetscBool shearCoords; /* Flag for shear transform */
16: PetscBool nonaffineCoords; /* Flag for non-affine transform */
17: /* Element definition */
18: PetscInt qorder; /* Order of the quadrature */
19: PetscInt numComponents; /* Number of field components */
20: PetscFE fe; /* The finite element */
21: /* Testing space */
22: PetscInt porder; /* Order of polynomials to test */
23: PetscBool convergence; /* Test for order of convergence */
24: PetscBool convRefine; /* Test for convergence using refinement, otherwise use coarsening */
25: PetscBool constraints; /* Test local constraints */
26: PetscBool tree; /* Test tree routines */
27: PetscBool testFEjacobian; /* Test finite element Jacobian assembly */
28: PetscBool testFVgrad; /* Test finite difference gradient routine */
29: PetscBool testInjector; /* Test finite element injection routines */
30: PetscInt treeCell; /* Cell to refine in tree test */
31: PetscReal constants[3]; /* Constant values for each dimension */
32: } AppCtx;
34: /* u = 1 */
35: PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
36: {
37: AppCtx *user = (AppCtx *)ctx;
38: PetscInt d;
39: for (d = 0; d < dim; ++d) u[d] = user->constants[d];
40: return PETSC_SUCCESS;
41: }
42: PetscErrorCode constantDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
43: {
44: PetscInt d;
45: for (d = 0; d < dim; ++d) u[d] = 0.0;
46: return PETSC_SUCCESS;
47: }
49: /* u = x */
50: PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
51: {
52: PetscInt d;
53: for (d = 0; d < dim; ++d) u[d] = coords[d];
54: return PETSC_SUCCESS;
55: }
56: PetscErrorCode linearDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
57: {
58: PetscInt d, e;
59: for (d = 0; d < dim; ++d) {
60: u[d] = 0.0;
61: for (e = 0; e < dim; ++e) u[d] += (d == e ? 1.0 : 0.0) * n[e];
62: }
63: return PETSC_SUCCESS;
64: }
66: /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */
67: PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
68: {
69: if (dim > 2) {
70: u[0] = coords[0] * coords[1];
71: u[1] = coords[1] * coords[2];
72: u[2] = coords[2] * coords[0];
73: } else if (dim > 1) {
74: u[0] = coords[0] * coords[0];
75: u[1] = coords[0] * coords[1];
76: } else if (dim > 0) {
77: u[0] = coords[0] * coords[0];
78: }
79: return PETSC_SUCCESS;
80: }
81: PetscErrorCode quadraticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
82: {
83: if (dim > 2) {
84: u[0] = coords[1] * n[0] + coords[0] * n[1];
85: u[1] = coords[2] * n[1] + coords[1] * n[2];
86: u[2] = coords[2] * n[0] + coords[0] * n[2];
87: } else if (dim > 1) {
88: u[0] = 2.0 * coords[0] * n[0];
89: u[1] = coords[1] * n[0] + coords[0] * n[1];
90: } else if (dim > 0) {
91: u[0] = 2.0 * coords[0] * n[0];
92: }
93: return PETSC_SUCCESS;
94: }
96: /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
97: PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
98: {
99: if (dim > 2) {
100: u[0] = coords[0] * coords[0] * coords[1];
101: u[1] = coords[1] * coords[1] * coords[2];
102: u[2] = coords[2] * coords[2] * coords[0];
103: } else if (dim > 1) {
104: u[0] = coords[0] * coords[0] * coords[0];
105: u[1] = coords[0] * coords[0] * coords[1];
106: } else if (dim > 0) {
107: u[0] = coords[0] * coords[0] * coords[0];
108: }
109: return PETSC_SUCCESS;
110: }
111: PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
112: {
113: if (dim > 2) {
114: u[0] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1];
115: u[1] = 2.0 * coords[1] * coords[2] * n[1] + coords[1] * coords[1] * n[2];
116: u[2] = 2.0 * coords[2] * coords[0] * n[2] + coords[2] * coords[2] * n[0];
117: } else if (dim > 1) {
118: u[0] = 3.0 * coords[0] * coords[0] * n[0];
119: u[1] = 2.0 * coords[0] * coords[1] * n[0] + coords[0] * coords[0] * n[1];
120: } else if (dim > 0) {
121: u[0] = 3.0 * coords[0] * coords[0] * n[0];
122: }
123: return PETSC_SUCCESS;
124: }
126: /* u = tanh(x) */
127: PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
128: {
129: PetscInt d;
130: for (d = 0; d < dim; ++d) u[d] = PetscTanhReal(coords[d] - 0.5);
131: return PETSC_SUCCESS;
132: }
133: PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
134: {
135: PetscInt d;
136: for (d = 0; d < dim; ++d) u[d] = 1.0 / PetscSqr(PetscCoshReal(coords[d] - 0.5)) * n[d];
137: return PETSC_SUCCESS;
138: }
140: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
141: {
142: PetscInt n = 3;
144: PetscFunctionBeginUser;
145: options->useDA = PETSC_FALSE;
146: options->shearCoords = PETSC_FALSE;
147: options->nonaffineCoords = PETSC_FALSE;
148: options->qorder = 0;
149: options->numComponents = PETSC_DEFAULT;
150: options->porder = 0;
151: options->convergence = PETSC_FALSE;
152: options->convRefine = PETSC_TRUE;
153: options->constraints = PETSC_FALSE;
154: options->tree = PETSC_FALSE;
155: options->treeCell = 0;
156: options->testFEjacobian = PETSC_FALSE;
157: options->testFVgrad = PETSC_FALSE;
158: options->testInjector = PETSC_FALSE;
159: options->constants[0] = 1.0;
160: options->constants[1] = 2.0;
161: options->constants[2] = 3.0;
163: PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
164: PetscCall(PetscOptionsBool("-use_da", "Flag for DMDA mesh", "ex3.c", options->useDA, &options->useDA, NULL));
165: PetscCall(PetscOptionsBool("-shear_coords", "Transform coordinates with a shear", "ex3.c", options->shearCoords, &options->shearCoords, NULL));
166: PetscCall(PetscOptionsBool("-non_affine_coords", "Transform coordinates with a non-affine transform", "ex3.c", options->nonaffineCoords, &options->nonaffineCoords, NULL));
167: PetscCall(PetscOptionsBoundedInt("-qorder", "The quadrature order", "ex3.c", options->qorder, &options->qorder, NULL, 0));
168: PetscCall(PetscOptionsBoundedInt("-num_comp", "The number of field components", "ex3.c", options->numComponents, &options->numComponents, NULL, PETSC_DEFAULT));
169: PetscCall(PetscOptionsBoundedInt("-porder", "The order of polynomials to test", "ex3.c", options->porder, &options->porder, NULL, 0));
170: PetscCall(PetscOptionsBool("-convergence", "Check the convergence rate", "ex3.c", options->convergence, &options->convergence, NULL));
171: PetscCall(PetscOptionsBool("-conv_refine", "Use refinement for the convergence rate", "ex3.c", options->convRefine, &options->convRefine, NULL));
172: PetscCall(PetscOptionsBool("-constraints", "Test local constraints (serial only)", "ex3.c", options->constraints, &options->constraints, NULL));
173: PetscCall(PetscOptionsBool("-tree", "Test tree routines", "ex3.c", options->tree, &options->tree, NULL));
174: PetscCall(PetscOptionsBoundedInt("-tree_cell", "cell to refine in tree test", "ex3.c", options->treeCell, &options->treeCell, NULL, 0));
175: PetscCall(PetscOptionsBool("-test_fe_jacobian", "Test finite element Jacobian assembly", "ex3.c", options->testFEjacobian, &options->testFEjacobian, NULL));
176: PetscCall(PetscOptionsBool("-test_fv_grad", "Test finite volume gradient reconstruction", "ex3.c", options->testFVgrad, &options->testFVgrad, NULL));
177: PetscCall(PetscOptionsBool("-test_injector", "Test finite element injection", "ex3.c", options->testInjector, &options->testInjector, NULL));
178: PetscCall(PetscOptionsRealArray("-constants", "Set the constant values", "ex3.c", options->constants, &n, NULL));
179: PetscOptionsEnd();
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: static PetscErrorCode TransformCoordinates(DM dm, AppCtx *user)
184: {
185: PetscSection coordSection;
186: Vec coordinates;
187: PetscScalar *coords;
188: PetscInt vStart, vEnd, v;
190: PetscFunctionBeginUser;
191: if (user->nonaffineCoords) {
192: /* x' = r^(1/p) (x/r), y' = r^(1/p) (y/r), z' = z */
193: PetscCall(DMGetCoordinateSection(dm, &coordSection));
194: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
195: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
196: PetscCall(VecGetArray(coordinates, &coords));
197: for (v = vStart; v < vEnd; ++v) {
198: PetscInt dof, off;
199: PetscReal p = 4.0, r;
201: PetscCall(PetscSectionGetDof(coordSection, v, &dof));
202: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
203: switch (dof) {
204: case 2:
205: r = PetscSqr(PetscRealPart(coords[off + 0])) + PetscSqr(PetscRealPart(coords[off + 1]));
206: coords[off + 0] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 0];
207: coords[off + 1] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 1];
208: break;
209: case 3:
210: r = PetscSqr(PetscRealPart(coords[off + 0])) + PetscSqr(PetscRealPart(coords[off + 1]));
211: coords[off + 0] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 0];
212: coords[off + 1] = r == 0.0 ? 0.0 : PetscPowReal(r, (1 - p) / (2 * p)) * coords[off + 1];
213: coords[off + 2] = coords[off + 2];
214: break;
215: }
216: }
217: PetscCall(VecRestoreArray(coordinates, &coords));
218: }
219: if (user->shearCoords) {
220: /* x' = x + m y + m z, y' = y + m z, z' = z */
221: PetscCall(DMGetCoordinateSection(dm, &coordSection));
222: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
223: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
224: PetscCall(VecGetArray(coordinates, &coords));
225: for (v = vStart; v < vEnd; ++v) {
226: PetscInt dof, off;
227: PetscReal m = 1.0;
229: PetscCall(PetscSectionGetDof(coordSection, v, &dof));
230: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
231: switch (dof) {
232: case 2:
233: coords[off + 0] = coords[off + 0] + m * coords[off + 1];
234: coords[off + 1] = coords[off + 1];
235: break;
236: case 3:
237: coords[off + 0] = coords[off + 0] + m * coords[off + 1] + m * coords[off + 2];
238: coords[off + 1] = coords[off + 1] + m * coords[off + 2];
239: coords[off + 2] = coords[off + 2];
240: break;
241: }
242: }
243: PetscCall(VecRestoreArray(coordinates, &coords));
244: }
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
249: {
250: PetscInt dim = 2;
251: PetscBool simplex;
253: PetscFunctionBeginUser;
254: if (user->useDA) {
255: switch (dim) {
256: case 2:
257: PetscCall(DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 2, 2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, dm));
258: PetscCall(DMSetFromOptions(*dm));
259: PetscCall(DMSetUp(*dm));
260: PetscCall(DMDASetVertexCoordinates(*dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
261: break;
262: default:
263: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot create structured mesh of dimension %" PetscInt_FMT, dim);
264: }
265: PetscCall(PetscObjectSetName((PetscObject)*dm, "Hexahedral Mesh"));
266: } else {
267: PetscCall(DMCreate(comm, dm));
268: PetscCall(DMSetType(*dm, DMPLEX));
269: PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
270: PetscCall(DMSetFromOptions(*dm));
272: PetscCall(DMGetDimension(*dm, &dim));
273: PetscCall(DMPlexIsSimplex(*dm, &simplex));
274: PetscCallMPI(MPI_Bcast(&simplex, 1, MPIU_BOOL, 0, comm));
275: if (user->tree) {
276: DM refTree, ncdm = NULL;
278: PetscCall(DMPlexCreateDefaultReferenceTree(comm, dim, simplex, &refTree));
279: PetscCall(DMViewFromOptions(refTree, NULL, "-reftree_dm_view"));
280: PetscCall(DMPlexSetReferenceTree(*dm, refTree));
281: PetscCall(DMDestroy(&refTree));
282: PetscCall(DMPlexTreeRefineCell(*dm, user->treeCell, &ncdm));
283: if (ncdm) {
284: PetscCall(DMDestroy(dm));
285: *dm = ncdm;
286: PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_FALSE));
287: }
288: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "tree_"));
289: PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
290: PetscCall(DMSetFromOptions(*dm));
291: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
292: } else {
293: PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_TRUE));
294: }
295: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "dist_"));
296: PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE));
297: PetscCall(DMSetFromOptions(*dm));
298: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, NULL));
299: if (simplex) PetscCall(PetscObjectSetName((PetscObject)*dm, "Simplicial Mesh"));
300: else PetscCall(PetscObjectSetName((PetscObject)*dm, "Hexahedral Mesh"));
301: }
302: PetscCall(DMSetFromOptions(*dm));
303: PetscCall(TransformCoordinates(*dm, user));
304: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: static void simple_mass(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
309: {
310: PetscInt d, e;
311: for (d = 0, e = 0; d < dim; d++, e += dim + 1) g0[e] = 1.;
312: }
314: /* < \nabla v, 1/2(\nabla u + {\nabla u}^T) > */
315: static void symmetric_gradient_inner_product(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar C[])
316: {
317: PetscInt compI, compJ, d, e;
319: for (compI = 0; compI < dim; ++compI) {
320: for (compJ = 0; compJ < dim; ++compJ) {
321: for (d = 0; d < dim; ++d) {
322: for (e = 0; e < dim; e++) {
323: if (d == e && d == compI && d == compJ) {
324: C[((compI * dim + compJ) * dim + d) * dim + e] = 1.0;
325: } else if ((d == compJ && e == compI) || (d == e && compI == compJ)) {
326: C[((compI * dim + compJ) * dim + d) * dim + e] = 0.5;
327: } else {
328: C[((compI * dim + compJ) * dim + d) * dim + e] = 0.0;
329: }
330: }
331: }
332: }
333: }
334: }
336: static PetscErrorCode SetupSection(DM dm, AppCtx *user)
337: {
338: PetscFunctionBeginUser;
339: if (user->constraints) {
340: /* test local constraints */
341: DM coordDM;
342: PetscInt fStart, fEnd, f, vStart, vEnd, v;
343: PetscInt edgesx = 2, vertsx;
344: PetscInt edgesy = 2, vertsy;
345: PetscMPIInt size;
346: PetscInt numConst;
347: PetscSection aSec;
348: PetscInt *anchors;
349: PetscInt offset;
350: IS aIS;
351: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
353: PetscCallMPI(MPI_Comm_size(comm, &size));
354: PetscCheck(size <= 1, comm, PETSC_ERR_SUP, "Local constraint test can only be performed in serial");
356: /* we are going to test constraints by using them to enforce periodicity
357: * in one direction, and comparing to the existing method of enforcing
358: * periodicity */
360: /* first create the coordinate section so that it does not clone the
361: * constraints */
362: PetscCall(DMGetCoordinateDM(dm, &coordDM));
364: /* create the constrained-to-anchor section */
365: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
366: PetscCall(DMPlexGetDepthStratum(dm, 1, &fStart, &fEnd));
367: PetscCall(PetscSectionCreate(PETSC_COMM_SELF, &aSec));
368: PetscCall(PetscSectionSetChart(aSec, PetscMin(fStart, vStart), PetscMax(fEnd, vEnd)));
370: /* define the constraints */
371: PetscCall(PetscOptionsGetInt(NULL, NULL, "-da_grid_x", &edgesx, NULL));
372: PetscCall(PetscOptionsGetInt(NULL, NULL, "-da_grid_y", &edgesy, NULL));
373: vertsx = edgesx + 1;
374: vertsy = edgesy + 1;
375: numConst = vertsy + edgesy;
376: PetscCall(PetscMalloc1(numConst, &anchors));
377: offset = 0;
378: for (v = vStart + edgesx; v < vEnd; v += vertsx) {
379: PetscCall(PetscSectionSetDof(aSec, v, 1));
380: anchors[offset++] = v - edgesx;
381: }
382: for (f = fStart + edgesx * vertsy + edgesx * edgesy; f < fEnd; f++) {
383: PetscCall(PetscSectionSetDof(aSec, f, 1));
384: anchors[offset++] = f - edgesx * edgesy;
385: }
386: PetscCall(PetscSectionSetUp(aSec));
387: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numConst, anchors, PETSC_OWN_POINTER, &aIS));
389: PetscCall(DMPlexSetAnchors(dm, aSec, aIS));
390: PetscCall(PetscSectionDestroy(&aSec));
391: PetscCall(ISDestroy(&aIS));
392: }
393: PetscCall(DMSetNumFields(dm, 1));
394: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)user->fe));
395: PetscCall(DMCreateDS(dm));
396: if (user->constraints) {
397: /* test getting local constraint matrix that matches section */
398: PetscSection aSec;
399: IS aIS;
401: PetscCall(DMPlexGetAnchors(dm, &aSec, &aIS));
402: if (aSec) {
403: PetscDS ds;
404: PetscSection cSec, section;
405: PetscInt cStart, cEnd, c, numComp;
406: Mat cMat, mass;
407: Vec local;
408: const PetscInt *anchors;
410: PetscCall(DMGetLocalSection(dm, §ion));
411: /* this creates the matrix and preallocates the matrix structure: we
412: * just have to fill in the values */
413: PetscCall(DMGetDefaultConstraints(dm, &cSec, &cMat, NULL));
414: PetscCall(PetscSectionGetChart(cSec, &cStart, &cEnd));
415: PetscCall(ISGetIndices(aIS, &anchors));
416: PetscCall(PetscFEGetNumComponents(user->fe, &numComp));
417: for (c = cStart; c < cEnd; c++) {
418: PetscInt cDof;
420: /* is this point constrained? (does it have an anchor?) */
421: PetscCall(PetscSectionGetDof(aSec, c, &cDof));
422: if (cDof) {
423: PetscInt cOff, a, aDof, aOff, j;
424: PetscCheck(cDof == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Found %" PetscInt_FMT " anchor points: should be just one", cDof);
426: /* find the anchor point */
427: PetscCall(PetscSectionGetOffset(aSec, c, &cOff));
428: a = anchors[cOff];
430: /* find the constrained dofs (row in constraint matrix) */
431: PetscCall(PetscSectionGetDof(cSec, c, &cDof));
432: PetscCall(PetscSectionGetOffset(cSec, c, &cOff));
434: /* find the anchor dofs (column in constraint matrix) */
435: PetscCall(PetscSectionGetDof(section, a, &aDof));
436: PetscCall(PetscSectionGetOffset(section, a, &aOff));
438: PetscCheck(cDof == aDof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point and anchor have different number of dofs: %" PetscInt_FMT ", %" PetscInt_FMT, cDof, aDof);
439: PetscCheck(cDof % numComp == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Point dofs not divisible by field components: %" PetscInt_FMT ", %" PetscInt_FMT, cDof, numComp);
441: /* put in a simple equality constraint */
442: for (j = 0; j < cDof; j++) PetscCall(MatSetValue(cMat, cOff + j, aOff + j, 1., INSERT_VALUES));
443: }
444: }
445: PetscCall(MatAssemblyBegin(cMat, MAT_FINAL_ASSEMBLY));
446: PetscCall(MatAssemblyEnd(cMat, MAT_FINAL_ASSEMBLY));
447: PetscCall(ISRestoreIndices(aIS, &anchors));
449: /* Now that we have constructed the constraint matrix, any FE matrix
450: * that we construct will apply the constraints during construction */
452: PetscCall(DMCreateMatrix(dm, &mass));
453: /* get a dummy local variable to serve as the solution */
454: PetscCall(DMGetLocalVector(dm, &local));
455: PetscCall(DMGetDS(dm, &ds));
456: /* set the jacobian to be the mass matrix */
457: PetscCall(PetscDSSetJacobian(ds, 0, 0, simple_mass, NULL, NULL, NULL));
458: /* build the mass matrix */
459: PetscCall(DMPlexSNESComputeJacobianFEM(dm, local, mass, mass, NULL));
460: PetscCall(MatView(mass, PETSC_VIEWER_STDOUT_WORLD));
461: PetscCall(MatDestroy(&mass));
462: PetscCall(DMRestoreLocalVector(dm, &local));
463: }
464: }
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: static PetscErrorCode TestFEJacobian(DM dm, AppCtx *user)
469: {
470: PetscFunctionBeginUser;
471: if (!user->useDA) {
472: Vec local;
473: const Vec *vecs;
474: Mat E;
475: MatNullSpace sp;
476: PetscBool isNullSpace, hasConst;
477: PetscInt dim, n, i;
478: Vec res = NULL, localX, localRes;
479: PetscDS ds;
481: PetscCall(DMGetDimension(dm, &dim));
482: PetscCheck(user->numComponents == dim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "The number of components %" PetscInt_FMT " must be equal to the dimension %" PetscInt_FMT " for this test", user->numComponents, dim);
483: PetscCall(DMGetDS(dm, &ds));
484: PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, symmetric_gradient_inner_product));
485: PetscCall(DMCreateMatrix(dm, &E));
486: PetscCall(DMGetLocalVector(dm, &local));
487: PetscCall(DMPlexSNESComputeJacobianFEM(dm, local, E, E, NULL));
488: PetscCall(DMPlexCreateRigidBody(dm, 0, &sp));
489: PetscCall(MatNullSpaceGetVecs(sp, &hasConst, &n, &vecs));
490: if (n) PetscCall(VecDuplicate(vecs[0], &res));
491: PetscCall(DMCreateLocalVector(dm, &localX));
492: PetscCall(DMCreateLocalVector(dm, &localRes));
493: for (i = 0; i < n; i++) { /* also test via matrix-free Jacobian application */
494: PetscReal resNorm;
496: PetscCall(VecSet(localRes, 0.));
497: PetscCall(VecSet(localX, 0.));
498: PetscCall(VecSet(local, 0.));
499: PetscCall(VecSet(res, 0.));
500: PetscCall(DMGlobalToLocalBegin(dm, vecs[i], INSERT_VALUES, localX));
501: PetscCall(DMGlobalToLocalEnd(dm, vecs[i], INSERT_VALUES, localX));
502: PetscCall(DMSNESComputeJacobianAction(dm, local, localX, localRes, NULL));
503: PetscCall(DMLocalToGlobalBegin(dm, localRes, ADD_VALUES, res));
504: PetscCall(DMLocalToGlobalEnd(dm, localRes, ADD_VALUES, res));
505: PetscCall(VecNorm(res, NORM_2, &resNorm));
506: if (resNorm > PETSC_SMALL) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Symmetric gradient action null space vector %" PetscInt_FMT " residual: %E\n", i, (double)resNorm));
507: }
508: PetscCall(VecDestroy(&localRes));
509: PetscCall(VecDestroy(&localX));
510: PetscCall(VecDestroy(&res));
511: PetscCall(MatNullSpaceTest(sp, E, &isNullSpace));
512: if (isNullSpace) {
513: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Symmetric gradient null space: PASS\n"));
514: } else {
515: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Symmetric gradient null space: FAIL\n"));
516: }
517: PetscCall(MatNullSpaceDestroy(&sp));
518: PetscCall(MatDestroy(&E));
519: PetscCall(DMRestoreLocalVector(dm, &local));
520: }
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: static PetscErrorCode TestInjector(DM dm, AppCtx *user)
525: {
526: DM refTree;
527: PetscMPIInt rank;
529: PetscFunctionBegin;
530: PetscCall(DMPlexGetReferenceTree(dm, &refTree));
531: if (refTree) {
532: Mat inj;
534: PetscCall(DMPlexComputeInjectorReferenceTree(refTree, &inj));
535: PetscCall(PetscObjectSetName((PetscObject)inj, "Reference Tree Injector"));
536: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
537: if (rank == 0) PetscCall(MatView(inj, PETSC_VIEWER_STDOUT_SELF));
538: PetscCall(MatDestroy(&inj));
539: }
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: static PetscErrorCode TestFVGrad(DM dm, AppCtx *user)
544: {
545: MPI_Comm comm;
546: DM dmRedist, dmfv, dmgrad, dmCell, refTree;
547: PetscFV fv;
548: PetscInt dim, nvecs, v, cStart, cEnd, cEndInterior;
549: PetscMPIInt size;
550: Vec cellgeom, grad, locGrad;
551: const PetscScalar *cgeom;
552: PetscReal allVecMaxDiff = 0., fvTol = 100. * PETSC_MACHINE_EPSILON;
554: PetscFunctionBeginUser;
555: comm = PetscObjectComm((PetscObject)dm);
556: PetscCall(DMGetDimension(dm, &dim));
557: /* duplicate DM, give dup. a FV discretization */
558: PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
559: PetscCallMPI(MPI_Comm_size(comm, &size));
560: dmRedist = NULL;
561: if (size > 1) PetscCall(DMPlexDistributeOverlap(dm, 1, NULL, &dmRedist));
562: if (!dmRedist) {
563: dmRedist = dm;
564: PetscCall(PetscObjectReference((PetscObject)dmRedist));
565: }
566: PetscCall(PetscFVCreate(comm, &fv));
567: PetscCall(PetscFVSetType(fv, PETSCFVLEASTSQUARES));
568: PetscCall(PetscFVSetNumComponents(fv, user->numComponents));
569: PetscCall(PetscFVSetSpatialDimension(fv, dim));
570: PetscCall(PetscFVSetFromOptions(fv));
571: PetscCall(PetscFVSetUp(fv));
572: {
573: PetscSF pointSF;
574: DMLabel label;
576: PetscCall(DMCreateLabel(dmRedist, "Face Sets"));
577: PetscCall(DMGetLabel(dmRedist, "Face Sets", &label));
578: PetscCall(DMGetPointSF(dmRedist, &pointSF));
579: PetscCall(PetscObjectReference((PetscObject)pointSF));
580: PetscCall(DMSetPointSF(dmRedist, NULL));
581: PetscCall(DMPlexMarkBoundaryFaces(dmRedist, 1, label));
582: PetscCall(DMSetPointSF(dmRedist, pointSF));
583: PetscCall(PetscSFDestroy(&pointSF));
584: }
585: PetscCall(DMPlexConstructGhostCells(dmRedist, NULL, NULL, &dmfv));
586: PetscCall(DMDestroy(&dmRedist));
587: PetscCall(DMSetNumFields(dmfv, 1));
588: PetscCall(DMSetField(dmfv, 0, NULL, (PetscObject)fv));
589: PetscCall(DMCreateDS(dmfv));
590: PetscCall(DMPlexGetReferenceTree(dm, &refTree));
591: if (refTree) PetscCall(DMCopyDisc(dmfv, refTree));
592: PetscCall(DMPlexGetGradientDM(dmfv, fv, &dmgrad));
593: PetscCall(DMPlexGetHeightStratum(dmfv, 0, &cStart, &cEnd));
594: nvecs = dim * (dim + 1) / 2;
595: PetscCall(DMPlexGetGeometryFVM(dmfv, NULL, &cellgeom, NULL));
596: PetscCall(VecGetDM(cellgeom, &dmCell));
597: PetscCall(VecGetArrayRead(cellgeom, &cgeom));
598: PetscCall(DMGetGlobalVector(dmgrad, &grad));
599: PetscCall(DMGetLocalVector(dmgrad, &locGrad));
600: PetscCall(DMPlexGetCellTypeStratum(dmgrad, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
601: cEndInterior = (cEndInterior < 0) ? cEnd : cEndInterior;
602: for (v = 0; v < nvecs; v++) {
603: Vec locX;
604: PetscInt c;
605: PetscScalar trueGrad[3][3] = {{0.}};
606: const PetscScalar *gradArray;
607: PetscReal maxDiff, maxDiffGlob;
609: PetscCall(DMGetLocalVector(dmfv, &locX));
610: /* get the local projection of the rigid body mode */
611: for (c = cStart; c < cEnd; c++) {
612: PetscFVCellGeom *cg;
613: PetscScalar cx[3] = {0., 0., 0.};
615: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
616: if (v < dim) {
617: cx[v] = 1.;
618: } else {
619: PetscInt w = v - dim;
621: cx[(w + 1) % dim] = cg->centroid[(w + 2) % dim];
622: cx[(w + 2) % dim] = -cg->centroid[(w + 1) % dim];
623: }
624: PetscCall(DMPlexVecSetClosure(dmfv, NULL, locX, c, cx, INSERT_ALL_VALUES));
625: }
626: /* TODO: this isn't in any header */
627: PetscCall(DMPlexReconstructGradientsFVM(dmfv, locX, grad));
628: PetscCall(DMGlobalToLocalBegin(dmgrad, grad, INSERT_VALUES, locGrad));
629: PetscCall(DMGlobalToLocalEnd(dmgrad, grad, INSERT_VALUES, locGrad));
630: PetscCall(VecGetArrayRead(locGrad, &gradArray));
631: /* compare computed gradient to exact gradient */
632: if (v >= dim) {
633: PetscInt w = v - dim;
635: trueGrad[(w + 1) % dim][(w + 2) % dim] = 1.;
636: trueGrad[(w + 2) % dim][(w + 1) % dim] = -1.;
637: }
638: maxDiff = 0.;
639: for (c = cStart; c < cEndInterior; c++) {
640: PetscScalar *compGrad;
641: PetscInt i, j, k;
642: PetscReal FrobDiff = 0.;
644: PetscCall(DMPlexPointLocalRead(dmgrad, c, gradArray, &compGrad));
646: for (i = 0, k = 0; i < dim; i++) {
647: for (j = 0; j < dim; j++, k++) {
648: PetscScalar diff = compGrad[k] - trueGrad[i][j];
649: FrobDiff += PetscRealPart(diff * PetscConj(diff));
650: }
651: }
652: FrobDiff = PetscSqrtReal(FrobDiff);
653: maxDiff = PetscMax(maxDiff, FrobDiff);
654: }
655: PetscCall(MPIU_Allreduce(&maxDiff, &maxDiffGlob, 1, MPIU_REAL, MPIU_MAX, comm));
656: allVecMaxDiff = PetscMax(allVecMaxDiff, maxDiffGlob);
657: PetscCall(VecRestoreArrayRead(locGrad, &gradArray));
658: PetscCall(DMRestoreLocalVector(dmfv, &locX));
659: }
660: if (allVecMaxDiff < fvTol) {
661: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Finite volume gradient reconstruction: PASS\n"));
662: } else {
663: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Finite volume gradient reconstruction: FAIL at tolerance %g with max difference %g\n", (double)fvTol, (double)allVecMaxDiff));
664: }
665: PetscCall(DMRestoreLocalVector(dmgrad, &locGrad));
666: PetscCall(DMRestoreGlobalVector(dmgrad, &grad));
667: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
668: PetscCall(DMDestroy(&dmfv));
669: PetscCall(PetscFVDestroy(&fv));
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
674: {
675: Vec u;
676: PetscReal n[3] = {1.0, 1.0, 1.0};
678: PetscFunctionBeginUser;
679: PetscCall(DMGetGlobalVector(dm, &u));
680: /* Project function into FE function space */
681: PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u));
682: PetscCall(VecViewFromOptions(u, NULL, "-projection_view"));
683: /* Compare approximation to exact in L_2 */
684: PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error));
685: PetscCall(DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer));
686: PetscCall(DMRestoreGlobalVector(dm, &u));
687: PetscFunctionReturn(PETSC_SUCCESS);
688: }
690: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
691: {
692: PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
693: PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
694: void *exactCtxs[3];
695: MPI_Comm comm;
696: PetscReal error, errorDer, tol = PETSC_SMALL;
698: PetscFunctionBeginUser;
699: exactCtxs[0] = user;
700: exactCtxs[1] = user;
701: exactCtxs[2] = user;
702: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
703: /* Setup functions to approximate */
704: switch (order) {
705: case 0:
706: exactFuncs[0] = constant;
707: exactFuncDers[0] = constantDer;
708: break;
709: case 1:
710: exactFuncs[0] = linear;
711: exactFuncDers[0] = linearDer;
712: break;
713: case 2:
714: exactFuncs[0] = quadratic;
715: exactFuncDers[0] = quadraticDer;
716: break;
717: case 3:
718: exactFuncs[0] = cubic;
719: exactFuncDers[0] = cubicDer;
720: break;
721: default:
722: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for order %" PetscInt_FMT, order);
723: }
724: PetscCall(ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
725: /* Report result */
726: if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
727: else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
728: if (errorDer > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer));
729: else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
730: PetscFunctionReturn(PETSC_SUCCESS);
731: }
733: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
734: {
735: PetscErrorCode (*exactFuncs[1])(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
736: PetscErrorCode (*exactFuncDers[1])(PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx);
737: PetscReal n[3] = {1.0, 1.0, 1.0};
738: void *exactCtxs[3];
739: DM rdm, idm, fdm;
740: Mat Interp;
741: Vec iu, fu, scaling;
742: MPI_Comm comm;
743: PetscInt dim;
744: PetscReal error, errorDer, tol = PETSC_SMALL;
746: PetscFunctionBeginUser;
747: exactCtxs[0] = user;
748: exactCtxs[1] = user;
749: exactCtxs[2] = user;
750: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
751: PetscCall(DMGetDimension(dm, &dim));
752: PetscCall(DMRefine(dm, comm, &rdm));
753: PetscCall(DMSetCoarseDM(rdm, dm));
754: PetscCall(DMPlexSetRegularRefinement(rdm, user->convRefine));
755: if (user->tree) {
756: DM refTree;
757: PetscCall(DMPlexGetReferenceTree(dm, &refTree));
758: PetscCall(DMPlexSetReferenceTree(rdm, refTree));
759: }
760: if (user->useDA) PetscCall(DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
761: PetscCall(SetupSection(rdm, user));
762: /* Setup functions to approximate */
763: switch (order) {
764: case 0:
765: exactFuncs[0] = constant;
766: exactFuncDers[0] = constantDer;
767: break;
768: case 1:
769: exactFuncs[0] = linear;
770: exactFuncDers[0] = linearDer;
771: break;
772: case 2:
773: exactFuncs[0] = quadratic;
774: exactFuncDers[0] = quadraticDer;
775: break;
776: case 3:
777: exactFuncs[0] = cubic;
778: exactFuncDers[0] = cubicDer;
779: break;
780: default:
781: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %" PetscInt_FMT " order %" PetscInt_FMT, dim, order);
782: }
783: idm = checkRestrict ? rdm : dm;
784: fdm = checkRestrict ? dm : rdm;
785: PetscCall(DMGetGlobalVector(idm, &iu));
786: PetscCall(DMGetGlobalVector(fdm, &fu));
787: PetscCall(DMSetApplicationContext(dm, user));
788: PetscCall(DMSetApplicationContext(rdm, user));
789: PetscCall(DMCreateInterpolation(dm, rdm, &Interp, &scaling));
790: /* Project function into initial FE function space */
791: PetscCall(DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu));
792: /* Interpolate function into final FE function space */
793: if (checkRestrict) {
794: PetscCall(MatRestrict(Interp, iu, fu));
795: PetscCall(VecPointwiseMult(fu, scaling, fu));
796: } else PetscCall(MatInterpolate(Interp, iu, fu));
797: /* Compare approximation to exact in L_2 */
798: PetscCall(DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error));
799: PetscCall(DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer));
800: /* Report result */
801: if (error > tol) PetscCall(PetscPrintf(comm, "Interpolation tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
802: else PetscCall(PetscPrintf(comm, "Interpolation tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
803: if (errorDer > tol) PetscCall(PetscPrintf(comm, "Interpolation tests FAIL for order %" PetscInt_FMT " derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer));
804: else PetscCall(PetscPrintf(comm, "Interpolation tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
805: PetscCall(DMRestoreGlobalVector(idm, &iu));
806: PetscCall(DMRestoreGlobalVector(fdm, &fu));
807: PetscCall(MatDestroy(&Interp));
808: PetscCall(VecDestroy(&scaling));
809: PetscCall(DMDestroy(&rdm));
810: PetscFunctionReturn(PETSC_SUCCESS);
811: }
813: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
814: {
815: DM odm = dm, rdm = NULL, cdm = NULL;
816: PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {trig};
817: PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx) = {trigDer};
818: void *exactCtxs[3];
819: PetscInt r, c, cStart, cEnd;
820: PetscReal errorOld, errorDerOld, error, errorDer, rel, len, lenOld;
821: double p;
823: PetscFunctionBeginUser;
824: if (!user->convergence) PetscFunctionReturn(PETSC_SUCCESS);
825: exactCtxs[0] = user;
826: exactCtxs[1] = user;
827: exactCtxs[2] = user;
828: PetscCall(PetscObjectReference((PetscObject)odm));
829: if (!user->convRefine) {
830: for (r = 0; r < Nr; ++r) {
831: PetscCall(DMRefine(odm, PetscObjectComm((PetscObject)dm), &rdm));
832: PetscCall(DMDestroy(&odm));
833: odm = rdm;
834: }
835: PetscCall(SetupSection(odm, user));
836: }
837: PetscCall(ComputeError(odm, exactFuncs, exactFuncDers, exactCtxs, &errorOld, &errorDerOld, user));
838: if (user->convRefine) {
839: for (r = 0; r < Nr; ++r) {
840: PetscCall(DMRefine(odm, PetscObjectComm((PetscObject)dm), &rdm));
841: if (user->useDA) PetscCall(DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
842: PetscCall(SetupSection(rdm, user));
843: PetscCall(ComputeError(rdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
844: p = PetscLog2Real(errorOld / error);
845: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Function convergence rate at refinement %" PetscInt_FMT ": %.2f\n", r, (double)p));
846: p = PetscLog2Real(errorDerOld / errorDer);
847: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Derivative convergence rate at refinement %" PetscInt_FMT ": %.2f\n", r, (double)p));
848: PetscCall(DMDestroy(&odm));
849: odm = rdm;
850: errorOld = error;
851: errorDerOld = errorDer;
852: }
853: } else {
854: /* PetscCall(ComputeLongestEdge(dm, &lenOld)); */
855: PetscCall(DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd));
856: lenOld = cEnd - cStart;
857: for (c = 0; c < Nr; ++c) {
858: PetscCall(DMCoarsen(odm, PetscObjectComm((PetscObject)dm), &cdm));
859: if (user->useDA) PetscCall(DMDASetVertexCoordinates(cdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
860: PetscCall(SetupSection(cdm, user));
861: PetscCall(ComputeError(cdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
862: /* PetscCall(ComputeLongestEdge(cdm, &len)); */
863: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
864: len = cEnd - cStart;
865: rel = error / errorOld;
866: p = PetscLogReal(rel) / PetscLogReal(lenOld / len);
867: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Function convergence rate at coarsening %" PetscInt_FMT ": %.2f\n", c, (double)p));
868: rel = errorDer / errorDerOld;
869: p = PetscLogReal(rel) / PetscLogReal(lenOld / len);
870: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Derivative convergence rate at coarsening %" PetscInt_FMT ": %.2f\n", c, (double)p));
871: PetscCall(DMDestroy(&odm));
872: odm = cdm;
873: errorOld = error;
874: errorDerOld = errorDer;
875: lenOld = len;
876: }
877: }
878: PetscCall(DMDestroy(&odm));
879: PetscFunctionReturn(PETSC_SUCCESS);
880: }
882: int main(int argc, char **argv)
883: {
884: DM dm;
885: AppCtx user; /* user-defined work context */
886: PetscInt dim = 2;
887: PetscBool simplex = PETSC_FALSE;
889: PetscFunctionBeginUser;
890: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
891: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
892: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
893: if (!user.useDA) {
894: PetscCall(DMGetDimension(dm, &dim));
895: PetscCall(DMPlexIsSimplex(dm, &simplex));
896: }
897: PetscCall(DMPlexMetricSetFromOptions(dm));
898: user.numComponents = user.numComponents < 0 ? dim : user.numComponents;
899: PetscCall(PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.numComponents, simplex, NULL, user.qorder, &user.fe));
900: PetscCall(SetupSection(dm, &user));
901: if (user.testFEjacobian) PetscCall(TestFEJacobian(dm, &user));
902: if (user.testFVgrad) PetscCall(TestFVGrad(dm, &user));
903: if (user.testInjector) PetscCall(TestInjector(dm, &user));
904: PetscCall(CheckFunctions(dm, user.porder, &user));
905: {
906: PetscDualSpace dsp;
907: PetscInt k;
909: PetscCall(PetscFEGetDualSpace(user.fe, &dsp));
910: PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
911: if (dim == 2 && simplex == PETSC_TRUE && user.tree == PETSC_FALSE && k == 0) {
912: PetscCall(CheckInterpolation(dm, PETSC_FALSE, user.porder, &user));
913: PetscCall(CheckInterpolation(dm, PETSC_TRUE, user.porder, &user));
914: }
915: }
916: PetscCall(CheckConvergence(dm, 3, &user));
917: PetscCall(PetscFEDestroy(&user.fe));
918: PetscCall(DMDestroy(&dm));
919: PetscCall(PetscFinalize());
920: return 0;
921: }
923: /*TEST
925: test:
926: suffix: 1
927: requires: triangle
929: # 2D P_1 on a triangle
930: test:
931: suffix: p1_2d_0
932: requires: triangle
933: args: -petscspace_degree 1 -qorder 1 -convergence
934: test:
935: suffix: p1_2d_1
936: requires: triangle
937: args: -petscspace_degree 1 -qorder 1 -porder 1
938: test:
939: suffix: p1_2d_2
940: requires: triangle
941: args: -petscspace_degree 1 -qorder 1 -porder 2
942: test:
943: suffix: p1_2d_3
944: requires: triangle mmg
945: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
946: test:
947: suffix: p1_2d_4
948: requires: triangle mmg
949: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
950: test:
951: suffix: p1_2d_5
952: requires: triangle mmg
953: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
955: # 3D P_1 on a tetrahedron
956: test:
957: suffix: p1_3d_0
958: requires: ctetgen
959: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -convergence
960: test:
961: suffix: p1_3d_1
962: requires: ctetgen
963: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 1
964: test:
965: suffix: p1_3d_2
966: requires: ctetgen
967: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 2
968: test:
969: suffix: p1_3d_3
970: requires: ctetgen mmg
971: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
972: test:
973: suffix: p1_3d_4
974: requires: ctetgen mmg
975: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
976: test:
977: suffix: p1_3d_5
978: requires: ctetgen mmg
979: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
981: # 2D P_2 on a triangle
982: test:
983: suffix: p2_2d_0
984: requires: triangle
985: args: -petscspace_degree 2 -qorder 2 -convergence
986: test:
987: suffix: p2_2d_1
988: requires: triangle
989: args: -petscspace_degree 2 -qorder 2 -porder 1
990: test:
991: suffix: p2_2d_2
992: requires: triangle
993: args: -petscspace_degree 2 -qorder 2 -porder 2
994: test:
995: suffix: p2_2d_3
996: requires: triangle mmg
997: args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
998: test:
999: suffix: p2_2d_4
1000: requires: triangle mmg
1001: args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1002: test:
1003: suffix: p2_2d_5
1004: requires: triangle mmg
1005: args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0
1007: # 3D P_2 on a tetrahedron
1008: test:
1009: suffix: p2_3d_0
1010: requires: ctetgen
1011: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -convergence
1012: test:
1013: suffix: p2_3d_1
1014: requires: ctetgen
1015: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 1
1016: test:
1017: suffix: p2_3d_2
1018: requires: ctetgen
1019: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 2
1020: test:
1021: suffix: p2_3d_3
1022: requires: ctetgen mmg
1023: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1024: test:
1025: suffix: p2_3d_4
1026: requires: ctetgen mmg
1027: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1028: test:
1029: suffix: p2_3d_5
1030: requires: ctetgen mmg
1031: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0
1033: # 2D Q_1 on a quadrilaterial DA
1034: test:
1035: suffix: q1_2d_da_0
1036: requires: broken
1037: args: -use_da 1 -petscspace_degree 1 -qorder 1 -convergence
1038: test:
1039: suffix: q1_2d_da_1
1040: requires: broken
1041: args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 1
1042: test:
1043: suffix: q1_2d_da_2
1044: requires: broken
1045: args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 2
1047: # 2D Q_1 on a quadrilaterial Plex
1048: test:
1049: suffix: q1_2d_plex_0
1050: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -convergence
1051: test:
1052: suffix: q1_2d_plex_1
1053: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1
1054: test:
1055: suffix: q1_2d_plex_2
1056: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2
1057: test:
1058: suffix: q1_2d_plex_3
1059: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1 -shear_coords
1060: test:
1061: suffix: q1_2d_plex_4
1062: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2 -shear_coords
1063: test:
1064: suffix: q1_2d_plex_5
1065: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 0 -non_affine_coords -convergence
1066: test:
1067: suffix: q1_2d_plex_6
1068: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 1 -non_affine_coords -convergence
1069: test:
1070: suffix: q1_2d_plex_7
1071: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 2 -non_affine_coords -convergence
1073: # 2D Q_2 on a quadrilaterial
1074: test:
1075: suffix: q2_2d_plex_0
1076: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1077: test:
1078: suffix: q2_2d_plex_1
1079: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1080: test:
1081: suffix: q2_2d_plex_2
1082: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1083: test:
1084: suffix: q2_2d_plex_3
1085: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1086: test:
1087: suffix: q2_2d_plex_4
1088: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1089: test:
1090: suffix: q2_2d_plex_5
1091: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 0 -non_affine_coords -convergence
1092: test:
1093: suffix: q2_2d_plex_6
1094: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 1 -non_affine_coords -convergence
1095: test:
1096: suffix: q2_2d_plex_7
1097: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 2 -non_affine_coords -convergence
1099: # 2D P_3 on a triangle
1100: test:
1101: suffix: p3_2d_0
1102: requires: triangle !single
1103: args: -petscspace_degree 3 -qorder 3 -convergence
1104: test:
1105: suffix: p3_2d_1
1106: requires: triangle !single
1107: args: -petscspace_degree 3 -qorder 3 -porder 1
1108: test:
1109: suffix: p3_2d_2
1110: requires: triangle !single
1111: args: -petscspace_degree 3 -qorder 3 -porder 2
1112: test:
1113: suffix: p3_2d_3
1114: requires: triangle !single
1115: args: -petscspace_degree 3 -qorder 3 -porder 3
1116: test:
1117: suffix: p3_2d_4
1118: requires: triangle mmg
1119: args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -convergence -conv_refine 0
1120: test:
1121: suffix: p3_2d_5
1122: requires: triangle mmg
1123: args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 1 -conv_refine 0
1124: test:
1125: suffix: p3_2d_6
1126: requires: triangle mmg
1127: args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 3 -conv_refine 0
1129: # 2D Q_3 on a quadrilaterial
1130: test:
1131: suffix: q3_2d_0
1132: requires: !single
1133: args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -convergence
1134: test:
1135: suffix: q3_2d_1
1136: requires: !single
1137: args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 1
1138: test:
1139: suffix: q3_2d_2
1140: requires: !single
1141: args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 2
1142: test:
1143: suffix: q3_2d_3
1144: requires: !single
1145: args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 3
1147: # 2D P_1disc on a triangle/quadrilateral
1148: test:
1149: suffix: p1d_2d_0
1150: requires: triangle
1151: args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1152: test:
1153: suffix: p1d_2d_1
1154: requires: triangle
1155: args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1156: test:
1157: suffix: p1d_2d_2
1158: requires: triangle
1159: args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1160: test:
1161: suffix: p1d_2d_3
1162: requires: triangle
1163: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1164: filter: sed -e "s/convergence rate at refinement 0: 2/convergence rate at refinement 0: 1.9/g"
1165: test:
1166: suffix: p1d_2d_4
1167: requires: triangle
1168: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1169: test:
1170: suffix: p1d_2d_5
1171: requires: triangle
1172: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1174: # 2D BDM_1 on a triangle
1175: test:
1176: suffix: bdm1_2d_0
1177: requires: triangle
1178: args: -petscspace_degree 1 -petscdualspace_type bdm \
1179: -num_comp 2 -qorder 1 -convergence
1180: test:
1181: suffix: bdm1_2d_1
1182: requires: triangle
1183: args: -petscspace_degree 1 -petscdualspace_type bdm \
1184: -num_comp 2 -qorder 1 -porder 1
1185: test:
1186: suffix: bdm1_2d_2
1187: requires: triangle
1188: args: -petscspace_degree 1 -petscdualspace_type bdm \
1189: -num_comp 2 -qorder 1 -porder 2
1191: # 2D BDM_1 on a quadrilateral
1192: test:
1193: suffix: bdm1q_2d_0
1194: requires: triangle
1195: args: -petscspace_degree 1 -petscdualspace_type bdm \
1196: -petscdualspace_lagrange_tensor 1 \
1197: -dm_plex_simplex 0 -num_comp 2 -qorder 1 -convergence
1198: test:
1199: suffix: bdm1q_2d_1
1200: requires: triangle
1201: args: -petscspace_degree 1 -petscdualspace_type bdm \
1202: -petscdualspace_lagrange_tensor 1 \
1203: -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 1
1204: test:
1205: suffix: bdm1q_2d_2
1206: requires: triangle
1207: args: -petscspace_degree 1 -petscdualspace_type bdm \
1208: -petscdualspace_lagrange_tensor 1 \
1209: -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 2
1211: # Test high order quadrature
1212: test:
1213: suffix: p1_quad_2
1214: requires: triangle
1215: args: -petscspace_degree 1 -qorder 2 -porder 1
1216: test:
1217: suffix: p1_quad_5
1218: requires: triangle
1219: args: -petscspace_degree 1 -qorder 5 -porder 1
1220: test:
1221: suffix: p2_quad_3
1222: requires: triangle
1223: args: -petscspace_degree 2 -qorder 3 -porder 2
1224: test:
1225: suffix: p2_quad_5
1226: requires: triangle
1227: args: -petscspace_degree 2 -qorder 5 -porder 2
1228: test:
1229: suffix: q1_quad_2
1230: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 2 -porder 1
1231: test:
1232: suffix: q1_quad_5
1233: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 5 -porder 1
1234: test:
1235: suffix: q2_quad_3
1236: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 3 -porder 1
1237: test:
1238: suffix: q2_quad_5
1239: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 5 -porder 1
1241: # Nonconforming tests
1242: test:
1243: suffix: constraints
1244: args: -dm_coord_space 0 -dm_plex_simplex 0 -petscspace_type tensor -petscspace_degree 1 -qorder 0 -constraints
1245: test:
1246: suffix: nonconforming_tensor_2
1247: nsize: 4
1248: args: -dist_dm_distribute -test_fe_jacobian -test_injector -petscpartitioner_type simple -tree -dm_plex_simplex 0 -dm_plex_max_projection_height 1 -petscspace_type tensor -petscspace_degree 2 -qorder 2 -dm_view ascii::ASCII_INFO_DETAIL
1249: test:
1250: suffix: nonconforming_tensor_3
1251: nsize: 4
1252: args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_max_projection_height 2 -petscspace_type tensor -petscspace_degree 1 -qorder 1 -dm_view ascii::ASCII_INFO_DETAIL
1253: test:
1254: suffix: nonconforming_tensor_2_fv
1255: nsize: 4
1256: args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_simplex 0 -num_comp 2
1257: test:
1258: suffix: nonconforming_tensor_3_fv
1259: nsize: 4
1260: args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -num_comp 3
1261: test:
1262: suffix: nonconforming_tensor_2_hi
1263: requires: !single
1264: nsize: 4
1265: args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_simplex 0 -dm_plex_max_projection_height 1 -petscspace_type tensor -petscspace_degree 4 -qorder 4
1266: test:
1267: suffix: nonconforming_tensor_3_hi
1268: requires: !single skip
1269: nsize: 4
1270: args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -dm_plex_max_projection_height 2 -petscspace_type tensor -petscspace_degree 4 -qorder 4
1271: test:
1272: suffix: nonconforming_simplex_2
1273: requires: triangle
1274: nsize: 4
1275: args: -dist_dm_distribute -test_fe_jacobian -test_injector -petscpartitioner_type simple -tree -dm_plex_max_projection_height 1 -petscspace_degree 2 -qorder 2 -dm_view ascii::ASCII_INFO_DETAIL
1276: test:
1277: suffix: nonconforming_simplex_2_hi
1278: requires: triangle !single
1279: nsize: 4
1280: args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_max_projection_height 1 -petscspace_degree 4 -qorder 4
1281: test:
1282: suffix: nonconforming_simplex_2_fv
1283: requires: triangle
1284: nsize: 4
1285: args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -num_comp 2
1286: test:
1287: suffix: nonconforming_simplex_3
1288: requires: ctetgen
1289: nsize: 4
1290: args: -dist_dm_distribute -test_fe_jacobian -test_injector -petscpartitioner_type simple -tree -dm_plex_dim 3 -dm_plex_max_projection_height 2 -petscspace_degree 2 -qorder 2 -dm_view ascii::ASCII_INFO_DETAIL
1291: test:
1292: suffix: nonconforming_simplex_3_hi
1293: requires: ctetgen skip
1294: nsize: 4
1295: args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_dim 3 -dm_plex_max_projection_height 2 -petscspace_degree 4 -qorder 4
1296: test:
1297: suffix: nonconforming_simplex_3_fv
1298: requires: ctetgen
1299: nsize: 4
1300: args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_dim 3 -num_comp 3
1302: # 3D WXY on a triangular prism
1303: test:
1304: suffix: wxy_0
1305: args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism -qorder 2 -porder 0 \
1306: -petscspace_type sum \
1307: -petscspace_variables 3 \
1308: -petscspace_components 3 \
1309: -petscspace_sum_spaces 2 \
1310: -petscspace_sum_concatenate false \
1311: -sumcomp_0_petscspace_variables 3 \
1312: -sumcomp_0_petscspace_components 3 \
1313: -sumcomp_0_petscspace_degree 1 \
1314: -sumcomp_1_petscspace_variables 3 \
1315: -sumcomp_1_petscspace_components 3 \
1316: -sumcomp_1_petscspace_type wxy \
1317: -petscdualspace_refcell triangular_prism \
1318: -petscdualspace_form_degree 0 \
1319: -petscdualspace_order 1 \
1320: -petscdualspace_components 3
1322: TEST*/
1324: /*
1325: # 2D Q_2 on a quadrilaterial Plex
1326: test:
1327: suffix: q2_2d_plex_0
1328: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1329: test:
1330: suffix: q2_2d_plex_1
1331: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1332: test:
1333: suffix: q2_2d_plex_2
1334: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1335: test:
1336: suffix: q2_2d_plex_3
1337: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1338: test:
1339: suffix: q2_2d_plex_4
1340: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1341: test:
1342: suffix: q2_2d_plex_5
1343: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 0 -non_affine_coords
1344: test:
1345: suffix: q2_2d_plex_6
1346: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 1 -non_affine_coords
1347: test:
1348: suffix: q2_2d_plex_7
1349: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 2 -non_affine_coords
1351: test:
1352: suffix: p1d_2d_6
1353: requires: mmg
1354: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1355: test:
1356: suffix: p1d_2d_7
1357: requires: mmg
1358: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1359: test:
1360: suffix: p1d_2d_8
1361: requires: mmg
1362: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1363: */