Actual source code: ex3.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: static char help[] = "Check that a DM can accurately represent and interpolate functions of a given polynomial order\n\n";

  3: #include <petscdmplex.h>
  4: #include <petscdm.h>
  5: #include <petscdmda.h>
  6: #include <petscfe.h>
  7: #include <petscds.h>
  8: #include <petscksp.h>
  9: #include <petscsnes.h>

 11: typedef struct {
 12:   PetscInt  debug;             /* The debugging level */
 13:   /* Domain and mesh definition */
 14:   PetscInt  dim;               /* The topological mesh dimension */
 15:   PetscBool simplex;           /* Flag for simplex or tensor product mesh */
 16:   PetscBool useDA;             /* Flag DMDA tensor product mesh */
 17:   PetscBool interpolate;       /* Generate intermediate mesh elements */
 18:   PetscReal refinementLimit;   /* The largest allowable cell volume */
 19:   /* Element definition */
 20:   PetscInt  qorder;            /* Order of the quadrature */
 21:   PetscInt  numComponents;     /* Number of field components */
 22:   PetscFE   fe;                /* The finite element */
 23:   /* Testing space */
 24:   PetscInt  porder;            /* Order of polynomials to test */
 25:   PetscBool convergence;       /* Test for order of convergence */
 26:   PetscBool convRefine;        /* Test for convergence using refinement, otherwise use coarsening */
 27:   PetscBool constraints;       /* Test local constraints */
 28:   PetscBool tree;              /* Test tree routines */
 29:   PetscBool testFEjacobian;    /* Test finite element Jacobian assembly */
 30:   PetscBool testFVgrad;        /* Test finite difference gradient routine */
 31:   PetscBool testInjector;      /* Test finite element injection routines */
 32:   PetscInt  treeCell;          /* Cell to refine in tree test */
 33:   PetscReal constants[3];      /* Constant values for each dimension */
 34: } AppCtx;

 36: /* u = 1 */
 37: PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 38: {
 39:   AppCtx *user = (AppCtx *) ctx;
 40:   PetscInt d;
 41:   for (d = 0; d < user->dim; ++d) u[d] = user->constants[d];
 42:   return 0;
 43: }
 44: PetscErrorCode constantDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 45: {
 46:   AppCtx *user = (AppCtx *) ctx;
 47:   PetscInt d;
 48:   for (d = 0; d < user->dim; ++d) u[d] = 0.0;
 49:   return 0;
 50: }

 52: /* u = x */
 53: PetscErrorCode linear(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 54: {
 55:   AppCtx *user = (AppCtx *) ctx;
 56:   PetscInt d;
 57:   for (d = 0; d < user->dim; ++d) u[d] = coords[d];
 58:   return 0;
 59: }
 60: PetscErrorCode linearDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 61: {
 62:   AppCtx *user = (AppCtx *) ctx;
 63:   PetscInt d, e;
 64:   for (d = 0; d < user->dim; ++d) {
 65:     u[d] = 0.0;
 66:     for (e = 0; e < user->dim; ++e) u[d] += (d == e ? 1.0 : 0.0) * n[e];
 67:   }
 68:   return 0;
 69: }

 71: /* u = x^2 or u = (x^2, xy) or u = (xy, yz, zx) */
 72: PetscErrorCode quadratic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 73: {
 74:   AppCtx *user = (AppCtx *) ctx;
 75:   if (user->dim > 2)      {u[0] = coords[0]*coords[1]; u[1] = coords[1]*coords[2]; u[2] = coords[2]*coords[0];}
 76:   else if (user->dim > 1) {u[0] = coords[0]*coords[0]; u[1] = coords[0]*coords[1];}
 77:   else if (user->dim > 0) {u[0] = coords[0]*coords[0];}
 78:   return 0;
 79: }
 80: PetscErrorCode quadraticDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 81: {
 82:   AppCtx *user = (AppCtx *) ctx;
 83:   if (user->dim > 2)      {u[0] = coords[1]*n[0] + coords[0]*n[1]; u[1] = coords[2]*n[1] + coords[1]*n[2]; u[2] = coords[2]*n[0] + coords[0]*n[2];}
 84:   else if (user->dim > 1) {u[0] = 2.0*coords[0]*n[0]; u[1] = coords[1]*n[0] + coords[0]*n[1];}
 85:   else if (user->dim > 0) {u[0] = 2.0*coords[0]*n[0];}
 86:   return 0;
 87: }

 89: /* u = x^3 or u = (x^3, x^2y) or u = (x^2y, y^2z, z^2x) */
 90: PetscErrorCode cubic(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
 91: {
 92:   AppCtx *user = (AppCtx *) ctx;
 93:   if (user->dim > 2)      {u[0] = coords[0]*coords[0]*coords[1]; u[1] = coords[1]*coords[1]*coords[2]; u[2] = coords[2]*coords[2]*coords[0];}
 94:   else if (user->dim > 1) {u[0] = coords[0]*coords[0]*coords[0]; u[1] = coords[0]*coords[0]*coords[1];}
 95:   else if (user->dim > 0) {u[0] = coords[0]*coords[0]*coords[0];}
 96:   return 0;
 97: }
 98: PetscErrorCode cubicDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
 99: {
100:   AppCtx *user = (AppCtx *) ctx;
101:   if (user->dim > 2)      {u[0] = 2.0*coords[0]*coords[1]*n[0] + coords[0]*coords[0]*n[1]; u[1] = 2.0*coords[1]*coords[2]*n[1] + coords[1]*coords[1]*n[2]; u[2] = 2.0*coords[2]*coords[0]*n[2] + coords[2]*coords[2]*n[0];}
102:   else if (user->dim > 1) {u[0] = 3.0*coords[0]*coords[0]*n[0]; u[1] = 2.0*coords[0]*coords[1]*n[0] + coords[0]*coords[0]*n[1];}
103:   else if (user->dim > 0) {u[0] = 3.0*coords[0]*coords[0]*n[0];}
104:   return 0;
105: }

107: /* u = tanh(x) */
108: PetscErrorCode trig(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nf, PetscScalar *u, void *ctx)
109: {
110:   AppCtx *user = (AppCtx *) ctx;
111:   PetscInt d;
112:   for (d = 0; d < user->dim; ++d) u[d] = tanh(coords[d] - 0.5);
113:   return 0;
114: }
115: PetscErrorCode trigDer(PetscInt dim, PetscReal time, const PetscReal coords[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx)
116: {
117:   AppCtx *user = (AppCtx *) ctx;
118:   PetscInt d;
119:   for (d = 0; d < user->dim; ++d) u[d] = 1.0/PetscSqr(cosh(coords[d] - 0.5)) * n[d];
120:   return 0;
121: }

125: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
126: {

130:   options->debug           = 0;
131:   options->dim             = 2;
132:   options->simplex         = PETSC_TRUE;
133:   options->useDA           = PETSC_TRUE;
134:   options->interpolate     = PETSC_TRUE;
135:   options->refinementLimit = 0.0;
136:   options->qorder          = 0;
137:   options->numComponents   = 1;
138:   options->porder          = 0;
139:   options->convergence     = PETSC_FALSE;
140:   options->convRefine      = PETSC_TRUE;
141:   options->constraints     = PETSC_FALSE;
142:   options->tree            = PETSC_FALSE;
143:   options->treeCell        = 0;
144:   options->testFEjacobian  = PETSC_FALSE;
145:   options->testFVgrad      = PETSC_FALSE;
146:   options->testInjector    = PETSC_FALSE;

148:   PetscOptionsBegin(comm, "", "Projection Test Options", "DMPlex");
149:   PetscOptionsInt("-debug", "The debugging level", "ex3.c", options->debug, &options->debug, NULL);
150:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex3.c", options->dim, &options->dim, NULL);
151:   PetscOptionsBool("-simplex", "Flag for simplices or hexhedra", "ex3.c", options->simplex, &options->simplex, NULL);
152:   PetscOptionsBool("-use_da", "Flag for DMDA mesh", "ex3.c", options->useDA, &options->useDA, NULL);
153:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex3.c", options->interpolate, &options->interpolate, NULL);
154:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex3.c", options->refinementLimit, &options->refinementLimit, NULL);
155:   PetscOptionsInt("-qorder", "The quadrature order", "ex3.c", options->qorder, &options->qorder, NULL);
156:   PetscOptionsInt("-num_comp", "The number of field components", "ex3.c", options->numComponents, &options->numComponents, NULL);
157:   PetscOptionsInt("-porder", "The order of polynomials to test", "ex3.c", options->porder, &options->porder, NULL);
158:   PetscOptionsBool("-convergence", "Check the convergence rate", "ex3.c", options->convergence, &options->convergence, NULL);
159:   PetscOptionsBool("-conv_refine", "Use refinement for the convergence rate", "ex3.c", options->convRefine, &options->convRefine, NULL);
160:   PetscOptionsBool("-constraints", "Test local constraints (serial only)", "ex3.c", options->constraints, &options->constraints, NULL);
161:   PetscOptionsBool("-tree", "Test tree routines", "ex3.c", options->tree, &options->tree, NULL);
162:   PetscOptionsInt("-tree_cell", "cell to refine in tree test", "ex3.c", options->treeCell, &options->treeCell, NULL);
163:   PetscOptionsBool("-test_fe_jacobian", "Test finite element Jacobian assembly", "ex3.c", options->testFEjacobian, &options->testFEjacobian, NULL);
164:   PetscOptionsBool("-test_fv_grad", "Test finite volume gradient reconstruction", "ex3.c", options->testFVgrad, &options->testFVgrad, NULL);
165:   PetscOptionsBool("-test_injector","Test finite element injection", "ex3.c", options->testInjector, &options->testInjector,NULL);
166:   PetscOptionsEnd();

168:   return(0);
169: }

173: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
174: {
175:   PetscInt       dim             = user->dim;
176:   PetscBool      interpolate     = user->interpolate;
177:   PetscReal      refinementLimit = user->refinementLimit;
178:   PetscBool      isPlex;

182:   if (user->simplex) {
183:     DM refinedMesh     = NULL;

185:     DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
186:     /* Refine mesh using a volume constraint */
187:     DMPlexSetRefinementLimit(*dm, refinementLimit);
188:     DMRefine(*dm, comm, &refinedMesh);
189:     if (refinedMesh) {
190:       DMDestroy(dm);
191:       *dm  = refinedMesh;
192:     }
193:     DMPlexSetRefinementUniform(*dm, PETSC_TRUE);
194:   } else {
195:     if (user->constraints || user->tree || !user->useDA) {
196:       PetscInt cells[3] = {2, 2, 2};

198:       PetscOptionsGetInt(NULL,NULL,"-da_grid_x",&cells[0],NULL);
199:       PetscOptionsGetInt(NULL,NULL,"-da_grid_y",&cells[1],NULL);
200:       PetscOptionsGetInt(NULL,NULL,"-da_grid_z",&cells[2],NULL);
201:       DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);
202:     } else {
203:       switch (user->dim) {
204:       case 2:
205:         DMDACreate2d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, -2, -2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, dm);
206:         DMDASetVertexCoordinates(*dm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
207:         break;
208:       default:
209:         SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot create structured mesh of dimension %d", dim);
210:       }
211:       PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");
212:     }
213:   }
214:   PetscObjectTypeCompare((PetscObject)*dm,DMPLEX,&isPlex);
215:   if (isPlex) {
216:     DM distributedMesh = NULL;
217:     if (user->tree) {
218:       DM refTree;
219:       DM ncdm = NULL;

221:       DMPlexCreateDefaultReferenceTree(comm,user->dim,user->simplex,&refTree);
222:       DMPlexSetReferenceTree(*dm,refTree);
223:       DMDestroy(&refTree);
224:       DMPlexTreeRefineCell(*dm,user->treeCell,&ncdm);
225:       if (ncdm) {
226:         DMDestroy(dm);
227:         *dm = ncdm;
228:         DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
229:       }
230:     }
231:     else {
232:       DMPlexSetRefinementUniform(*dm, PETSC_TRUE);
233:     }
234:     /* Distribute mesh over processes */
235:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
236:     if (distributedMesh) {
237:       DMDestroy(dm);
238:       *dm  = distributedMesh;
239:     }
240:     if (user->simplex) {
241:       PetscObjectSetName((PetscObject) *dm, "Simplicial Mesh");
242:     }
243:     else {
244:       PetscObjectSetName((PetscObject) *dm, "Hexahedral Mesh");
245:     }
246:   }
247:   DMSetFromOptions(*dm);
248:   DMViewFromOptions(*dm,NULL,"-dm_view");
249:   return(0);
250: }

254: static void simple_mass(PetscInt dim, PetscInt Nf, PetscInt NfAux,
255:                         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
256:                         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
257:                         PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])
258: {
259:   PetscInt d, e;
260:   for (d = 0, e = 0; d < dim; d++, e+=dim+1) {
261:     g0[e] = 1.;
262:   }
263: }

265: /* < \nabla v, 1/2(\nabla u + {\nabla u}^T) > */
268: static void symmetric_gradient_inner_product(PetscInt dim, PetscInt Nf, PetscInt NfAux,
269:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
270:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
271:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar C[])
272: {
273:   PetscInt compI, compJ, d, e;

275:   for (compI = 0; compI < dim; ++compI) {
276:     for (compJ = 0; compJ < dim; ++compJ) {
277:       for (d = 0; d < dim; ++d) {
278:         for (e = 0; e < dim; e++) {
279:           if (d == e && d == compI && d == compJ) {
280:             C[((compI*dim+compJ)*dim+d)*dim+e] = 1.0;
281:           }
282:           else if ((d == compJ && e == compI) || (d == e && compI == compJ)) {
283:             C[((compI*dim+compJ)*dim+d)*dim+e] = 0.5;
284:           }
285:           else {
286:             C[((compI*dim+compJ)*dim+d)*dim+e] = 0.0;
287:           }
288:         }
289:       }
290:     }
291:   }
292: }

296: static PetscErrorCode SetupSection(DM dm, AppCtx *user)
297: {
298:   PetscBool      isPlex;

302:   if (!user->simplex && user->constraints) {
303:     /* test local constraints */
304:     DM            coordDM;
305:     PetscInt      fStart, fEnd, f, vStart, vEnd, v;
306:     PetscInt      edgesx = 2, vertsx;
307:     PetscInt      edgesy = 2, vertsy;
308:     PetscMPIInt   size;
309:     PetscInt      numConst;
310:     PetscSection  aSec;
311:     PetscInt     *anchors;
312:     PetscInt      offset;
313:     IS            aIS;
314:     MPI_Comm comm = PetscObjectComm((PetscObject)dm);

316:     MPI_Comm_size(comm,&size);
317:     if (size > 1) SETERRQ(comm,PETSC_ERR_SUP,"Local constraint test can only be performed in serial");

319:     /* we are going to test constraints by using them to enforce periodicity
320:      * in one direction, and comparing to the existing method of enforcing
321:      * periodicity */

323:     /* first create the coordinate section so that it does not clone the
324:      * constraints */
325:     DMGetCoordinateDM(dm,&coordDM);

327:     /* create the constrained-to-anchor section */
328:     DMPlexGetDepthStratum(dm,0,&vStart,&vEnd);
329:     DMPlexGetDepthStratum(dm,1,&fStart,&fEnd);
330:     PetscSectionCreate(PETSC_COMM_SELF,&aSec);
331:     PetscSectionSetChart(aSec,PetscMin(fStart,vStart),PetscMax(fEnd,vEnd));

333:     /* define the constraints */
334:     PetscOptionsGetInt(NULL,NULL,"-da_grid_x",&edgesx,NULL);
335:     PetscOptionsGetInt(NULL,NULL,"-da_grid_y",&edgesy,NULL);
336:     vertsx = edgesx + 1;
337:     vertsy = edgesy + 1;
338:     numConst = vertsy + edgesy;
339:     PetscMalloc1(numConst,&anchors);
340:     offset = 0;
341:     for (v = vStart + edgesx; v < vEnd; v+= vertsx) {
342:       PetscSectionSetDof(aSec,v,1);
343:       anchors[offset++] = v - edgesx;
344:     }
345:     for (f = fStart + edgesx * vertsy + edgesx * edgesy; f < fEnd; f++) {
346:       PetscSectionSetDof(aSec,f,1);
347:       anchors[offset++] = f - edgesx * edgesy;
348:     }
349:     PetscSectionSetUp(aSec);
350:     ISCreateGeneral(PETSC_COMM_SELF,numConst,anchors,PETSC_OWN_POINTER,&aIS);

352:     DMPlexSetAnchors(dm,aSec,aIS);
353:     PetscSectionDestroy(&aSec);
354:     ISDestroy(&aIS);
355:   }
356:   DMSetNumFields(dm, 1);
357:   DMSetField(dm, 0, (PetscObject) user->fe);
358:   PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
359:   if (!isPlex) {
360:       PetscSection    section;
361:       const PetscInt *numDof;
362:       PetscInt        numComp;

364:       PetscFEGetNumComponents(user->fe, &numComp);
365:       PetscFEGetNumDof(user->fe, &numDof);
366:       DMDACreateSection(dm, &numComp, numDof, NULL, &section);
367:       DMSetDefaultSection(dm, section);
368:       PetscSectionDestroy(&section);
369:   }
370:   if (!user->simplex && user->constraints) {
371:     /* test getting local constraint matrix that matches section */
372:     PetscSection aSec;
373:     IS           aIS;

375:     DMPlexGetAnchors(dm,&aSec,&aIS);
376:     if (aSec) {
377:       PetscDS         ds;
378:       PetscSection    cSec, section;
379:       PetscInt        cStart, cEnd, c, numComp;
380:       Mat             cMat, mass;
381:       Vec             local;
382:       const PetscInt *anchors;

384:       DMGetDefaultSection(dm,&section);
385:       /* this creates the matrix and preallocates the matrix structure: we
386:        * just have to fill in the values */
387:       DMGetDefaultConstraints(dm,&cSec,&cMat);
388:       PetscSectionGetChart(cSec,&cStart,&cEnd);
389:       ISGetIndices(aIS,&anchors);
390:       PetscFEGetNumComponents(user->fe, &numComp);
391:       for (c = cStart; c < cEnd; c++) {
392:         PetscInt cDof;

394:         /* is this point constrained? (does it have an anchor?) */
395:         PetscSectionGetDof(aSec,c,&cDof);
396:         if (cDof) {
397:           PetscInt cOff, a, aDof, aOff, j;
398:           if (cDof != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Found %d anchor points: should be just one",cDof);

400:           /* find the anchor point */
401:           PetscSectionGetOffset(aSec,c,&cOff);
402:           a    = anchors[cOff];

404:           /* find the constrained dofs (row in constraint matrix) */
405:           PetscSectionGetDof(cSec,c,&cDof);
406:           PetscSectionGetOffset(cSec,c,&cOff);

408:           /* find the anchor dofs (column in constraint matrix) */
409:           PetscSectionGetDof(section,a,&aDof);
410:           PetscSectionGetOffset(section,a,&aOff);

412:           if (cDof != aDof) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point and anchor have different number of dofs: %d, %d\n",cDof,aDof);
413:           if (cDof % numComp) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Point dofs not divisible by field components: %d, %d\n",cDof,numComp);

415:           /* put in a simple equality constraint */
416:           for (j = 0; j < cDof; j++) {
417:             MatSetValue(cMat,cOff+j,aOff+j,1.,INSERT_VALUES);
418:           }
419:         }
420:       }
421:       MatAssemblyBegin(cMat,MAT_FINAL_ASSEMBLY);
422:       MatAssemblyEnd(cMat,MAT_FINAL_ASSEMBLY);
423:       ISRestoreIndices(aIS,&anchors);

425:       /* Now that we have constructed the constraint matrix, any FE matrix
426:        * that we construct will apply the constraints during construction */

428:       DMCreateMatrix(dm,&mass);
429:       /* get a dummy local variable to serve as the solution */
430:       DMGetLocalVector(dm,&local);
431:       DMGetDS(dm,&ds);
432:       /* set the jacobian to be the mass matrix */
433:       PetscDSSetJacobian(ds, 0, 0, simple_mass, NULL,  NULL, NULL);
434:       /* build the mass matrix */
435:       DMPlexSNESComputeJacobianFEM(dm,local,mass,mass,NULL);
436:       MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
437:       MatDestroy(&mass);
438:       DMRestoreLocalVector(dm,&local);
439: #if 0
440:       {
441:         /* compare this to periodicity with DMDA: this is broken right now
442:          * because DMCreateMatrix() doesn't respect the default section that I
443:          * set */
444:         DM              dmda;
445:         PetscSection    section;
446:         const PetscInt *numDof;
447:         PetscInt        numComp;

449:                                                               /* periodic x */
450:         DMDACreate2d(PetscObjectComm((PetscObject)dm), DM_BOUNDARY_PERIODIC, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, -2, -2, PETSC_DETERMINE, PETSC_DETERMINE, 1, 1, NULL, NULL, &dmda);
451:         DMDASetVertexCoordinates(dmda, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);


454:         PetscFEGetNumComponents(user->fe, &numComp);
455:         PetscFEGetNumDof(user->fe, &numDof);
456:         DMDACreateSection(dmda, &numComp, numDof, NULL, &section);
457:         DMSetDefaultSection(dmda, section);
458:         PetscSectionDestroy(&section);
459:         DMCreateMatrix(dmda,&mass);
460:         /* there isn't a DMDA equivalent of DMPlexSNESComputeJacobianFEM()
461:          * right now, but we can at least verify the nonzero structure */
462:         MatView(mass,PETSC_VIEWER_STDOUT_WORLD);
463:         MatDestroy(&mass);
464:         DMDestroy(&dmda);
465:       }
466: #endif
467:     }
468:   }
469:   return(0);
470: }

474: static PetscErrorCode TestFEJacobian(DM dm, AppCtx *user)
475: {
476:   PetscBool      isPlex;

480:   PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
481:   if (isPlex) {
482:     Vec          local;
483:     Mat          E;
484:     MatNullSpace sp;
485:     PetscBool    isNullSpace;
486:     PetscDS ds;

488:     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);
489:     DMGetDS(dm,&ds);
490:     PetscDSSetJacobian(ds,0,0,NULL,NULL,NULL,symmetric_gradient_inner_product);
491:     DMCreateMatrix(dm,&E);
492:     DMGetLocalVector(dm,&local);
493:     DMPlexSNESComputeJacobianFEM(dm,local,E,E,NULL);
494:     DMPlexCreateRigidBody(dm,&sp);
495:     MatNullSpaceTest(sp,E,&isNullSpace);
496:     if (isNullSpace) {
497:       PetscPrintf(PetscObjectComm((PetscObject)dm),"Symmetric gradient null space: PASS\n");
498:     }
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: }

511: static PetscErrorCode TestInjector(DM dm, AppCtx *user)
512: {
513:   DM             refTree;
514:   PetscMPIInt    rank;

518:   DMPlexGetReferenceTree(dm,&refTree);
519:   if (refTree) {
520:     Mat inj;

522:     DMPlexComputeInjectorReferenceTree(refTree,&inj);
523:     PetscObjectSetName((PetscObject)inj,"Reference Tree Injector");
524:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
525:     if (!rank) {
526:       MatView(inj,PETSC_VIEWER_STDOUT_SELF);
527:     }
528:     MatDestroy(&inj);
529:   }
530:   return(0);
531: }

535: static PetscErrorCode TestFVGrad(DM dm, AppCtx *user)
536: {
537:   MPI_Comm comm;
538:   DM dmRedist, dmfv, dmgrad, dmCell, refTree;
539:   PetscFV fv;
540:   PetscInt nvecs, v, cStart, cEnd, cEndInterior;
541:   PetscMPIInt size;
542:   Vec cellgeom, grad, locGrad;
543:   const PetscScalar *cgeom;
544:   PetscReal allVecMaxDiff = 0., fvTol = 100. * PETSC_MACHINE_EPSILON;

548:   comm = PetscObjectComm((PetscObject)dm);
549:   /* duplicate DM, give dup. a FV discretization */
550:   DMPlexSetAdjacencyUseCone(dm,PETSC_TRUE);
551:   DMPlexSetAdjacencyUseClosure(dm,PETSC_FALSE);
552:   MPI_Comm_size(comm,&size);
553:   dmRedist = NULL;
554:   if (size > 1) {
555:     DMPlexDistributeOverlap(dm,1,NULL,&dmRedist);
556:   }
557:   if (!dmRedist) {
558:     dmRedist = dm;
559:     PetscObjectReference((PetscObject)dmRedist);
560:   }
561:   PetscFVCreate(comm,&fv);
562:   PetscFVSetType(fv,PETSCFVLEASTSQUARES);
563:   PetscFVSetNumComponents(fv,user->numComponents);
564:   PetscFVSetSpatialDimension(fv,user->dim);
565:   PetscFVSetFromOptions(fv);
566:   PetscFVSetUp(fv);
567:   DMPlexConstructGhostCells(dmRedist,NULL,NULL,&dmfv);
568:   DMDestroy(&dmRedist);
569:   DMSetNumFields(dmfv,1);
570:   DMSetField(dmfv, 0, (PetscObject) fv);
571:   DMPlexGetReferenceTree(dm,&refTree);
572:   if (refTree) {
573:     PetscDS ds;
574:     DMGetDS(dmfv,&ds);
575:     DMSetDS(refTree,ds);
576:   }
577:   DMPlexSNESGetGradientDM(dmfv, fv, &dmgrad);
578:   DMPlexGetHeightStratum(dmfv,0,&cStart,&cEnd);
579:   nvecs = user->dim * (user->dim+1) / 2;
580:   DMPlexSNESGetGeometryFVM(dmfv,NULL,&cellgeom,NULL);
581:   VecGetDM(cellgeom,&dmCell);
582:   VecGetArrayRead(cellgeom,&cgeom);
583:   DMGetGlobalVector(dmgrad,&grad);
584:   DMGetLocalVector(dmgrad,&locGrad);
585:   DMPlexGetHybridBounds(dmgrad,&cEndInterior,NULL,NULL,NULL);
586:   cEndInterior = (cEndInterior < 0) ? cEnd: cEndInterior;
587:   for (v = 0; v < nvecs; v++) {
588:     Vec locX;
589:     PetscInt c;
590:     PetscScalar trueGrad[3][3] = {{0.}};
591:     const PetscScalar *gradArray;
592:     PetscReal maxDiff, maxDiffGlob;

594:     DMGetLocalVector(dmfv,&locX);
595:     /* get the local projection of the rigid body mode */
596:     for (c = cStart; c < cEnd; c++) {
597:       PetscFVCellGeom *cg;
598:       PetscScalar cx[3] = {0.,0.,0.};

600:       DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
601:       if (v < user->dim) {
602:         cx[v] = 1.;
603:       }
604:       else {
605:         PetscInt w = v - user->dim;

607:         cx[(w + 1) % user->dim] =  cg->centroid[(w + 2) % user->dim];
608:         cx[(w + 2) % user->dim] = -cg->centroid[(w + 1) % user->dim];
609:       }
610:       DMPlexVecSetClosure(dmfv,NULL,locX,c,cx,INSERT_ALL_VALUES);
611:     }
612:     /* TODO: this isn't in any header */
613:     DMPlexReconstructGradientsFVM(dmfv,locX,grad);
614:     DMGlobalToLocalBegin(dmgrad,grad,INSERT_VALUES,locGrad);
615:     DMGlobalToLocalEnd(dmgrad,grad,INSERT_VALUES,locGrad);
616:     VecGetArrayRead(locGrad,&gradArray);
617:     /* compare computed gradient to exact gradient */
618:     if (v >= user->dim) {
619:       PetscInt w = v - user->dim;

621:       trueGrad[(w + 1) % user->dim][(w + 2) % user->dim] =  1.;
622:       trueGrad[(w + 2) % user->dim][(w + 1) % user->dim] = -1.;
623:     }
624:     maxDiff = 0.;
625:     for (c = cStart; c < cEndInterior; c++) {
626:       PetscScalar *compGrad;
627:       PetscInt i, j, k;
628:       PetscReal FrobDiff = 0.;

630:       DMPlexPointLocalRead(dmgrad, c, gradArray, &compGrad);

632:       for (i = 0, k = 0; i < user->dim; i++) {
633:         for (j = 0; j < user->dim; j++, k++) {
634:           PetscScalar diff = compGrad[k] - trueGrad[i][j];
635:           FrobDiff += PetscRealPart(diff * PetscConj(diff));
636:         }
637:       }
638:       FrobDiff = PetscSqrtReal(FrobDiff);
639:       maxDiff  = PetscMax(maxDiff,FrobDiff);
640:     }
641:     MPI_Allreduce(&maxDiff,&maxDiffGlob,1,MPIU_REAL,MPIU_MAX,comm);
642:     allVecMaxDiff = PetscMax(allVecMaxDiff,maxDiffGlob);
643:     VecRestoreArrayRead(locGrad,&gradArray);
644:     DMRestoreLocalVector(dmfv,&locX);
645:   }
646:   if (allVecMaxDiff < fvTol) {
647:     PetscPrintf(PetscObjectComm((PetscObject)dm),"Finite volume gradient reconstruction: PASS\n");
648:   }
649:   else {
650:     PetscPrintf(PetscObjectComm((PetscObject)dm),"Finite volume gradient reconstruction: FAIL at tolerance %g with max difference %g\n",fvTol,allVecMaxDiff);
651:   }
652:   DMRestoreLocalVector(dmgrad,&locGrad);
653:   DMRestoreGlobalVector(dmgrad,&grad);
654:   VecRestoreArrayRead(cellgeom,&cgeom);
655:   DMDestroy(&dmfv);
656:   PetscFVDestroy(&fv);
657:   return(0);
658: }

662: static PetscErrorCode ComputeError(DM dm, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *),
663:                                    PetscErrorCode (**exactFuncDers)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *),
664:                                    void **exactCtxs, PetscReal *error, PetscReal *errorDer, AppCtx *user)
665: {
666:   Vec            u;
667:   PetscReal      n[3] = {1.0, 1.0, 1.0};

671:   DMGetGlobalVector(dm, &u);
672:   /* Project function into FE function space */
673:   DMProjectFunction(dm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, u);
674:   /* Compare approximation to exact in L_2 */
675:   DMComputeL2Diff(dm, 0.0, exactFuncs, exactCtxs, u, error);
676:   DMComputeL2GradientDiff(dm, 0.0, exactFuncDers, exactCtxs, u, n, errorDer);
677:   DMRestoreGlobalVector(dm, &u);
678:   return(0);
679: }

