Actual source code: plexfem.c

petsc-master 2018-05-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: static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx)
 59: {
 60:   const PetscInt eps[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
 61:   PetscInt *ctxInt  = (PetscInt *) ctx;
 62:   PetscInt  dim2    = ctxInt[0];
 63:   PetscInt  d       = ctxInt[1];
 64:   PetscInt  i, j, k = dim > 2 ? d - dim : d;

 67:   if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2);
 68:   for (i = 0; i < dim; i++) mode[i] = 0.;
 69:   if (d < dim) {
 70:     mode[d] = 1.; /* Translation along axis d */
 71:   } else {
 72:     for (i = 0; i < dim; i++) {
 73:       for (j = 0; j < dim; j++) {
 74:         mode[j] += eps[i][j][k]*X[i]; /* Rotation about axis d */
 75:       }
 76:     }
 77:   }
 78:   return(0);
 79: }

 81: /*@C
 82:   DMPlexCreateRigidBody - For the default global section, create rigid body modes by function space interpolation

 84:   Collective on DM

 86:   Input Arguments:
 87: . dm - the DM

 89:   Output Argument:
 90: . sp - the null space

 92:   Note: This is necessary to provide a suitable coarse space for algebraic multigrid

 94:   Level: advanced

 96: .seealso: MatNullSpaceCreate(), PCGAMG
 97: @*/
 98: PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
 99: {
100:   MPI_Comm       comm;
101:   Vec            mode[6];
102:   PetscSection   section, globalSection;
103:   PetscInt       dim, dimEmbed, n, m, d, i, j;

107:   PetscObjectGetComm((PetscObject)dm,&comm);
108:   DMGetDimension(dm, &dim);
109:   DMGetCoordinateDim(dm, &dimEmbed);
110:   if (dim == 1) {
111:     MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
112:     return(0);
113:   }
114:   DMGetDefaultSection(dm, &section);
115:   DMGetDefaultGlobalSection(dm, &globalSection);
116:   PetscSectionGetConstrainedStorageSize(globalSection, &n);
117:   m    = (dim*(dim+1))/2;
118:   VecCreate(comm, &mode[0]);
119:   VecSetSizes(mode[0], n, PETSC_DETERMINE);
120:   VecSetUp(mode[0]);
121:   for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
122:   for (d = 0; d < m; d++) {
123:     PetscInt         ctx[2];
124:     PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
125:     void            *voidctx = (void *) (&ctx[0]);

127:     ctx[0] = dimEmbed;
128:     ctx[1] = d;
129:     DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);
130:   }
131:   for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
132:   /* Orthonormalize system */
133:   for (i = dim; i < m; ++i) {
134:     PetscScalar dots[6];

136:     VecMDot(mode[i], i, mode, dots);
137:     for (j = 0; j < i; ++j) dots[j] *= -1.0;
138:     VecMAXPY(mode[i], i, dots, mode);
139:     VecNormalize(mode[i], NULL);
140:   }
141:   MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
142:   for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
143:   return(0);
144: }

146: /*@C
147:   DMPlexCreateRigidBodies - For the default global section, create rigid body modes by function space interpolation

149:   Collective on DM

151:   Input Arguments:
152: + dm    - the DM
153: . nb    - The number of bodies
154: . label - The DMLabel marking each domain
155: . nids  - The number of ids per body
156: - ids   - An array of the label ids in sequence for each domain

158:   Output Argument:
159: . sp - the null space

161:   Note: This is necessary to provide a suitable coarse space for algebraic multigrid

163:   Level: advanced

165: .seealso: MatNullSpaceCreate()
166: @*/
167: PetscErrorCode DMPlexCreateRigidBodies(DM dm, PetscInt nb, DMLabel label, const PetscInt nids[], const PetscInt ids[], MatNullSpace *sp)
168: {
169:   MPI_Comm       comm;
170:   PetscSection   section, globalSection;
171:   Vec           *mode;
172:   PetscScalar   *dots;
173:   PetscInt       dim, dimEmbed, n, m, b, d, i, j, off;

177:   PetscObjectGetComm((PetscObject)dm,&comm);
178:   DMGetDimension(dm, &dim);
179:   DMGetCoordinateDim(dm, &dimEmbed);
180:   DMGetDefaultSection(dm, &section);
181:   DMGetDefaultGlobalSection(dm, &globalSection);
182:   PetscSectionGetConstrainedStorageSize(globalSection, &n);
183:   m    = nb * (dim*(dim+1))/2;
184:   PetscMalloc2(m, &mode, m, &dots);
185:   VecCreate(comm, &mode[0]);
186:   VecSetSizes(mode[0], n, PETSC_DETERMINE);
187:   VecSetUp(mode[0]);
188:   for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
189:   for (b = 0, off = 0; b < nb; ++b) {
190:     for (d = 0; d < m/nb; ++d) {
191:       PetscInt         ctx[2];
192:       PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
193:       void            *voidctx = (void *) (&ctx[0]);

195:       ctx[0] = dimEmbed;
196:       ctx[1] = d;
197:       DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);
198:       off   += nids[b];
199:     }
200:   }
201:   for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
202:   /* Orthonormalize system */
203:   for (i = 0; i < m; ++i) {
204:     VecMDot(mode[i], i, mode, dots);
205:     for (j = 0; j < i; ++j) dots[j] *= -1.0;
206:     VecMAXPY(mode[i], i, dots, mode);
207:     VecNormalize(mode[i], NULL);
208:   }
209:   MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
210:   for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
211:   PetscFree2(mode, dots);
212:   return(0);
213: }

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

224:   Input Parameters:
225: + dm - the DMPlex object
226: - height - the maximum projection height >= 0

228:   Level: advanced

230: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
231: @*/
232: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
233: {
234:   DM_Plex *plex = (DM_Plex *) dm->data;

238:   plex->maxProjectionHeight = height;
239:   return(0);
240: }

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

246:   Input Parameters:
247: . dm - the DMPlex object

249:   Output Parameters:
250: . height - the maximum projection height

252:   Level: intermediate

254: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
255: @*/
256: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
257: {
258:   DM_Plex *plex = (DM_Plex *) dm->data;

262:   *height = plex->maxProjectionHeight;
263:   return(0);
264: }

266: /*@C
267:   DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector

269:   Input Parameters:
270: + dm     - The DM, with a PetscDS that matches the problem being constrained
271: . time   - The time
272: . field  - The field to constrain
273: . Nc     - The number of constrained field components, or 0 for all components
274: . comps  - An array of constrained component numbers, or NULL for all components
275: . label  - The DMLabel defining constrained points
276: . numids - The number of DMLabel ids for constrained points
277: . ids    - An array of ids for constrained points
278: . func   - A pointwise function giving boundary values
279: - ctx    - An optional user context for bcFunc

281:   Output Parameter:
282: . locX   - A local vector to receives the boundary values

284:   Level: developer

286: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
287: @*/
288: 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)
289: {
290:   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
291:   void            **ctxs;
292:   PetscInt          numFields;
293:   PetscErrorCode    ierr;

296:   DMGetNumFields(dm, &numFields);
297:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
298:   funcs[field] = func;
299:   ctxs[field]  = ctx;
300:   DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);
301:   PetscFree2(funcs,ctxs);
302:   return(0);
303: }

