Actual source code: plexfem.c

petsc-master 2017-07-25
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petscsf.h>

  4:  #include <petsc/private/petscfeimpl.h>
  5:  #include <petsc/private/petscfvimpl.h>

  7: /*@
  8:   DMPlexGetScale - Get the scale for the specified fundamental unit

 10:   Not collective

 12:   Input Arguments:
 13: + dm   - the DM
 14: - unit - The SI unit

 16:   Output Argument:
 17: . scale - The value used to scale all quantities with this unit

 19:   Level: advanced

 21: .seealso: DMPlexSetScale(), PetscUnit
 22: @*/
 23: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
 24: {
 25:   DM_Plex *mesh = (DM_Plex*) dm->data;

 30:   *scale = mesh->scale[unit];
 31:   return(0);
 32: }

 34: /*@
 35:   DMPlexSetScale - Set the scale for the specified fundamental unit

 37:   Not collective

 39:   Input Arguments:
 40: + dm   - the DM
 41: . unit - The SI unit
 42: - scale - The value used to scale all quantities with this unit

 44:   Level: advanced

 46: .seealso: DMPlexGetScale(), PetscUnit
 47: @*/
 48: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
 49: {
 50:   DM_Plex *mesh = (DM_Plex*) dm->data;

 54:   mesh->scale[unit] = scale;
 55:   return(0);
 56: }

 58: PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
 59: {
 60:   switch (i) {
 61:   case 0:
 62:     switch (j) {
 63:     case 0: return 0;
 64:     case 1:
 65:       switch (k) {
 66:       case 0: return 0;
 67:       case 1: return 0;
 68:       case 2: return 1;
 69:       }
 70:     case 2:
 71:       switch (k) {
 72:       case 0: return 0;
 73:       case 1: return -1;
 74:       case 2: return 0;
 75:       }
 76:     }
 77:   case 1:
 78:     switch (j) {
 79:     case 0:
 80:       switch (k) {
 81:       case 0: return 0;
 82:       case 1: return 0;
 83:       case 2: return -1;
 84:       }
 85:     case 1: return 0;
 86:     case 2:
 87:       switch (k) {
 88:       case 0: return 1;
 89:       case 1: return 0;
 90:       case 2: return 0;
 91:       }
 92:     }
 93:   case 2:
 94:     switch (j) {
 95:     case 0:
 96:       switch (k) {
 97:       case 0: return 0;
 98:       case 1: return 1;
 99:       case 2: return 0;
100:       }
101:     case 1:
102:       switch (k) {
103:       case 0: return -1;
104:       case 1: return 0;
105:       case 2: return 0;
106:       }
107:     case 2: return 0;
108:     }
109:   }
110:   return 0;
111: }

113: static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx)
114: {
115:   PetscInt *ctxInt  = (PetscInt *) ctx;
116:   PetscInt  dim2    = ctxInt[0];
117:   PetscInt  d       = ctxInt[1];
118:   PetscInt  i, j, k = dim > 2 ? d - dim : d;

121:   if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2);
122:   for (i = 0; i < dim; i++) mode[i] = 0.;
123:   if (d < dim) {
124:     mode[d] = 1.;
125:   } else {
126:     for (i = 0; i < dim; i++) {
127:       for (j = 0; j < dim; j++) {
128:         mode[j] += epsilon(i, j, k)*X[i];
129:       }
130:     }
131:   }
132:   return(0);
133: }

135: /*@C
136:   DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates

138:   Collective on DM

140:   Input Arguments:
141: . dm - the DM

143:   Output Argument:
144: . sp - the null space

146:   Note: This is necessary to take account of Dirichlet conditions on the displacements

148:   Level: advanced

150: .seealso: MatNullSpaceCreate()
151: @*/
152: PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
153: {
154:   MPI_Comm       comm;
155:   Vec            mode[6];
156:   PetscSection   section, globalSection;
157:   PetscInt       dim, dimEmbed, n, m, d, i, j;

161:   PetscObjectGetComm((PetscObject)dm,&comm);
162:   DMGetDimension(dm, &dim);
163:   DMGetCoordinateDim(dm, &dimEmbed);
164:   if (dim == 1) {
165:     MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
166:     return(0);
167:   }
168:   DMGetDefaultSection(dm, &section);
169:   DMGetDefaultGlobalSection(dm, &globalSection);
170:   PetscSectionGetConstrainedStorageSize(globalSection, &n);
171:   m    = (dim*(dim+1))/2;
172:   VecCreate(comm, &mode[0]);
173:   VecSetSizes(mode[0], n, PETSC_DETERMINE);
174:   VecSetUp(mode[0]);
175:   for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
176:   for (d = 0; d < m; d++) {
177:     PetscInt         ctx[2];
178:     PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
179:     void            *voidctx = (void *) (&ctx[0]);

181:     ctx[0] = dimEmbed;
182:     ctx[1] = d;
183:     DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);
184:   }
185:   for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
186:   /* Orthonormalize system */
187:   for (i = dim; i < m; ++i) {
188:     PetscScalar dots[6];

190:     VecMDot(mode[i], i, mode, dots);
191:     for (j = 0; j < i; ++j) dots[j] *= -1.0;
192:     VecMAXPY(mode[i], i, dots, mode);
193:     VecNormalize(mode[i], NULL);
194:   }
195:   MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
196:   for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
197:   return(0);
198: }

200: /*@
201:   DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs
202:   are computed by associating the basis function with one of the mesh points in its transitively-closed support, and
203:   evaluating the dual space basis of that point.  A basis function is associated with the point in its
204:   transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum
205:   projection height, which is set with this function.  By default, the maximum projection height is zero, which means
206:   that only mesh cells are used to project basis functions.  A height of one, for example, evaluates a cell-interior
207:   basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face.

209:   Input Parameters:
210: + dm - the DMPlex object
211: - height - the maximum projection height >= 0

213:   Level: advanced

215: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
216: @*/
217: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
218: {
219:   DM_Plex *plex = (DM_Plex *) dm->data;

223:   plex->maxProjectionHeight = height;
224:   return(0);
225: }

227: /*@
228:   DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in
229:   DMPlexProjectXXXLocal() functions.

231:   Input Parameters:
232: . dm - the DMPlex object

234:   Output Parameters:
235: . height - the maximum projection height

237:   Level: intermediate

239: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
240: @*/
241: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
242: {
243:   DM_Plex *plex = (DM_Plex *) dm->data;

247:   *height = plex->maxProjectionHeight;
248:   return(0);
249: }

