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