305: /*@C
306:   DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector

308:   Input Parameters:
309: + dm     - The DM, with a PetscDS that matches the problem being constrained
310: . time   - The time
311: . locU   - A local vector with the input solution values
312: . field  - The field to constrain
313: . Nc     - The number of constrained field components, or 0 for all components
314: . comps  - An array of constrained component numbers, or NULL for all components
315: . label  - The DMLabel defining constrained points
316: . numids - The number of DMLabel ids for constrained points
317: . ids    - An array of ids for constrained points
318: . func   - A pointwise function giving boundary values
319: - ctx    - An optional user context for bcFunc

321:   Output Parameter:
322: . locX   - A local vector to receives the boundary values

324:   Level: developer

326: .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary()
327: @*/
328: PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
329:                                                         void (*func)(PetscInt, PetscInt, PetscInt,
330:                                                                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
331:                                                                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
332:                                                                      PetscReal, const PetscReal[], PetscInt, const PetscScalar[],
333:                                                                      PetscScalar[]),
334:                                                         void *ctx, Vec locX)
335: {
336:   void (**funcs)(PetscInt, PetscInt, PetscInt,
337:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
338:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
339:                  PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
340:   void            **ctxs;
341:   PetscInt          numFields;
342:   PetscErrorCode    ierr;

345:   DMGetNumFields(dm, &numFields);
346:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
347:   funcs[field] = func;
348:   ctxs[field]  = ctx;
349:   DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
350:   PetscFree2(funcs,ctxs);
351:   return(0);
352: }

354: /*@C
355:   DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector

357:   Input Parameters:
358: + dm     - The DM, with a PetscDS that matches the problem being constrained
359: . time   - The time
360: . faceGeometry - A vector with the FVM face geometry information
361: . cellGeometry - A vector with the FVM cell geometry information
362: . Grad         - A vector with the FVM cell gradient information
363: . field  - The field to constrain
364: . Nc     - The number of constrained field components, or 0 for all components
365: . comps  - An array of constrained component numbers, or NULL for all components
366: . label  - The DMLabel defining constrained points
367: . numids - The number of DMLabel ids for constrained points
368: . ids    - An array of ids for constrained points
369: . func   - A pointwise function giving boundary values
370: - ctx    - An optional user context for bcFunc

372:   Output Parameter:
373: . locX   - A local vector to receives the boundary values

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

377:   Level: developer

379: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
380: @*/
381: 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[],
382:                                                  PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
383: {
384:   PetscDS            prob;
385:   PetscSF            sf;
386:   DM                 dmFace, dmCell, dmGrad;
387:   const PetscScalar *facegeom, *cellgeom = NULL, *grad;
388:   const PetscInt    *leaves;
389:   PetscScalar       *x, *fx;
390:   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
391:   PetscErrorCode     ierr, ierru = 0;

394:   DMGetPointSF(dm, &sf);
395:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
396:   nleaves = PetscMax(0, nleaves);
397:   DMGetDimension(dm, &dim);
398:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
399:   DMGetDS(dm, &prob);
400:   VecGetDM(faceGeometry, &dmFace);
401:   VecGetArrayRead(faceGeometry, &facegeom);
402:   if (cellGeometry) {
403:     VecGetDM(cellGeometry, &dmCell);
404:     VecGetArrayRead(cellGeometry, &cellgeom);
405:   }
406:   if (Grad) {
407:     PetscFV fv;

409:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
410:     VecGetDM(Grad, &dmGrad);
411:     VecGetArrayRead(Grad, &grad);
412:     PetscFVGetNumComponents(fv, &pdim);
413:     DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);
414:   }
415:   VecGetArray(locX, &x);
416:   for (i = 0; i < numids; ++i) {
417:     IS              faceIS;
418:     const PetscInt *faces;
419:     PetscInt        numFaces, f;

421:     DMLabelGetStratumIS(label, ids[i], &faceIS);
422:     if (!faceIS) continue; /* No points with that id on this process */
423:     ISGetLocalSize(faceIS, &numFaces);
424:     ISGetIndices(faceIS, &faces);
425:     for (f = 0; f < numFaces; ++f) {
426:       const PetscInt         face = faces[f], *cells;
427:       PetscFVFaceGeom        *fg;

429:       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
430:       PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
431:       if (loc >= 0) continue;
432:       DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
433:       DMPlexGetSupport(dm, face, &cells);
434:       if (Grad) {
435:         PetscFVCellGeom       *cg;
436:         PetscScalar           *cx, *cgrad;
437:         PetscScalar           *xG;
438:         PetscReal              dx[3];
439:         PetscInt               d;

441:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
442:         DMPlexPointLocalRead(dm, cells[0], x, &cx);
443:         DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
444:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
445:         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
446:         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
447:         ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
448:         if (ierru) {
449:           ISRestoreIndices(faceIS, &faces);
450:           ISDestroy(&faceIS);
451:           goto cleanup;
452:         }
453:       } else {
454:         PetscScalar       *xI;
455:         PetscScalar       *xG;

457:         DMPlexPointLocalRead(dm, cells[0], x, &xI);
458:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
459:         ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
460:         if (ierru) {
461:           ISRestoreIndices(faceIS, &faces);
462:           ISDestroy(&faceIS);
463:           goto cleanup;
464:         }
465:       }
466:     }
467:     ISRestoreIndices(faceIS, &faces);
468:     ISDestroy(&faceIS);
469:   }
470:   cleanup:
471:   VecRestoreArray(locX, &x);
472:   if (Grad) {
473:     DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);
474:     VecRestoreArrayRead(Grad, &grad);
475:   }
476:   if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
477:   VecRestoreArrayRead(faceGeometry, &facegeom);
478:   CHKERRQ(ierru);
479:   return(0);
480: }

482: PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
483: {
484:   PetscInt       numBd, b;

488:   PetscDSGetNumBoundary(dm->prob, &numBd);
489:   for (b = 0; b < numBd; ++b) {
490:     DMBoundaryConditionType type;
491:     const char             *labelname;
492:     DMLabel                 label;
493:     PetscInt                field, Nc;
494:     const PetscInt         *comps;
495:     PetscObject             obj;
496:     PetscClassId            id;
497:     void                    (*func)(void);
498:     PetscInt                numids;
499:     const PetscInt         *ids;
500:     void                   *ctx;

502:     DMGetBoundary(dm, b, &type, NULL, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);
503:     if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
504:     DMGetLabel(dm, labelname, &label);
505:     DMGetField(dm, field, &obj);
506:     PetscObjectGetClassId(obj, &id);
507:     if (id == PETSCFE_CLASSID) {
508:       switch (type) {
509:         /* for FEM, there is no insertion to be done for non-essential boundary conditions */
510:       case DM_BC_ESSENTIAL:
511:         DMPlexLabelAddCells(dm,label);
512:         DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
513:         DMPlexLabelClearCells(dm,label);
514:         break;
515:       case DM_BC_ESSENTIAL_FIELD:
516:         DMPlexLabelAddCells(dm,label);
517:         DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
518:                                                         (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
519:                                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
520:                                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);
521:         DMPlexLabelClearCells(dm,label);
522:         break;
523:       default: break;
524:       }
525:     } else if (id == PETSCFV_CLASSID) {
526:       if (!faceGeomFVM) continue;
527:       DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids,
528:                                                (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
529:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
530:   }
531:   return(0);
532: }

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

537:   Input Parameters:
538: + dm - The DM
539: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
540: . time - The time
541: . faceGeomFVM - Face geometry data for FV discretizations
542: . cellGeomFVM - Cell geometry data for FV discretizations
543: - gradFVM - Gradient reconstruction data for FV discretizations

545:   Output Parameters:
546: . locX - Solution updated with boundary values

548:   Level: developer

550: .seealso: DMProjectFunctionLabelLocal()
551: @*/
552: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
553: {

562:   PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));
563:   return(0);
564: }

566: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
567: {
568:   Vec              localX;
569:   PetscErrorCode   ierr;

572:   DMGetLocalVector(dm, &localX);
573:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);
574:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
575:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
576:   DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);
577:   DMRestoreLocalVector(dm, &localX);
578:   return(0);
579: }

581: /*@C
582:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

584:   Input Parameters:
585: + dm     - The DM
586: . time   - The time
587: . funcs  - The functions to evaluate for each field component
588: . ctxs   - Optional array of contexts to pass to each function, or NULL.
589: - localX - The coefficient vector u_h, a local vector

591:   Output Parameter:
592: . diff - The diff ||u - u_h||_2

594:   Level: developer

596: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
597: @*/
598: PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff)
599: {
600:   const PetscInt   debug = 0;
601:   PetscSection     section;
602:   PetscQuadrature  quad;
603:   PetscScalar     *funcVal, *interpolant;
604:   PetscReal       *coords, *detJ, *J;
605:   PetscReal        localDiff = 0.0;
606:   const PetscReal *quadWeights;
607:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, cEndInterior, c, field, fieldOffset;
608:   PetscErrorCode   ierr;

611:   DMGetDimension(dm, &dim);
612:   DMGetCoordinateDim(dm, &coordDim);
613:   DMGetDefaultSection(dm, &section);
614:   PetscSectionGetNumFields(section, &numFields);
615:   for (field = 0; field < numFields; ++field) {
616:     PetscObject  obj;
617:     PetscClassId id;
618:     PetscInt     Nc;

620:     DMGetField(dm, field, &obj);
621:     PetscObjectGetClassId(obj, &id);
622:     if (id == PETSCFE_CLASSID) {
623:       PetscFE fe = (PetscFE) obj;

625:       PetscFEGetQuadrature(fe, &quad);
626:       PetscFEGetNumComponents(fe, &Nc);
627:     } else if (id == PETSCFV_CLASSID) {
628:       PetscFV fv = (PetscFV) obj;

630:       PetscFVGetQuadrature(fv, &quad);
631:       PetscFVGetNumComponents(fv, &Nc);
632:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
633:     numComponents += Nc;
634:   }
635:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
636:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
637:   PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
638:   DMPlexGetVTKCellHeight(dm, &cellHeight);
639:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
640:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
641:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
642:   for (c = cStart; c < cEnd; ++c) {
643:     PetscScalar *x = NULL;
644:     PetscReal    elemDiff = 0.0;
645:     PetscInt     qc = 0;

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

650:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
651:       PetscObject  obj;
652:       PetscClassId id;
653:       void * const ctx = ctxs ? ctxs[field] : NULL;
654:       PetscInt     Nb, Nc, q, fc;

656:       DMGetField(dm, field, &obj);
657:       PetscObjectGetClassId(obj, &id);
658:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
659:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
660:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
661:       if (debug) {
662:         char title[1024];
663:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
664:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
665:       }
666:       for (q = 0; q < Nq; ++q) {
667:         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);
668:         (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx);
669:         if (ierr) {
670:           PetscErrorCode ierr2;
671:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
672:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
673:           ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
674: 
675:         }
676:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
677:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
678:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
679:         for (fc = 0; fc < Nc; ++fc) {
680:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
681:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
682:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
683:         }
684:       }
685:       fieldOffset += Nb;
686:       qc += Nc;
687:     }
688:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
689:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
690:     localDiff += elemDiff;
691:   }
692:   PetscFree5(funcVal,interpolant,coords,detJ,J);
693:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
694:   *diff = PetscSqrtReal(*diff);
695:   return(0);
696: }