251: /*@C
252:   DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector

254:   Input Parameters:
255: + dm     - The DM, with a PetscDS that matches the problem being constrained
256: . time   - The time
257: . field  - The field to constrain
258: . Nc     - The number of constrained field components, or 0 for all components
259: . comps  - An array of constrained component numbers, or NULL for all components
260: . label  - The DMLabel defining constrained points
261: . numids - The number of DMLabel ids for constrained points
262: . ids    - An array of ids for constrained points
263: . func   - A pointwise function giving boundary values
264: - ctx    - An optional user context for bcFunc

266:   Output Parameter:
267: . locX   - A local vector to receives the boundary values

269:   Level: developer

271: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
272: @*/
273: PetscErrorCode DMPlexInsertBoundaryValuesEssential(DM dm, PetscReal time, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX)
274: {
275:   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
276:   void            **ctxs;
277:   PetscInt          numFields;
278:   PetscErrorCode    ierr;

281:   DMGetNumFields(dm, &numFields);
282:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
283:   funcs[field] = func;
284:   ctxs[field]  = ctx;
285:   DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);
286:   PetscFree2(funcs,ctxs);
287:   return(0);
288: }

290: /*@C
291:   DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector

293:   Input Parameters:
294: + dm     - The DM, with a PetscDS that matches the problem being constrained
295: . time   - The time
296: . locU   - A local vector with the input solution values
297: . field  - The field to constrain
298: . Nc     - The number of constrained field components, or 0 for all components
299: . comps  - An array of constrained component numbers, or NULL for all components
300: . label  - The DMLabel defining constrained points
301: . numids - The number of DMLabel ids for constrained points
302: . ids    - An array of ids for constrained points
303: . func   - A pointwise function giving boundary values
304: - ctx    - An optional user context for bcFunc

306:   Output Parameter:
307: . locX   - A local vector to receives the boundary values

309:   Level: developer

311: .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary()
312: @*/
313: PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
314:                                                         void (*func)(PetscInt, PetscInt, PetscInt,
315:                                                                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
316:                                                                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
317:                                                                      PetscReal, const PetscReal[], PetscInt, const PetscScalar[],
318:                                                                      PetscScalar[]),
319:                                                         void *ctx, Vec locX)
320: {
321:   void (**funcs)(PetscInt, PetscInt, PetscInt,
322:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
323:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
324:                  PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
325:   void            **ctxs;
326:   PetscInt          numFields;
327:   PetscErrorCode    ierr;

330:   DMGetNumFields(dm, &numFields);
331:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
332:   funcs[field] = func;
333:   ctxs[field]  = ctx;
334:   DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
335:   PetscFree2(funcs,ctxs);
336:   return(0);
337: }

339: /*@C
340:   DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector

342:   Input Parameters:
343: + dm     - The DM, with a PetscDS that matches the problem being constrained
344: . time   - The time
345: . faceGeometry - A vector with the FVM face geometry information
346: . cellGeometry - A vector with the FVM cell geometry information
347: . Grad         - A vector with the FVM cell gradient information
348: . field  - The field to constrain
349: . Nc     - The number of constrained field components, or 0 for all components
350: . comps  - An array of constrained component numbers, or NULL for all components
351: . label  - The DMLabel defining constrained points
352: . numids - The number of DMLabel ids for constrained points
353: . ids    - An array of ids for constrained points
354: . func   - A pointwise function giving boundary values
355: - ctx    - An optional user context for bcFunc

357:   Output Parameter:
358: . locX   - A local vector to receives the boundary values

360:   Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary()

362:   Level: developer

364: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
365: @*/
366: PetscErrorCode DMPlexInsertBoundaryValuesRiemann(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
367:                                                  PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
368: {
369:   PetscDS            prob;
370:   PetscSF            sf;
371:   DM                 dmFace, dmCell, dmGrad;
372:   const PetscScalar *facegeom, *cellgeom = NULL, *grad;
373:   const PetscInt    *leaves;
374:   PetscScalar       *x, *fx;
375:   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
376:   PetscErrorCode     ierr, ierru = 0;

379:   DMGetPointSF(dm, &sf);
380:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
381:   nleaves = PetscMax(0, nleaves);
382:   DMGetDimension(dm, &dim);
383:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
384:   DMGetDS(dm, &prob);
385:   VecGetDM(faceGeometry, &dmFace);
386:   VecGetArrayRead(faceGeometry, &facegeom);
387:   if (cellGeometry) {
388:     VecGetDM(cellGeometry, &dmCell);
389:     VecGetArrayRead(cellGeometry, &cellgeom);
390:   }
391:   if (Grad) {
392:     PetscFV fv;

394:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
395:     VecGetDM(Grad, &dmGrad);
396:     VecGetArrayRead(Grad, &grad);
397:     PetscFVGetNumComponents(fv, &pdim);
398:     DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
399:   }
400:   VecGetArray(locX, &x);
401:   for (i = 0; i < numids; ++i) {
402:     IS              faceIS;
403:     const PetscInt *faces;
404:     PetscInt        numFaces, f;

406:     DMLabelGetStratumIS(label, ids[i], &faceIS);
407:     if (!faceIS) continue; /* No points with that id on this process */
408:     ISGetLocalSize(faceIS, &numFaces);
409:     ISGetIndices(faceIS, &faces);
410:     for (f = 0; f < numFaces; ++f) {
411:       const PetscInt         face = faces[f], *cells;
412:       PetscFVFaceGeom        *fg;

414:       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
415:       PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
416:       if (loc >= 0) continue;
417:       DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
418:       DMPlexGetSupport(dm, face, &cells);
419:       if (Grad) {
420:         PetscFVCellGeom       *cg;
421:         PetscScalar           *cx, *cgrad;
422:         PetscScalar           *xG;
423:         PetscReal              dx[3];
424:         PetscInt               d;

426:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
427:         DMPlexPointLocalRead(dm, cells[0], x, &cx);
428:         DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
429:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
430:         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
431:         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
432:         ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
433:         if (ierru) {
434:           ISRestoreIndices(faceIS, &faces);
435:           ISDestroy(&faceIS);
436:           goto cleanup;
437:         }
438:       } else {
439:         PetscScalar       *xI;
440:         PetscScalar       *xG;

442:         DMPlexPointLocalRead(dm, cells[0], x, &xI);
443:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
444:         ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
445:         if (ierru) {
446:           ISRestoreIndices(faceIS, &faces);
447:           ISDestroy(&faceIS);
448:           goto cleanup;
449:         }
450:       }
451:     }
452:     ISRestoreIndices(faceIS, &faces);
453:     ISDestroy(&faceIS);
454:   }
455:   cleanup:
456:   VecRestoreArray(locX, &x);
457:   if (Grad) {
458:     DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);
459:     VecRestoreArrayRead(Grad, &grad);
460:   }
461:   if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
462:   VecRestoreArrayRead(faceGeometry, &facegeom);
463:   CHKERRQ(ierru);
464:   return(0);
465: }

467: PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
468: {
469:   PetscInt       numBd, b;

473:   PetscDSGetNumBoundary(dm->prob, &numBd);
474:   for (b = 0; b < numBd; ++b) {
475:     DMBoundaryConditionType type;
476:     const char             *labelname;
477:     DMLabel                 label;
478:     PetscInt                field, Nc;
479:     const PetscInt         *comps;
480:     PetscObject             obj;
481:     PetscClassId            id;
482:     void                    (*func)(void);
483:     PetscInt                numids;
484:     const PetscInt         *ids;
485:     void                   *ctx;

487:     DMGetBoundary(dm, b, &type, NULL, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);
488:     if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
489:     DMGetLabel(dm, labelname, &label);
490:     DMGetField(dm, field, &obj);
491:     PetscObjectGetClassId(obj, &id);
492:     if (id == PETSCFE_CLASSID) {
493:       switch (type) {
494:         /* for FEM, there is no insertion to be done for non-essential boundary conditions */
495:       case DM_BC_ESSENTIAL:
496:         DMPlexLabelAddCells(dm,label);
497:         DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
498:         DMPlexLabelClearCells(dm,label);
499:         break;
500:       case DM_BC_ESSENTIAL_FIELD:
501:         DMPlexLabelAddCells(dm,label);
502:         DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
503:                                                         (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
504:                                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
505:                                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);
506:         DMPlexLabelClearCells(dm,label);
507:         break;
508:       default: break;
509:       }
510:     } else if (id == PETSCFV_CLASSID) {
511:       if (!faceGeomFVM) continue;
512:       DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids,
513:                                                (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
514:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
515:   }
516:   return(0);
517: }

519: /*@
520:   DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector

522:   Input Parameters:
523: + dm - The DM
524: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
525: . time - The time
526: . faceGeomFVM - Face geometry data for FV discretizations
527: . cellGeomFVM - Cell geometry data for FV discretizations
528: - gradFVM - Gradient reconstruction data for FV discretizations

530:   Output Parameters:
531: . locX - Solution updated with boundary values

533:   Level: developer

535: .seealso: DMProjectFunctionLabelLocal()
536: @*/
537: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
538: {

547:   PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));
548:   return(0);
549: }

551: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
552: {
553:   const PetscInt   debug = 0;
554:   PetscSection     section;
555:   PetscQuadrature  quad;
556:   Vec              localX;
557:   PetscScalar     *funcVal, *interpolant;
558:   PetscReal       *coords, *detJ, *J;
559:   PetscReal        localDiff = 0.0;
560:   const PetscReal *quadWeights;
561:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
562:   PetscErrorCode   ierr;

565:   DMGetDimension(dm, &dim);
566:   DMGetCoordinateDim(dm, &coordDim);
567:   DMGetDefaultSection(dm, &section);
568:   PetscSectionGetNumFields(section, &numFields);
569:   DMGetLocalVector(dm, &localX);
570:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
571:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
572:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
573:   for (field = 0; field < numFields; ++field) {
574:     PetscObject  obj;
575:     PetscClassId id;
576:     PetscInt     Nc;

578:     DMGetField(dm, field, &obj);
579:     PetscObjectGetClassId(obj, &id);
580:     if (id == PETSCFE_CLASSID) {
581:       PetscFE fe = (PetscFE) obj;

583:       PetscFEGetQuadrature(fe, &quad);
584:       PetscFEGetNumComponents(fe, &Nc);
585:     } else if (id == PETSCFV_CLASSID) {
586:       PetscFV fv = (PetscFV) obj;

588:       PetscFVGetQuadrature(fv, &quad);
589:       PetscFVGetNumComponents(fv, &Nc);
590:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
591:     numComponents += Nc;
592:   }
593:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
594:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
595:   PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
596:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
597:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
598:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
599:   for (c = cStart; c < cEnd; ++c) {
600:     PetscScalar *x = NULL;
601:     PetscReal    elemDiff = 0.0;
602:     PetscInt     qc = 0;

604:     DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
605:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

607:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
608:       PetscObject  obj;
609:       PetscClassId id;
610:       void * const ctx = ctxs ? ctxs[field] : NULL;
611:       PetscInt     Nb, Nc, q, fc;

613:       DMGetField(dm, field, &obj);
614:       PetscObjectGetClassId(obj, &id);
615:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
616:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
617:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
618:       if (debug) {
619:         char title[1024];
620:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
621:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
622:       }
623:       for (q = 0; q < Nq; ++q) {
624:         if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, point %D", detJ[q], c, q);
625:         (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx);
626:         if (ierr) {
627:           PetscErrorCode ierr2;
628:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
629:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
630:           ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
631: 
632:         }
633:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
634:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
635:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
636:         for (fc = 0; fc < Nc; ++fc) {
637:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
638:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
639:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
640:         }
641:       }
642:       fieldOffset += Nb;
643:       qc += Nc;
644:     }
645:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
646:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
647:     localDiff += elemDiff;
648:   }
649:   PetscFree5(funcVal,interpolant,coords,detJ,J);
650:   DMRestoreLocalVector(dm, &localX);
651:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
652:   *diff = PetscSqrtReal(*diff);
653:   return(0);
654: }

656: PetscErrorCode DMComputeL2GradientDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
657: {
658:   const PetscInt   debug = 0;
659:   PetscSection     section;
660:   PetscQuadrature  quad;
661:   Vec              localX;
662:   PetscScalar     *funcVal, *interpolant;
663:   const PetscReal *quadPoints, *quadWeights;
664:   PetscReal       *coords, *realSpaceDer, *J, *invJ, *detJ;
665:   PetscReal        localDiff = 0.0;
666:   PetscInt         dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset;
667:   PetscErrorCode   ierr;

670:   DMGetDimension(dm, &dim);
671:   DMGetCoordinateDim(dm, &coordDim);
672:   DMGetDefaultSection(dm, &section);
673:   PetscSectionGetNumFields(section, &numFields);
674:   DMGetLocalVector(dm, &localX);
675:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
676:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
677:   for (field = 0; field < numFields; ++field) {
678:     PetscFE  fe;
679:     PetscInt Nc;

681:     DMGetField(dm, field, (PetscObject *) &fe);
682:     PetscFEGetQuadrature(fe, &quad);
683:     PetscFEGetNumComponents(fe, &Nc);
684:     numComponents += Nc;
685:   }
686:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
687:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
688:   /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
689:   PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);
690:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
691:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
692:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
693:   for (c = cStart; c < cEnd; ++c) {
694:     PetscScalar *x = NULL;
695:     PetscReal    elemDiff = 0.0;
696:     PetscInt     qc = 0;

698:     DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);
699:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

701:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
702:       PetscFE          fe;
703:       void * const     ctx = ctxs ? ctxs[field] : NULL;
704:       PetscReal       *basisDer;
705:       PetscInt         Nb, Nc, q, fc;

707:       DMGetField(dm, field, (PetscObject *) &fe);
708:       PetscFEGetDimension(fe, &Nb);
709:       PetscFEGetNumComponents(fe, &Nc);
710:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
711:       if (debug) {
712:         char title[1024];
713:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
714:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
715:       }
716:       for (q = 0; q < Nq; ++q) {
717:         if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", detJ[q], c, q);
718:         (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx);
719:         if (ierr) {
720:           PetscErrorCode ierr2;
721:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
722:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
723:           ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2);
724: 
725:         }
726:         PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);
727:         for (fc = 0; fc < Nc; ++fc) {
728:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
729:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
730:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
731:         }
732:       }
733:       fieldOffset += Nb;
734:       qc          += Nc;
735:     }
736:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
737:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
738:     localDiff += elemDiff;
739:   }
740:   PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);
741:   DMRestoreLocalVector(dm, &localX);
742:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
743:   *diff = PetscSqrtReal(*diff);
744:   return(0);
745: }

