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