698: 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)
699: {
700:   const PetscInt   debug = 0;
701:   PetscSection     section;
702:   PetscQuadrature  quad;
703:   Vec              localX;
704:   PetscScalar     *funcVal, *interpolant;
705:   const PetscReal *quadPoints, *quadWeights;
706:   PetscReal       *coords, *realSpaceDer, *J, *invJ, *detJ;
707:   PetscReal        localDiff = 0.0;
708:   PetscInt         dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset;
709:   PetscErrorCode   ierr;

712:   DMGetDimension(dm, &dim);
713:   DMGetCoordinateDim(dm, &coordDim);
714:   DMGetDefaultSection(dm, &section);
715:   PetscSectionGetNumFields(section, &numFields);
716:   DMGetLocalVector(dm, &localX);
717:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
718:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
719:   for (field = 0; field < numFields; ++field) {
720:     PetscFE  fe;
721:     PetscInt Nc;

723:     DMGetField(dm, field, (PetscObject *) &fe);
724:     PetscFEGetQuadrature(fe, &quad);
725:     PetscFEGetNumComponents(fe, &Nc);
726:     numComponents += Nc;
727:   }
728:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
729:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
730:   /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
731:   PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);
732:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
733:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
734:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
735:   for (c = cStart; c < cEnd; ++c) {
736:     PetscScalar *x = NULL;
737:     PetscReal    elemDiff = 0.0;
738:     PetscInt     qc = 0;

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

743:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
744:       PetscFE          fe;
745:       void * const     ctx = ctxs ? ctxs[field] : NULL;
746:       PetscReal       *basisDer;
747:       PetscInt         Nb, Nc, q, fc;

749:       DMGetField(dm, field, (PetscObject *) &fe);
750:       PetscFEGetDimension(fe, &Nb);
751:       PetscFEGetNumComponents(fe, &Nc);
752:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
753:       if (debug) {
754:         char title[1024];
755:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
756:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
757:       }
758:       for (q = 0; q < Nq; ++q) {
759:         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);
760:         (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx);
761:         if (ierr) {
762:           PetscErrorCode ierr2;
763:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
764:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
765:           ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2);
766: 
767:         }
768:         PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);
769:         for (fc = 0; fc < Nc; ++fc) {
770:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
771:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
772:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
773:         }
774:       }
775:       fieldOffset += Nb;
776:       qc          += Nc;
777:     }
778:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
779:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
780:     localDiff += elemDiff;
781:   }
782:   PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);
783:   DMRestoreLocalVector(dm, &localX);
784:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
785:   *diff = PetscSqrtReal(*diff);
786:   return(0);
787: }

789: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
790: {
791:   const PetscInt   debug = 0;
792:   PetscSection     section;
793:   PetscQuadrature  quad;
794:   Vec              localX;
795:   PetscScalar     *funcVal, *interpolant;
796:   PetscReal       *coords, *detJ, *J;
797:   PetscReal       *localDiff;
798:   const PetscReal *quadPoints, *quadWeights;
799:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
800:   PetscErrorCode   ierr;

803:   DMGetDimension(dm, &dim);
804:   DMGetCoordinateDim(dm, &coordDim);
805:   DMGetDefaultSection(dm, &section);
806:   PetscSectionGetNumFields(section, &numFields);
807:   DMGetLocalVector(dm, &localX);
808:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
809:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
810:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
811:   for (field = 0; field < numFields; ++field) {
812:     PetscObject  obj;
813:     PetscClassId id;
814:     PetscInt     Nc;

816:     DMGetField(dm, field, &obj);
817:     PetscObjectGetClassId(obj, &id);
818:     if (id == PETSCFE_CLASSID) {
819:       PetscFE fe = (PetscFE) obj;

821:       PetscFEGetQuadrature(fe, &quad);
822:       PetscFEGetNumComponents(fe, &Nc);
823:     } else if (id == PETSCFV_CLASSID) {
824:       PetscFV fv = (PetscFV) obj;

826:       PetscFVGetQuadrature(fv, &quad);
827:       PetscFVGetNumComponents(fv, &Nc);
828:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
829:     numComponents += Nc;
830:   }
831:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
832:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
833:   PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
834:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
835:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
836:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
837:   for (c = cStart; c < cEnd; ++c) {
838:     PetscScalar *x = NULL;
839:     PetscInt     qc = 0;

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

844:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
845:       PetscObject  obj;
846:       PetscClassId id;
847:       void * const ctx = ctxs ? ctxs[field] : NULL;
848:       PetscInt     Nb, Nc, q, fc;

850:       PetscReal       elemDiff = 0.0;

852:       DMGetField(dm, field, &obj);
853:       PetscObjectGetClassId(obj, &id);
854:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
855:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
856:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
857:       if (debug) {
858:         char title[1024];
859:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
860:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
861:       }
862:       for (q = 0; q < Nq; ++q) {
863:         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);
864:         (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx);
865:         if (ierr) {
866:           PetscErrorCode ierr2;
867:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
868:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
869:           ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
870: 
871:         }
872:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
873:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
874:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
875:         for (fc = 0; fc < Nc; ++fc) {
876:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
877:           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]);}
878:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
879:         }
880:       }
881:       fieldOffset += Nb;
882:       qc          += Nc;
883:       localDiff[field] += elemDiff;
884:     }
885:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
886:   }
887:   DMRestoreLocalVector(dm, &localX);
888:   MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
889:   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
890:   PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);
891:   return(0);
892: }

894: /*@C
895:   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.

897:   Input Parameters:
898: + dm    - The DM
899: . time  - The time
900: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
901: . ctxs  - Optional array of contexts to pass to each function, or NULL.
902: - X     - The coefficient vector u_h

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

907:   Level: developer

909: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
910: @*/
911: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
912: {
913:   PetscSection     section;
914:   PetscQuadrature  quad;
915:   Vec              localX;
916:   PetscScalar     *funcVal, *interpolant;
917:   PetscReal       *coords, *detJ, *J;
918:   const PetscReal *quadPoints, *quadWeights;
919:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
920:   PetscErrorCode   ierr;

923:   VecSet(D, 0.0);
924:   DMGetDimension(dm, &dim);
925:   DMGetCoordinateDim(dm, &coordDim);
926:   DMGetDefaultSection(dm, &section);
927:   PetscSectionGetNumFields(section, &numFields);
928:   DMGetLocalVector(dm, &localX);
929:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
930:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
931:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
932:   for (field = 0; field < numFields; ++field) {
933:     PetscObject  obj;
934:     PetscClassId id;
935:     PetscInt     Nc;

937:     DMGetField(dm, field, &obj);
938:     PetscObjectGetClassId(obj, &id);
939:     if (id == PETSCFE_CLASSID) {
940:       PetscFE fe = (PetscFE) obj;

942:       PetscFEGetQuadrature(fe, &quad);
943:       PetscFEGetNumComponents(fe, &Nc);
944:     } else if (id == PETSCFV_CLASSID) {
945:       PetscFV fv = (PetscFV) obj;

947:       PetscFVGetQuadrature(fv, &quad);
948:       PetscFVGetNumComponents(fv, &Nc);
949:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
950:     numComponents += Nc;
951:   }
952:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
953:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
954:   PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
955:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
956:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
957:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
958:   for (c = cStart; c < cEnd; ++c) {
959:     PetscScalar *x = NULL;
960:     PetscScalar  elemDiff = 0.0;
961:     PetscInt     qc = 0;

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

966:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
967:       PetscObject  obj;
968:       PetscClassId id;
969:       void * const ctx = ctxs ? ctxs[field] : NULL;
970:       PetscInt     Nb, Nc, q, fc;

972:       DMGetField(dm, field, &obj);
973:       PetscObjectGetClassId(obj, &id);
974:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
975:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
976:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
977:       if (funcs[field]) {
978:         for (q = 0; q < Nq; ++q) {
979:           if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)detJ[q], c, q);
980:           (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx);
981:           if (ierr) {
982:             PetscErrorCode ierr2;
983:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
984:             ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
985:             ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
986: 
987:           }
988:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
989:           else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
990:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
991:           for (fc = 0; fc < Nc; ++fc) {
992:             const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
993:             elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
994:           }
995:         }
996:       }
997:       fieldOffset += Nb;
998:       qc          += Nc;
999:     }
1000:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1001:     VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1002:   }
1003:   PetscFree5(funcVal,interpolant,coords,detJ,J);
1004:   DMRestoreLocalVector(dm, &localX);
1005:   VecSqrtAbs(D);
1006:   return(0);
1007: }

