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