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