1009: /*@C
1010:   DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec.

1012:   Input Parameters:
1013: + dm - The DM
1014: - LocX  - The coefficient vector u_h

1016:   Output Parameter:
1017: . locC - A Vec which holds the Clement interpolant of the gradient

1019:   Notes:
1020:     Add citation to (Clement, 1975) and definition of the interpolant
1021:   \nabla u_h(v_i) = \sum_{T_i \in support(v_i)} |T_i| \nabla u_h(T_i) / \sum_{T_i \in support(v_i)} |T_i| where |T_i| is the cell volume

1023:   Level: developer

1025: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1026: @*/
1027: PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC)
1028: {
1029:   DM_Plex         *mesh  = (DM_Plex *) dm->data;
1030:   PetscInt         debug = mesh->printFEM;
1031:   DM               dmC;
1032:   PetscSection     section;
1033:   PetscQuadrature  quad;
1034:   PetscScalar     *interpolant, *gradsum;
1035:   PetscReal       *coords, *detJ, *J, *invJ;
1036:   const PetscReal *quadPoints, *quadWeights;
1037:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset;
1038:   PetscErrorCode   ierr;

1041:   VecGetDM(locC, &dmC);
1042:   VecSet(locC, 0.0);
1043:   DMGetDimension(dm, &dim);
1044:   DMGetCoordinateDim(dm, &coordDim);
1045:   DMGetDefaultSection(dm, &section);
1046:   PetscSectionGetNumFields(section, &numFields);
1047:   for (field = 0; field < numFields; ++field) {
1048:     PetscObject  obj;
1049:     PetscClassId id;
1050:     PetscInt     Nc;

1052:     DMGetField(dm, field, &obj);
1053:     PetscObjectGetClassId(obj, &id);
1054:     if (id == PETSCFE_CLASSID) {
1055:       PetscFE fe = (PetscFE) obj;

1057:       PetscFEGetQuadrature(fe, &quad);
1058:       PetscFEGetNumComponents(fe, &Nc);
1059:     } else if (id == PETSCFV_CLASSID) {
1060:       PetscFV fv = (PetscFV) obj;

1062:       PetscFVGetQuadrature(fv, &quad);
1063:       PetscFVGetNumComponents(fv, &Nc);
1064:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1065:     numComponents += Nc;
1066:   }
1067:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1068:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1069:   PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ);
1070:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1071:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1072:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1073:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1074:   for (v = vStart; v < vEnd; ++v) {
1075:     PetscScalar volsum = 0.0;
1076:     PetscInt   *star = NULL;
1077:     PetscInt    starSize, st, d, fc;

1079:     PetscMemzero(gradsum, coordDim*numComponents * sizeof(PetscScalar));
1080:     DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1081:     for (st = 0; st < starSize*2; st += 2) {
1082:       const PetscInt cell = star[st];
1083:       PetscScalar   *grad = &gradsum[coordDim*numComponents];
1084:       PetscScalar   *x    = NULL;
1085:       PetscReal      vol  = 0.0;

1087:       if ((cell < cStart) || (cell >= cEnd)) continue;
1088:       DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, J, invJ, detJ);
1089:       DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
1090:       for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1091:         PetscObject  obj;
1092:         PetscClassId id;
1093:         PetscInt     Nb, Nc, q, qc = 0;

1095:         PetscMemzero(grad, coordDim*numComponents * sizeof(PetscScalar));
1096:         DMGetField(dm, field, &obj);
1097:         PetscObjectGetClassId(obj, &id);
1098:         if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1099:         else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1100:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1101:         for (q = 0; q < Nq; ++q) {
1102:           if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)detJ[q], cell, q);
1103:           if (ierr) {
1104:             PetscErrorCode ierr2;
1105:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
1106:             ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
1107:             ierr2 = PetscFree4(interpolant,coords,detJ,J);CHKERRQ(ierr2);
1108: 
1109:           }
1110:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], coordDim, invJ, NULL, q, interpolant);}
1111:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1112:           for (fc = 0; fc < Nc; ++fc) {
1113:             const PetscReal wt = quadWeights[q*qNc+qc+fc];

1115:             for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*detJ[q];
1116:           }
1117:           vol += quadWeights[q*qNc]*detJ[q];
1118:         }
1119:         fieldOffset += Nb;
1120:         qc          += Nc;
1121:       }
1122:       DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
1123:       for (fc = 0; fc < numComponents; ++fc) {
1124:         for (d = 0; d < coordDim; ++d) {
1125:           gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
1126:         }
1127:       }
1128:       volsum += vol;
1129:       if (debug) {
1130:         PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
1131:         for (fc = 0; fc < numComponents; ++fc) {
1132:           for (d = 0; d < coordDim; ++d) {
1133:             if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
1134:             PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
1135:           }
1136:         }
1137:         PetscPrintf(PETSC_COMM_SELF, "]\n");
1138:       }
1139:     }
1140:     for (fc = 0; fc < numComponents; ++fc) {
1141:       for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum;
1142:     }
1143:     DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1144:     DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
1145:   }
1146:   PetscFree6(gradsum,interpolant,coords,detJ,J,invJ);
1147:   return(0);
1148: }

