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