Actual source code: ex3.c
petsc-3.10.5 2019-03-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, (PetscObject) user->fe);
423: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
424: if (!isPlex) {
425: PetscSection section;
426: const PetscInt *numDof;
427: PetscInt numComp;
429: PetscFEGetNumComponents(user->fe, &numComp);
430: PetscFEGetNumDof(user->fe, &numDof);
431: DMDACreateSection(dm, &numComp, numDof, NULL, §ion);
432: DMSetSection(dm, section);
433: PetscSectionDestroy(§ion);
434: }
435: if (!user->simplex && user->constraints) {
436: /* test getting local constraint matrix that matches section */
437: PetscSection aSec;
438: IS aIS;
440: DMPlexGetAnchors(dm,&aSec,&aIS);
441: if (aSec) {
442: PetscDS ds;
443: PetscSection cSec, section;
444: PetscInt cStart, cEnd, c, numComp;
445: Mat cMat, mass;
446: Vec local;
447: const PetscInt *anchors;
449: DMGetSection(dm,§ion);
450: /* this creates the matrix and preallocates the matrix structure: we
451: * just have to fill in the values */
452: DMGetDefaultConstraints(dm,&cSec,&cMat);
453: PetscSectionGetChart(cSec,&cStart,&cEnd);
454: ISGetIndices(aIS,&anchors);
455: PetscFEGetNumComponents(user->fe, &numComp);
456: for (c = cStart; c < cEnd; c++) {
457: PetscInt cDof;
459: /* is this point constrained? (does it have an anchor?) */
460: PetscSectionGetDof(aSec,c,&cDof);
461: if (cDof) {
462: PetscInt cOff, a, aDof, aOff, j;
463: if (cDof != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Found %d anchor points: should be just one",cDof);
465: /* find the anchor point */
466: PetscSectionGetOffset(aSec,c,&cOff);
467: a = anchors[cOff];
469: /* find the constrained dofs (row in constraint matrix) */
470: PetscSectionGetDof(cSec,c,&cDof);
471: PetscSectionGetOffset(cSec,c,&cOff);
473: /* find the anchor dofs (column in constraint matrix) */
474: PetscSectionGetDof(section,a,&aDof);
475: PetscSectionGetOffset(section,a,&aOff);
477: if (cDof != aDof) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point and anchor have different number of dofs: %d, %d\n",cDof,aDof);
478: if (cDof % numComp) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point dofs not divisible by field components: %d, %d\n",cDof,numComp);
480: /* put in a simple equality constraint */
481: for (j = 0; j < cDof; j++) {
482: MatSetValue(cMat,cOff+j,aOff+j,1.,INSERT_VALUES);
483: }
484: }
485: }
486: MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
487: MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
488: ISRestoreIndices(aIS,&anchors);
490: /* Now that we have constructed the constraint matrix, any FE matrix
491: * that we construct will apply the constraints during construction */
493: DMCreateMatrix(dm,&mass);
494: /* get a dummy local variable to serve as the solution */
495: DMGetLocalVector(dm,&local);
496: DMGetDS(dm,&ds);
497: /* set the jacobian to be the mass matrix */
498: PetscDSSetJacobian(ds, 0, 0, simple_mass, NULL, NULL, NULL);
499: /* build the mass matrix */
500: DMPlexSNESComputeJacobianFEM(dm,local,mass,mass,NULL);
501: MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
502: MatDestroy(&mass);
503: DMRestoreLocalVector(dm,&local);
504: #if 0
505: {
506: /* compare this to periodicity with DMDA: this is broken right now
507: * because DMCreateMatrix() doesn't respect the default section that I
508: * set */
509: DM dmda;
510: PetscSection section;
511: const PetscInt *numDof;
512: PetscInt numComp;
514: /* periodic x */
515: DMDACreate2d(PetscObjectComm((PetscObject)dm), DM_BOUNDARY_PERIODIC, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, -2, -2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dmda);
516: DMSetFromOptions(dmda);
517: DMSetUp(dmda);
518: DMDASetVertexCoordinates(dmda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
521: PetscFEGetNumComponents(user->fe, &numComp);
522: PetscFEGetNumDof(user->fe, &numDof);
523: DMDACreateSection(dmda, &numComp, numDof, NULL, §ion);
524: DMSetSection(dmda, section);
525: PetscSectionDestroy(§ion);
526: DMCreateMatrix(dmda,&mass);
527: /* there isn't a DMDA equivalent of DMPlexSNESComputeJacobianFEM()
528: * right now, but we can at least verify the nonzero structure */
529: MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
530: MatDestroy(&mass);
531: DMDestroy(&dmda);
532: }
533: #endif
534: }
535: }
536: return(0);
537: }
539: static PetscErrorCode TestFEJacobian(DM dm, AppCtx *user)
540: {
541: PetscBool isPlex;
545: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
546: if (isPlex) {
547: Vec local;
548: const Vec *vecs;
549: Mat E;
550: MatNullSpace sp;
551: PetscBool isNullSpace, hasConst;
552: PetscInt n, i;
553: Vec res, localX, localRes;
554: PetscDS ds;
556: 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);
557: DMGetDS(dm,&ds);
558: PetscDSSetJacobian(ds,0,0,NULL,NULL,NULL,symmetric_gradient_inner_product);
559: DMCreateMatrix(dm,&E);
560: DMGetLocalVector(dm,&local);
561: DMPlexSNESComputeJacobianFEM(dm,local,E,E,NULL);
562: DMPlexCreateRigidBody(dm,&sp);
563: MatNullSpaceGetVecs(sp,&hasConst,&n,&vecs);
564: VecDuplicate(vecs[0],&res);
565: DMCreateLocalVector(dm,&localX);
566: DMCreateLocalVector(dm,&localRes);
567: for (i = 0; i < n; i++) { /* also test via matrix-free Jacobian application */
568: PetscReal resNorm;
570: VecSet(localRes,0.);
571: VecSet(localX,0.);
572: VecSet(local,0.);
573: VecSet(res,0.);
574: DMGlobalToLocalBegin(dm,vecs[i],INSERT_VALUES,localX);
575: DMGlobalToLocalEnd(dm,vecs[i],INSERT_VALUES,localX);
576: DMPlexComputeJacobianAction(dm,NULL,0,0,local,NULL,localX,localRes,NULL);
577: DMLocalToGlobalBegin(dm,localRes,ADD_VALUES,res);
578: DMLocalToGlobalEnd(dm,localRes,ADD_VALUES,res);
579: VecNorm(res,NORM_2,&resNorm);
580: if (resNorm > PETSC_SMALL) {
581: PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient action null space vector %D residual: %E\n",i,resNorm);
582: }
583: }
584: VecDestroy(&localRes);
585: VecDestroy(&localX);
586: VecDestroy(&res);
587: MatNullSpaceTest(sp,E,&isNullSpace);
588: if (isNullSpace) {
589: PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient null space: PASS\n");
590: } else {
591: PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient null space: FAIL\n");
592: }
593: MatNullSpaceDestroy(&sp);
594: MatDestroy(&E);
595: DMRestoreLocalVector(dm,&local);
596: }
597: return(0);
598: }
600: static PetscErrorCode TestInjector(DM dm, AppCtx *user)
601: {
602: DM refTree;
603: PetscMPIInt rank;
607: DMPlexGetReferenceTree(dm,&refTree);
608: if (refTree) {
609: Mat inj;
611: DMPlexComputeInjectorReferenceTree(refTree,&inj);
612: PetscObjectSetName((PetscObject)inj,"Reference Tree Injector");
613: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
614: if (!rank) {
615: MatView(inj,PETSC_VIEWER_STDOUT_SELF);
616: }
617: MatDestroy(&inj);
618: }
619: return(0);
620: }
622: static PetscErrorCode TestFVGrad(DM dm, AppCtx *user)
623: {
624: MPI_Comm comm;
625: DM dmRedist, dmfv, dmgrad, dmCell, refTree;
626: PetscFV fv;
627: PetscInt nvecs, v, cStart, cEnd, cEndInterior;
628: PetscMPIInt size;
629: Vec cellgeom, grad, locGrad;
630: const PetscScalar *cgeom;
631: PetscReal allVecMaxDiff = 0., fvTol = 100. * PETSC_MACHINE_EPSILON;
632: PetscErrorCode ierr;
635: comm = PetscObjectComm((PetscObject)dm);
636: /* duplicate DM, give dup. a FV discretization */
637: DMPlexSetAdjacencyUseCone(dm,PETSC_TRUE);
638: DMPlexSetAdjacencyUseClosure(dm,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, (PetscObject) fv);
658: DMPlexGetReferenceTree(dm,&refTree);
659: if (refTree) {
660: PetscDS ds;
661: DMGetDS(dmfv,&ds);
662: DMSetDS(refTree,ds);
663: }
664: DMPlexSNESGetGradientDM(dmfv, fv, &dmgrad);
665: DMPlexGetHeightStratum(dmfv,0,&cStart,&cEnd);
666: nvecs = user->dim * (user->dim+1) / 2;
667: DMPlexSNESGetGeometryFVM(dmfv,NULL,&cellgeom,NULL);
668: VecGetDM(cellgeom,&dmCell);
669: VecGetArrayRead(cellgeom,&cgeom);
670: DMGetGlobalVector(dmgrad,&grad);
671: DMGetLocalVector(dmgrad,&locGrad);
672: DMPlexGetHybridBounds(dmgrad,&cEndInterior,NULL,NULL,NULL);
673: cEndInterior = (cEndInterior < 0) ? cEnd: cEndInterior;
674: for (v = 0; v < nvecs; v++) {
675: Vec locX;
676: PetscInt c;
677: PetscScalar trueGrad[3][3] = {{0.}};
678: const PetscScalar *gradArray;
679: PetscReal maxDiff, maxDiffGlob;
681: DMGetLocalVector(dmfv,&locX);
682: /* get the local projection of the rigid body mode */
683: for (c = cStart; c < cEnd; c++) {
684: PetscFVCellGeom *cg;
685: PetscScalar cx[3] = {0.,0.,0.};
687: DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
688: if (v < user->dim) {
689: cx[v] = 1.;
690: } else {
691: PetscInt w = v - user->dim;
693: cx[(w + 1) % user->dim] = cg->centroid[(w + 2) % user->dim];
694: cx[(w + 2) % user->dim] = -cg->centroid[(w + 1) % user->dim];
695: }
696: DMPlexVecSetClosure(dmfv,NULL,locX,c,cx,INSERT_ALL_VALUES);
697: }
698: /* TODO: this isn't in any header */
699: DMPlexReconstructGradientsFVM(dmfv,locX,grad);
700: DMGlobalToLocalBegin(dmgrad,grad,INSERT_VALUES,locGrad);
701: DMGlobalToLocalEnd(dmgrad,grad,INSERT_VALUES,locGrad);
702: VecGetArrayRead(locGrad,&gradArray);
703: /* compare computed gradient to exact gradient */
704: if (v >= user->dim) {
705: PetscInt w = v - user->dim;
707: trueGrad[(w + 1) % user->dim][(w + 2) % user->dim] = 1.;
708: trueGrad[(w + 2) % user->dim][(w + 1) % user->dim] = -1.;
709: }
710: maxDiff = 0.;
711: for (c = cStart; c < cEndInterior; c++) {
712: PetscScalar *compGrad;
713: PetscInt i, j, k;
714: PetscReal FrobDiff = 0.;
716: DMPlexPointLocalRead(dmgrad, c, gradArray, &compGrad);
718: for (i = 0, k = 0; i < user->dim; i++) {
719: for (j = 0; j < user->dim; j++, k++) {
720: PetscScalar diff = compGrad[k] - trueGrad[i][j];
721: FrobDiff += PetscRealPart(diff * PetscConj(diff));
722: }
723: }
724: FrobDiff = PetscSqrtReal(FrobDiff);
725: maxDiff = PetscMax(maxDiff,FrobDiff);
726: }
727: MPI_Allreduce(&maxDiff,&maxDiffGlob,1,MPIU_REAL,MPIU_MAX,comm);
728: allVecMaxDiff = PetscMax(allVecMaxDiff,maxDiffGlob);
729: VecRestoreArrayRead(locGrad,&gradArray);
730: DMRestoreLocalVector(dmfv,&locX);
731: }
732: if (allVecMaxDiff < fvTol) {
733: PetscPrintf(PetscObjectComm((PetscObject)dm),"Finite volume gradient reconstruction: PASS\n");
734: } else {
735: PetscPrintf(PetscObjectComm((PetscObject)dm),"Finite volume gradient reconstruction: FAIL at tolerance %g with max difference %g\n",fvTol,allVecMaxDiff);
736: }
737: DMRestoreLocalVector(dmgrad,&locGrad);
738: DMRestoreGlobalVector(dmgrad,&grad);
739: VecRestoreArrayRead(cellgeom,&cgeom);
740: DMDestroy(&dmfv);
741: PetscFVDestroy(&fv);
742: return(0);
743: }
745: static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *),
746: PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *),
747: void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
748: {
749: Vec u;
750: PetscReal n[3] = {1.0, 1.0, 1.0};
754: DMGetGlobalVector(dm, &u);
755: /* Project function into FE function space */
756: DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u);
757: VecViewFromOptions(u, NULL, "-projection_view");
758: /* Compare approximation to exact in L_2 */
759: DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error);
760: DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer);
761: DMRestoreGlobalVector(dm, &u);
762: return(0);
763: }
765: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
766: {
767: PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
768: PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
769: void *exactCtxs[3];
770: MPI_Comm comm;
771: PetscReal error, errorDer, tol = PETSC_SMALL;
772: PetscErrorCode ierr;
775: exactCtxs[0] = user;
776: exactCtxs[1] = user;
777: exactCtxs[2] = user;
778: user->constants[0] = 1.0;
779: user->constants[1] = 2.0;
780: user->constants[2] = 3.0;
781: PetscObjectGetComm((PetscObject)dm, &comm);
782: /* Setup functions to approximate */
783: switch (order) {
784: case 0:
785: exactFuncs[0] = constant;
786: exactFuncDers[0] = constantDer;
787: break;
788: case 1:
789: exactFuncs[0] = linear;
790: exactFuncDers[0] = linearDer;
791: break;
792: case 2:
793: exactFuncs[0] = quadratic;
794: exactFuncDers[0] = quadraticDer;
795: break;
796: case 3:
797: exactFuncs[0] = cubic;
798: exactFuncDers[0] = cubicDer;
799: break;
800: default:
801: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %d order %d", user->dim, order);
802: }
803: ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
804: /* Report result */
805: if (error > tol) {PetscPrintf(comm, "Function tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol,(double) error);}
806: else {PetscPrintf(comm, "Function tests pass for order %D at tolerance %g\n", order, (double)tol);}
807: if (errorDer > tol) {PetscPrintf(comm, "Function tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);}
808: else {PetscPrintf(comm, "Function tests pass for order %D derivatives at tolerance %g\n", order, (double)tol);}
809: return(0);
810: }
812: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
813: {
814: PetscErrorCode (*exactFuncs[1]) (PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
815: PetscErrorCode (*exactFuncDers[1]) (PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx);
816: PetscReal n[3] = {1.0, 1.0, 1.0};
817: void *exactCtxs[3];
818: DM rdm, idm, fdm;
819: Mat Interp;
820: Vec iu, fu, scaling;
821: MPI_Comm comm;
822: PetscInt dim = user->dim;
823: PetscReal error, errorDer, tol = PETSC_SMALL;
824: PetscBool isPlex, isDA;
825: PetscErrorCode ierr;
828: exactCtxs[0] = user;
829: exactCtxs[1] = user;
830: exactCtxs[2] = user;
831: user->constants[0] = 1.0;
832: user->constants[1] = 2.0;
833: user->constants[2] = 3.0;
834: PetscObjectGetComm((PetscObject)dm,&comm);
835: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
836: PetscObjectTypeCompare((PetscObject) dm, DMDA, &isDA);
837: DMRefine(dm, comm, &rdm);
838: DMSetCoarseDM(rdm, dm);
839: DMPlexSetRegularRefinement(rdm, user->convRefine);
840: if (user->tree && isPlex) {
841: DM refTree;
842: DMPlexGetReferenceTree(dm,&refTree);
843: DMPlexSetReferenceTree(rdm,refTree);
844: }
845: if (!user->simplex && !user->constraints) {DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
846: SetupSection(rdm, user);
847: /* Setup functions to approximate */
848: switch (order) {
849: case 0:
850: exactFuncs[0] = constant;
851: exactFuncDers[0] = constantDer;
852: break;
853: case 1:
854: exactFuncs[0] = linear;
855: exactFuncDers[0] = linearDer;
856: break;
857: case 2:
858: exactFuncs[0] = quadratic;
859: exactFuncDers[0] = quadraticDer;
860: break;
861: case 3:
862: exactFuncs[0] = cubic;
863: exactFuncDers[0] = cubicDer;
864: break;
865: default:
866: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %D order %D", dim, order);
867: }
868: idm = checkRestrict ? rdm : dm;
869: fdm = checkRestrict ? dm : rdm;
870: DMGetGlobalVector(idm, &iu);
871: DMGetGlobalVector(fdm, &fu);
872: DMSetApplicationContext(dm, user);
873: DMSetApplicationContext(rdm, user);
874: DMCreateInterpolation(dm, rdm, &Interp, &scaling);
875: /* Project function into initial FE function space */
876: DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu);
877: /* Interpolate function into final FE function space */
878: if (checkRestrict) {MatRestrict(Interp, iu, fu);VecPointwiseMult(fu, scaling, fu);}
879: else {MatInterpolate(Interp, iu, fu);}
880: /* Compare approximation to exact in L_2 */
881: DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error);
882: DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer);
883: /* Report result */
884: if (error > tol) {PetscPrintf(comm, "Interpolation tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol, (double)error);}
885: else {PetscPrintf(comm, "Interpolation tests pass for order %D at tolerance %g\n", order, (double)tol);}
886: if (errorDer > tol) {PetscPrintf(comm, "Interpolation tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);}
887: else {PetscPrintf(comm, "Interpolation tests pass for order %D derivatives at tolerance %g\n", order, (double)tol);}
888: DMRestoreGlobalVector(idm, &iu);
889: DMRestoreGlobalVector(fdm, &fu);
890: MatDestroy(&Interp);
891: VecDestroy(&scaling);
892: DMDestroy(&rdm);
893: return(0);
894: }
896: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
897: {
898: DM odm = dm, rdm = NULL, cdm = NULL;
899: PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {trig};
900: PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx) = {trigDer};
901: void *exactCtxs[3];
902: PetscInt r, c, cStart, cEnd;
903: PetscReal errorOld, errorDerOld, error, errorDer, rel, len, lenOld;
904: double p;
905: PetscErrorCode ierr;
908: if (!user->convergence) return(0);
909: exactCtxs[0] = user;
910: exactCtxs[1] = user;
911: exactCtxs[2] = user;
912: PetscObjectReference((PetscObject) odm);
913: if (!user->convRefine) {
914: for (r = 0; r < Nr; ++r) {
915: DMRefine(odm, PetscObjectComm((PetscObject) dm), &rdm);
916: DMDestroy(&odm);
917: odm = rdm;
918: }
919: SetupSection(odm, user);
920: }
921: ComputeError(odm, exactFuncs, exactFuncDers, exactCtxs, &errorOld, &errorDerOld, user);
922: if (user->convRefine) {
923: for (r = 0; r < Nr; ++r) {
924: DMRefine(odm, PetscObjectComm((PetscObject) dm), &rdm);
925: if (!user->simplex && user->useDA) {DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
926: SetupSection(rdm, user);
927: ComputeError(rdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
928: p = PetscLog2Real(errorOld/error);
929: PetscPrintf(PetscObjectComm((PetscObject) dm), "Function convergence rate at refinement %D: %.2g\n", r, (double)p);
930: p = PetscLog2Real(errorDerOld/errorDer);
931: PetscPrintf(PetscObjectComm((PetscObject) dm), "Derivative convergence rate at refinement %D: %.2g\n", r, (double)p);
932: DMDestroy(&odm);
933: odm = rdm;
934: errorOld = error;
935: errorDerOld = errorDer;
936: }
937: } else {
938: /* ComputeLongestEdge(dm, &lenOld); */
939: DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd);
940: lenOld = cEnd - cStart;
941: for (c = 0; c < Nr; ++c) {
942: DMCoarsen(odm, PetscObjectComm((PetscObject) dm), &cdm);
943: if (!user->simplex && user->useDA) {DMDASetVertexCoordinates(cdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
944: SetupSection(cdm, user);
945: ComputeError(cdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
946: /* ComputeLongestEdge(cdm, &len); */
947: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
948: len = cEnd - cStart;
949: rel = error/errorOld;
950: p = PetscLogReal(rel) / PetscLogReal(lenOld / len);
951: PetscPrintf(PetscObjectComm((PetscObject) dm), "Function convergence rate at coarsening %D: %.2g\n", c, (double)p);
952: rel = errorDer/errorDerOld;
953: p = PetscLogReal(rel) / PetscLogReal(lenOld / len);
954: PetscPrintf(PetscObjectComm((PetscObject) dm), "Derivative convergence rate at coarsening %D: %.2g\n", c, (double)p);
955: DMDestroy(&odm);
956: odm = cdm;
957: errorOld = error;
958: errorDerOld = errorDer;
959: lenOld = len;
960: }
961: }
962: DMDestroy(&odm);
963: return(0);
964: }
966: int main(int argc, char **argv)
967: {
968: DM dm;
969: AppCtx user; /* user-defined work context */
972: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
973: ProcessOptions(PETSC_COMM_WORLD, &user);
974: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
975: PetscFECreateDefault(PETSC_COMM_WORLD, user.dim, user.numComponents, user.simplex, NULL, user.qorder, &user.fe);
976: SetupSection(dm, &user);
977: if (user.testFEjacobian) {TestFEJacobian(dm, &user);}
978: if (user.testFVgrad) {TestFVGrad(dm, &user);}
979: if (user.testInjector) {TestInjector(dm, &user);}
980: CheckFunctions(dm, user.porder, &user);
981: if (user.dim == 2 && user.simplex == PETSC_TRUE && user.tree == PETSC_FALSE) {
982: CheckInterpolation(dm, PETSC_FALSE, user.porder, &user);
983: CheckInterpolation(dm, PETSC_TRUE, user.porder, &user);
984: }
985: CheckConvergence(dm, 3, &user);
986: PetscFEDestroy(&user.fe);
987: DMDestroy(&dm);
988: PetscFinalize();
989: return ierr;
990: }
992: /*TEST
994: test:
995: suffix: 1
996: requires: triangle
998: # 2D P_1 on a triangle
999: test:
1000: suffix: p1_2d_0
1001: requires: triangle
1002: args: -petscspace_degree 1 -qorder 1 -convergence
1003: test:
1004: suffix: p1_2d_1
1005: requires: triangle
1006: args: -petscspace_degree 1 -qorder 1 -porder 1
1007: test:
1008: suffix: p1_2d_2
1009: requires: triangle
1010: args: -petscspace_degree 1 -qorder 1 -porder 2
1011: test:
1012: suffix: p1_2d_3
1013: requires: triangle pragmatic
1014: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1015: test:
1016: suffix: p1_2d_4
1017: requires: triangle pragmatic
1018: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1019: test:
1020: suffix: p1_2d_5
1021: requires: triangle pragmatic
1022: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1024: # 3D P_1 on a tetrahedron
1025: test:
1026: suffix: p1_3d_0
1027: requires: ctetgen
1028: args: -dim 3 -petscspace_degree 1 -qorder 1 -convergence
1029: test:
1030: suffix: p1_3d_1
1031: requires: ctetgen
1032: args: -dim 3 -petscspace_degree 1 -qorder 1 -porder 1
1033: test:
1034: suffix: p1_3d_2
1035: requires: ctetgen
1036: args: -dim 3 -petscspace_degree 1 -qorder 1 -porder 2
1037: test:
1038: suffix: p1_3d_3
1039: requires: ctetgen pragmatic
1040: args: -dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1041: test:
1042: suffix: p1_3d_4
1043: requires: ctetgen pragmatic
1044: args: -dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1045: test:
1046: suffix: p1_3d_5
1047: requires: ctetgen pragmatic
1048: args: -dim 3 -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1050: # 2D P_2 on a triangle
1051: test:
1052: suffix: p2_2d_0
1053: requires: triangle
1054: args: -petscspace_degree 2 -qorder 2 -convergence
1055: test:
1056: suffix: p2_2d_1
1057: requires: triangle
1058: args: -petscspace_degree 2 -qorder 2 -porder 1
1059: test:
1060: suffix: p2_2d_2
1061: requires: triangle
1062: args: -petscspace_degree 2 -qorder 2 -porder 2
1063: test:
1064: suffix: p2_2d_3
1065: requires: triangle pragmatic
1066: args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1067: test:
1068: suffix: p2_2d_4
1069: requires: triangle pragmatic
1070: args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1071: test:
1072: suffix: p2_2d_5
1073: requires: triangle pragmatic
1074: args: -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0
1076: # 3D P_2 on a tetrahedron
1077: test:
1078: suffix: p2_3d_0
1079: requires: ctetgen
1080: args: -dim 3 -petscspace_degree 2 -qorder 2 -convergence
1081: test:
1082: suffix: p2_3d_1
1083: requires: ctetgen
1084: args: -dim 3 -petscspace_degree 2 -qorder 2 -porder 1
1085: test:
1086: suffix: p2_3d_2
1087: requires: ctetgen
1088: args: -dim 3 -petscspace_degree 2 -qorder 2 -porder 2
1089: test:
1090: suffix: p2_3d_3
1091: requires: ctetgen pragmatic
1092: args: -dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -convergence -conv_refine 0
1093: test:
1094: suffix: p2_3d_4
1095: requires: ctetgen pragmatic
1096: args: -dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 1 -conv_refine 0
1097: test:
1098: suffix: p2_3d_5
1099: requires: ctetgen pragmatic
1100: args: -dim 3 -petscspace_degree 2 -qorder 2 -dm_plex_hash_location -porder 2 -conv_refine 0
1102: # 2D Q_1 on a quadrilaterial DA
1103: test:
1104: suffix: q1_2d_da_0
1105: requires: mpi_type_get_envelope broken
1106: args: -simplex 0 -petscspace_degree 1 -qorder 1 -convergence
1107: test:
1108: suffix: q1_2d_da_1
1109: requires: mpi_type_get_envelope broken
1110: args: -simplex 0 -petscspace_degree 1 -qorder 1 -porder 1
1111: test:
1112: suffix: q1_2d_da_2
1113: requires: mpi_type_get_envelope broken
1114: args: -simplex 0 -petscspace_degree 1 -qorder 1 -porder 2
1116: # 2D Q_1 on a quadrilaterial Plex
1117: test:
1118: suffix: q1_2d_plex_0
1119: args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 1 -convergence
1120: test:
1121: suffix: q1_2d_plex_1
1122: args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 1 -porder 1
1123: test:
1124: suffix: q1_2d_plex_2
1125: args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 1 -porder 2
1126: test:
1127: suffix: q1_2d_plex_3
1128: args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 1 -porder 1 -shear_coords
1129: test:
1130: suffix: q1_2d_plex_4
1131: args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 1 -porder 2 -shear_coords
1132: test:
1133: suffix: q1_2d_plex_5
1134: args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 0 -non_affine_coords -convergence
1135: test:
1136: suffix: q1_2d_plex_6
1137: args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 1 -non_affine_coords -convergence
1138: test:
1139: suffix: q1_2d_plex_7
1140: args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscspace_type tensor -qorder 1 -porder 2 -non_affine_coords -convergence
1142: # 2D Q_2 on a quadrilaterial
1143: test:
1144: suffix: q2_2d_plex_0
1145: requires: mpi_type_get_envelope
1146: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1147: test:
1148: suffix: q2_2d_plex_1
1149: requires: mpi_type_get_envelope
1150: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1151: test:
1152: suffix: q2_2d_plex_2
1153: requires: mpi_type_get_envelope
1154: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1155: test:
1156: suffix: q2_2d_plex_3
1157: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1158: test:
1159: suffix: q2_2d_plex_4
1160: requires: mpi_type_get_envelope
1161: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1162: test:
1163: suffix: q2_2d_plex_5
1164: requires: mpi_type_get_envelope
1165: args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 0 -non_affine_coords -convergence
1166: test:
1167: suffix: q2_2d_plex_6
1168: requires: mpi_type_get_envelope
1169: args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 1 -non_affine_coords -convergence
1170: test:
1171: suffix: q2_2d_plex_7
1172: requires: mpi_type_get_envelope
1173: args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_type tensor -qorder 2 -porder 2 -non_affine_coords -convergence
1176: # 2D P_3 on a triangle
1177: test:
1178: suffix: p3_2d_0
1179: requires: triangle !single
1180: args: -petscspace_degree 3 -qorder 3 -convergence
1181: test:
1182: suffix: p3_2d_1
1183: requires: triangle !single
1184: args: -petscspace_degree 3 -qorder 3 -porder 1
1185: test:
1186: suffix: p3_2d_2
1187: requires: triangle !single
1188: args: -petscspace_degree 3 -qorder 3 -porder 2
1189: test:
1190: suffix: p3_2d_3
1191: requires: triangle !single
1192: args: -petscspace_degree 3 -qorder 3 -porder 3
1193: test:
1194: suffix: p3_2d_4
1195: requires: triangle pragmatic
1196: args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -convergence -conv_refine 0
1197: test:
1198: suffix: p3_2d_5
1199: requires: triangle pragmatic
1200: args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 1 -conv_refine 0
1201: test:
1202: suffix: p3_2d_6
1203: requires: triangle pragmatic
1204: args: -petscspace_degree 3 -qorder 3 -dm_plex_hash_location -porder 3 -conv_refine 0
1206: # 2D Q_3 on a quadrilaterial
1207: test:
1208: suffix: q3_2d_0
1209: requires: mpi_type_get_envelope !single
1210: args: -use_da 0 -simplex 0 -petscspace_degree 3 -qorder 3 -convergence
1211: test:
1212: suffix: q3_2d_1
1213: requires: mpi_type_get_envelope !single
1214: args: -use_da 0 -simplex 0 -petscspace_degree 3 -qorder 3 -porder 1
1215: test:
1216: suffix: q3_2d_2
1217: requires: mpi_type_get_envelope !single
1218: args: -use_da 0 -simplex 0 -petscspace_degree 3 -qorder 3 -porder 2
1219: test:
1220: suffix: q3_2d_3
1221: requires: mpi_type_get_envelope !single
1222: args: -use_da 0 -simplex 0 -petscspace_degree 3 -qorder 3 -porder 3
1224: # 2D P_1disc on a triangle/quadrilateral
1225: test:
1226: suffix: p1d_2d_0
1227: requires: triangle
1228: args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1229: test:
1230: suffix: p1d_2d_1
1231: requires: triangle
1232: args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1233: test:
1234: suffix: p1d_2d_2
1235: requires: triangle
1236: args: -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1237: test:
1238: suffix: p1d_2d_3
1239: requires: triangle
1240: args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -convergence
1241: filter: sed -e "s/convergence rate at refinement 0: 2/convergence rate at refinement 0: 1.9/g"
1242: test:
1243: suffix: p1d_2d_4
1244: requires: triangle
1245: args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 1
1246: test:
1247: suffix: p1d_2d_5
1248: requires: triangle
1249: args: -use_da 0 -simplex 0 -petscspace_degree 1 -petscdualspace_lagrange_continuity 0 -qorder 1 -porder 2
1251: # Test high order quadrature
1252: test:
1253: suffix: p1_quad_2
1254: requires: triangle
1255: args: -petscspace_degree 1 -qorder 2 -porder 1
1256: test:
1257: suffix: p1_quad_5
1258: requires: triangle
1259: args: -petscspace_degree 1 -qorder 5 -porder 1
1260: test:
1261: suffix: p2_quad_3
1262: requires: triangle
1263: args: -petscspace_degree 2 -qorder 3 -porder 2
1264: test:
1265: suffix: p2_quad_5
1266: requires: triangle
1267: args: -petscspace_degree 2 -qorder 5 -porder 2
1268: test:
1269: suffix: q1_quad_2
1270: requires: mpi_type_get_envelope
1271: args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 2 -porder 1
1272: test:
1273: suffix: q1_quad_5
1274: requires: mpi_type_get_envelope
1275: args: -use_da 0 -simplex 0 -petscspace_degree 1 -qorder 5 -porder 1
1276: test:
1277: suffix: q2_quad_3
1278: requires: mpi_type_get_envelope
1279: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 3 -porder 1
1280: test:
1281: suffix: q2_quad_5
1282: requires: mpi_type_get_envelope
1283: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 5 -porder 1
1286: # Nonconforming tests
1287: test:
1288: suffix: constraints
1289: args: -simplex 0 -petscspace_type tensor -petscspace_degree 1 -qorder 0 -constraints
1290: test:
1291: suffix: nonconforming_tensor_2
1292: nsize: 4
1293: 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
1294: test:
1295: suffix: nonconforming_tensor_3
1296: nsize: 4
1297: 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
1298: test:
1299: suffix: nonconforming_tensor_2_fv
1300: nsize: 4
1301: args: -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -simplex 0 -dim 2 -num_comp 2
1302: test:
1303: suffix: nonconforming_tensor_3_fv
1304: nsize: 4
1305: args: -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -simplex 0 -dim 3 -num_comp 3
1306: test:
1307: suffix: nonconforming_tensor_2_hi
1308: requires: !single
1309: nsize: 4
1310: 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
1311: test:
1312: suffix: nonconforming_tensor_3_hi
1313: requires: !single skip
1314: nsize: 4
1315: 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
1316: test:
1317: suffix: nonconforming_simplex_2
1318: requires: triangle
1319: nsize: 4
1320: 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
1321: test:
1322: suffix: nonconforming_simplex_2_hi
1323: requires: triangle !single
1324: nsize: 4
1325: args: -test_fe_jacobian -petscpartitioner_type simple -tree -simplex 1 -dim 2 -dm_plex_max_projection_height 1 -petscspace_degree 4 -qorder 4
1326: test:
1327: suffix: nonconforming_simplex_2_fv
1328: requires: triangle
1329: nsize: 4
1330: args: -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -simplex 1 -dim 2 -num_comp 2
1331: test:
1332: suffix: nonconforming_simplex_3
1333: requires: ctetgen
1334: nsize: 4
1335: 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
1336: test:
1337: suffix: nonconforming_simplex_3_hi
1338: requires: ctetgen skip
1339: nsize: 4
1340: args: -test_fe_jacobian -petscpartitioner_type simple -tree -simplex 1 -dim 3 -dm_plex_max_projection_height 2 -petscspace_degree 4 -qorder 4
1341: test:
1342: suffix: nonconforming_simplex_3_fv
1343: requires: ctetgen
1344: nsize: 4
1345: args: -test_fv_grad -test_injector -petsclimiter_type none -petscpartitioner_type simple -tree -simplex 1 -dim 3 -num_comp 3
1347: TEST*/
1349: /*
1350: # 2D Q_2 on a quadrilaterial Plex
1351: test:
1352: suffix: q2_2d_plex_0
1353: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -convergence
1354: test:
1355: suffix: q2_2d_plex_1
1356: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 1
1357: test:
1358: suffix: q2_2d_plex_2
1359: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 2
1360: test:
1361: suffix: q2_2d_plex_3
1362: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 1 -shear_coords
1363: test:
1364: suffix: q2_2d_plex_4
1365: args: -use_da 0 -simplex 0 -petscspace_degree 2 -qorder 2 -porder 2 -shear_coords
1366: test:
1367: suffix: q2_2d_plex_5
1368: args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 0 -non_affine_coords
1369: test:
1370: suffix: q2_2d_plex_6
1371: args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 1 -non_affine_coords
1372: test:
1373: suffix: q2_2d_plex_7
1374: args: -use_da 0 -simplex 0 -petscspace_degree 2 -petscspace_poly_tensor 1 -qorder 2 -porder 2 -non_affine_coords
1376: test:
1377: suffix: p1d_2d_6
1378: requires: pragmatic
1379: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -convergence -conv_refine 0
1380: test:
1381: suffix: p1d_2d_7
1382: requires: pragmatic
1383: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 1 -conv_refine 0
1384: test:
1385: suffix: p1d_2d_8
1386: requires: pragmatic
1387: args: -petscspace_degree 1 -qorder 1 -dm_plex_hash_location -porder 2 -conv_refine 0
1388: */