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