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