747: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
748: {
749:   const PetscInt   debug = 0;
750:   PetscSection     section;
751:   PetscQuadrature  quad;
752:   Vec              localX;
753:   PetscScalar     *funcVal, *interpolant;
754:   PetscReal       *coords, *detJ, *J;
755:   PetscReal       *localDiff;
756:   const PetscReal *quadPoints, *quadWeights;
757:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
758:   PetscErrorCode   ierr;

761:   DMGetDimension(dm, &dim);
762:   DMGetCoordinateDim(dm, &coordDim);
763:   DMGetDefaultSection(dm, &section);
764:   PetscSectionGetNumFields(section, &numFields);
765:   DMGetLocalVector(dm, &localX);
766:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
767:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
768:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
769:   for (field = 0; field < numFields; ++field) {
770:     PetscObject  obj;
771:     PetscClassId id;
772:     PetscInt     Nc;

774:     DMGetField(dm, field, &obj);
775:     PetscObjectGetClassId(obj, &id);
776:     if (id == PETSCFE_CLASSID) {
777:       PetscFE fe = (PetscFE) obj;

779:       PetscFEGetQuadrature(fe, &quad);
780:       PetscFEGetNumComponents(fe, &Nc);
781:     } else if (id == PETSCFV_CLASSID) {
782:       PetscFV fv = (PetscFV) obj;

784:       PetscFVGetQuadrature(fv, &quad);
785:       PetscFVGetNumComponents(fv, &Nc);
786:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
787:     numComponents += Nc;
788:   }
789:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
790:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
791:   PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
792:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
793:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
794:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
795:   for (c = cStart; c < cEnd; ++c) {
796:     PetscScalar *x = NULL;
797:     PetscInt     qc = 0;

799:     DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
800:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

802:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
803:       PetscObject  obj;
804:       PetscClassId id;
805:       void * const ctx = ctxs ? ctxs[field] : NULL;
806:       PetscInt     Nb, Nc, q, fc;

808:       PetscReal       elemDiff = 0.0;

810:       DMGetField(dm, field, &obj);
811:       PetscObjectGetClassId(obj, &id);
812:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
813:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
814:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
815:       if (debug) {
816:         char title[1024];
817:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
818:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
819:       }
820:       for (q = 0; q < Nq; ++q) {
821:         if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature point %D", detJ, c, q);
822:         (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx);
823:         if (ierr) {
824:           PetscErrorCode ierr2;
825:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
826:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
827:           ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
828: 
829:         }
830:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
831:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
832:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
833:         for (fc = 0; fc < Nc; ++fc) {
834:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
835:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d point %g %g %g diff %g\n", c, field, coordDim > 0 ? coords[0] : 0., coordDim > 1 ? coords[1] : 0., coordDim > 2 ? coords[2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
836:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
837:         }
838:       }
839:       fieldOffset += Nb;
840:       qc          += Nc;
841:       localDiff[field] += elemDiff;
842:     }
843:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
844:   }
845:   DMRestoreLocalVector(dm, &localX);
846:   MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
847:   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
848:   PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);
849:   return(0);
850: }