1150: static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user)
1151: {
1152:   DM                 dmAux = NULL;
1153:   PetscDS            prob,    probAux = NULL;
1154:   PetscSection       section, sectionAux;
1155:   Vec                locX,    locA;
1156:   PetscInt           dim, numCells = cEnd - cStart, c, f;
1157:   PetscBool          useFVM = PETSC_FALSE;
1158:   /* DS */
1159:   PetscInt           Nf,    totDim,    *uOff, *uOff_x, numConstants;
1160:   PetscInt           NfAux, totDimAux, *aOff;
1161:   PetscScalar       *u, *a;
1162:   const PetscScalar *constants;
1163:   /* Geometry */
1164:   PetscFEGeom       *cgeomFEM;
1165:   DM                 dmGrad;
1166:   PetscQuadrature    affineQuad = NULL;
1167:   Vec                cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
1168:   PetscFVCellGeom   *cgeomFVM;
1169:   const PetscScalar *lgrad;
1170:   PetscBool          isAffine;
1171:   DMField            coordField;
1172:   IS                 cellIS;
1173:   PetscErrorCode     ierr;

1176:   DMGetDS(dm, &prob);
1177:   DMGetDimension(dm, &dim);
1178:   DMGetDefaultSection(dm, &section);
1179:   PetscSectionGetNumFields(section, &Nf);
1180:   /* Determine which discretizations we have */
1181:   for (f = 0; f < Nf; ++f) {
1182:     PetscObject  obj;
1183:     PetscClassId id;

1185:     PetscDSGetDiscretization(prob, f, &obj);
1186:     PetscObjectGetClassId(obj, &id);
1187:     if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE;
1188:   }
1189:   /* Get local solution with boundary values */
1190:   DMGetLocalVector(dm, &locX);
1191:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1192:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1193:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1194:   /* Read DS information */
1195:   PetscDSGetTotalDimension(prob, &totDim);
1196:   PetscDSGetComponentOffsets(prob, &uOff);
1197:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1198:   ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);
1199:   PetscDSGetConstants(prob, &numConstants, &constants);
1200:   /* Read Auxiliary DS information */
1201:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1202:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1203:   if (dmAux) {
1204:     DMGetDS(dmAux, &probAux);
1205:     PetscDSGetNumFields(probAux, &NfAux);
1206:     DMGetDefaultSection(dmAux, &sectionAux);
1207:     PetscDSGetTotalDimension(probAux, &totDimAux);
1208:     PetscDSGetComponentOffsets(probAux, &aOff);
1209:   }
1210:   /* Allocate data  arrays */
1211:   PetscCalloc1(numCells*totDim, &u);
1212:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1213:   /* Read out geometry */
1214:   DMGetCoordinateField(dm,&coordField);
1215:   DMFieldGetFEInvariance(coordField,cellIS,NULL,&isAffine,NULL);
1216:   if (isAffine) {
1217:     DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
1218:     if (affineQuad) {
1219:       DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);
1220:     }
1221:   }
1222:   if (useFVM) {
1223:     PetscFV   fv = NULL;
1224:     Vec       grad;
1225:     PetscInt  fStart, fEnd;
1226:     PetscBool compGrad;

1228:     for (f = 0; f < Nf; ++f) {
1229:       PetscObject  obj;
1230:       PetscClassId id;

1232:       PetscDSGetDiscretization(prob, f, &obj);
1233:       PetscObjectGetClassId(obj, &id);
1234:       if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;}
1235:     }
1236:     PetscFVGetComputeGradients(fv, &compGrad);
1237:     PetscFVSetComputeGradients(fv, PETSC_TRUE);
1238:     DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1239:     DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1240:     PetscFVSetComputeGradients(fv, compGrad);
1241:     VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1242:     /* Reconstruct and limit cell gradients */
1243:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1244:     DMGetGlobalVector(dmGrad, &grad);
1245:     DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1246:     /* Communicate gradient values */
1247:     DMGetLocalVector(dmGrad, &locGrad);
1248:     DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1249:     DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1250:     DMRestoreGlobalVector(dmGrad, &grad);
1251:     /* Handle non-essential (e.g. outflow) boundary values */
1252:     DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1253:     VecGetArrayRead(locGrad, &lgrad);
1254:   }
1255:   /* Read out data from inputs */
1256:   for (c = cStart; c < cEnd; ++c) {
1257:     PetscScalar *x = NULL;
1258:     PetscInt     i;

1260:     DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1261:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1262:     DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1263:     if (dmAux) {
1264:       DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1265:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1266:       DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1267:     }
1268:   }
1269:   /* Do integration for each field */
1270:   for (f = 0; f < Nf; ++f) {
1271:     PetscObject  obj;
1272:     PetscClassId id;
1273:     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1275:     PetscDSGetDiscretization(prob, f, &obj);
1276:     PetscObjectGetClassId(obj, &id);
1277:     if (id == PETSCFE_CLASSID) {
1278:       PetscFE         fe = (PetscFE) obj;
1279:       PetscQuadrature q;
1280:       PetscFEGeom     *chunkGeom = NULL;
1281:       PetscInt        Nq, Nb;

1283:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1284:       PetscFEGetQuadrature(fe, &q);
1285:       PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1286:       PetscFEGetDimension(fe, &Nb);
1287:       blockSize = Nb*Nq;
1288:       batchSize = numBlocks * blockSize;
1289:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1290:       numChunks = numCells / (numBatches*batchSize);
1291:       Ne        = numChunks*numBatches*batchSize;
1292:       Nr        = numCells % (numBatches*batchSize);
1293:       offset    = numCells - Nr;
1294:       if (!affineQuad) {
1295:         DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);
1296:       }
1297:       PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
1298:       PetscFEIntegrate(fe, prob, f, Ne, chunkGeom, u, probAux, a, cintegral);
1299:       PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);
1300:       PetscFEIntegrate(fe, prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);
1301:       PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);
1302:       if (!affineQuad) {
1303:         PetscFEGeomDestroy(&cgeomFEM);
1304:       }
1305:     } else if (id == PETSCFV_CLASSID) {
1306:       PetscInt       foff;
1307:       PetscPointFunc obj_func;
1308:       PetscScalar    lint;

1310:       PetscDSGetObjective(prob, f, &obj_func);
1311:       PetscDSGetFieldOffset(prob, f, &foff);
1312:       if (obj_func) {
1313:         for (c = 0; c < numCells; ++c) {
1314:           PetscScalar *u_x;

1316:           DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1317:           obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, &a[totDimAux*c], NULL, NULL, 0.0, cgeomFVM[c].centroid, numConstants, constants, &lint);
1318:           cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1319:         }
1320:       }
1321:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1322:   }
1323:   /* Cleanup data arrays */
1324:   if (useFVM) {
1325:     VecRestoreArrayRead(locGrad, &lgrad);
1326:     VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1327:     DMRestoreLocalVector(dmGrad, &locGrad);
1328:     VecDestroy(&faceGeometryFVM);
1329:     VecDestroy(&cellGeometryFVM);
1330:     DMDestroy(&dmGrad);
1331:   }
1332:   if (dmAux) {PetscFree(a);}
1333:   PetscFree(u);
1334:   /* Cleanup */
1335:   if (affineQuad) {
1336:     PetscFEGeomDestroy(&cgeomFEM);
1337:   }
1338:   PetscQuadratureDestroy(&affineQuad);
1339:   ISDestroy(&cellIS);
1340:   DMRestoreLocalVector(dm, &locX);
1341:   return(0);
1342: }

1344: /*@
1345:   DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user

1347:   Input Parameters:
1348: + dm - The mesh
1349: . X  - Global input vector
1350: - user - The user context

1352:   Output Parameter:
1353: . integral - Integral for each field

1355:   Level: developer

1357: .seealso: DMPlexComputeResidualFEM()
1358: @*/
1359: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
1360: {
1361:   DM_Plex       *mesh = (DM_Plex *) dm->data;
1362:   PetscScalar   *cintegral, *lintegral;
1363:   PetscInt       Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;

1370:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1371:   DMGetNumFields(dm, &Nf);
1372:   DMPlexGetVTKCellHeight(dm, &cellHeight);
1373:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1374:   DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1375:   cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1376:   /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1377:   PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);
1378:   DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1379:   /* Sum up values */
1380:   for (cell = cStart; cell < cEnd; ++cell) {
1381:     const PetscInt c = cell - cStart;

1383:     if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1384:     for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f];
1385:   }
1386:   MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1387:   if (mesh->printFEM) {
1388:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");
1389:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));}
1390:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1391:   }
1392:   PetscFree2(lintegral, cintegral);
1393:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1394:   return(0);
1395: }

1397: /*@
1398:   DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user

1400:   Input Parameters:
1401: + dm - The mesh
1402: . X  - Global input vector
1403: - user - The user context

1405:   Output Parameter:
1406: . integral - Cellwise integrals for each field

1408:   Level: developer

1410: .seealso: DMPlexComputeResidualFEM()
1411: @*/
1412: PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user)
1413: {
1414:   DM_Plex       *mesh = (DM_Plex *) dm->data;
1415:   DM             dmF;
1416:   PetscSection   sectionF;
1417:   PetscScalar   *cintegral, *af;
1418:   PetscInt       Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;

1425:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1426:   DMGetNumFields(dm, &Nf);
1427:   DMPlexGetVTKCellHeight(dm, &cellHeight);
1428:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1429:   DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1430:   cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1431:   /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1432:   PetscCalloc1((cEnd-cStart)*Nf, &cintegral);
1433:   DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1434:   /* Put values in F*/
1435:   VecGetDM(F, &dmF);
1436:   DMGetDefaultSection(dmF, &sectionF);
1437:   VecGetArray(F, &af);
1438:   for (cell = cStart; cell < cEnd; ++cell) {
1439:     const PetscInt c = cell - cStart;
1440:     PetscInt       dof, off;

1442:     if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1443:     PetscSectionGetDof(sectionF, cell, &dof);
1444:     PetscSectionGetOffset(sectionF, cell, &off);
1445:     if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf);
1446:     for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f];
1447:   }
1448:   VecRestoreArray(F, &af);
1449:   PetscFree(cintegral);
1450:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1451:   return(0);
1452: }

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

1457:   Input Parameters:
1458: + dmf  - The fine mesh
1459: . dmc  - The coarse mesh
1460: - user - The user context

1462:   Output Parameter:
1463: . In  - The interpolation matrix

1465:   Level: developer

