Actual source code: ex3.c
petsc-3.6.1 2015-08-06
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>
10: typedef struct {
11: PetscInt debug; /* The debugging level */
12: /* Domain and mesh definition */
13: PetscInt dim; /* The topological mesh dimension */
14: PetscBool simplex; /* Flag for simplex or tensor product mesh */
15: PetscBool useDA; /* Flag DMDA tensor product mesh */
16: PetscBool interpolate; /* Generate intermediate mesh elements */
17: PetscReal refinementLimit; /* The largest allowable cell volume */
18: /* Element definition */
19: PetscInt qorder; /* Order of the quadrature */
20: PetscInt numComponents; /* Number of field components */
21: PetscFE fe; /* The finite element */
22: /* Testing space */
23: PetscInt porder; /* Order of polynomials to test */
24: PetscBool convergence; /* Test for order of convergence */
25: PetscBool constraints; /* Test local constraints */
26: PetscBool tree; /* Test tree routines */
27: PetscInt treeCell; /* Cell to refine in tree test */
28: PetscReal constants[3]; /* Constant values for each dimension */
29: } AppCtx;
31: /* u = 1 */
32: PetscErrorCode constant(PetscInt dim, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
33: {
34: AppCtx *user = (AppCtx *) ctx;
35: PetscInt d;
36: for (d = 0; d < user->dim; ++d) u[d] = user->constants[d];
37: return 0;
38: }
39: PetscErrorCode constantDer(PetscInt dim, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
40: {
41: AppCtx *user = (AppCtx *) ctx;
42: PetscInt d;
43: for (d = 0; d < user->dim; ++d) u[d] = 0.0;
44: return 0;
45: }
47: /* u = x */
48: PetscErrorCode linear(PetscInt dim, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
49: {
50: AppCtx *user = (AppCtx *) ctx;
51: PetscInt d;
52: for (d = 0; d < user->dim; ++d) u[d] = coords[d];
53: return 0;
54: }
55: PetscErrorCode linearDer(PetscInt dim, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
56: {
57: AppCtx *user = (AppCtx *) ctx;
58: PetscInt d, e;
59: for (d = 0; d < user->dim; ++d) {
60: u[d] = 0.0;
61: for (e = 0; e < user->dim; ++e) u[d] += (d == e ? 1.0 : 0.0) * n[e];
62: }
63: return 0;
64: }
66: /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */
67: PetscErrorCode quadratic(PetscInt dim, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
68: {
69: AppCtx *user = (AppCtx *) ctx;
70: if (user->dim > 2) {u[0] = coords[0]*coords[1]; u[1] = coords[1]*coords[2]; u[2] = coords[2]*coords[0];}
71: else if (user->dim > 1) {u[0] = coords[0]*coords[0]; u[1] = coords[0]*coords[1];}
72: else if (user->dim > 0) {u[0] = coords[0]*coords[0];}
73: return 0;
74: }
75: PetscErrorCode quadraticDer(PetscInt dim, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
76: {
77: AppCtx *user = (AppCtx *) ctx;
78: if (user->dim > 2) {u[0] = coords[1]*n[0] + coords[0]*n[1]; u[1] = coords[2]*n[1] + coords[1]*n[2]; u[2] = coords[2]*n[0] + coords[0]*n[2];}
79: else if (user->dim > 1) {u[0] = 2.0*coords[0]*n[0]; u[1] = coords[1]*n[0] + coords[0]*n[1];}
80: else if (user->dim > 0) {u[0] = 2.0*coords[0]*n[0];}
81: return 0;
82: }
84: /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
85: PetscErrorCode cubic(PetscInt dim, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
86: {
87: AppCtx *user = (AppCtx *) ctx;
88: if (user->dim > 2) {u[0] = coords[0]*coords[0]*coords[1]; u[1] = coords[1]*coords[1]*coords[2]; u[2] = coords[2]*coords[2]*coords[0];}
89: else if (user->dim > 1) {u[0] = coords[0]*coords[0]*coords[0]; u[1] = coords[0]*coords[0]*coords[1];}
90: else if (user->dim > 0) {u[0] = coords[0]*coords[0]*coords[0];}
91: return 0;
92: }
93: PetscErrorCode cubicDer(PetscInt dim, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
94: {
95: AppCtx *user = (AppCtx *) ctx;
96: if (user->dim > 2) {u[0] = 2.0*coords[0]*coords[1]*n[0] + coords[0]*coords[0]*n[1]; u[1] = 2.0*coords[1]*coords[2]*n[1] + coords[1]*coords[1]*n[2]; u[2] = 2.0*coords[2]*coords[0]*n[2] + coords[2]*coords[2]*n[0];}
97: else if (user->dim > 1) {u[0] = 3.0*coords[0]*coords[0]*n[0]; u[1] = 2.0*coords[0]*coords[1]*n[0] + coords[0]*coords[0]*n[1];}
98: else if (user->dim > 0) {u[0] = 3.0*coords[0]*coords[0]*n[0];}
99: return 0;
100: }
102: /* u = sin(x) */
103: PetscErrorCode trig(PetscInt dim, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
104: {
105: AppCtx *user = (AppCtx *) ctx;
106: PetscInt d;
107: for (d = 0; d < user->dim; ++d) u[d] = tanh(coords[d] - 0.5);
108: return 0;
109: }
110: PetscErrorCode trigDer(PetscInt dim, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
111: {
112: AppCtx *user = (AppCtx *) ctx;
113: PetscInt d;
114: for (d = 0; d < user->dim; ++d) u[d] = 1.0/PetscSqr(cosh(coords[d] - 0.5)) * n[d];
115: return 0;
116: }
120: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
121: {
125: options->debug = 0;
126: options->dim = 2;
127: options->simplex = PETSC_TRUE;
128: options->useDA = PETSC_TRUE;
129: options->interpolate = PETSC_TRUE;
130: options->refinementLimit = 0.0;
131: options->qorder = 0;
132: options->numComponents = 1;
133: options->porder = 0;
134: options->convergence = PETSC_FALSE;
135: options->constraints = PETSC_FALSE;
136: options->tree = PETSC_FALSE;
137: options->treeCell = 0;
139: PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
140: PetscOptionsInt("-debug", "The debugging level", "ex3.c", options->debug, &options->debug, NULL);
141: PetscOptionsInt("-dim", "The topological mesh dimension", "ex3.c", options->dim, &options->dim, NULL);
142: PetscOptionsBool("-simplex", "Flag for simplices or hexhedra", "ex3.c", options->simplex, &options->simplex, NULL);
143: PetscOptionsBool("-use_da", "Flag for DMDA mesh", "ex3.c", options->useDA, &options->useDA, NULL);
144: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex3.c", options->interpolate, &options->interpolate, NULL);
145: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex3.c", options->refinementLimit, &options->refinementLimit, NULL);
146: PetscOptionsInt("-qorder", "The quadrature order", "ex3.c", options->qorder, &options->qorder, NULL);
147: PetscOptionsInt("-num_comp", "The number of field components", "ex3.c", options->numComponents, &options->numComponents, NULL);
148: PetscOptionsInt("-porder", "The order of polynomials to test", "ex3.c", options->porder, &options->porder, NULL);
149: PetscOptionsBool("-convergence", "Check the convergence rate", "ex3.c", options->convergence, &options->convergence, NULL);
150: PetscOptionsBool("-constraints", "Test local constraints (serial only)", "ex3.c", options->constraints, &options->constraints, NULL);
151: PetscOptionsBool("-tree", "Test tree routines", "ex3.c", options->tree, &options->tree, NULL);
152: PetscOptionsInt("-tree_cell", "cell to refine in tree test", "ex3.c", options->treeCell, &options->treeCell, NULL);
153: PetscOptionsEnd();
155: return(0);
156: }
160: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
161: {
162: PetscInt dim = user->dim;
163: PetscBool interpolate = user->interpolate;
164: PetscReal refinementLimit = user->refinementLimit;
165: PetscBool isPlex;
169: if (user->simplex) {
170: DM refinedMesh = NULL;
172: DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
173: /* Refine mesh using a volume constraint */
174: DMPlexSetRefinementLimit(*dm, refinementLimit);
175: DMRefine(*dm, comm, &refinedMesh);
176: if (refinedMesh) {
177: DMDestroy(dm);
178: *dm = refinedMesh;
179: }
180: DMPlexSetRefinementUniform(*dm, PETSC_TRUE);
181: } else {
182: if (user->constraints || user->tree || !user->useDA) {
183: PetscInt cells[3] = {2, 2, 2};
185: PetscOptionsGetInt(NULL,"-da_grid_x",&cells[0],NULL);
186: PetscOptionsGetInt(NULL,"-da_grid_y",&cells[1],NULL);
187: PetscOptionsGetInt(NULL,"-da_grid_z",&cells[2],NULL);
188: DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);
189: } else {
190: switch (user->dim) {
191: case 2:
192: DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, -2, -2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, dm);
193: DMDASetVertexCoordinates(*dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
194: break;
195: default:
196: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot create structured mesh of dimension %d", dim);
197: }
198: PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");
199: }
200: }
201: PetscObjectTypeCompare((PetscObject)*dm,DMPLEX,&isPlex);
202: if (isPlex) {
203: DM distributedMesh = NULL;
204: if (user->tree) {
205: DM refTree;
206: DM ncdm = NULL;
208: DMPlexCreateDefaultReferenceTree(comm,user->dim,user->simplex,&refTree);
209: DMPlexSetReferenceTree(*dm,refTree);
210: DMDestroy(&refTree);
211: DMPlexTreeRefineCell(*dm,user->treeCell,&ncdm);
212: if (ncdm) {
213: DMDestroy(dm);
214: *dm = ncdm;
215: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
216: }
217: }
218: else {
219: DMPlexSetRefinementUniform(*dm, PETSC_TRUE);
220: }
221: /* Distribute mesh over processes */
222: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
223: if (distributedMesh) {
224: DMDestroy(dm);
225: *dm = distributedMesh;
226: }
227: if (user->simplex) {
228: PetscObjectSetName((PetscObject) *dm, "Simplicial Mesh");
229: }
230: else {
231: PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");
232: }
233: }
234: DMSetFromOptions(*dm);
235: DMViewFromOptions(*dm,NULL,"-dm_view");
236: return(0);
237: }
241: static void simple_mass(PetscInt dim, PetscInt Nf, PetscInt NfAux,
242: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
243: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
244: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
245: {
246: PetscInt d, e;
247: for (d = 0, e = 0; d < dim; d++, e+=dim+1) {
248: g0[e] = 1.;
249: }
250: }
252: /* < \nabla v, 1/2(\nabla u + {\nabla u}^T) > */
255: static void symmetric_gradient_inner_product(PetscInt dim, PetscInt Nf, PetscInt NfAux,
256: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
257: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
258: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar C[])
259: {
260: PetscInt compI, compJ, d, e;
262: for (compI = 0; compI < dim; ++compI) {
263: for (compJ = 0; compJ < dim; ++compJ) {
264: for (d = 0; d < dim; ++d) {
265: for (e = 0; e < dim; e++) {
266: if (d == e && d == compI && d == compJ) {
267: C[((compI*dim+compJ)*dim+d)*dim+e] = 1.0;
268: }
269: else if ((d == compJ && e == compI) || (d == e && compI == compJ)) {
270: C[((compI*dim+compJ)*dim+d)*dim+e] = 0.5;
271: }
272: else {
273: C[((compI*dim+compJ)*dim+d)*dim+e] = 0.0;
274: }
275: }
276: }
277: }
278: }
279: }
283: static PetscErrorCode SetupSection(DM dm, AppCtx *user)
284: {
285: PetscBool isPlex;
289: if (!user->simplex && user->constraints) {
290: /* test local constraints */
291: DM coordDM;
292: PetscInt fStart, fEnd, f, vStart, vEnd, v;
293: PetscInt edgesx = 2, vertsx;
294: PetscInt edgesy = 2, vertsy;
295: PetscMPIInt size;
296: PetscInt numConst;
297: PetscSection aSec;
298: PetscInt *anchors;
299: PetscInt offset;
300: IS aIS;
301: MPI_Comm comm = PetscObjectComm((PetscObject)dm);
303: MPI_Comm_size(comm,&size);
304: if (size > 1) SETERRQ(comm,PETSC_ERR_SUP,"Local constraint test can only be performed in serial");
306: /* we are going to test constraints by using them to enforce periodicity
307: * in one direction, and comparing to the existing method of enforcing
308: * periodicity */
310: /* first create the coordinate section so that it does not clone the
311: * constraints */
312: DMGetCoordinateDM(dm,&coordDM);
314: /* create the constrained-to-anchor section */
315: DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
316: DMPlexGetDepthStratum(dm,1,&fStart,&fEnd);
317: PetscSectionCreate(PETSC_COMM_SELF,&aSec);
318: PetscSectionSetChart(aSec,PetscMin(fStart,vStart),PetscMax(fEnd,vEnd));
320: /* define the constraints */
321: PetscOptionsGetInt(NULL,"-da_grid_x",&edgesx,NULL);
322: PetscOptionsGetInt(NULL,"-da_grid_y",&edgesy,NULL);
323: vertsx = edgesx + 1;
324: vertsy = edgesy + 1;
325: numConst = vertsy + edgesy;
326: PetscMalloc1(numConst,&anchors);
327: offset = 0;
328: for (v = vStart + edgesx; v < vEnd; v+= vertsx) {
329: PetscSectionSetDof(aSec,v,1);
330: anchors[offset++] = v - edgesx;
331: }
332: for (f = fStart + edgesx * vertsy + edgesx * edgesy; f < fEnd; f++) {
333: PetscSectionSetDof(aSec,f,1);
334: anchors[offset++] = f - edgesx * edgesy;
335: }
336: PetscSectionSetUp(aSec);
337: ISCreateGeneral(PETSC_COMM_SELF,numConst,anchors,PETSC_OWN_POINTER,&aIS);
339: DMPlexSetAnchors(dm,aSec,aIS);
340: PetscSectionDestroy(&aSec);
341: ISDestroy(&aIS);
342: }
343: DMSetNumFields(dm, 1);
344: DMSetField(dm, 0, (PetscObject) user->fe);
345: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
346: if (!isPlex) {
347: PetscSection section;
348: const PetscInt *numDof;
349: PetscInt numComp;
351: PetscFEGetNumComponents(user->fe, &numComp);
352: PetscFEGetNumDof(user->fe, &numDof);
353: DMDACreateSection(dm, &numComp, numDof, NULL, §ion);
354: DMSetDefaultSection(dm, section);
355: PetscSectionDestroy(§ion);
356: }
357: if (!user->simplex && user->constraints) {
358: /* test getting local constraint matrix that matches section */
359: PetscSection aSec;
360: IS aIS;
362: DMPlexGetAnchors(dm,&aSec,&aIS);
363: if (aSec) {
364: PetscDS ds;
365: PetscSection cSec, section;
366: PetscInt cStart, cEnd, c, numComp;
367: Mat cMat, mass;
368: Vec local;
369: const PetscInt *anchors;
371: DMGetDefaultSection(dm,§ion);
372: /* this creates the matrix and preallocates the matrix structure: we
373: * just have to fill in the values */
374: DMGetDefaultConstraints(dm,&cSec,&cMat);
375: PetscSectionGetChart(cSec,&cStart,&cEnd);
376: ISGetIndices(aIS,&anchors);
377: PetscFEGetNumComponents(user->fe, &numComp);
378: for (c = cStart; c < cEnd; c++) {
379: PetscInt cDof;
381: /* is this point constrained? (does it have an anchor?) */
382: PetscSectionGetDof(aSec,c,&cDof);
383: if (cDof) {
384: PetscInt cOff, a, aDof, aOff, j;
385: if (cDof != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Found %d anchor points: should be just one",cDof);
387: /* find the anchor point */
388: PetscSectionGetOffset(aSec,c,&cOff);
389: a = anchors[cOff];
391: /* find the constrained dofs (row in constraint matrix) */
392: PetscSectionGetDof(cSec,c,&cDof);
393: PetscSectionGetOffset(cSec,c,&cOff);
395: /* find the anchor dofs (column in constraint matrix) */
396: PetscSectionGetDof(section,a,&aDof);
397: PetscSectionGetOffset(section,a,&aOff);
399: if (cDof != aDof) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point and anchor have different number of dofs: %d, %d\n",cDof,aDof);
400: if (cDof % numComp) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point dofs not divisible by field components: %d, %d\n",cDof,numComp);
402: /* put in a simple equality constraint */
403: for (j = 0; j < cDof; j++) {
404: MatSetValue(cMat,cOff+j,aOff+j,1.,INSERT_VALUES);
405: }
406: }
407: }
408: MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
409: MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
410: ISRestoreIndices(aIS,&anchors);
412: /* Now that we have constructed the constraint matrix, any FE matrix
413: * that we construct will apply the constraints during construction */
415: DMCreateMatrix(dm,&mass);
416: /* get a dummy local variable to serve as the solution */
417: DMGetLocalVector(dm,&local);
418: DMGetDS(dm,&ds);
419: /* set the jacobian to be the mass matrix */
420: PetscDSSetJacobian(ds, 0, 0, simple_mass, NULL, NULL, NULL);
421: /* build the mass matrix */
422: DMPlexSNESComputeJacobianFEM(dm,local,mass,mass,NULL);
423: MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
424: MatDestroy(&mass);
425: DMRestoreLocalVector(dm,&local);
426: #if 0
427: {
428: /* compare this to periodicity with DMDA: this is broken right now
429: * because DMCreateMatrix() doesn't respect the default section that I
430: * set */
431: DM dmda;
432: PetscSection section;
433: const PetscInt *numDof;
434: PetscInt numComp;
436: /* periodic x */
437: DMDACreate2d(PetscObjectComm((PetscObject)dm), DM_BOUNDARY_PERIODIC, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, -2, -2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dmda);
438: DMDASetVertexCoordinates(dmda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
441: PetscFEGetNumComponents(user->fe, &numComp);
442: PetscFEGetNumDof(user->fe, &numDof);
443: DMDACreateSection(dmda, &numComp, numDof, NULL, §ion);
444: DMSetDefaultSection(dmda, section);
445: PetscSectionDestroy(§ion);
446: DMCreateMatrix(dmda,&mass);
447: /* there isn't a DMDA equivalent of DMPlexSNESComputeJacobianFEM()
448: * right now, but we can at least verify the nonzero structure */
449: MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
450: MatDestroy(&mass);
451: DMDestroy(&dmda);
452: }
453: #endif
454: }
455: }
456: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
457: if (user->tree && isPlex) {
458: Vec local;
459: Mat E;
460: MatNullSpace sp;
461: PetscBool isNullSpace;
462: PetscDS ds;
464: if (user->numComponents != user->dim) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "The number of components %d must be equal to the dimension %d for this test", user->numComponents, user->dim);
465: DMGetDS(dm,&ds);
466: PetscDSSetJacobian(ds,0,0,NULL,NULL,NULL,symmetric_gradient_inner_product);
467: DMCreateMatrix(dm,&E);
468: DMGetLocalVector(dm,&local);
469: DMPlexSNESComputeJacobianFEM(dm,local,E,E,NULL);
470: DMPlexCreateRigidBody(dm,&sp);
471: MatNullSpaceTest(sp,E,&isNullSpace);
472: MatNullSpaceDestroy(&sp);
473: MatDestroy(&E);
474: DMRestoreLocalVector(dm,&local);
475: }
476: return(0);
477: }
481: static PetscErrorCode ComputeError_Plex(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, const PetscReal[], PetscInt, PetscScalar *, void *),
482: PetscErrorCode (**exactFuncDers)(PetscInt, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *),
483: void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
484: {
485: Vec u;
486: PetscReal n[3] = {1.0, 1.0, 1.0};
490: DMGetGlobalVector(dm, &u);
491: /* Project function into FE function space */
492: DMPlexProjectFunction(dm, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u);
493: /* Compare approximation to exact in L_2 */
494: DMPlexComputeL2Diff(dm, exactFuncs, exactCtxs, u, error);
495: DMPlexComputeL2GradientDiff(dm, exactFuncDers, exactCtxs, u, n, errorDer);
496: DMRestoreGlobalVector(dm, &u);
497: return(0);
498: }
502: static PetscErrorCode ComputeError_DA(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, const PetscReal[], PetscInt, PetscScalar *, void *),
503: PetscErrorCode (**exactFuncDers)(PetscInt, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *),
504: void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
505: {
506: Vec u;
507: PetscReal n[3] = {1.0, 1.0, 1.0};
511: DMGetGlobalVector(dm, &u);
512: /* Project function into FE function space */
513: DMDAProjectFunction(dm, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u);
514: /* Compare approximation to exact in L_2 */
515: DMDAComputeL2Diff(dm, exactFuncs, exactCtxs, u, error);
516: DMDAComputeL2GradientDiff(dm, exactFuncDers, exactCtxs, u, n, errorDer);
517: DMRestoreGlobalVector(dm, &u);
518: return(0);
519: }
523: static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, const PetscReal[], PetscInt, PetscScalar *, void *),
524: PetscErrorCode (**exactFuncDers)(PetscInt, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *),
525: void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
526: {
527: PetscBool isPlex, isDA;
531: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
532: PetscObjectTypeCompare((PetscObject) dm, DMDA, &isDA);
533: if (isPlex) {
534: ComputeError_Plex(dm, exactFuncs, exactFuncDers, exactCtxs, error, errorDer, user);
535: } else if (isDA) {
536: ComputeError_DA(dm, exactFuncs, exactFuncDers, exactCtxs, error, errorDer, user);
537: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "No FEM projection routine for this type of DM");
538: return(0);
539: }
543: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
544: {
545: PetscErrorCode (*exactFuncs[1]) (PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
546: PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
547: void *exactCtxs[3] = {user, user, user};
548: MPI_Comm comm;
549: PetscReal error, errorDer, tol = 1.0e-10;
550: PetscErrorCode ierr;
553: user->constants[0] = 1.0;
554: user->constants[1] = 2.0;
555: user->constants[2] = 3.0;
556: PetscObjectGetComm((PetscObject)dm, &comm);
557: /* Setup functions to approximate */
558: switch (order) {
559: case 0:
560: exactFuncs[0] = constant;
561: exactFuncDers[0] = constantDer;
562: break;
563: case 1:
564: exactFuncs[0] = linear;
565: exactFuncDers[0] = linearDer;
566: break;
567: case 2:
568: exactFuncs[0] = quadratic;
569: exactFuncDers[0] = quadraticDer;
570: break;
571: case 3:
572: exactFuncs[0] = cubic;
573: exactFuncDers[0] = cubicDer;
574: break;
575: default:
576: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %d order %d", user->dim, order);
577: }
578: ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
579: /* Report result */
580: if (error > tol) {PetscPrintf(comm, "Function tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol,(double) error);}
581: else {PetscPrintf(comm, "Function tests pass for order %D at tolerance %g\n", order, (double)tol);}
582: if (errorDer > tol) {PetscPrintf(comm, "Function tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);}
583: else {PetscPrintf(comm, "Function tests pass for order %D derivatives at tolerance %g\n", order, (double)tol);}
584: return(0);
585: }
589: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
590: {
591: PetscErrorCode (*exactFuncs[1]) (PetscInt, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
592: PetscErrorCode (*exactFuncDers[1]) (PetscInt, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx);
593: PetscReal n[3] = {1.0, 1.0, 1.0};
594: void *exactCtxs[3] = {user, user, user};
595: DM rdm, idm, fdm;
596: Mat Interp;
597: Vec iu, fu, scaling;
598: MPI_Comm comm;
599: PetscInt dim = user->dim;
600: PetscReal error, errorDer, tol = 1.0e-10;
601: PetscBool isPlex, isDA;
602: PetscErrorCode ierr;
605: user->constants[0] = 1.0;
606: user->constants[1] = 2.0;
607: user->constants[2] = 3.0;
608: PetscObjectGetComm((PetscObject)dm,&comm);
609: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
610: PetscObjectTypeCompare((PetscObject) dm, DMDA, &isDA);
611: DMRefine(dm, comm, &rdm);
612: if (user->tree && isPlex) {
613: DM refTree;
614: DMPlexGetReferenceTree(dm,&refTree);
615: DMPlexSetReferenceTree(rdm,refTree);
616: }
617: if (!user->simplex && !user->constraints) {DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
618: SetupSection(rdm, user);
619: /* Setup functions to approximate */
620: switch (order) {
621: case 0:
622: exactFuncs[0] = constant;
623: exactFuncDers[0] = constantDer;
624: break;
625: case 1:
626: exactFuncs[0] = linear;
627: exactFuncDers[0] = linearDer;
628: break;
629: case 2:
630: exactFuncs[0] = quadratic;
631: exactFuncDers[0] = quadraticDer;
632: break;
633: default:
634: SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %D order %D", dim, order);
635: }
636: idm = checkRestrict ? rdm : dm;
637: fdm = checkRestrict ? dm : rdm;
638: DMGetGlobalVector(idm, &iu);
639: DMGetGlobalVector(fdm, &fu);
640: DMSetApplicationContext(dm, user);
641: DMSetApplicationContext(rdm, user);
642: DMCreateInterpolation(dm, rdm, &Interp, &scaling);
643: /* Project function into initial FE function space */
644: if (isPlex) {
645: DMPlexProjectFunction(idm, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu);
646: } else if (isDA) {
647: DMDAProjectFunction(idm, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu);
648: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "No FEM projection routine for this type of DM");
649: /* Interpolate function into final FE function space */
650: if (checkRestrict) {MatRestrict(Interp, iu, fu);VecPointwiseMult(fu, scaling, fu);}
651: else {MatInterpolate(Interp, iu, fu);}
652: /* Compare approximation to exact in L_2 */
653: if (isPlex) {
654: DMPlexComputeL2Diff(fdm, exactFuncs, exactCtxs, fu, &error);
655: DMPlexComputeL2GradientDiff(fdm, exactFuncDers, exactCtxs, fu, n, &errorDer);
656: } else if (isDA) {
657: DMDAComputeL2Diff(fdm, exactFuncs, exactCtxs, fu, &error);
658: DMDAComputeL2GradientDiff(dm, exactFuncDers, exactCtxs, fu, n, &errorDer);
659: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "No FEM L_2 difference routine for this type of DM");
660: /* Report result */
661: if (error > tol) {PetscPrintf(comm, "Interpolation tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol, (double)error);}
662: else {PetscPrintf(comm, "Interpolation tests pass for order %D at tolerance %g\n", order, (double)tol);}
663: if (errorDer > tol) {PetscPrintf(comm, "Interpolation tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);}
664: else {PetscPrintf(comm, "Interpolation tests pass for order %D derivatives at tolerance %g\n", order, (double)tol);}
665: DMRestoreGlobalVector(idm, &iu);
666: DMRestoreGlobalVector(fdm, &fu);
667: MatDestroy(&Interp);
668: VecDestroy(&scaling);
669: DMDestroy(&rdm);
670: return(0);
671: }
675: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
676: {
677: DM odm = dm, rdm = NULL;
678: PetscErrorCode (*exactFuncs[1]) (PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {trig};
679: PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx) = {trigDer};
680: void *exactCtxs[3] = {user, user, user};
681: PetscInt r;
682: PetscReal errorOld, errorDerOld, error, errorDer;
683: double p;
684: PetscErrorCode ierr;
687: if (!user->convergence) return(0);
688: PetscObjectReference((PetscObject) odm);
689: ComputeError(odm, exactFuncs, exactFuncDers, exactCtxs, &errorOld, &errorDerOld, user);
690: for (r = 0; r < Nr; ++r) {
691: DMRefine(odm, PetscObjectComm((PetscObject) dm), &rdm);
692: if (!user->simplex) {DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
693: SetupSection(rdm, user);
694: ComputeError(rdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
695: p = PetscLog2Real(errorOld/error);
696: PetscPrintf(PetscObjectComm((PetscObject) dm), "Function convergence rate at refinement %D: %.2g\n", r, (double)p);
697: p = PetscLog2Real(errorDerOld/errorDer);
698: PetscPrintf(PetscObjectComm((PetscObject) dm), "Derivative convergence rate at refinement %D: %.2g\n", r, (double)p);
699: DMDestroy(&odm);
700: odm = rdm;
701: errorOld = error;
702: errorDerOld = errorDer;
703: }
704: DMDestroy(&odm);
705: return(0);
706: }
710: int main(int argc, char **argv)
711: {
712: DM dm;
713: AppCtx user; /* user-defined work context */
716: PetscInitialize(&argc, &argv, NULL, help);
717: ProcessOptions(PETSC_COMM_WORLD, &user);
718: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
719: PetscFECreateDefault(dm, user.dim, user.numComponents, user.simplex, NULL, user.qorder, &user.fe);
720: SetupSection(dm, &user);
721: CheckFunctions(dm, user.porder, &user);
722: if (user.dim == 2 && user.simplex == PETSC_TRUE && user.tree == PETSC_FALSE) {
723: CheckInterpolation(dm, PETSC_FALSE, user.porder, &user);
724: CheckInterpolation(dm, PETSC_TRUE, user.porder, &user);
725: }
726: CheckConvergence(dm, 3, &user);
727: PetscFEDestroy(&user.fe);
728: DMDestroy(&dm);
729: PetscFinalize();
730: return 0;
731: }