852: /*@C
853:   DMPlexComputeL2DiffVec - This function computes the cellwise L_2 difference between a function u and an FEM interpolant solution u_h, and stores it in a Vec.

855:   Input Parameters:
856: + dm    - The DM
857: . time  - The time
858: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
859: . ctxs  - Optional array of contexts to pass to each function, or NULL.
860: - X     - The coefficient vector u_h

862:   Output Parameter:
863: . D - A Vec which holds the difference ||u - u_h||_2 for each cell

865:   Level: developer

867: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
868: @*/
869: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
870: {
871:   PetscSection     section;
872:   PetscQuadrature  quad;
873:   Vec              localX;
874:   PetscScalar     *funcVal, *interpolant;
875:   PetscReal       *coords, *detJ, *J;
876:   const PetscReal *quadPoints, *quadWeights;
877:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
878:   PetscErrorCode   ierr;

881:   VecSet(D, 0.0);
882:   DMGetDimension(dm, &dim);
883:   DMGetCoordinateDim(dm, &coordDim);
884:   DMGetDefaultSection(dm, &section);
885:   PetscSectionGetNumFields(section, &numFields);
886:   DMGetLocalVector(dm, &localX);
887:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
888:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
889:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
890:   for (field = 0; field < numFields; ++field) {
891:     PetscObject  obj;
892:     PetscClassId id;
893:     PetscInt     Nc;

895:     DMGetField(dm, field, &obj);
896:     PetscObjectGetClassId(obj, &id);
897:     if (id == PETSCFE_CLASSID) {
898:       PetscFE fe = (PetscFE) obj;

900:       PetscFEGetQuadrature(fe, &quad);
901:       PetscFEGetNumComponents(fe, &Nc);
902:     } else if (id == PETSCFV_CLASSID) {
903:       PetscFV fv = (PetscFV) obj;

905:       PetscFVGetQuadrature(fv, &quad);
906:       PetscFVGetNumComponents(fv, &Nc);
907:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
908:     numComponents += Nc;
909:   }
910:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
911:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
912:   PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
913:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
914:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
915:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
916:   for (c = cStart; c < cEnd; ++c) {
917:     PetscScalar *x = NULL;
918:     PetscScalar  elemDiff = 0.0;
919:     PetscInt     qc = 0;

921:     DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
922:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

924:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
925:       PetscObject  obj;
926:       PetscClassId id;
927:       void * const ctx = ctxs ? ctxs[field] : NULL;
928:       PetscInt     Nb, Nc, q, fc;

930:       DMGetField(dm, field, &obj);
931:       PetscObjectGetClassId(obj, &id);
932:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
933:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
934:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
935:       if (funcs[field]) {
936:         for (q = 0; q < Nq; ++q) {
937:           if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", detJ[q], c, q);
938:           (*funcs[field])(dim, time, &coords[q*coordDim], Nc, funcVal, ctx);
939:           if (ierr) {
940:             PetscErrorCode ierr2;
941:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
942:             ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);
943:             ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
944: 
945:           }
946:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
947:           else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
948:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
949:           for (fc = 0; fc < Nc; ++fc) {
950:             const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
951:             elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
952:           }
953:         }
954:       }
955:       fieldOffset += Nb;
956:       qc          += Nc;
957:     }
958:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
959:     VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
960:   }
961:   PetscFree5(funcVal,interpolant,coords,detJ,J);
962:   DMRestoreLocalVector(dm, &localX);
963:   VecSqrtAbs(D);
964:   return(0);
965: }