1467: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1468: @*/
1469: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1470: {
1471:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1472:   const char       *name  = "Interpolator";
1473:   PetscDS           prob;
1474:   PetscFE          *feRef;
1475:   PetscFV          *fvRef;
1476:   PetscSection      fsection, fglobalSection;
1477:   PetscSection      csection, cglobalSection;
1478:   PetscScalar      *elemMat;
1479:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1480:   PetscInt          cTotDim, rTotDim = 0;
1481:   PetscErrorCode    ierr;

1484:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1485:   DMGetDimension(dmf, &dim);
1486:   DMGetDefaultSection(dmf, &fsection);
1487:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1488:   DMGetDefaultSection(dmc, &csection);
1489:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1490:   PetscSectionGetNumFields(fsection, &Nf);
1491:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1492:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1493:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1494:   DMGetDS(dmf, &prob);
1495:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1496:   for (f = 0; f < Nf; ++f) {
1497:     PetscObject  obj;
1498:     PetscClassId id;
1499:     PetscInt     rNb = 0, Nc = 0;

1501:     PetscDSGetDiscretization(prob, f, &obj);
1502:     PetscObjectGetClassId(obj, &id);
1503:     if (id == PETSCFE_CLASSID) {
1504:       PetscFE fe = (PetscFE) obj;

1506:       PetscFERefine(fe, &feRef[f]);
1507:       PetscFEGetDimension(feRef[f], &rNb);
1508:       PetscFEGetNumComponents(fe, &Nc);
1509:     } else if (id == PETSCFV_CLASSID) {
1510:       PetscFV        fv = (PetscFV) obj;
1511:       PetscDualSpace Q;

1513:       PetscFVRefine(fv, &fvRef[f]);
1514:       PetscFVGetDualSpace(fvRef[f], &Q);
1515:       PetscDualSpaceGetDimension(Q, &rNb);
1516:       PetscFVGetNumComponents(fv, &Nc);
1517:     }
1518:     rTotDim += rNb;
1519:   }
1520:   PetscDSGetTotalDimension(prob, &cTotDim);
1521:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1522:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1523:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1524:     PetscDualSpace   Qref;
1525:     PetscQuadrature  f;
1526:     const PetscReal *qpoints, *qweights;
1527:     PetscReal       *points;
1528:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1530:     /* Compose points from all dual basis functionals */
1531:     if (feRef[fieldI]) {
1532:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
1533:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
1534:     } else {
1535:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1536:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1537:     }
1538:     PetscDualSpaceGetDimension(Qref, &fpdim);
1539:     for (i = 0; i < fpdim; ++i) {
1540:       PetscDualSpaceGetFunctional(Qref, i, &f);
1541:       PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
1542:       npoints += Np;
1543:     }
1544:     PetscMalloc1(npoints*dim,&points);
1545:     for (i = 0, k = 0; i < fpdim; ++i) {
1546:       PetscDualSpaceGetFunctional(Qref, i, &f);
1547:       PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1548:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1549:     }

1551:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1552:       PetscObject  obj;
1553:       PetscClassId id;
1554:       PetscReal   *B;
1555:       PetscInt     NcJ = 0, cpdim = 0, j, qNc;

1557:       PetscDSGetDiscretization(prob, fieldJ, &obj);
1558:       PetscObjectGetClassId(obj, &id);
1559:       if (id == PETSCFE_CLASSID) {
1560:         PetscFE fe = (PetscFE) obj;

1562:         /* Evaluate basis at points */
1563:         PetscFEGetNumComponents(fe, &NcJ);
1564:         PetscFEGetDimension(fe, &cpdim);
1565:         /* For now, fields only interpolate themselves */
1566:         if (fieldI == fieldJ) {
1567:           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);
1568:           PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1569:           for (i = 0, k = 0; i < fpdim; ++i) {
1570:             PetscDualSpaceGetFunctional(Qref, i, &f);
1571:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1572:             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);
1573:             for (p = 0; p < Np; ++p, ++k) {
1574:               for (j = 0; j < cpdim; ++j) {
1575:                 /*
1576:                    cTotDim:            Total columns in element interpolation matrix, sum of number of dual basis functionals in each field
1577:                    offsetI, offsetJ:   Offsets into the larger element interpolation matrix for different fields
1578:                    fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals
1579:                    qNC, Nc, Ncj, c:    Number of components in this field
1580:                    Np, p:              Number of quad points in the fine grid functional i
1581:                    k:                  i*Np + p, overall point number for the interpolation
1582:                 */
1583:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
1584:               }
1585:             }
1586:           }
1587:           PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1588:         }
1589:       } else if (id == PETSCFV_CLASSID) {
1590:         PetscFV        fv = (PetscFV) obj;

1592:         /* Evaluate constant function at points */
1593:         PetscFVGetNumComponents(fv, &NcJ);
1594:         cpdim = 1;
1595:         /* For now, fields only interpolate themselves */
1596:         if (fieldI == fieldJ) {
1597:           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);
1598:           for (i = 0, k = 0; i < fpdim; ++i) {
1599:             PetscDualSpaceGetFunctional(Qref, i, &f);
1600:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1601:             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);
1602:             for (p = 0; p < Np; ++p, ++k) {
1603:               for (j = 0; j < cpdim; ++j) {
1604:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c];
1605:               }
1606:             }
1607:           }
1608:         }
1609:       }
1610:       offsetJ += cpdim;
1611:     }
1612:     offsetI += fpdim;
1613:     PetscFree(points);
1614:   }
1615:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1616:   /* Preallocate matrix */
1617:   {
1618:     Mat          preallocator;
1619:     PetscScalar *vals;
1620:     PetscInt    *cellCIndices, *cellFIndices;
1621:     PetscInt     locRows, locCols, cell;

1623:     MatGetLocalSize(In, &locRows, &locCols);
1624:     MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1625:     MatSetType(preallocator, MATPREALLOCATOR);
1626:     MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1627:     MatSetUp(preallocator);
1628:     PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1629:     for (cell = cStart; cell < cEnd; ++cell) {
1630:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1631:       MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1632:     }
1633:     PetscFree3(vals,cellCIndices,cellFIndices);
1634:     MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1635:     MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1636:     MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1637:     MatDestroy(&preallocator);
1638:   }
1639:   /* Fill matrix */
1640:   MatZeroEntries(In);
1641:   for (c = cStart; c < cEnd; ++c) {
1642:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1643:   }
1644:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1645:   PetscFree2(feRef,fvRef);
1646:   PetscFree(elemMat);
1647:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1648:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1649:   if (mesh->printFEM) {
1650:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1651:     MatChop(In, 1.0e-10);
1652:     MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1653:   }
1654:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1655:   return(0);
1656: }

1658: PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user)
1659: {
1660:   SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness");
1661: }

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

1666:   Input Parameters:
1667: + dmf  - The fine mesh
1668: . dmc  - The coarse mesh
1669: - user - The user context

1671:   Output Parameter:
1672: . In  - The interpolation matrix

1674:   Level: developer

