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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, PetscCtx 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, MPI_C_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 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: PetscScalar trueGrad[3][3] = {{0.}};
657: const PetscScalar *gradArray;
658: PetscReal maxDiff;
660: PetscCall(DMGetLocalVector(dmfv, &locX));
661: /* get the local projection of the rigid body mode */
662: for (PetscInt c = cStart; c < cEnd; c++) {
663: PetscFVCellGeom *cg;
664: PetscScalar cx[3] = {0., 0., 0.};
666: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
667: if (v < dim) {
668: cx[v] = 1.;
669: } else {
670: PetscInt w = v - dim;
672: cx[(w + 1) % dim] = cg->centroid[(w + 2) % dim];
673: cx[(w + 2) % dim] = -cg->centroid[(w + 1) % dim];
674: }
675: PetscCall(DMPlexVecSetClosure(dmfv, NULL, locX, c, cx, INSERT_ALL_VALUES));
676: }
677: /* TODO: this isn't in any header */
678: PetscCall(DMPlexReconstructGradientsFVM(dmfv, locX, grad));
679: PetscCall(DMGlobalToLocalBegin(dmgrad, grad, INSERT_VALUES, locGrad));
680: PetscCall(DMGlobalToLocalEnd(dmgrad, grad, INSERT_VALUES, locGrad));
681: PetscCall(VecGetArrayRead(locGrad, &gradArray));
682: /* compare computed gradient to exact gradient */
683: if (v >= dim) {
684: PetscInt w = v - dim;
686: trueGrad[(w + 1) % dim][(w + 2) % dim] = 1.;
687: trueGrad[(w + 2) % dim][(w + 1) % dim] = -1.;
688: }
689: maxDiff = 0.;
690: for (PetscInt c = cStart; c < cEndInterior; c++) {
691: PetscScalar *compGrad;
692: PetscInt i, j, k;
693: PetscReal FrobDiff = 0.;
695: PetscCall(DMPlexPointLocalRead(dmgrad, c, gradArray, &compGrad));
697: for (i = 0, k = 0; i < dim; i++) {
698: for (j = 0; j < dim; j++, k++) {
699: PetscScalar diff = compGrad[k] - trueGrad[i][j];
700: FrobDiff += PetscRealPart(diff * PetscConj(diff));
701: }
702: }
703: FrobDiff = PetscSqrtReal(FrobDiff);
704: maxDiff = PetscMax(maxDiff, FrobDiff);
705: }
706: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &maxDiff, 1, MPIU_REAL, MPIU_MAX, comm));
707: allVecMaxDiff = PetscMax(allVecMaxDiff, maxDiff);
708: PetscCall(VecRestoreArrayRead(locGrad, &gradArray));
709: PetscCall(DMRestoreLocalVector(dmfv, &locX));
710: }
711: if (allVecMaxDiff < fvTol) {
712: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Finite volume gradient reconstruction: PASS\n"));
713: } else {
714: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Finite volume gradient reconstruction: FAIL at tolerance %g with max difference %g\n", (double)fvTol, (double)allVecMaxDiff));
715: }
716: PetscCall(DMRestoreLocalVector(dmgrad, &locGrad));
717: PetscCall(DMRestoreGlobalVector(dmgrad, &grad));
718: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
719: PetscCall(DMDestroy(&dmfv));
720: PetscCall(PetscFVDestroy(&fv));
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: 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)
725: {
726: Vec u;
727: PetscReal n[3] = {1.0, 1.0, 1.0};
729: PetscFunctionBeginUser;
730: PetscCall(DMGetGlobalVector(dm, &u));
731: /* Project function into FE function space */
732: PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u));
733: PetscCall(VecViewFromOptions(u, NULL, "-projection_view"));
734: /* Compare approximation to exact in L_2 */
735: PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error));
736: PetscCall(DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer));
737: PetscCall(DMRestoreGlobalVector(dm, &u));
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
742: {
743: PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
744: PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
745: void *exactCtxs[3];
746: MPI_Comm comm;
747: PetscReal error, errorDer, tol = PETSC_SMALL;
749: PetscFunctionBeginUser;
750: exactCtxs[0] = user;
751: exactCtxs[1] = user;
752: exactCtxs[2] = user;
753: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
754: /* Setup functions to approximate */
755: switch (order) {
756: case 0:
757: if (user->RT) {
758: exactFuncs[0] = rt0;
759: exactFuncDers[0] = rt0Der;
760: } else {
761: exactFuncs[0] = constant;
762: exactFuncDers[0] = constantDer;
763: }
764: break;
765: case 1:
766: if (user->RT) {
767: exactFuncs[0] = rt1;
768: exactFuncDers[0] = rt1Der;
769: } else {
770: exactFuncs[0] = linear;
771: exactFuncDers[0] = linearDer;
772: }
773: break;
774: case 2:
775: exactFuncs[0] = quadratic;
776: exactFuncDers[0] = quadraticDer;
777: break;
778: case 3:
779: exactFuncs[0] = cubic;
780: exactFuncDers[0] = cubicDer;
781: break;
782: default:
783: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for order %" PetscInt_FMT, order);
784: }
785: PetscCall(ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
786: /* Report result */
787: if (error > tol) PetscCall(PetscPrintf(comm, "Function tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
788: else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
789: 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));
790: else PetscCall(PetscPrintf(comm, "Function tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
795: {
796: PetscErrorCode (*exactFuncs[1])(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, PetscCtx ctx);
797: PetscErrorCode (*exactFuncDers[1])(PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, PetscCtx ctx);
798: PetscReal n[3] = {1.0, 1.0, 1.0};
799: void *exactCtxs[3];
800: DM rdm, idm, fdm;
801: Mat Interp;
802: Vec iu, fu, scaling;
803: MPI_Comm comm;
804: PetscInt dim;
805: PetscReal error, errorDer, tol = PETSC_SMALL;
807: PetscFunctionBeginUser;
808: exactCtxs[0] = user;
809: exactCtxs[1] = user;
810: exactCtxs[2] = user;
811: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
812: PetscCall(DMGetDimension(dm, &dim));
813: PetscCall(DMRefine(dm, comm, &rdm));
814: PetscCall(DMSetCoarseDM(rdm, dm));
815: PetscCall(DMGetCoordinatesLocalSetUp(rdm));
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(DMGetCoordinatesLocalSetUp(fdm));
861: PetscCall(DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error));
862: PetscCall(DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer));
863: /* Report result */
864: if (error > tol) PetscCall(PetscPrintf(comm, "Interpolation tests FAIL for order %" PetscInt_FMT " at tolerance %g error %g\n", order, (double)tol, (double)error));
865: else PetscCall(PetscPrintf(comm, "Interpolation tests pass for order %" PetscInt_FMT " at tolerance %g\n", order, (double)tol));
866: 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));
867: else PetscCall(PetscPrintf(comm, "Interpolation tests pass for order %" PetscInt_FMT " derivatives at tolerance %g\n", order, (double)tol));
868: PetscCall(DMRestoreGlobalVector(idm, &iu));
869: PetscCall(DMRestoreGlobalVector(fdm, &fu));
870: PetscCall(MatDestroy(&Interp));
871: PetscCall(VecDestroy(&scaling));
872: PetscCall(DMDestroy(&rdm));
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
877: {
878: DM odm = dm, rdm = NULL, cdm = NULL;
879: PetscErrorCode (*exactFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx) = {trig};
880: PetscErrorCode (*exactFuncDers[1])(PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, PetscCtx ctx) = {trigDer};
881: void *exactCtxs[3];
882: PetscInt r, c, cStart, cEnd;
883: PetscReal errorOld, errorDerOld, error, errorDer, rel, len, lenOld;
884: double p;
886: PetscFunctionBeginUser;
887: if (!user->convergence) PetscFunctionReturn(PETSC_SUCCESS);
888: exactCtxs[0] = user;
889: exactCtxs[1] = user;
890: exactCtxs[2] = user;
891: PetscCall(PetscObjectReference((PetscObject)odm));
892: if (!user->convRefine) {
893: for (r = 0; r < Nr; ++r) {
894: PetscCall(DMRefine(odm, PetscObjectComm((PetscObject)dm), &rdm));
895: PetscCall(DMDestroy(&odm));
896: odm = rdm;
897: }
898: PetscCall(SetupSection(odm, user));
899: }
900: PetscCall(ComputeError(odm, exactFuncs, exactFuncDers, exactCtxs, &errorOld, &errorDerOld, user));
901: if (user->convRefine) {
902: for (r = 0; r < Nr; ++r) {
903: PetscCall(DMRefine(odm, PetscObjectComm((PetscObject)dm), &rdm));
904: if (user->useDA) PetscCall(DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
905: PetscCall(SetupSection(rdm, user));
906: PetscCall(ComputeError(rdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
907: p = PetscLog2Real(errorOld / error);
908: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Function convergence rate at refinement %" PetscInt_FMT ": %.2f\n", r, p));
909: p = PetscLog2Real(errorDerOld / errorDer);
910: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Derivative convergence rate at refinement %" PetscInt_FMT ": %.2f\n", r, p));
911: PetscCall(DMDestroy(&odm));
912: odm = rdm;
913: errorOld = error;
914: errorDerOld = errorDer;
915: }
916: } else {
917: /* PetscCall(ComputeLongestEdge(dm, &lenOld)); */
918: PetscCall(DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd));
919: lenOld = cEnd - cStart;
920: for (c = 0; c < Nr; ++c) {
921: PetscCall(DMCoarsen(odm, PetscObjectComm((PetscObject)dm), &cdm));
922: if (user->useDA) PetscCall(DMDASetVertexCoordinates(cdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
923: PetscCall(SetupSection(cdm, user));
924: PetscCall(ComputeError(cdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user));
925: /* PetscCall(ComputeLongestEdge(cdm, &len)); */
926: PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
927: len = cEnd - cStart;
928: rel = error / errorOld;
929: p = PetscLogReal(rel) / PetscLogReal(lenOld / len);
930: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Function convergence rate at coarsening %" PetscInt_FMT ": %.2f\n", c, p));
931: rel = errorDer / errorDerOld;
932: p = PetscLogReal(rel) / PetscLogReal(lenOld / len);
933: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Derivative convergence rate at coarsening %" PetscInt_FMT ": %.2f\n", c, p));
934: PetscCall(DMDestroy(&odm));
935: odm = cdm;
936: errorOld = error;
937: errorDerOld = errorDer;
938: lenOld = len;
939: }
940: }
941: PetscCall(DMDestroy(&odm));
942: PetscFunctionReturn(PETSC_SUCCESS);
943: }
945: int main(int argc, char **argv)
946: {
947: DM dm;
948: AppCtx user; /* user-defined work context */
949: PetscInt dim = 2;
950: PetscBool simplex = PETSC_FALSE;
952: PetscFunctionBeginUser;
953: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
954: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
955: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
956: if (!user.useDA) {
957: PetscCall(DMGetDimension(dm, &dim));
958: PetscCall(DMPlexIsSimplex(dm, &simplex));
959: }
960: PetscCall(DMPlexMetricSetFromOptions(dm));
961: user.numComponents = user.numComponents < 0 ? dim : user.numComponents;
962: PetscCall(PetscFECreateDefault(PETSC_COMM_WORLD, dim, user.numComponents, simplex, NULL, user.qorder, &user.fe));
963: PetscCall(SetupSection(dm, &user));
964: if (user.testFEjacobian) PetscCall(TestFEJacobian(dm, &user));
965: if (user.testFVgrad) PetscCall(TestFVGrad(dm, &user));
966: if (user.testInjector) PetscCall(TestInjector(dm, &user));
967: PetscCall(CheckFunctions(dm, user.porder, &user));
968: {
969: PetscDualSpace dsp;
970: PetscInt k;
972: PetscCall(PetscFEGetDualSpace(user.fe, &dsp));
973: PetscCall(PetscDualSpaceGetDeRahm(dsp, &k));
974: if (dim == 2 && user.constraints == PETSC_FALSE && user.tree == PETSC_FALSE && k == 0) {
975: PetscCall(CheckInterpolation(dm, PETSC_FALSE, user.porder, &user));
976: PetscCall(CheckInterpolation(dm, PETSC_TRUE, user.porder, &user));
977: }
978: }
979: PetscCall(CheckConvergence(dm, 3, &user));
980: PetscCall(PetscFEDestroy(&user.fe));
981: PetscCall(DMDestroy(&dm));
982: PetscCall(PetscFinalize());
983: return 0;
984: }
986: /*TEST
988: test:
989: suffix: 1
990: requires: triangle
992: # 2D P_0 on a triangle
993: test:
994: suffix: p0_2d_0
995: requires: triangle
996: args: -petscspace_degree 1 -qorder 1 -convergence
997: test:
998: suffix: p0_2d_1
999: requires: triangle
1000: args: -petscspace_degree 1 -qorder 1 -porder 1
1002: # 2D P_1 on a triangle
1003: test:
1004: suffix: p1_2d_0
1005: requires: triangle
1006: args: -petscspace_degree 1 -qorder 1 -convergence
1007: test:
1008: suffix: p1_2d_1
1009: requires: triangle
1010: args: -petscspace_degree 1 -qorder 1 -porder 1
1011: test:
1012: suffix: p1_2d_2
1013: requires: triangle
1014: args: -petscspace_degree 1 -qorder 1 -porder 2
1015: test:
1016: suffix: p1_2d_3
1017: requires: triangle mmg
1018: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1019: test:
1020: suffix: p1_2d_4
1021: requires: triangle mmg
1022: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1023: test:
1024: suffix: p1_2d_5
1025: requires: triangle mmg
1026: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1028: # 3D P_1 on a tetrahedron
1029: test:
1030: suffix: p1_3d_0
1031: requires: ctetgen
1032: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -convergence
1033: test:
1034: suffix: p1_3d_1
1035: requires: ctetgen
1036: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 1
1037: test:
1038: suffix: p1_3d_2
1039: requires: ctetgen
1040: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -porder 2
1041: test:
1042: suffix: p1_3d_3
1043: requires: ctetgen mmg
1044: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1045: test:
1046: suffix: p1_3d_4
1047: requires: ctetgen mmg
1048: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1049: test:
1050: suffix: p1_3d_5
1051: requires: ctetgen mmg
1052: args: -dm_plex_dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1054: # 2D P_2 on a triangle
1055: test:
1056: suffix: p2_2d_0
1057: requires: triangle
1058: args: -petscspace_degree 2 -qorder 2 -convergence
1059: test:
1060: suffix: p2_2d_1
1061: requires: triangle
1062: args: -petscspace_degree 2 -qorder 2 -porder 1
1063: test:
1064: suffix: p2_2d_2
1065: requires: triangle
1066: args: -petscspace_degree 2 -qorder 2 -porder 2
1067: test:
1068: suffix: p2_2d_3
1069: requires: triangle mmg
1070: args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1071: test:
1072: suffix: p2_2d_4
1073: requires: triangle mmg
1074: args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1075: test:
1076: suffix: p2_2d_5
1077: requires: triangle mmg
1078: args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0
1080: # 3D P_2 on a tetrahedron
1081: test:
1082: suffix: p2_3d_0
1083: requires: ctetgen
1084: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -convergence
1085: test:
1086: suffix: p2_3d_1
1087: requires: ctetgen
1088: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 1
1089: test:
1090: suffix: p2_3d_2
1091: requires: ctetgen
1092: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -porder 2
1093: test:
1094: suffix: p2_3d_3
1095: requires: ctetgen mmg
1096: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1097: test:
1098: suffix: p2_3d_4
1099: requires: ctetgen mmg
1100: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1101: test:
1102: suffix: p2_3d_5
1103: requires: ctetgen mmg
1104: args: -dm_plex_dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0
1106: # 2D Q_1 on a quadrilaterial DA
1107: test:
1108: suffix: q1_2d_da_0
1109: TODO: broken
1110: args: -use_da 1 -petscspace_degree 1 -qorder 1 -convergence
1111: test:
1112: suffix: q1_2d_da_1
1113: TODO: broken
1114: args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 1
1115: test:
1116: suffix: q1_2d_da_2
1117: TODO: broken
1118: args: -use_da 1 -petscspace_degree 1 -qorder 1 -porder 2
1120: # 2D P_0 on a quadrilaterial Plex
1121: test:
1122: suffix: p0_2d_plex_0
1123: args: -dm_plex_simplex 0 -petscspace_degree 0 -qorder 1 -convergence
1124: test:
1125: suffix: p0_2d_plex_1
1126: args: -dm_plex_simplex 0 -petscspace_degree 0 -qorder 1 -porder 1
1128: # 2D Q_1 on a quadrilaterial Plex
1129: test:
1130: suffix: q1_2d_plex_0
1131: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -convergence
1132: test:
1133: suffix: q1_2d_plex_1
1134: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1
1135: test:
1136: suffix: q1_2d_plex_2
1137: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2
1138: test:
1139: suffix: q1_2d_plex_3
1140: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 1 -shear_coords
1141: test:
1142: suffix: q1_2d_plex_4
1143: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 1 -porder 2 -shear_coords
1144: test:
1145: suffix: q1_2d_plex_5
1146: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 0 -non_affine_coords -convergence
1147: test:
1148: suffix: q1_2d_plex_6
1149: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 1 -non_affine_coords -convergence
1150: test:
1151: suffix: q1_2d_plex_7
1152: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 2 -non_affine_coords -convergence
1153: test:
1154: suffix: q1_2d_plex_8
1155: requires: triangle
1156: args: -dist_dm_refine 1 -dist_dm_plex_transform_type refine_tobox -petscspace_degree 1 -qorder 1 -convergence
1158: # 2D Q_2 on a quadrilaterial
1159: # The derivative interpolation test fails in single because we lose precision
1160: test:
1161: suffix: q2_2d_plex_0
1162: requires: !single
1163: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1164: test:
1165: suffix: q2_2d_plex_1
1166: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1167: test:
1168: suffix: q2_2d_plex_2
1169: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1170: test:
1171: suffix: q2_2d_plex_3
1172: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1173: test:
1174: suffix: q2_2d_plex_4
1175: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1176: # The derivative interpolation test fails in single because we lose precision
1177: test:
1178: suffix: q2_2d_plex_5
1179: requires: !single
1180: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 0 -non_affine_coords -convergence
1181: # The derivative interpolation test fails in single because we lose precision
1182: test:
1183: suffix: q2_2d_plex_6
1184: requires: !single
1185: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 1 -non_affine_coords -convergence
1186: test:
1187: suffix: q2_2d_plex_7
1188: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 2 -non_affine_coords -convergence
1190: # 2D P_3 on a triangle
1191: test:
1192: suffix: p3_2d_0
1193: requires: triangle !single
1194: args: -petscspace_degree 3 -qorder 3 -convergence
1195: test:
1196: suffix: p3_2d_1
1197: requires: triangle !single
1198: args: -petscspace_degree 3 -qorder 3 -porder 1
1199: test:
1200: suffix: p3_2d_2
1201: requires: triangle !single
1202: args: -petscspace_degree 3 -qorder 3 -porder 2
1203: test:
1204: suffix: p3_2d_3
1205: requires: triangle !single
1206: args: -petscspace_degree 3 -qorder 3 -porder 3
1207: test:
1208: suffix: p3_2d_4
1209: requires: triangle mmg
1210: args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -convergence -conv_refine 0
1211: test:
1212: suffix: p3_2d_5
1213: requires: triangle mmg
1214: args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 1 -conv_refine 0
1215: test:
1216: suffix: p3_2d_6
1217: requires: triangle mmg
1218: args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 3 -conv_refine 0
1220: # 2D Q_3 on a quadrilaterial
1221: test:
1222: suffix: q3_2d_0
1223: requires: !single
1224: args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -convergence
1225: test:
1226: suffix: q3_2d_1
1227: requires: !single
1228: args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 1
1229: test:
1230: suffix: q3_2d_2
1231: requires: !single
1232: args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 2
1233: test:
1234: suffix: q3_2d_3
1235: requires: !single
1236: args: -dm_plex_simplex 0 -petscspace_degree 3 -qorder 3 -porder 3
1238: # 2D P_1disc on a triangle/quadrilateral
1239: test:
1240: suffix: p1d_2d_0
1241: requires: triangle
1242: args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1243: test:
1244: suffix: p1d_2d_1
1245: requires: triangle
1246: args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1247: test:
1248: suffix: p1d_2d_2
1249: requires: triangle
1250: args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1251: test:
1252: suffix: p1d_2d_3
1253: requires: triangle
1254: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1255: filter: sed -e "s/convergence rate at refinement 0: 2/convergence rate at refinement 0: 1.9/g"
1256: test:
1257: suffix: p1d_2d_4
1258: requires: triangle
1259: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1260: test:
1261: suffix: p1d_2d_5
1262: requires: triangle
1263: args: -dm_plex_simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1265: # 2D BDM_1 on a triangle
1266: test:
1267: suffix: bdm1_2d_0
1268: requires: triangle
1269: args: -petscspace_degree 1 -petscdualspace_type bdm \
1270: -num_comp 2 -qorder 1 -convergence
1271: test:
1272: suffix: bdm1_2d_1
1273: requires: triangle
1274: args: -petscspace_degree 1 -petscdualspace_type bdm \
1275: -num_comp 2 -qorder 1 -porder 1
1276: test:
1277: suffix: bdm1_2d_2
1278: requires: triangle
1279: args: -petscspace_degree 1 -petscdualspace_type bdm \
1280: -num_comp 2 -qorder 1 -porder 2
1282: # 2D BDM_1 on a quadrilateral
1283: test:
1284: suffix: bdm1q_2d_0
1285: requires: triangle
1286: args: -petscspace_degree 1 -petscdualspace_type bdm \
1287: -petscdualspace_lagrange_tensor 1 \
1288: -dm_plex_simplex 0 -num_comp 2 -qorder 1 -convergence
1289: test:
1290: suffix: bdm1q_2d_1
1291: requires: triangle
1292: args: -petscspace_degree 1 -petscdualspace_type bdm \
1293: -petscdualspace_lagrange_tensor 1 \
1294: -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 1
1295: test:
1296: suffix: bdm1q_2d_2
1297: requires: triangle
1298: args: -petscspace_degree 1 -petscdualspace_type bdm \
1299: -petscdualspace_lagrange_tensor 1 \
1300: -dm_plex_simplex 0 -num_comp 2 -qorder 1 -porder 2
1302: # Test high order quadrature
1303: test:
1304: suffix: p1_quad_2
1305: requires: triangle
1306: args: -petscspace_degree 1 -qorder 2 -porder 1
1307: test:
1308: suffix: p1_quad_5
1309: requires: triangle
1310: args: -petscspace_degree 1 -qorder 5 -porder 1
1311: test:
1312: suffix: p2_quad_3
1313: requires: triangle
1314: args: -petscspace_degree 2 -qorder 3 -porder 2
1315: test:
1316: suffix: p2_quad_5
1317: requires: triangle
1318: args: -petscspace_degree 2 -qorder 5 -porder 2
1319: test:
1320: suffix: q1_quad_2
1321: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 2 -porder 1
1322: test:
1323: suffix: q1_quad_5
1324: args: -dm_plex_simplex 0 -petscspace_degree 1 -qorder 5 -porder 1
1325: test:
1326: suffix: q2_quad_3
1327: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 3 -porder 1
1328: test:
1329: suffix: q2_quad_5
1330: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 5 -porder 1
1332: # Nonconforming tests
1333: test:
1334: suffix: constraints
1335: args: -dm_coord_space 0 -dm_plex_simplex 0 -petscspace_type tensor -petscspace_degree 1 -qorder 0 -constraints
1336: test:
1337: suffix: nonconforming_tensor_2
1338: nsize: 4
1339: 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
1340: test:
1341: suffix: nonconforming_tensor_3
1342: nsize: 4
1343: 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
1344: test:
1345: suffix: nonconforming_tensor_2_fv
1346: nsize: 4
1347: args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_simplex 0 -num_comp 2
1348: test:
1349: suffix: nonconforming_tensor_3_fv
1350: nsize: 4
1351: 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
1352: test:
1353: suffix: nonconforming_tensor_2_hi
1354: requires: !single
1355: nsize: 4
1356: 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
1357: test:
1358: suffix: nonconforming_tensor_3_hi
1359: requires: !single skip
1360: nsize: 4
1361: 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
1362: test:
1363: suffix: nonconforming_simplex_2
1364: requires: triangle
1365: nsize: 4
1366: 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
1367: test:
1368: suffix: nonconforming_simplex_2_hi
1369: requires: triangle !single
1370: nsize: 4
1371: args: -dist_dm_distribute -test_fe_jacobian -petscpartitioner_type simple -tree -dm_plex_max_projection_height 1 -petscspace_degree 4 -qorder 4
1372: test:
1373: suffix: nonconforming_simplex_2_fv
1374: requires: triangle
1375: nsize: 4
1376: args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -num_comp 2
1377: test:
1378: suffix: nonconforming_simplex_3
1379: requires: ctetgen
1380: nsize: 4
1381: 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
1382: test:
1383: suffix: nonconforming_simplex_3_hi
1384: requires: ctetgen skip
1385: nsize: 4
1386: 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
1387: test:
1388: suffix: nonconforming_simplex_3_fv
1389: requires: ctetgen
1390: nsize: 4
1391: args: -dist_dm_distribute -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -dm_plex_dim 3 -num_comp 3
1393: # 3D WXY on a triangular prism
1394: test:
1395: suffix: wxy_0
1396: args: -dm_plex_reference_cell_domain -dm_plex_cell triangular_prism -qorder 2 -porder 0 \
1397: -petscspace_type sum \
1398: -petscspace_variables 3 \
1399: -petscspace_components 3 \
1400: -petscspace_sum_spaces 2 \
1401: -petscspace_sum_concatenate false \
1402: -sumcomp_0_petscspace_variables 3 \
1403: -sumcomp_0_petscspace_components 3 \
1404: -sumcomp_0_petscspace_degree 1 \
1405: -sumcomp_1_petscspace_variables 3 \
1406: -sumcomp_1_petscspace_components 3 \
1407: -sumcomp_1_petscspace_type wxy \
1408: -petscdualspace_refcell triangular_prism \
1409: -petscdualspace_form_degree 0 \
1410: -petscdualspace_order 1 \
1411: -petscdualspace_components 3
1413: # 2D RT_0 on a triangle
1414: test:
1415: suffix: rt0_2d_tri
1416: requires: triangle
1417: args: -qorder 1 -porder 0 -RT \
1418: -petscspace_type ptrimmed \
1419: -petscspace_components 2 \
1420: -petscspace_ptrimmed_form_degree -1 \
1421: -petscdualspace_order 1 \
1422: -petscdualspace_form_degree -1 \
1423: -petscdualspace_lagrange_trimmed true
1425: # 2D RT_0 on a quadrilateral
1426: test:
1427: suffix: rt0_2d_quad
1428: requires: triangle
1429: args: -dm_plex_simplex 0 -qorder 1 -porder 0 -RT \
1430: -petscspace_degree 1 \
1431: -petscspace_type sum \
1432: -petscspace_variables 2 \
1433: -petscspace_components 2 \
1434: -petscspace_sum_spaces 2 \
1435: -petscspace_sum_concatenate true \
1436: -sumcomp_0_petscspace_variables 2 \
1437: -sumcomp_0_petscspace_type tensor \
1438: -sumcomp_0_petscspace_tensor_spaces 2 \
1439: -sumcomp_0_petscspace_tensor_uniform false \
1440: -sumcomp_0_tensorcomp_0_petscspace_degree 1 \
1441: -sumcomp_0_tensorcomp_1_petscspace_degree 0 \
1442: -sumcomp_1_petscspace_variables 2 \
1443: -sumcomp_1_petscspace_type tensor \
1444: -sumcomp_1_petscspace_tensor_spaces 2 \
1445: -sumcomp_1_petscspace_tensor_uniform false \
1446: -sumcomp_1_tensorcomp_0_petscspace_degree 0 \
1447: -sumcomp_1_tensorcomp_1_petscspace_degree 1 \
1448: -petscdualspace_form_degree -1 \
1449: -petscdualspace_order 1 \
1450: -petscdualspace_lagrange_trimmed true
1452: TEST*/
1454: /*
1455: # 2D Q_2 on a quadrilaterial Plex
1456: test:
1457: suffix: q2_2d_plex_0
1458: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1459: test:
1460: suffix: q2_2d_plex_1
1461: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1462: test:
1463: suffix: q2_2d_plex_2
1464: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1465: test:
1466: suffix: q2_2d_plex_3
1467: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1468: test:
1469: suffix: q2_2d_plex_4
1470: args: -dm_plex_simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1471: test:
1472: suffix: q2_2d_plex_5
1473: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 0 -non_affine_coords
1474: test:
1475: suffix: q2_2d_plex_6
1476: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 1 -non_affine_coords
1477: test:
1478: suffix: q2_2d_plex_7
1479: args: -dm_plex_simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 2 -non_affine_coords
1481: test:
1482: suffix: p1d_2d_6
1483: requires: mmg
1484: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1485: test:
1486: suffix: p1d_2d_7
1487: requires: mmg
1488: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1489: test:
1490: suffix: p1d_2d_8
1491: requires: mmg
1492: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1493: */