967: /*@
968:   DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user

970:   Input Parameters:
971: + dm - The mesh
972: . X  - Local input vector
973: - user - The user context

975:   Output Parameter:
976: . integral - Local integral for each field

978:   Level: developer

980: .seealso: DMPlexComputeResidualFEM()
981: @*/
982: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
983: {
984:   DM_Plex           *mesh  = (DM_Plex *) dm->data;
985:   DM                 dmAux, dmGrad;
986:   Vec                localX, A, cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
987:   PetscDS            prob, probAux = NULL;
988:   PetscSection       section, sectionAux;
989:   PetscFV            fvm = NULL;
990:   PetscFECellGeom   *cgeomFEM;
991:   PetscFVFaceGeom   *fgeomFVM;
992:   PetscFVCellGeom   *cgeomFVM;
993:   PetscScalar       *u, *a = NULL;
994:   const PetscScalar *constants, *lgrad;
995:   PetscReal         *lintegral;
996:   PetscInt          *uOff, *uOff_x, *aOff = NULL;
997:   PetscInt           dim, numConstants, Nf, NfAux = 0, f, numCells, cStart, cEnd, cEndInterior, c;
998:   PetscInt           totDim, totDimAux;
999:   PetscBool          useFVM = PETSC_FALSE;
1000:   PetscErrorCode     ierr;

1003:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1004:   DMGetLocalVector(dm, &localX);
1005:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);
1006:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1007:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1008:   DMGetDimension(dm, &dim);
1009:   DMGetDefaultSection(dm, &section);
1010:   DMGetDS(dm, &prob);
1011:   PetscDSGetTotalDimension(prob, &totDim);
1012:   PetscDSGetComponentOffsets(prob, &uOff);
1013:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1014:   PetscDSGetConstants(prob, &numConstants, &constants);
1015:   PetscSectionGetNumFields(section, &Nf);
1016:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1017:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1018:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1019:   numCells = cEnd - cStart;
1020:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1021:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1022:   if (dmAux) {
1023:     DMGetDS(dmAux, &probAux);
1024:     PetscDSGetNumFields(probAux, &NfAux);
1025:     DMGetDefaultSection(dmAux, &sectionAux);
1026:     PetscDSGetTotalDimension(probAux, &totDimAux);
1027:     PetscDSGetComponentOffsets(probAux, &aOff);
1028:     PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);
1029:   }
1030:   for (f = 0; f < Nf; ++f) {
1031:     PetscObject  obj;
1032:     PetscClassId id;

1034:     PetscDSGetDiscretization(prob, f, &obj);
1035:     PetscObjectGetClassId(obj, &id);
1036:     if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1037:   }
1038:   if (useFVM) {
1039:     Vec       grad;
1040:     PetscInt  fStart, fEnd;
1041:     PetscBool compGrad;

1043:     PetscFVGetComputeGradients(fvm, &compGrad);
1044:     PetscFVSetComputeGradients(fvm, PETSC_TRUE);
1045:     DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1046:     DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1047:     PetscFVSetComputeGradients(fvm, compGrad);
1048:     VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1049:     VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1050:     /* Reconstruct and limit cell gradients */
1051:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1052:     DMGetGlobalVector(dmGrad, &grad);
1053:     DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);
1054:     /* Communicate gradient values */
1055:     DMGetLocalVector(dmGrad, &locGrad);
1056:     DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1057:     DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1058:     DMRestoreGlobalVector(dmGrad, &grad);
1059:     /* Handle non-essential (e.g. outflow) boundary values */
1060:     DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1061:     VecGetArrayRead(locGrad, &lgrad);
1062:   }
1063:   PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);
1064:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1065:   for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
1066:   for (c = cStart; c < cEnd; ++c) {
1067:     PetscScalar *x = NULL;
1068:     PetscInt     i;

1070:     DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);
1071:     if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c);
1072:     DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
1073:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1074:     DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
1075:     if (dmAux) {
1076:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1077:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1078:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1079:     }
1080:   }
1081:   for (f = 0; f < Nf; ++f) {
1082:     PetscObject  obj;
1083:     PetscClassId id;
1084:     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1086:     PetscDSGetDiscretization(prob, f, &obj);
1087:     PetscObjectGetClassId(obj, &id);
1088:     if (id == PETSCFE_CLASSID) {
1089:       PetscFE         fe = (PetscFE) obj;
1090:       PetscQuadrature q;
1091:       PetscInt        Nq, Nb;

1093:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1094:       PetscFEGetQuadrature(fe, &q);
1095:       PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1096:       PetscFEGetDimension(fe, &Nb);
1097:       blockSize = Nb*Nq;
1098:       batchSize = numBlocks * blockSize;
1099:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1100:       numChunks = numCells / (numBatches*batchSize);
1101:       Ne        = numChunks*numBatches*batchSize;
1102:       Nr        = numCells % (numBatches*batchSize);
1103:       offset    = numCells - Nr;
1104:       PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);
1105:       PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);
1106:     } else if (id == PETSCFV_CLASSID) {
1107:       /* PetscFV  fv = (PetscFV) obj; */
1108:       PetscInt       foff;
1109:       PetscPointFunc obj_func;
1110:       PetscScalar    lint;

1112:       PetscDSGetObjective(prob, f, &obj_func);
1113:       PetscDSGetFieldOffset(prob, f, &foff);
1114:       if (obj_func) {
1115:         for (c = 0; c < numCells; ++c) {
1116:           PetscScalar *u_x;

1118:           DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1119:           obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, a, NULL, NULL, 0.0, cgeomFVM[c].centroid, numConstants, constants, &lint);
1120:           lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1121:         }
1122:       }
1123:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1124:   }
1125:   if (useFVM) {
1126:     VecRestoreArrayRead(locGrad, &lgrad);
1127:     VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1128:     VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1129:     DMRestoreLocalVector(dmGrad, &locGrad);
1130:     VecDestroy(&faceGeometryFVM);
1131:     VecDestroy(&cellGeometryFVM);
1132:     DMDestroy(&dmGrad);
1133:   }
1134:   if (dmAux) {PetscFree(a);}
1135:   if (mesh->printFEM) {
1136:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1137:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);}
1138:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1139:   }
1140:   DMRestoreLocalVector(dm, &localX);
1141:   MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1142:   PetscFree3(lintegral,u,cgeomFEM);
1143:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1144:   return(0);
1145: }