683: static PetscErrorCode CheckFunctions(DM dm, PetscInt order, AppCtx *user)
684: {
685:   PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
686:   PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx);
687:   void            *exactCtxs[3];
688:   MPI_Comm         comm;
689:   PetscReal        error, errorDer, tol = PETSC_SMALL;
690:   PetscErrorCode   ierr;

693:   exactCtxs[0]       = user;
694:   exactCtxs[1]       = user;
695:   exactCtxs[2]       = user;
696:   user->constants[0] = 1.0;
697:   user->constants[1] = 2.0;
698:   user->constants[2] = 3.0;
699:   PetscObjectGetComm((PetscObject)dm, &comm);
700:   /* Setup functions to approximate */
701:   switch (order) {
702:   case 0:
703:     exactFuncs[0]    = constant;
704:     exactFuncDers[0] = constantDer;
705:     break;
706:   case 1:
707:     exactFuncs[0]    = linear;
708:     exactFuncDers[0] = linearDer;
709:     break;
710:   case 2:
711:     exactFuncs[0]    = quadratic;
712:     exactFuncDers[0] = quadraticDer;
713:     break;
714:   case 3:
715:     exactFuncs[0]    = cubic;
716:     exactFuncDers[0] = cubicDer;
717:     break;
718:   default:
719:     SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %d order %d", user->dim, order);
720:   }
721:   ComputeError(dm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
722:   /* Report result */
723:   if (error > tol)    {PetscPrintf(comm, "Function tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol,(double) error);}
724:   else                {PetscPrintf(comm, "Function tests pass for order %D at tolerance %g\n", order, (double)tol);}
725:   if (errorDer > tol) {PetscPrintf(comm, "Function tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);}
726:   else                {PetscPrintf(comm, "Function tests pass for order %D derivatives at tolerance %g\n", order, (double)tol);}
727:   return(0);
728: }

732: static PetscErrorCode CheckInterpolation(DM dm, PetscBool checkRestrict, PetscInt order, AppCtx *user)
733: {
734:   PetscErrorCode (*exactFuncs[1]) (PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
735:   PetscErrorCode (*exactFuncDers[1]) (PetscInt, PetscReal, const PetscReal x[], const PetscReal n[], PetscInt, PetscScalar *u, void *ctx);
736:   PetscReal       n[3]         = {1.0, 1.0, 1.0};
737:   void           *exactCtxs[3];
738:   DM              rdm, idm, fdm;
739:   Mat             Interp;
740:   Vec             iu, fu, scaling;
741:   MPI_Comm        comm;
742:   PetscInt        dim  = user->dim;
743:   PetscReal       error, errorDer, tol = 1.0e-10;
744:   PetscBool       isPlex, isDA;
745:   PetscErrorCode  ierr;

748:   exactCtxs[0]       = user;
749:   exactCtxs[1]       = user;
750:   exactCtxs[2]       = user;
751:   user->constants[0] = 1.0;
752:   user->constants[1] = 2.0;
753:   user->constants[2] = 3.0;
754:   PetscObjectGetComm((PetscObject)dm,&comm);
755:   PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
756:   PetscObjectTypeCompare((PetscObject) dm, DMDA,   &isDA);
757:   DMRefine(dm, comm, &rdm);
758:   DMSetCoarseDM(rdm, dm);
759:   DMPlexSetRegularRefinement(rdm, user->convRefine);
760:   if (user->tree && isPlex) {
761:     DM refTree;
762:     DMPlexGetReferenceTree(dm,&refTree);
763:     DMPlexSetReferenceTree(rdm,refTree);
764:   }
765:   if (!user->simplex && !user->constraints) {DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
766:   SetupSection(rdm, user);
767:   /* Setup functions to approximate */
768:   switch (order) {
769:   case 0:
770:     exactFuncs[0]    = constant;
771:     exactFuncDers[0] = constantDer;
772:     break;
773:   case 1:
774:     exactFuncs[0]    = linear;
775:     exactFuncDers[0] = linearDer;
776:     break;
777:   case 2:
778:     exactFuncs[0]    = quadratic;
779:     exactFuncDers[0] = quadraticDer;
780:     break;
781:   default:
782:     SETERRQ2(comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not determine functions to test for dimension %D order %D", dim, order);
783:   }
784:   idm  = checkRestrict ? rdm :  dm;
785:   fdm  = checkRestrict ?  dm : rdm;
786:   DMGetGlobalVector(idm, &iu);
787:   DMGetGlobalVector(fdm, &fu);
788:   DMSetApplicationContext(dm, user);
789:   DMSetApplicationContext(rdm, user);
790:   DMCreateInterpolation(dm, rdm, &Interp, &scaling);
791:   /* Project function into initial FE function space */
792:   DMProjectFunction(idm, 0.0, exactFuncs, exactCtxs, INSERT_ALL_VALUES, iu);
793:   /* Interpolate function into final FE function space */
794:   if (checkRestrict) {MatRestrict(Interp, iu, fu);VecPointwiseMult(fu, scaling, fu);}
795:   else               {MatInterpolate(Interp, iu, fu);}
796:   /* Compare approximation to exact in L_2 */
797:   DMComputeL2Diff(fdm, 0.0, exactFuncs, exactCtxs, fu, &error);
798:   DMComputeL2GradientDiff(fdm, 0.0, exactFuncDers, exactCtxs, fu, n, &errorDer);
799:   /* Report result */
800:   if (error > tol)    {PetscPrintf(comm, "Interpolation tests FAIL for order %D at tolerance %g error %g\n", order, (double)tol, (double)error);}
801:   else                {PetscPrintf(comm, "Interpolation tests pass for order %D at tolerance %g\n", order, (double)tol);}
802:   if (errorDer > tol) {PetscPrintf(comm, "Interpolation tests FAIL for order %D derivatives at tolerance %g error %g\n", order, (double)tol, (double)errorDer);}
803:   else                {PetscPrintf(comm, "Interpolation tests pass for order %D derivatives at tolerance %g\n", order, (double)tol);}
804:   DMRestoreGlobalVector(idm, &iu);
805:   DMRestoreGlobalVector(fdm, &fu);
806:   MatDestroy(&Interp);
807:   VecDestroy(&scaling);
808:   DMDestroy(&rdm);
809:   return(0);
810: }

814: static PetscErrorCode CheckConvergence(DM dm, PetscInt Nr, AppCtx *user)
815: {
816:   DM               odm = dm, rdm = NULL, cdm = NULL;
817:   PetscErrorCode (*exactFuncs[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {trig};
818:   PetscErrorCode (*exactFuncDers[1]) (PetscInt dim, PetscReal time, const PetscReal x[], const PetscReal n[], PetscInt Nf, PetscScalar *u, void *ctx) = {trigDer};
819:   void            *exactCtxs[3];
820:   PetscInt         r, c, cStart, cEnd;
821:   PetscReal        errorOld, errorDerOld, error, errorDer, rel, len, lenOld;
822:   double           p;
823:   PetscErrorCode   ierr;

826:   if (!user->convergence) return(0);
827:   exactCtxs[0] = user;
828:   exactCtxs[1] = user;
829:   exactCtxs[2] = user;
830:   PetscObjectReference((PetscObject) odm);
831:   if (!user->convRefine) {
832:     for (r = 0; r < Nr; ++r) {
833:       DMRefine(odm, PetscObjectComm((PetscObject) dm), &rdm);
834:       DMDestroy(&odm);
835:       odm  = rdm;
836:     }
837:     SetupSection(odm, user);
838:   }
839:   ComputeError(odm, exactFuncs, exactFuncDers, exactCtxs, &errorOld, &errorDerOld, user);
840:   if (user->convRefine) {
841:     for (r = 0; r < Nr; ++r) {
842:       DMRefine(odm, PetscObjectComm((PetscObject) dm), &rdm);
843:       if (!user->simplex) {DMDASetVertexCoordinates(rdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
844:       SetupSection(rdm, user);
845:       ComputeError(rdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
846:       p    = PetscLog2Real(errorOld/error);
847:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Function   convergence rate at refinement %D: %.2g\n", r, (double)p);
848:       p    = PetscLog2Real(errorDerOld/errorDer);
849:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Derivative convergence rate at refinement %D: %.2g\n", r, (double)p);
850:       DMDestroy(&odm);
851:       odm         = rdm;
852:       errorOld    = error;
853:       errorDerOld = errorDer;
854:     }
855:   } else {
856:     /* ComputeLongestEdge(dm, &lenOld); */
857:     DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd);
858:     lenOld = cEnd - cStart;
859:     for (c = 0; c < Nr; ++c) {
860:       DMCoarsen(odm, PetscObjectComm((PetscObject) dm), &cdm);
861:       if (!user->simplex) {DMDASetVertexCoordinates(cdm, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);}
862:       SetupSection(cdm, user);
863:       ComputeError(cdm, exactFuncs, exactFuncDers, exactCtxs, &error, &errorDer, user);
864:       /* ComputeLongestEdge(cdm, &len); */
865:       DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
866:       len  = cEnd - cStart;
867:       rel  = error/errorOld;
868:       p    = PetscLogReal(rel) / PetscLogReal(lenOld / len);
869:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Function   convergence rate at coarsening %D: %.2g\n", c, (double)p);
870:       rel  = errorDer/errorDerOld;
871:       p    = PetscLogReal(rel) / PetscLogReal(lenOld / len);
872:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Derivative convergence rate at coarsening %D: %.2g\n", c, (double)p);
873:       DMDestroy(&odm);
874:       odm         = cdm;
875:       errorOld    = error;
876:       errorDerOld = errorDer;
877:       lenOld      = len;
878:     }
879:   }
880:   DMDestroy(&odm);
881:   return(0);
882: }

886: int main(int argc, char **argv)
887: {
888:   DM             dm;
889:   AppCtx         user;                 /* user-defined work context */

892:   PetscInitialize(&argc, &argv, NULL, help);
893:   ProcessOptions(PETSC_COMM_WORLD, &user);
894:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
895:   PetscFECreateDefault(dm, user.dim, user.numComponents, user.simplex, NULL, user.qorder, &user.fe);
896:   SetupSection(dm, &user);
897:   if (user.testFEjacobian) {
898:     TestFEJacobian(dm, &user);
899:   }
900:   if (user.testFVgrad) {
901:     TestFVGrad(dm, &user);
902:   }
903:   if (user.testInjector) {
904:     TestInjector(dm, &user);
905:   }
906:   CheckFunctions(dm, user.porder, &user);
907:   if (user.dim == 2 && user.simplex == PETSC_TRUE && user.tree == PETSC_FALSE) {
908:     CheckInterpolation(dm, PETSC_FALSE, user.porder, &user);
909:     CheckInterpolation(dm, PETSC_TRUE,  user.porder, &user);
910:   }
911:   CheckConvergence(dm, 3, &user);
912:   PetscFEDestroy(&user.fe);
913:   DMDestroy(&dm);
914:   PetscFinalize();
915:   return 0;
916: }