1676: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1677: @*/
1678: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1679: {
1680:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
1681:   const char    *name = "Interpolator";
1682:   PetscDS        prob;
1683:   PetscSection   fsection, csection, globalFSection, globalCSection;
1684:   PetscHashJK    ht;
1685:   PetscLayout    rLayout;
1686:   PetscInt      *dnz, *onz;
1687:   PetscInt       locRows, rStart, rEnd;
1688:   PetscReal     *x, *v0, *J, *invJ, detJ;
1689:   PetscReal     *v0c, *Jc, *invJc, detJc;
1690:   PetscScalar   *elemMat;
1691:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

1695:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1696:   DMGetCoordinateDim(dmc, &dim);
1697:   DMGetDS(dmc, &prob);
1698:   PetscDSGetRefCoordArrays(prob, &x, NULL);
1699:   PetscDSGetNumFields(prob, &Nf);
1700:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1701:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1702:   DMGetDefaultSection(dmf, &fsection);
1703:   DMGetDefaultGlobalSection(dmf, &globalFSection);
1704:   DMGetDefaultSection(dmc, &csection);
1705:   DMGetDefaultGlobalSection(dmc, &globalCSection);
1706:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1707:   PetscDSGetTotalDimension(prob, &totDim);
1708:   PetscMalloc1(totDim, &elemMat);

1710:   MatGetLocalSize(In, &locRows, NULL);
1711:   PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1712:   PetscLayoutSetLocalSize(rLayout, locRows);
1713:   PetscLayoutSetBlockSize(rLayout, 1);
1714:   PetscLayoutSetUp(rLayout);
1715:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1716:   PetscLayoutDestroy(&rLayout);
1717:   PetscCalloc2(locRows,&dnz,locRows,&onz);
1718:   PetscHashJKCreate(&ht);
1719:   for (field = 0; field < Nf; ++field) {
1720:     PetscObject      obj;
1721:     PetscClassId     id;
1722:     PetscDualSpace   Q = NULL;
1723:     PetscQuadrature  f;
1724:     const PetscReal *qpoints;
1725:     PetscInt         Nc, Np, fpdim, i, d;

1727:     PetscDSGetDiscretization(prob, field, &obj);
1728:     PetscObjectGetClassId(obj, &id);
1729:     if (id == PETSCFE_CLASSID) {
1730:       PetscFE fe = (PetscFE) obj;

1732:       PetscFEGetDualSpace(fe, &Q);
1733:       PetscFEGetNumComponents(fe, &Nc);
1734:     } else if (id == PETSCFV_CLASSID) {
1735:       PetscFV fv = (PetscFV) obj;

1737:       PetscFVGetDualSpace(fv, &Q);
1738:       Nc   = 1;
1739:     }
1740:     PetscDualSpaceGetDimension(Q, &fpdim);
1741:     /* For each fine grid cell */
1742:     for (cell = cStart; cell < cEnd; ++cell) {
1743:       PetscInt *findices,   *cindices;
1744:       PetscInt  numFIndices, numCIndices;

1746:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1747:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1748:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
1749:       for (i = 0; i < fpdim; ++i) {
1750:         Vec             pointVec;
1751:         PetscScalar    *pV;
1752:         PetscSF         coarseCellSF = NULL;
1753:         const PetscSFNode *coarseCells;
1754:         PetscInt        numCoarseCells, q, c;

1756:         /* Get points from the dual basis functional quadrature */
1757:         PetscDualSpaceGetFunctional(Q, i, &f);
1758:         PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1759:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1760:         VecSetBlockSize(pointVec, dim);
1761:         VecGetArray(pointVec, &pV);
1762:         for (q = 0; q < Np; ++q) {
1763:           const PetscReal xi0[3] = {-1., -1., -1.};

1765:           /* Transform point to real space */
1766:           CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
1767:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1768:         }
1769:         VecRestoreArray(pointVec, &pV);
1770:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1771:         /* OPT: Pack all quad points from fine cell */
1772:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1773:         PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
1774:         /* Update preallocation info */
1775:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1776:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1777:         {
1778:           PetscHashJKKey  key;
1779:           PetscHashJKIter missing, iter;

1781:           key.j = findices[i];
1782:           if (key.j >= 0) {
1783:             /* Get indices for coarse elements */
1784:             for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1785:               DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1786:               for (c = 0; c < numCIndices; ++c) {
1787:                 key.k = cindices[c];
1788:                 if (key.k < 0) continue;
1789:                 PetscHashJKPut(ht, key, &missing, &iter);
1790:                 if (missing) {
1791:                   PetscHashJKSet(ht, iter, 1);
1792:                   if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1793:                   else                                     ++onz[key.j-rStart];
1794:                 }
1795:               }
1796:               DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1797:             }
1798:           }
1799:         }
1800:         PetscSFDestroy(&coarseCellSF);
1801:         VecDestroy(&pointVec);
1802:       }
1803:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1804:     }
1805:   }
1806:   PetscHashJKDestroy(&ht);
1807:   MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1808:   MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1809:   PetscFree2(dnz,onz);
1810:   for (field = 0; field < Nf; ++field) {
1811:     PetscObject      obj;
1812:     PetscClassId     id;
1813:     PetscDualSpace   Q = NULL;
1814:     PetscQuadrature  f;
1815:     const PetscReal *qpoints, *qweights;
1816:     PetscInt         Nc, qNc, Np, fpdim, i, d;

1818:     PetscDSGetDiscretization(prob, field, &obj);
1819:     PetscObjectGetClassId(obj, &id);
1820:     if (id == PETSCFE_CLASSID) {
1821:       PetscFE fe = (PetscFE) obj;

1823:       PetscFEGetDualSpace(fe, &Q);
1824:       PetscFEGetNumComponents(fe, &Nc);
1825:     } else if (id == PETSCFV_CLASSID) {
1826:       PetscFV fv = (PetscFV) obj;

1828:       PetscFVGetDualSpace(fv, &Q);
1829:       Nc   = 1;
1830:     } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %d",field);
1831:     PetscDualSpaceGetDimension(Q, &fpdim);
1832:     /* For each fine grid cell */
1833:     for (cell = cStart; cell < cEnd; ++cell) {
1834:       PetscInt *findices,   *cindices;
1835:       PetscInt  numFIndices, numCIndices;

1837:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1838:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1839:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
1840:       for (i = 0; i < fpdim; ++i) {
1841:         Vec             pointVec;
1842:         PetscScalar    *pV;
1843:         PetscSF         coarseCellSF = NULL;
1844:         const PetscSFNode *coarseCells;
1845:         PetscInt        numCoarseCells, cpdim, q, c, j;

1847:         /* Get points from the dual basis functional quadrature */
1848:         PetscDualSpaceGetFunctional(Q, i, &f);
1849:         PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
1850:         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);
1851:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1852:         VecSetBlockSize(pointVec, dim);
1853:         VecGetArray(pointVec, &pV);
1854:         for (q = 0; q < Np; ++q) {
1855:           const PetscReal xi0[3] = {-1., -1., -1.};

1857:           /* Transform point to real space */
1858:           CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
1859:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1860:         }
1861:         VecRestoreArray(pointVec, &pV);
1862:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1863:         /* OPT: Read this out from preallocation information */
1864:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1865:         /* Update preallocation info */
1866:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1867:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1868:         VecGetArray(pointVec, &pV);
1869:         for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1870:           PetscReal pVReal[3];
1871:           const PetscReal xi0[3] = {-1., -1., -1.};

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

1879:           if (id == PETSCFE_CLASSID) {
1880:             PetscFE    fe = (PetscFE) obj;
1881:             PetscReal *B;

1883:             /* Evaluate coarse basis on contained point */
1884:             PetscFEGetDimension(fe, &cpdim);
1885:             PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
1886:             PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
1887:             /* Get elemMat entries by multiplying by weight */
1888:             for (j = 0; j < cpdim; ++j) {
1889:               for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c];
1890:             }
1891:             PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
1892:           } else {
1893:             cpdim = 1;
1894:             for (j = 0; j < cpdim; ++j) {
1895:               for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
1896:             }
1897:           }
1898:           /* Update interpolator */
1899:           if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
1900:           if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
1901:           MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
1902:           DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1903:         }
1904:         VecRestoreArray(pointVec, &pV);
1905:         PetscSFDestroy(&coarseCellSF);
1906:         VecDestroy(&pointVec);
1907:       }
1908:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1909:     }
1910:   }
1911:   PetscFree3(v0,J,invJ);
1912:   PetscFree3(v0c,Jc,invJc);
1913:   PetscFree(elemMat);
1914:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1915:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1916:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1917:   return(0);
1918: }

1920: /*@
1921:   DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM.

1923:   Input Parameters:
1924: + dmf  - The fine mesh
1925: . dmc  - The coarse mesh
1926: - user - The user context

1928:   Output Parameter:
1929: . mass  - The mass matrix

1931:   Level: developer

1933: .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1934: @*/
1935: PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user)
1936: {
1937:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
1938:   const char    *name = "Mass Matrix";
1939:   PetscDS        prob;
1940:   PetscSection   fsection, csection, globalFSection, globalCSection;
1941:   PetscHashJK    ht;
1942:   PetscLayout    rLayout;
1943:   PetscInt      *dnz, *onz;
1944:   PetscInt       locRows, rStart, rEnd;
1945:   PetscReal     *x, *v0, *J, *invJ, detJ;
1946:   PetscReal     *v0c, *Jc, *invJc, detJc;
1947:   PetscScalar   *elemMat;
1948:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

1952:   DMGetCoordinateDim(dmc, &dim);
1953:   DMGetDS(dmc, &prob);
1954:   PetscDSGetRefCoordArrays(prob, &x, NULL);
1955:   PetscDSGetNumFields(prob, &Nf);
1956:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1957:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1958:   DMGetDefaultSection(dmf, &fsection);
1959:   DMGetDefaultGlobalSection(dmf, &globalFSection);
1960:   DMGetDefaultSection(dmc, &csection);
1961:   DMGetDefaultGlobalSection(dmc, &globalCSection);
1962:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1963:   PetscDSGetTotalDimension(prob, &totDim);
1964:   PetscMalloc1(totDim, &elemMat);

1966:   MatGetLocalSize(mass, &locRows, NULL);
1967:   PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);
1968:   PetscLayoutSetLocalSize(rLayout, locRows);
1969:   PetscLayoutSetBlockSize(rLayout, 1);
1970:   PetscLayoutSetUp(rLayout);
1971:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1972:   PetscLayoutDestroy(&rLayout);
1973:   PetscCalloc2(locRows,&dnz,locRows,&onz);
1974:   PetscHashJKCreate(&ht);
1975:   for (field = 0; field < Nf; ++field) {
1976:     PetscObject      obj;
1977:     PetscClassId     id;
1978:     PetscQuadrature  quad;
1979:     const PetscReal *qpoints;
1980:     PetscInt         Nq, Nc, i, d;

1982:     PetscDSGetDiscretization(prob, field, &obj);
1983:     PetscObjectGetClassId(obj, &id);
1984:     if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);}
1985:     else                       {PetscFVGetQuadrature((PetscFV) obj, &quad);}
1986:     PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);
1987:     /* For each fine grid cell */
1988:     for (cell = cStart; cell < cEnd; ++cell) {
1989:       Vec                pointVec;
1990:       PetscScalar       *pV;
1991:       PetscSF            coarseCellSF = NULL;
1992:       const PetscSFNode *coarseCells;
1993:       PetscInt           numCoarseCells, q, c;
1994:       PetscInt          *findices,   *cindices;
1995:       PetscInt           numFIndices, numCIndices;

1997:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1998:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1999:       /* Get points from the quadrature */
2000:       VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2001:       VecSetBlockSize(pointVec, dim);
2002:       VecGetArray(pointVec, &pV);
2003:       for (q = 0; q < Nq; ++q) {
2004:         const PetscReal xi0[3] = {-1., -1., -1.};

2006:         /* Transform point to real space */
2007:         CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2008:         for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2009:       }
2010:       VecRestoreArray(pointVec, &pV);
2011:       /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2012:       DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2013:       PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2014:       /* Update preallocation info */
2015:       PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2016:       if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2017:       {
2018:         PetscHashJKKey  key;
2019:         PetscHashJKIter missing, iter;

2021:         for (i = 0; i < numFIndices; ++i) {
2022:           key.j = findices[i];
2023:           if (key.j >= 0) {
2024:             /* Get indices for coarse elements */
2025:             for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2026:               DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2027:               for (c = 0; c < numCIndices; ++c) {
2028:                 key.k = cindices[c];
2029:                 if (key.k < 0) continue;
2030:                 PetscHashJKPut(ht, key, &missing, &iter);
2031:                 if (missing) {
2032:                   PetscHashJKSet(ht, iter, 1);
2033:                   if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
2034:                   else                                     ++onz[key.j-rStart];
2035:                 }
2036:               }
2037:               DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2038:             }
2039:           }
2040:         }
2041:       }
2042:       PetscSFDestroy(&coarseCellSF);
2043:       VecDestroy(&pointVec);
2044:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2045:     }
2046:   }
2047:   PetscHashJKDestroy(&ht);
2048:   MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
2049:   MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2050:   PetscFree2(dnz,onz);
2051:   for (field = 0; field < Nf; ++field) {
2052:     PetscObject      obj;
2053:     PetscClassId     id;
2054:     PetscQuadrature  quad;
2055:     PetscReal       *Bfine;
2056:     const PetscReal *qpoints, *qweights;
2057:     PetscInt         Nq, Nc, i, d;

2059:     PetscDSGetDiscretization(prob, field, &obj);
2060:     PetscObjectGetClassId(obj, &id);
2061:     if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);}
2062:     else                       {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2063:     PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);
2064:     /* For each fine grid cell */
2065:     for (cell = cStart; cell < cEnd; ++cell) {
2066:       Vec                pointVec;
2067:       PetscScalar       *pV;
2068:       PetscSF            coarseCellSF = NULL;
2069:       const PetscSFNode *coarseCells;
2070:       PetscInt           numCoarseCells, cpdim, q, c, j;
2071:       PetscInt          *findices,   *cindices;
2072:       PetscInt           numFIndices, numCIndices;

2074:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2075:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2076:       /* Get points from the quadrature */
2077:       VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2078:       VecSetBlockSize(pointVec, dim);
2079:       VecGetArray(pointVec, &pV);
2080:       for (q = 0; q < Nq; ++q) {
2081:         const PetscReal xi0[3] = {-1., -1., -1.};

2083:         /* Transform point to real space */
2084:         CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2085:         for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2086:       }
2087:       VecRestoreArray(pointVec, &pV);
2088:       /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2089:       DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2090:       /* Update matrix */
2091:       PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2092:       if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2093:       VecGetArray(pointVec, &pV);
2094:       for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2095:         PetscReal pVReal[3];
2096:         const PetscReal xi0[3] = {-1., -1., -1.};


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

2105:         if (id == PETSCFE_CLASSID) {
2106:           PetscFE    fe = (PetscFE) obj;
2107:           PetscReal *B;

2109:           /* Evaluate coarse basis on contained point */
2110:           PetscFEGetDimension(fe, &cpdim);
2111:           PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2112:           /* Get elemMat entries by multiplying by weight */
2113:           for (i = 0; i < numFIndices; ++i) {
2114:             PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2115:             for (j = 0; j < cpdim; ++j) {
2116:               for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ;
2117:             }
2118:             /* Update interpolator */
2119:             if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2120:             if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2121:             MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2122:           }
2123:           PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2124:         } else {
2125:           cpdim = 1;
2126:           for (i = 0; i < numFIndices; ++i) {
2127:             PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2128:             for (j = 0; j < cpdim; ++j) {
2129:               for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ;
2130:             }
2131:             /* Update interpolator */
2132:             if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2133:             PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);
2134:             if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2135:             MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2136:           }
2137:         }
2138:         DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2139:       }
2140:       VecRestoreArray(pointVec, &pV);
2141:       PetscSFDestroy(&coarseCellSF);
2142:       VecDestroy(&pointVec);
2143:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2144:     }
2145:   }
2146:   PetscFree3(v0,J,invJ);
2147:   PetscFree3(v0c,Jc,invJc);
2148:   PetscFree(elemMat);
2149:   MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
2150:   MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
2151:   return(0);
2152: }

