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