1147: /*@
1148:   DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.

1150:   Input Parameters:
1151: + dmf  - The fine mesh
1152: . dmc  - The coarse mesh
1153: - user - The user context

1155:   Output Parameter:
1156: . In  - The interpolation matrix

1158:   Level: developer

1160: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1161: @*/
1162: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1163: {
1164:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1165:   const char       *name  = "Interpolator";
1166:   PetscDS           prob;
1167:   PetscFE          *feRef;
1168:   PetscFV          *fvRef;
1169:   PetscSection      fsection, fglobalSection;
1170:   PetscSection      csection, cglobalSection;
1171:   PetscScalar      *elemMat;
1172:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1173:   PetscInt          cTotDim, rTotDim = 0;
1174:   PetscErrorCode    ierr;

1177:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1178:   DMGetDimension(dmf, &dim);
1179:   DMGetDefaultSection(dmf, &fsection);
1180:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1181:   DMGetDefaultSection(dmc, &csection);
1182:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1183:   PetscSectionGetNumFields(fsection, &Nf);
1184:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1185:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1186:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1187:   DMGetDS(dmf, &prob);
1188:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1189:   for (f = 0; f < Nf; ++f) {
1190:     PetscObject  obj;
1191:     PetscClassId id;
1192:     PetscInt     rNb = 0, Nc = 0;

1194:     PetscDSGetDiscretization(prob, f, &obj);
1195:     PetscObjectGetClassId(obj, &id);
1196:     if (id == PETSCFE_CLASSID) {
1197:       PetscFE fe = (PetscFE) obj;

1199:       PetscFERefine(fe, &feRef[f]);
1200:       PetscFEGetDimension(feRef[f], &rNb);
1201:       PetscFEGetNumComponents(fe, &Nc);
1202:     } else if (id == PETSCFV_CLASSID) {
1203:       PetscFV        fv = (PetscFV) obj;
1204:       PetscDualSpace Q;

1206:       PetscFVRefine(fv, &fvRef[f]);
1207:       PetscFVGetDualSpace(fvRef[f], &Q);
1208:       PetscDualSpaceGetDimension(Q, &rNb);
1209:       PetscFVGetNumComponents(fv, &Nc);
1210:     }
1211:     rTotDim += rNb;
1212:   }
1213:   PetscDSGetTotalDimension(prob, &cTotDim);
1214:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1215:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1216:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1217:     PetscDualSpace   Qref;
1218:     PetscQuadrature  f;
1219:     const PetscReal *qpoints, *qweights;
1220:     PetscReal       *points;
1221:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1223:     /* Compose points from all dual basis functionals */
1224:     if (feRef[fieldI]) {
1225:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
1226:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
1227:     } else {
1228:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1229:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1230:     }
1231:     PetscDualSpaceGetDimension(Qref, &fpdim);
1232:     for (i = 0; i < fpdim; ++i) {
1233:       PetscDualSpaceGetFunctional(Qref, i, &f);
1234:       PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
1235:       npoints += Np;
1236:     }
1237:     PetscMalloc1(npoints*dim,&points);
1238:     for (i = 0, k = 0; i < fpdim; ++i) {
1239:       PetscDualSpaceGetFunctional(Qref, i, &f);
1240:       PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1241:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1242:     }

1244:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1245:       PetscObject  obj;
1246:       PetscClassId id;
1247:       PetscReal   *B;
1248:       PetscInt     NcJ = 0, cpdim = 0, j, qNc;

1250:       PetscDSGetDiscretization(prob, fieldJ, &obj);
1251:       PetscObjectGetClassId(obj, &id);
1252:       if (id == PETSCFE_CLASSID) {
1253:         PetscFE fe = (PetscFE) obj;

1255:         /* Evaluate basis at points */
1256:         PetscFEGetNumComponents(fe, &NcJ);
1257:         PetscFEGetDimension(fe, &cpdim);
1258:         /* For now, fields only interpolate themselves */
1259:         if (fieldI == fieldJ) {
1260:           if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %D does not match coarse field %D", Nc, NcJ);
1261:           PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1262:           for (i = 0, k = 0; i < fpdim; ++i) {
1263:             PetscDualSpaceGetFunctional(Qref, i, &f);
1264:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1265:             if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ);
1266:             for (p = 0; p < Np; ++p, ++k) {
1267:               for (j = 0; j < cpdim; ++j) {
1268:                 /* NOTE: This is not quite right, unless fpdim == number of fine grid functional quad points */
1269:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
1270:               }
1271:             }
1272:           }
1273:           PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1274:         }
1275:       } else if (id == PETSCFV_CLASSID) {
1276:         PetscFV        fv = (PetscFV) obj;

1278:         /* Evaluate constant function at points */
1279:         PetscFVGetNumComponents(fv, &NcJ);
1280:         cpdim = 1;
1281:         /* For now, fields only interpolate themselves */
1282:         if (fieldI == fieldJ) {
1283:           if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ);
1284:           for (i = 0, k = 0; i < fpdim; ++i) {
1285:             PetscDualSpaceGetFunctional(Qref, i, &f);
1286:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1287:             if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ);
1288:             for (p = 0; p < Np; ++p, ++k) {
1289:               for (j = 0; j < cpdim; ++j) {
1290:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c];
1291:               }
1292:             }
1293:           }
1294:         }
1295:       }
1296:       offsetJ += cpdim*NcJ;
1297:     }
1298:     offsetI += fpdim*Nc;
1299:     PetscFree(points);
1300:   }
1301:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1302:   /* Preallocate matrix */
1303:   {
1304:     Mat          preallocator;
1305:     PetscScalar *vals;
1306:     PetscInt    *cellCIndices, *cellFIndices;
1307:     PetscInt     locRows, locCols, cell;

1309:     MatGetLocalSize(In, &locRows, &locCols);
1310:     MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1311:     MatSetType(preallocator, MATPREALLOCATOR);
1312:     MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1313:     MatSetUp(preallocator);
1314:     PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1315:     for (cell = cStart; cell < cEnd; ++cell) {
1316:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1317:       MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1318:     }
1319:     PetscFree3(vals,cellCIndices,cellFIndices);
1320:     MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1321:     MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1322:     MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1323:     MatDestroy(&preallocator);
1324:   }
1325:   /* Fill matrix */
1326:   MatZeroEntries(In);
1327:   for (c = cStart; c < cEnd; ++c) {
1328:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1329:   }
1330:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1331:   PetscFree2(feRef,fvRef);
1332:   PetscFree(elemMat);
1333:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1334:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1335:   if (mesh->printFEM) {
1336:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1337:     MatChop(In, 1.0e-10);
1338:     MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1339:   }
1340:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1341:   return(0);
1342: }