2154: /*@
2155:   DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns

2157:   Input Parameters:
2158: + dmc  - The coarse mesh
2159: - dmf  - The fine mesh
2160: - user - The user context

2162:   Output Parameter:
2163: . sc   - The mapping

2165:   Level: developer

2167: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2168: @*/
2169: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
2170: {
2171:   PetscDS        prob;
2172:   PetscFE       *feRef;
2173:   PetscFV       *fvRef;
2174:   Vec            fv, cv;
2175:   IS             fis, cis;
2176:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
2177:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
2178:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;

2182:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2183:   DMGetDimension(dmf, &dim);
2184:   DMGetDefaultSection(dmf, &fsection);
2185:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
2186:   DMGetDefaultSection(dmc, &csection);
2187:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
2188:   PetscSectionGetNumFields(fsection, &Nf);
2189:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
2190:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
2191:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2192:   DMGetDS(dmc, &prob);
2193:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
2194:   for (f = 0; f < Nf; ++f) {
2195:     PetscObject  obj;
2196:     PetscClassId id;
2197:     PetscInt     fNb = 0, Nc = 0;

2199:     PetscDSGetDiscretization(prob, f, &obj);
2200:     PetscObjectGetClassId(obj, &id);
2201:     if (id == PETSCFE_CLASSID) {
2202:       PetscFE fe = (PetscFE) obj;

2204:       PetscFERefine(fe, &feRef[f]);
2205:       PetscFEGetDimension(feRef[f], &fNb);
2206:       PetscFEGetNumComponents(fe, &Nc);
2207:     } else if (id == PETSCFV_CLASSID) {
2208:       PetscFV        fv = (PetscFV) obj;
2209:       PetscDualSpace Q;

2211:       PetscFVRefine(fv, &fvRef[f]);
2212:       PetscFVGetDualSpace(fvRef[f], &Q);
2213:       PetscDualSpaceGetDimension(Q, &fNb);
2214:       PetscFVGetNumComponents(fv, &Nc);
2215:     }
2216:     fTotDim += fNb;
2217:   }
2218:   PetscDSGetTotalDimension(prob, &cTotDim);
2219:   PetscMalloc1(cTotDim,&cmap);
2220:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
2221:     PetscFE        feC;
2222:     PetscFV        fvC;
2223:     PetscDualSpace QF, QC;
2224:     PetscInt       order = -1, NcF, NcC, fpdim, cpdim;

2226:     if (feRef[field]) {
2227:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
2228:       PetscFEGetNumComponents(feC, &NcC);
2229:       PetscFEGetNumComponents(feRef[field], &NcF);
2230:       PetscFEGetDualSpace(feRef[field], &QF);
2231:       PetscDualSpaceGetOrder(QF, &order);
2232:       PetscDualSpaceGetDimension(QF, &fpdim);
2233:       PetscFEGetDualSpace(feC, &QC);
2234:       PetscDualSpaceGetDimension(QC, &cpdim);
2235:     } else {
2236:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
2237:       PetscFVGetNumComponents(fvC, &NcC);
2238:       PetscFVGetNumComponents(fvRef[field], &NcF);
2239:       PetscFVGetDualSpace(fvRef[field], &QF);
2240:       PetscDualSpaceGetDimension(QF, &fpdim);
2241:       PetscFVGetDualSpace(fvC, &QC);
2242:       PetscDualSpaceGetDimension(QC, &cpdim);
2243:     }
2244:     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);
2245:     for (c = 0; c < cpdim; ++c) {
2246:       PetscQuadrature  cfunc;
2247:       const PetscReal *cqpoints, *cqweights;
2248:       PetscInt         NqcC, NpC;
2249:       PetscBool        found = PETSC_FALSE;

2251:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
2252:       PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);
2253:       if (NqcC != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcC, NcC);
2254:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
2255:       for (f = 0; f < fpdim; ++f) {
2256:         PetscQuadrature  ffunc;
2257:         const PetscReal *fqpoints, *fqweights;
2258:         PetscReal        sum = 0.0;
2259:         PetscInt         NqcF, NpF;

2261:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
2262:         PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);
2263:         if (NqcF != NcF) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcF, NcF);
2264:         if (NpC != NpF) continue;
2265:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
2266:         if (sum > 1.0e-9) continue;
2267:         for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]);
2268:         if (sum < 1.0e-9) continue;
2269:         cmap[offsetC+c] = offsetF+f;
2270:         found = PETSC_TRUE;
2271:         break;
2272:       }
2273:       if (!found) {
2274:         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
2275:         if (fvRef[field] || (feRef[field] && order == 0)) {
2276:           cmap[offsetC+c] = offsetF+0;
2277:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
2278:       }
2279:     }
2280:     offsetC += cpdim;
2281:     offsetF += fpdim;
2282:   }
2283:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
2284:   PetscFree2(feRef,fvRef);

2286:   DMGetGlobalVector(dmf, &fv);
2287:   DMGetGlobalVector(dmc, &cv);
2288:   VecGetOwnershipRange(cv, &startC, &endC);
2289:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
2290:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
2291:   PetscMalloc1(m,&cindices);
2292:   PetscMalloc1(m,&findices);
2293:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
2294:   for (c = cStart; c < cEnd; ++c) {
2295:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
2296:     for (d = 0; d < cTotDim; ++d) {
2297:       if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
2298:       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]]);
2299:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
2300:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
2301:     }
2302:   }
2303:   PetscFree(cmap);
2304:   PetscFree2(cellCIndices,cellFIndices);

2306:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
2307:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
2308:   VecScatterCreate(cv, cis, fv, fis, sc);
2309:   ISDestroy(&cis);
2310:   ISDestroy(&fis);
2311:   DMRestoreGlobalVector(dmf, &fv);
2312:   DMRestoreGlobalVector(dmc, &cv);
2313:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2314:   return(0);
2315: }