1344: /*@
1345:   DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM.

1347:   Input Parameters:
1348: + dmf  - The fine mesh
1349: . dmc  - The coarse mesh
1350: - user - The user context

1352:   Output Parameter:
1353: . In  - The interpolation matrix

1355:   Level: developer

1357: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1358: @*/
1359: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1360: {
1361:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
1362:   const char    *name = "Interpolator";
1363:   PetscDS        prob;
1364:   PetscSection   fsection, csection, globalFSection, globalCSection;
1365:   PetscHashJK    ht;
1366:   PetscLayout    rLayout;
1367:   PetscInt      *dnz, *onz;
1368:   PetscInt       locRows, rStart, rEnd;
1369:   PetscReal     *x, *v0, *J, *invJ, detJ;
1370:   PetscReal     *v0c, *Jc, *invJc, detJc;
1371:   PetscScalar   *elemMat;
1372:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

1376:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1377:   DMGetCoordinateDim(dmc, &dim);
1378:   DMGetDS(dmc, &prob);
1379:   PetscDSGetRefCoordArrays(prob, &x, NULL);
1380:   PetscDSGetNumFields(prob, &Nf);
1381:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1382:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1383:   DMGetDefaultSection(dmf, &fsection);
1384:   DMGetDefaultGlobalSection(dmf, &globalFSection);
1385:   DMGetDefaultSection(dmc, &csection);
1386:   DMGetDefaultGlobalSection(dmc, &globalCSection);
1387:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1388:   PetscDSGetTotalDimension(prob, &totDim);
1389:   PetscMalloc1(totDim, &elemMat);

1391:   MatGetLocalSize(In, &locRows, NULL);
1392:   PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1393:   PetscLayoutSetLocalSize(rLayout, locRows);
1394:   PetscLayoutSetBlockSize(rLayout, 1);
1395:   PetscLayoutSetUp(rLayout);
1396:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1397:   PetscLayoutDestroy(&rLayout);
1398:   PetscCalloc2(locRows,&dnz,locRows,&onz);
1399:   PetscHashJKCreate(&ht);
1400:   for (field = 0; field < Nf; ++field) {
1401:     PetscObject      obj;
1402:     PetscClassId     id;
1403:     PetscDualSpace   Q = NULL;
1404:     PetscQuadrature  f;
1405:     const PetscReal *qpoints;
1406:     PetscInt         Nc, Np, fpdim, i, d;

1408:     PetscDSGetDiscretization(prob, field, &obj);
1409:     PetscObjectGetClassId(obj, &id);
1410:     if (id == PETSCFE_CLASSID) {
1411:       PetscFE fe = (PetscFE) obj;

1413:       PetscFEGetDualSpace(fe, &Q);
1414:       PetscFEGetNumComponents(fe, &Nc);
1415:     } else if (id == PETSCFV_CLASSID) {
1416:       PetscFV fv = (PetscFV) obj;

1418:       PetscFVGetDualSpace(fv, &Q);
1419:       Nc   = 1;
1420:     }
1421:     PetscDualSpaceGetDimension(Q, &fpdim);
1422:     /* For each fine grid cell */
1423:     for (cell = cStart; cell < cEnd; ++cell) {
1424:       PetscInt *findices,   *cindices;
1425:       PetscInt  numFIndices, numCIndices;

1427:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1428:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1429:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
1430:       for (i = 0; i < fpdim; ++i) {
1431:         Vec             pointVec;
1432:         PetscScalar    *pV;
1433:         PetscSF         coarseCellSF = NULL;
1434:         const PetscSFNode *coarseCells;
1435:         PetscInt        numCoarseCells, q, c;

1437:         /* Get points from the dual basis functional quadrature */
1438:         PetscDualSpaceGetFunctional(Q, i, &f);
1439:         PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1440:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1441:         VecSetBlockSize(pointVec, dim);
1442:         VecGetArray(pointVec, &pV);
1443:         for (q = 0; q < Np; ++q) {
1444:           /* Transform point to real space */
1445:           CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1446:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1447:         }
1448:         VecRestoreArray(pointVec, &pV);
1449:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1450:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1451:         PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
1452:         /* Update preallocation info */
1453:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1454:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1455:         {
1456:           PetscHashJKKey  key;
1457:           PetscHashJKIter missing, iter;

1459:           key.j = findices[i];
1460:           if (key.j >= 0) {
1461:             /* Get indices for coarse elements */
1462:             for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1463:               DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1464:               for (c = 0; c < numCIndices; ++c) {
1465:                 key.k = cindices[c];
1466:                 if (key.k < 0) continue;
1467:                 PetscHashJKPut(ht, key, &missing, &iter);
1468:                 if (missing) {
1469:                   PetscHashJKSet(ht, iter, 1);
1470:                   if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1471:                   else                                     ++onz[key.j-rStart];
1472:                 }
1473:               }
1474:               DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1475:             }
1476:           }
1477:         }
1478:         PetscSFDestroy(&coarseCellSF);
1479:         VecDestroy(&pointVec);
1480:       }
1481:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1482:     }
1483:   }
1484:   PetscHashJKDestroy(&ht);
1485:   MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1486:   MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1487:   PetscFree2(dnz,onz);
1488:   for (field = 0; field < Nf; ++field) {
1489:     PetscObject      obj;
1490:     PetscClassId     id;
1491:     PetscDualSpace   Q = NULL;
1492:     PetscQuadrature  f;
1493:     const PetscReal *qpoints, *qweights;
1494:     PetscInt         Nc, qNc, Np, fpdim, i, d;

1496:     PetscDSGetDiscretization(prob, field, &obj);
1497:     PetscObjectGetClassId(obj, &id);
1498:     if (id == PETSCFE_CLASSID) {
1499:       PetscFE fe = (PetscFE) obj;

1501:       PetscFEGetDualSpace(fe, &Q);
1502:       PetscFEGetNumComponents(fe, &Nc);
1503:     } else if (id == PETSCFV_CLASSID) {
1504:       PetscFV fv = (PetscFV) obj;

1506:       PetscFVGetDualSpace(fv, &Q);
1507:       Nc   = 1;
1508:     }
1509:     PetscDualSpaceGetDimension(Q, &fpdim);
1510:     /* For each fine grid cell */
1511:     for (cell = cStart; cell < cEnd; ++cell) {
1512:       PetscInt *findices,   *cindices;
1513:       PetscInt  numFIndices, numCIndices;

1515:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1516:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1517:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
1518:       for (i = 0; i < fpdim; ++i) {
1519:         Vec             pointVec;
1520:         PetscScalar    *pV;
1521:         PetscSF         coarseCellSF = NULL;
1522:         const PetscSFNode *coarseCells;
1523:         PetscInt        numCoarseCells, cpdim, q, c, j;

1525:         /* Get points from the dual basis functional quadrature */
1526:         PetscDualSpaceGetFunctional(Q, i, &f);
1527:         PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
1528:         if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, Nc);
1529:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1530:         VecSetBlockSize(pointVec, dim);
1531:         VecGetArray(pointVec, &pV);
1532:         for (q = 0; q < Np; ++q) {
1533:           /* Transform point to real space */
1534:           CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1535:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1536:         }
1537:         VecRestoreArray(pointVec, &pV);
1538:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1539:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1540:         /* Update preallocation info */
1541:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1542:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1543:         VecGetArray(pointVec, &pV);
1544:         for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1545:           PetscReal pVReal[3];

1547:           DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1548:           /* Transform points from real space to coarse reference space */
1549:           DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
1550:           for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
1551:           CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x);

1553:           if (id == PETSCFE_CLASSID) {
1554:             PetscFE    fe = (PetscFE) obj;
1555:             PetscReal *B;

1557:             /* Evaluate coarse basis on contained point */
1558:             PetscFEGetDimension(fe, &cpdim);
1559:             PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
1560:             PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
1561:             /* Get elemMat entries by multiplying by weight */
1562:             for (j = 0; j < cpdim; ++j) {
1563:               for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c];
1564:             }
1565:             PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
1566:           } else {
1567:             cpdim = 1;
1568:             for (j = 0; j < cpdim; ++j) {
1569:               for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
1570:             }
1571:           }
1572:           /* Update interpolator */
1573:           if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
1574:           if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
1575:           MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
1576:           DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1577:         }
1578:         VecRestoreArray(pointVec, &pV);
1579:         PetscSFDestroy(&coarseCellSF);
1580:         VecDestroy(&pointVec);
1581:       }
1582:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1583:     }
1584:   }
1585:   PetscFree3(v0,J,invJ);
1586:   PetscFree3(v0c,Jc,invJc);
1587:   PetscFree(elemMat);
1588:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1589:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1590:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1591:   return(0);
1592: }

1594: /*@
1595:   DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns

1597:   Input Parameters:
1598: + dmc  - The coarse mesh
1599: - dmf  - The fine mesh
1600: - user - The user context

1602:   Output Parameter:
1603: . sc   - The mapping

1605:   Level: developer

1607: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1608: @*/
1609: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1610: {
1611:   PetscDS        prob;
1612:   PetscFE       *feRef;
1613:   PetscFV       *fvRef;
1614:   Vec            fv, cv;
1615:   IS             fis, cis;
1616:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1617:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1618:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;

1622:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1623:   DMGetDimension(dmf, &dim);
1624:   DMGetDefaultSection(dmf, &fsection);
1625:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1626:   DMGetDefaultSection(dmc, &csection);
1627:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1628:   PetscSectionGetNumFields(fsection, &Nf);
1629:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1630:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1631:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1632:   DMGetDS(dmc, &prob);
1633:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1634:   for (f = 0; f < Nf; ++f) {
1635:     PetscObject  obj;
1636:     PetscClassId id;
1637:     PetscInt     fNb = 0, Nc = 0;

1639:     PetscDSGetDiscretization(prob, f, &obj);
1640:     PetscObjectGetClassId(obj, &id);
1641:     if (id == PETSCFE_CLASSID) {
1642:       PetscFE fe = (PetscFE) obj;

1644:       PetscFERefine(fe, &feRef[f]);
1645:       PetscFEGetDimension(feRef[f], &fNb);
1646:       PetscFEGetNumComponents(fe, &Nc);
1647:     } else if (id == PETSCFV_CLASSID) {
1648:       PetscFV        fv = (PetscFV) obj;
1649:       PetscDualSpace Q;

1651:       PetscFVRefine(fv, &fvRef[f]);
1652:       PetscFVGetDualSpace(fvRef[f], &Q);
1653:       PetscDualSpaceGetDimension(Q, &fNb);
1654:       PetscFVGetNumComponents(fv, &Nc);
1655:     }
1656:     fTotDim += fNb*Nc;
1657:   }
1658:   PetscDSGetTotalDimension(prob, &cTotDim);
1659:   PetscMalloc1(cTotDim,&cmap);
1660:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1661:     PetscFE        feC;
1662:     PetscFV        fvC;
1663:     PetscDualSpace QF, QC;
1664:     PetscInt       NcF, NcC, fpdim, cpdim;

1666:     if (feRef[field]) {
1667:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1668:       PetscFEGetNumComponents(feC, &NcC);
1669:       PetscFEGetNumComponents(feRef[field], &NcF);
1670:       PetscFEGetDualSpace(feRef[field], &QF);
1671:       PetscDualSpaceGetDimension(QF, &fpdim);
1672:       PetscFEGetDualSpace(feC, &QC);
1673:       PetscDualSpaceGetDimension(QC, &cpdim);
1674:     } else {
1675:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
1676:       PetscFVGetNumComponents(fvC, &NcC);
1677:       PetscFVGetNumComponents(fvRef[field], &NcF);
1678:       PetscFVGetDualSpace(fvRef[field], &QF);
1679:       PetscDualSpaceGetDimension(QF, &fpdim);
1680:       PetscFVGetDualSpace(fvC, &QC);
1681:       PetscDualSpaceGetDimension(QC, &cpdim);
1682:     }
1683:     if (NcF != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", NcF, NcC);
1684:     for (c = 0; c < cpdim; ++c) {
1685:       PetscQuadrature  cfunc;
1686:       const PetscReal *cqpoints;
1687:       PetscInt         NpC;
1688:       PetscBool        found = PETSC_FALSE;

1690:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
1691:       PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);
1692:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1693:       for (f = 0; f < fpdim; ++f) {
1694:         PetscQuadrature  ffunc;
1695:         const PetscReal *fqpoints;
1696:         PetscReal        sum = 0.0;
1697:         PetscInt         NpF, comp;

1699:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
1700:         PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);
1701:         if (NpC != NpF) continue;
1702:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1703:         if (sum > 1.0e-9) continue;
1704:         for (comp = 0; comp < NcC; ++comp) {
1705:           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1706:         }
1707:         found = PETSC_TRUE;
1708:         break;
1709:       }
1710:       if (!found) {
1711:         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1712:         if (fvRef[field]) {
1713:           PetscInt comp;
1714:           for (comp = 0; comp < NcC; ++comp) {
1715:             cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1716:           }
1717:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1718:       }
1719:     }
1720:     offsetC += cpdim*NcC;
1721:     offsetF += fpdim*NcF;
1722:   }
1723:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
1724:   PetscFree2(feRef,fvRef);

1726:   DMGetGlobalVector(dmf, &fv);
1727:   DMGetGlobalVector(dmc, &cv);
1728:   VecGetOwnershipRange(cv, &startC, &endC);
1729:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1730:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1731:   PetscMalloc1(m,&cindices);
1732:   PetscMalloc1(m,&findices);
1733:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1734:   for (c = cStart; c < cEnd; ++c) {
1735:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1736:     for (d = 0; d < cTotDim; ++d) {
1737:       if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
1738:       if ((findices[cellCIndices[d]-startC] >= 0) && (findices[cellCIndices[d]-startC] != cellFIndices[cmap[d]])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coarse dof %d maps to both %d and %d", cindices[cellCIndices[d]-startC], findices[cellCIndices[d]-startC], cellFIndices[cmap[d]]);
1739:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1740:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1741:     }
1742:   }
1743:   PetscFree(cmap);
1744:   PetscFree2(cellCIndices,cellFIndices);

1746:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1747:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1748:   VecScatterCreate(cv, cis, fv, fis, sc);
1749:   ISDestroy(&cis);
1750:   ISDestroy(&fis);
1751:   DMRestoreGlobalVector(dmf, &fv);
1752:   DMRestoreGlobalVector(dmc, &cv);
1753:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1754:   return(0);
1755: }