Actual source code: plexfem.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petscsf.h>

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

  8: static PetscErrorCode DMPlexConvertPlex(DM dm, DM *plex, PetscBool copy)
  9: {
 10:   PetscBool      isPlex;

 14:   PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
 15:   if (isPlex) {
 16:     *plex = dm;
 17:     PetscObjectReference((PetscObject) dm);
 18:   } else {
 19:     PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
 20:     if (!*plex) {
 21:       DMConvert(dm,DMPLEX,plex);
 22:       PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
 23:       if (copy) {
 24:         const char *comps[3] = {"A", "dmAux"};
 25:         PetscObject obj;
 26:         PetscInt    i;

 28:         {
 29:           /* Run the subdomain hook (this will copy the DMSNES/DMTS) */
 30:           DMSubDomainHookLink link;
 31:           for (link = dm->subdomainhook; link; link = link->next) {
 32:             if (link->ddhook) {(*link->ddhook)(dm, *plex, link->ctx);}
 33:           }
 34:         }
 35:         for (i = 0; i < 3; i++) {
 36:           PetscObjectQuery((PetscObject) dm, comps[i], &obj);
 37:           PetscObjectCompose((PetscObject) *plex, comps[i], obj);
 38:         }
 39:       }
 40:     } else {
 41:       PetscObjectReference((PetscObject) *plex);
 42:     }
 43:   }
 44:   return(0);
 45: }

 47: static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom (void *ctx)
 48: {
 49:   PetscFEGeom *geom = (PetscFEGeom *) ctx;

 53:   PetscFEGeomDestroy(&geom);
 54:   return(0);
 55: }

 57: static PetscErrorCode DMPlexGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
 58: {
 59:   char            composeStr[33] = {0};
 60:   PetscObjectId   id;
 61:   PetscContainer  container;
 62:   PetscErrorCode  ierr;

 65:   PetscObjectGetId((PetscObject)quad,&id);
 66:   PetscSNPrintf(composeStr, 32, "DMPlexGetFEGeom_%x\n", id);
 67:   PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
 68:   if (container) {
 69:     PetscContainerGetPointer(container, (void **) geom);
 70:   } else {
 71:     DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
 72:     PetscContainerCreate(PETSC_COMM_SELF,&container);
 73:     PetscContainerSetPointer(container, (void *) *geom);
 74:     PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
 75:     PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
 76:     PetscContainerDestroy(&container);
 77:   }
 78:   return(0);
 79: }

 81: static PetscErrorCode DMPlexRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
 82: {
 84:   *geom = NULL;
 85:   return(0);
 86: }

 88: /*@
 89:   DMPlexGetScale - Get the scale for the specified fundamental unit

 91:   Not collective

 93:   Input Arguments:
 94: + dm   - the DM
 95: - unit - The SI unit

 97:   Output Argument:
 98: . scale - The value used to scale all quantities with this unit

100:   Level: advanced

102: .seealso: DMPlexSetScale(), PetscUnit
103: @*/
104: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
105: {
106:   DM_Plex *mesh = (DM_Plex*) dm->data;

111:   *scale = mesh->scale[unit];
112:   return(0);
113: }

115: /*@
116:   DMPlexSetScale - Set the scale for the specified fundamental unit

118:   Not collective

120:   Input Arguments:
121: + dm   - the DM
122: . unit - The SI unit
123: - scale - The value used to scale all quantities with this unit

125:   Level: advanced

127: .seealso: DMPlexGetScale(), PetscUnit
128: @*/
129: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
130: {
131:   DM_Plex *mesh = (DM_Plex*) dm->data;

135:   mesh->scale[unit] = scale;
136:   return(0);
137: }

139: static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx)
140: {
141:   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}}};
142:   PetscInt *ctxInt  = (PetscInt *) ctx;
143:   PetscInt  dim2    = ctxInt[0];
144:   PetscInt  d       = ctxInt[1];
145:   PetscInt  i, j, k = dim > 2 ? d - dim : d;

148:   if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2);
149:   for (i = 0; i < dim; i++) mode[i] = 0.;
150:   if (d < dim) {
151:     mode[d] = 1.; /* Translation along axis d */
152:   } else {
153:     for (i = 0; i < dim; i++) {
154:       for (j = 0; j < dim; j++) {
155:         mode[j] += eps[i][j][k]*X[i]; /* Rotation about axis d */
156:       }
157:     }
158:   }
159:   return(0);
160: }

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

165:   Collective on DM

167:   Input Arguments:
168: . dm - the DM

170:   Output Argument:
171: . sp - the null space

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

175:   Level: advanced

177: .seealso: MatNullSpaceCreate(), PCGAMG
178: @*/
179: PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
180: {
181:   MPI_Comm       comm;
182:   Vec            mode[6];
183:   PetscSection   section, globalSection;
184:   PetscInt       dim, dimEmbed, n, m, mmin, d, i, j;

188:   PetscObjectGetComm((PetscObject)dm,&comm);
189:   DMGetDimension(dm, &dim);
190:   DMGetCoordinateDim(dm, &dimEmbed);
191:   if (dim == 1) {
192:     MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
193:     return(0);
194:   }
195:   DMGetSection(dm, &section);
196:   DMGetGlobalSection(dm, &globalSection);
197:   PetscSectionGetConstrainedStorageSize(globalSection, &n);
198:   m    = (dim*(dim+1))/2;
199:   VecCreate(comm, &mode[0]);
200:   VecSetSizes(mode[0], n, PETSC_DETERMINE);
201:   VecSetUp(mode[0]);
202:   VecGetSize(mode[0], &n);
203:   mmin = PetscMin(m, n);
204:   for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
205:   for (d = 0; d < m; d++) {
206:     PetscInt         ctx[2];
207:     PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
208:     void            *voidctx = (void *) (&ctx[0]);

210:     ctx[0] = dimEmbed;
211:     ctx[1] = d;
212:     DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);
213:   }
214:   for (i = 0; i < PetscMin(dim, mmin); ++i) {VecNormalize(mode[i], NULL);}
215:   /* Orthonormalize system */
216:   for (i = dim; i < mmin; ++i) {
217:     PetscScalar dots[6];

219:     VecMDot(mode[i], i, mode, dots);
220:     for (j = 0; j < i; ++j) dots[j] *= -1.0;
221:     VecMAXPY(mode[i], i, dots, mode);
222:     VecNormalize(mode[i], NULL);
223:   }
224:   MatNullSpaceCreate(comm, PETSC_FALSE, mmin, mode, sp);
225:   for (i = 0; i < m; ++i) {VecDestroy(&mode[i]);}
226:   return(0);
227: }

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

232:   Collective on DM

234:   Input Arguments:
235: + dm    - the DM
236: . nb    - The number of bodies
237: . label - The DMLabel marking each domain
238: . nids  - The number of ids per body
239: - ids   - An array of the label ids in sequence for each domain

241:   Output Argument:
242: . sp - the null space

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

246:   Level: advanced

248: .seealso: MatNullSpaceCreate()
249: @*/
250: PetscErrorCode DMPlexCreateRigidBodies(DM dm, PetscInt nb, DMLabel label, const PetscInt nids[], const PetscInt ids[], MatNullSpace *sp)
251: {
252:   MPI_Comm       comm;
253:   PetscSection   section, globalSection;
254:   Vec           *mode;
255:   PetscScalar   *dots;
256:   PetscInt       dim, dimEmbed, n, m, b, d, i, j, off;

260:   PetscObjectGetComm((PetscObject)dm,&comm);
261:   DMGetDimension(dm, &dim);
262:   DMGetCoordinateDim(dm, &dimEmbed);
263:   DMGetSection(dm, &section);
264:   DMGetGlobalSection(dm, &globalSection);
265:   PetscSectionGetConstrainedStorageSize(globalSection, &n);
266:   m    = nb * (dim*(dim+1))/2;
267:   PetscMalloc2(m, &mode, m, &dots);
268:   VecCreate(comm, &mode[0]);
269:   VecSetSizes(mode[0], n, PETSC_DETERMINE);
270:   VecSetUp(mode[0]);
271:   for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
272:   for (b = 0, off = 0; b < nb; ++b) {
273:     for (d = 0; d < m/nb; ++d) {
274:       PetscInt         ctx[2];
275:       PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
276:       void            *voidctx = (void *) (&ctx[0]);

278:       ctx[0] = dimEmbed;
279:       ctx[1] = d;
280:       DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);
281:       off   += nids[b];
282:     }
283:   }
284:   for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
285:   /* Orthonormalize system */
286:   for (i = 0; i < m; ++i) {
287:     VecMDot(mode[i], i, mode, dots);
288:     for (j = 0; j < i; ++j) dots[j] *= -1.0;
289:     VecMAXPY(mode[i], i, dots, mode);
290:     VecNormalize(mode[i], NULL);
291:   }
292:   MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
293:   for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
294:   PetscFree2(mode, dots);
295:   return(0);
296: }

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

307:   Input Parameters:
308: + dm - the DMPlex object
309: - height - the maximum projection height >= 0

311:   Level: advanced

313: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
314: @*/
315: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
316: {
317:   DM_Plex *plex = (DM_Plex *) dm->data;

321:   plex->maxProjectionHeight = height;
322:   return(0);
323: }

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

329:   Input Parameters:
330: . dm - the DMPlex object

332:   Output Parameters:
333: . height - the maximum projection height

335:   Level: intermediate

337: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
338: @*/
339: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
340: {
341:   DM_Plex *plex = (DM_Plex *) dm->data;

345:   *height = plex->maxProjectionHeight;
346:   return(0);
347: }

349: /*@C
350:   DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector

352:   Input Parameters:
353: + dm     - The DM, with a PetscDS that matches the problem being constrained
354: . time   - The time
355: . field  - The field to constrain
356: . Nc     - The number of constrained field components, or 0 for all components
357: . comps  - An array of constrained component numbers, or NULL for all components
358: . label  - The DMLabel defining constrained points
359: . numids - The number of DMLabel ids for constrained points
360: . ids    - An array of ids for constrained points
361: . func   - A pointwise function giving boundary values
362: - ctx    - An optional user context for bcFunc

364:   Output Parameter:
365: . locX   - A local vector to receives the boundary values

367:   Level: developer

369: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
370: @*/
371: 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)
372: {
373:   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
374:   void            **ctxs;
375:   PetscInt          numFields;
376:   PetscErrorCode    ierr;

379:   DMGetNumFields(dm, &numFields);
380:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
381:   funcs[field] = func;
382:   ctxs[field]  = ctx;
383:   DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);
384:   PetscFree2(funcs,ctxs);
385:   return(0);
386: }

388: /*@C
389:   DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector

391:   Input Parameters:
392: + dm     - The DM, with a PetscDS that matches the problem being constrained
393: . time   - The time
394: . locU   - A local vector with the input solution values
395: . field  - The field to constrain
396: . Nc     - The number of constrained field components, or 0 for all components
397: . comps  - An array of constrained component numbers, or NULL for all components
398: . label  - The DMLabel defining constrained points
399: . numids - The number of DMLabel ids for constrained points
400: . ids    - An array of ids for constrained points
401: . func   - A pointwise function giving boundary values
402: - ctx    - An optional user context for bcFunc

404:   Output Parameter:
405: . locX   - A local vector to receives the boundary values

407:   Level: developer

409: .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary()
410: @*/
411: PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
412:                                                         void (*func)(PetscInt, PetscInt, PetscInt,
413:                                                                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
414:                                                                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
415:                                                                      PetscReal, const PetscReal[], PetscInt, const PetscScalar[],
416:                                                                      PetscScalar[]),
417:                                                         void *ctx, Vec locX)
418: {
419:   void (**funcs)(PetscInt, PetscInt, PetscInt,
420:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
421:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
422:                  PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
423:   void            **ctxs;
424:   PetscInt          numFields;
425:   PetscErrorCode    ierr;

428:   DMGetNumFields(dm, &numFields);
429:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
430:   funcs[field] = func;
431:   ctxs[field]  = ctx;
432:   DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
433:   PetscFree2(funcs,ctxs);
434:   return(0);
435: }

437: /*@C
438:   DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector

440:   Input Parameters:
441: + dm     - The DM, with a PetscDS that matches the problem being constrained
442: . time   - The time
443: . faceGeometry - A vector with the FVM face geometry information
444: . cellGeometry - A vector with the FVM cell geometry information
445: . Grad         - A vector with the FVM cell gradient information
446: . field  - The field to constrain
447: . Nc     - The number of constrained field components, or 0 for all components
448: . comps  - An array of constrained component numbers, or NULL for all components
449: . label  - The DMLabel defining constrained points
450: . numids - The number of DMLabel ids for constrained points
451: . ids    - An array of ids for constrained points
452: . func   - A pointwise function giving boundary values
453: - ctx    - An optional user context for bcFunc

455:   Output Parameter:
456: . locX   - A local vector to receives the boundary values

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

460:   Level: developer

462: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
463: @*/
464: 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[],
465:                                                  PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
466: {
467:   PetscDS            prob;
468:   PetscSF            sf;
469:   DM                 dmFace, dmCell, dmGrad;
470:   const PetscScalar *facegeom, *cellgeom = NULL, *grad;
471:   const PetscInt    *leaves;
472:   PetscScalar       *x, *fx;
473:   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
474:   PetscErrorCode     ierr, ierru = 0;

477:   DMGetPointSF(dm, &sf);
478:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
479:   nleaves = PetscMax(0, nleaves);
480:   DMGetDimension(dm, &dim);
481:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
482:   DMGetDS(dm, &prob);
483:   VecGetDM(faceGeometry, &dmFace);
484:   VecGetArrayRead(faceGeometry, &facegeom);
485:   if (cellGeometry) {
486:     VecGetDM(cellGeometry, &dmCell);
487:     VecGetArrayRead(cellGeometry, &cellgeom);
488:   }
489:   if (Grad) {
490:     PetscFV fv;

492:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
493:     VecGetDM(Grad, &dmGrad);
494:     VecGetArrayRead(Grad, &grad);
495:     PetscFVGetNumComponents(fv, &pdim);
496:     DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);
497:   }
498:   VecGetArray(locX, &x);
499:   for (i = 0; i < numids; ++i) {
500:     IS              faceIS;
501:     const PetscInt *faces;
502:     PetscInt        numFaces, f;

504:     DMLabelGetStratumIS(label, ids[i], &faceIS);
505:     if (!faceIS) continue; /* No points with that id on this process */
506:     ISGetLocalSize(faceIS, &numFaces);
507:     ISGetIndices(faceIS, &faces);
508:     for (f = 0; f < numFaces; ++f) {
509:       const PetscInt         face = faces[f], *cells;
510:       PetscFVFaceGeom        *fg;

512:       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
513:       PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
514:       if (loc >= 0) continue;
515:       DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
516:       DMPlexGetSupport(dm, face, &cells);
517:       if (Grad) {
518:         PetscFVCellGeom       *cg;
519:         PetscScalar           *cx, *cgrad;
520:         PetscScalar           *xG;
521:         PetscReal              dx[3];
522:         PetscInt               d;

524:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
525:         DMPlexPointLocalRead(dm, cells[0], x, &cx);
526:         DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
527:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
528:         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
529:         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
530:         ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
531:         if (ierru) {
532:           ISRestoreIndices(faceIS, &faces);
533:           ISDestroy(&faceIS);
534:           goto cleanup;
535:         }
536:       } else {
537:         PetscScalar       *xI;
538:         PetscScalar       *xG;

540:         DMPlexPointLocalRead(dm, cells[0], x, &xI);
541:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
542:         ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
543:         if (ierru) {
544:           ISRestoreIndices(faceIS, &faces);
545:           ISDestroy(&faceIS);
546:           goto cleanup;
547:         }
548:       }
549:     }
550:     ISRestoreIndices(faceIS, &faces);
551:     ISDestroy(&faceIS);
552:   }
553:   cleanup:
554:   VecRestoreArray(locX, &x);
555:   if (Grad) {
556:     DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);
557:     VecRestoreArrayRead(Grad, &grad);
558:   }
559:   if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
560:   VecRestoreArrayRead(faceGeometry, &facegeom);
561:   CHKERRQ(ierru);
562:   return(0);
563: }

565: PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
566: {
567:   PetscDS        prob;
568:   PetscInt       numBd, b;

572:   DMGetDS(dm, &prob);
573:   PetscDSGetNumBoundary(prob, &numBd);
574:   for (b = 0; b < numBd; ++b) {
575:     DMBoundaryConditionType type;
576:     const char             *labelname;
577:     DMLabel                 label;
578:     PetscInt                field, Nc;
579:     const PetscInt         *comps;
580:     PetscObject             obj;
581:     PetscClassId            id;
582:     void                    (*func)(void);
583:     PetscInt                numids;
584:     const PetscInt         *ids;
585:     void                   *ctx;

587:     DMGetBoundary(dm, b, &type, NULL, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);
588:     if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
589:     DMGetLabel(dm, labelname, &label);
590:     DMGetField(dm, field, NULL, &obj);
591:     PetscObjectGetClassId(obj, &id);
592:     if (id == PETSCFE_CLASSID) {
593:       switch (type) {
594:         /* for FEM, there is no insertion to be done for non-essential boundary conditions */
595:       case DM_BC_ESSENTIAL:
596:         DMPlexLabelAddCells(dm,label);
597:         DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
598:         DMPlexLabelClearCells(dm,label);
599:         break;
600:       case DM_BC_ESSENTIAL_FIELD:
601:         DMPlexLabelAddCells(dm,label);
602:         DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
603:                                                         (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
604:                                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
605:                                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);
606:         DMPlexLabelClearCells(dm,label);
607:         break;
608:       default: break;
609:       }
610:     } else if (id == PETSCFV_CLASSID) {
611:       if (!faceGeomFVM) continue;
612:       DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids,
613:                                                (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
614:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
615:   }
616:   return(0);
617: }

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

622:   Input Parameters:
623: + dm - The DM
624: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
625: . time - The time
626: . faceGeomFVM - Face geometry data for FV discretizations
627: . cellGeomFVM - Cell geometry data for FV discretizations
628: - gradFVM - Gradient reconstruction data for FV discretizations

630:   Output Parameters:
631: . locX - Solution updated with boundary values

633:   Level: developer

635: .seealso: DMProjectFunctionLabelLocal()
636: @*/
637: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
638: {

647:   PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));
648:   return(0);
649: }

651: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
652: {
653:   Vec              localX;
654:   PetscErrorCode   ierr;

657:   DMGetLocalVector(dm, &localX);
658:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);
659:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
660:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
661:   DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);
662:   DMRestoreLocalVector(dm, &localX);
663:   return(0);
664: }

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

669:   Input Parameters:
670: + dm     - The DM
671: . time   - The time
672: . funcs  - The functions to evaluate for each field component
673: . ctxs   - Optional array of contexts to pass to each function, or NULL.
674: - localX - The coefficient vector u_h, a local vector

676:   Output Parameter:
677: . diff - The diff ||u - u_h||_2

679:   Level: developer

681: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
682: @*/
683: PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff)
684: {
685:   const PetscInt   debug = ((DM_Plex*)dm->data)->printL2;
686:   PetscSection     section;
687:   PetscQuadrature  quad;
688:   PetscScalar     *funcVal, *interpolant;
689:   PetscReal       *coords, *detJ, *J;
690:   PetscReal        localDiff = 0.0;
691:   const PetscReal *quadWeights;
692:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, cEndInterior, c, field, fieldOffset;
693:   PetscErrorCode   ierr;

696:   DMGetDimension(dm, &dim);
697:   DMGetCoordinateDim(dm, &coordDim);
698:   DMGetSection(dm, &section);
699:   PetscSectionGetNumFields(section, &numFields);
700:   for (field = 0; field < numFields; ++field) {
701:     PetscObject  obj;
702:     PetscClassId id;
703:     PetscInt     Nc;

705:     DMGetField(dm, field, NULL, &obj);
706:     PetscObjectGetClassId(obj, &id);
707:     if (id == PETSCFE_CLASSID) {
708:       PetscFE fe = (PetscFE) obj;

710:       PetscFEGetQuadrature(fe, &quad);
711:       PetscFEGetNumComponents(fe, &Nc);
712:     } else if (id == PETSCFV_CLASSID) {
713:       PetscFV fv = (PetscFV) obj;

715:       PetscFVGetQuadrature(fv, &quad);
716:       PetscFVGetNumComponents(fv, &Nc);
717:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
718:     numComponents += Nc;
719:   }
720:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
721:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
722:   PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
723:   DMPlexGetVTKCellHeight(dm, &cellHeight);
724:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
725:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
726:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
727:   for (c = cStart; c < cEnd; ++c) {
728:     PetscScalar *x = NULL;
729:     PetscReal    elemDiff = 0.0;
730:     PetscInt     qc = 0;

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

735:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
736:       PetscObject  obj;
737:       PetscClassId id;
738:       void * const ctx = ctxs ? ctxs[field] : NULL;
739:       PetscInt     Nb, Nc, q, fc;

741:       DMGetField(dm, field, NULL, &obj);
742:       PetscObjectGetClassId(obj, &id);
743:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
744:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
745:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
746:       if (debug) {
747:         char title[1024];
748:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
749:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
750:       }
751:       for (q = 0; q < Nq; ++q) {
752:         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);
753:         (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx);
754:         if (ierr) {
755:           PetscErrorCode ierr2;
756:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
757:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
758:           ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
759: 
760:         }
761:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
762:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
763:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
764:         for (fc = 0; fc < Nc; ++fc) {
765:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
766:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
767:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
768:         }
769:       }
770:       fieldOffset += Nb;
771:       qc += Nc;
772:     }
773:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
774:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
775:     localDiff += elemDiff;
776:   }
777:   PetscFree5(funcVal,interpolant,coords,detJ,J);
778:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
779:   *diff = PetscSqrtReal(*diff);
780:   return(0);
781: }

783: 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)
784: {
785:   const PetscInt   debug = ((DM_Plex*)dm->data)->printL2;
786:   PetscSection     section;
787:   PetscQuadrature  quad;
788:   Vec              localX;
789:   PetscScalar     *funcVal, *interpolant;
790:   const PetscReal *quadPoints, *quadWeights;
791:   PetscReal       *coords, *realSpaceDer, *J, *invJ, *detJ;
792:   PetscReal        localDiff = 0.0;
793:   PetscInt         dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset;
794:   PetscErrorCode   ierr;

797:   DMGetDimension(dm, &dim);
798:   DMGetCoordinateDim(dm, &coordDim);
799:   DMGetSection(dm, &section);
800:   PetscSectionGetNumFields(section, &numFields);
801:   DMGetLocalVector(dm, &localX);
802:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
803:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
804:   for (field = 0; field < numFields; ++field) {
805:     PetscFE  fe;
806:     PetscInt Nc;

808:     DMGetField(dm, field, NULL, (PetscObject *) &fe);
809:     PetscFEGetQuadrature(fe, &quad);
810:     PetscFEGetNumComponents(fe, &Nc);
811:     numComponents += Nc;
812:   }
813:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
814:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
815:   /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
816:   PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);
817:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
818:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
819:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
820:   for (c = cStart; c < cEnd; ++c) {
821:     PetscScalar *x = NULL;
822:     PetscReal    elemDiff = 0.0;
823:     PetscInt     qc = 0;

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

828:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
829:       PetscFE          fe;
830:       void * const     ctx = ctxs ? ctxs[field] : NULL;
831:       PetscReal       *basisDer;
832:       PetscInt         Nb, Nc, q, fc;

834:       DMGetField(dm, field, NULL, (PetscObject *) &fe);
835:       PetscFEGetDimension(fe, &Nb);
836:       PetscFEGetNumComponents(fe, &Nc);
837:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
838:       if (debug) {
839:         char title[1024];
840:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
841:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
842:       }
843:       for (q = 0; q < Nq; ++q) {
844:         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);
845:         (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx);
846:         if (ierr) {
847:           PetscErrorCode ierr2;
848:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
849:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
850:           ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2);
851: 
852:         }
853:         PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);
854:         for (fc = 0; fc < Nc; ++fc) {
855:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
856:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
857:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
858:         }
859:       }
860:       fieldOffset += Nb;
861:       qc          += Nc;
862:     }
863:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
864:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
865:     localDiff += elemDiff;
866:   }
867:   PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);
868:   DMRestoreLocalVector(dm, &localX);
869:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
870:   *diff = PetscSqrtReal(*diff);
871:   return(0);
872: }

874: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
875: {
876:   const PetscInt   debug = ((DM_Plex*)dm->data)->printL2;
877:   PetscSection     section;
878:   PetscQuadrature  quad;
879:   Vec              localX;
880:   PetscScalar     *funcVal, *interpolant;
881:   PetscReal       *coords, *detJ, *J;
882:   PetscReal       *localDiff;
883:   const PetscReal *quadPoints, *quadWeights;
884:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
885:   PetscErrorCode   ierr;

888:   DMGetDimension(dm, &dim);
889:   DMGetCoordinateDim(dm, &coordDim);
890:   DMGetSection(dm, &section);
891:   PetscSectionGetNumFields(section, &numFields);
892:   DMGetLocalVector(dm, &localX);
893:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
894:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
895:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
896:   for (field = 0; field < numFields; ++field) {
897:     PetscObject  obj;
898:     PetscClassId id;
899:     PetscInt     Nc;

901:     DMGetField(dm, field, NULL, &obj);
902:     PetscObjectGetClassId(obj, &id);
903:     if (id == PETSCFE_CLASSID) {
904:       PetscFE fe = (PetscFE) obj;

906:       PetscFEGetQuadrature(fe, &quad);
907:       PetscFEGetNumComponents(fe, &Nc);
908:     } else if (id == PETSCFV_CLASSID) {
909:       PetscFV fv = (PetscFV) obj;

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

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

929:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
930:       PetscObject  obj;
931:       PetscClassId id;
932:       void * const ctx = ctxs ? ctxs[field] : NULL;
933:       PetscInt     Nb, Nc, q, fc;

935:       PetscReal       elemDiff = 0.0;

937:       DMGetField(dm, field, NULL, &obj);
938:       PetscObjectGetClassId(obj, &id);
939:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
940:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
941:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
942:       if (debug) {
943:         char title[1024];
944:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
945:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
946:       }
947:       for (q = 0; q < Nq; ++q) {
948:         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);
949:         (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx);
950:         if (ierr) {
951:           PetscErrorCode ierr2;
952:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
953:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
954:           ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
955: 
956:         }
957:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
958:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
959:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
960:         for (fc = 0; fc < Nc; ++fc) {
961:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
962:           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]);}
963:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
964:         }
965:       }
966:       fieldOffset += Nb;
967:       qc          += Nc;
968:       localDiff[field] += elemDiff;
969:     }
970:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
971:   }
972:   DMRestoreLocalVector(dm, &localX);
973:   MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
974:   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
975:   PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);
976:   return(0);
977: }

979: /*@C
980:   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.

982:   Input Parameters:
983: + dm    - The DM
984: . time  - The time
985: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
986: . ctxs  - Optional array of contexts to pass to each function, or NULL.
987: - X     - The coefficient vector u_h

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

992:   Level: developer

994: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
995: @*/
996: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
997: {
998:   PetscSection     section;
999:   PetscQuadrature  quad;
1000:   Vec              localX;
1001:   PetscScalar     *funcVal, *interpolant;
1002:   PetscReal       *coords, *detJ, *J;
1003:   const PetscReal *quadPoints, *quadWeights;
1004:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1005:   PetscErrorCode   ierr;

1008:   VecSet(D, 0.0);
1009:   DMGetDimension(dm, &dim);
1010:   DMGetCoordinateDim(dm, &coordDim);
1011:   DMGetSection(dm, &section);
1012:   PetscSectionGetNumFields(section, &numFields);
1013:   DMGetLocalVector(dm, &localX);
1014:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1015:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1016:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1017:   for (field = 0; field < numFields; ++field) {
1018:     PetscObject  obj;
1019:     PetscClassId id;
1020:     PetscInt     Nc;

1022:     DMGetField(dm, field, NULL, &obj);
1023:     PetscObjectGetClassId(obj, &id);
1024:     if (id == PETSCFE_CLASSID) {
1025:       PetscFE fe = (PetscFE) obj;

1027:       PetscFEGetQuadrature(fe, &quad);
1028:       PetscFEGetNumComponents(fe, &Nc);
1029:     } else if (id == PETSCFV_CLASSID) {
1030:       PetscFV fv = (PetscFV) obj;

1032:       PetscFVGetQuadrature(fv, &quad);
1033:       PetscFVGetNumComponents(fv, &Nc);
1034:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1035:     numComponents += Nc;
1036:   }
1037:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1038:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1039:   PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
1040:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1041:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1042:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1043:   for (c = cStart; c < cEnd; ++c) {
1044:     PetscScalar *x = NULL;
1045:     PetscScalar  elemDiff = 0.0;
1046:     PetscInt     qc = 0;

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

1051:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1052:       PetscObject  obj;
1053:       PetscClassId id;
1054:       void * const ctx = ctxs ? ctxs[field] : NULL;
1055:       PetscInt     Nb, Nc, q, fc;

1057:       DMGetField(dm, field, NULL, &obj);
1058:       PetscObjectGetClassId(obj, &id);
1059:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1060:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1061:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1062:       if (funcs[field]) {
1063:         for (q = 0; q < Nq; ++q) {
1064:           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);
1065:           (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx);
1066:           if (ierr) {
1067:             PetscErrorCode ierr2;
1068:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1069:             ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
1070:             ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1071: 
1072:           }
1073:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
1074:           else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1075:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1076:           for (fc = 0; fc < Nc; ++fc) {
1077:             const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1078:             elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
1079:           }
1080:         }
1081:       }
1082:       fieldOffset += Nb;
1083:       qc          += Nc;
1084:     }
1085:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1086:     VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1087:   }
1088:   PetscFree5(funcVal,interpolant,coords,detJ,J);
1089:   DMRestoreLocalVector(dm, &localX);
1090:   VecSqrtAbs(D);
1091:   return(0);
1092: }

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

1097:   Input Parameters:
1098: + dm - The DM
1099: - LocX  - The coefficient vector u_h

1101:   Output Parameter:
1102: . locC - A Vec which holds the Clement interpolant of the gradient

1104:   Notes:
1105:     Add citation to (Clement, 1975) and definition of the interpolant
1106:   \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

1108:   Level: developer

1110: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1111: @*/
1112: PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC)
1113: {
1114:   DM_Plex         *mesh  = (DM_Plex *) dm->data;
1115:   PetscInt         debug = mesh->printFEM;
1116:   DM               dmC;
1117:   PetscSection     section;
1118:   PetscQuadrature  quad;
1119:   PetscScalar     *interpolant, *gradsum;
1120:   PetscReal       *coords, *detJ, *J, *invJ;
1121:   const PetscReal *quadPoints, *quadWeights;
1122:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset;
1123:   PetscErrorCode   ierr;

1126:   VecGetDM(locC, &dmC);
1127:   VecSet(locC, 0.0);
1128:   DMGetDimension(dm, &dim);
1129:   DMGetCoordinateDim(dm, &coordDim);
1130:   DMGetSection(dm, &section);
1131:   PetscSectionGetNumFields(section, &numFields);
1132:   for (field = 0; field < numFields; ++field) {
1133:     PetscObject  obj;
1134:     PetscClassId id;
1135:     PetscInt     Nc;

1137:     DMGetField(dm, field, NULL, &obj);
1138:     PetscObjectGetClassId(obj, &id);
1139:     if (id == PETSCFE_CLASSID) {
1140:       PetscFE fe = (PetscFE) obj;

1142:       PetscFEGetQuadrature(fe, &quad);
1143:       PetscFEGetNumComponents(fe, &Nc);
1144:     } else if (id == PETSCFV_CLASSID) {
1145:       PetscFV fv = (PetscFV) obj;

1147:       PetscFVGetQuadrature(fv, &quad);
1148:       PetscFVGetNumComponents(fv, &Nc);
1149:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1150:     numComponents += Nc;
1151:   }
1152:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1153:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1154:   PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ);
1155:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1156:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1157:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1158:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1159:   for (v = vStart; v < vEnd; ++v) {
1160:     PetscScalar volsum = 0.0;
1161:     PetscInt   *star = NULL;
1162:     PetscInt    starSize, st, d, fc;

1164:     PetscMemzero(gradsum, coordDim*numComponents * sizeof(PetscScalar));
1165:     DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1166:     for (st = 0; st < starSize*2; st += 2) {
1167:       const PetscInt cell = star[st];
1168:       PetscScalar   *grad = &gradsum[coordDim*numComponents];
1169:       PetscScalar   *x    = NULL;
1170:       PetscReal      vol  = 0.0;

1172:       if ((cell < cStart) || (cell >= cEnd)) continue;
1173:       DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, J, invJ, detJ);
1174:       DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
1175:       for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1176:         PetscObject  obj;
1177:         PetscClassId id;
1178:         PetscInt     Nb, Nc, q, qc = 0;

1180:         PetscMemzero(grad, coordDim*numComponents * sizeof(PetscScalar));
1181:         DMGetField(dm, field, NULL, &obj);
1182:         PetscObjectGetClassId(obj, &id);
1183:         if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1184:         else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1185:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1186:         for (q = 0; q < Nq; ++q) {
1187:           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);
1188:           if (ierr) {
1189:             PetscErrorCode ierr2;
1190:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
1191:             ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
1192:             ierr2 = PetscFree4(interpolant,coords,detJ,J);CHKERRQ(ierr2);
1193: 
1194:           }
1195:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], coordDim, invJ, NULL, q, interpolant);}
1196:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1197:           for (fc = 0; fc < Nc; ++fc) {
1198:             const PetscReal wt = quadWeights[q*qNc+qc+fc];

1200:             for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*detJ[q];
1201:           }
1202:           vol += quadWeights[q*qNc]*detJ[q];
1203:         }
1204:         fieldOffset += Nb;
1205:         qc          += Nc;
1206:       }
1207:       DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
1208:       for (fc = 0; fc < numComponents; ++fc) {
1209:         for (d = 0; d < coordDim; ++d) {
1210:           gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
1211:         }
1212:       }
1213:       volsum += vol;
1214:       if (debug) {
1215:         PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
1216:         for (fc = 0; fc < numComponents; ++fc) {
1217:           for (d = 0; d < coordDim; ++d) {
1218:             if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
1219:             PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
1220:           }
1221:         }
1222:         PetscPrintf(PETSC_COMM_SELF, "]\n");
1223:       }
1224:     }
1225:     for (fc = 0; fc < numComponents; ++fc) {
1226:       for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum;
1227:     }
1228:     DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1229:     DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
1230:   }
1231:   PetscFree6(gradsum,interpolant,coords,detJ,J,invJ);
1232:   return(0);
1233: }

1235: static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user)
1236: {
1237:   DM                 dmAux = NULL;
1238:   PetscDS            prob,    probAux = NULL;
1239:   PetscSection       section, sectionAux;
1240:   Vec                locX,    locA;
1241:   PetscInt           dim, numCells = cEnd - cStart, c, f;
1242:   PetscBool          useFVM = PETSC_FALSE;
1243:   /* DS */
1244:   PetscInt           Nf,    totDim,    *uOff, *uOff_x, numConstants;
1245:   PetscInt           NfAux, totDimAux, *aOff;
1246:   PetscScalar       *u, *a;
1247:   const PetscScalar *constants;
1248:   /* Geometry */
1249:   PetscFEGeom       *cgeomFEM;
1250:   DM                 dmGrad;
1251:   PetscQuadrature    affineQuad = NULL;
1252:   Vec                cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
1253:   PetscFVCellGeom   *cgeomFVM;
1254:   const PetscScalar *lgrad;
1255:   PetscInt           maxDegree;
1256:   DMField            coordField;
1257:   IS                 cellIS;
1258:   PetscErrorCode     ierr;

1261:   DMGetDS(dm, &prob);
1262:   DMGetDimension(dm, &dim);
1263:   DMGetSection(dm, &section);
1264:   PetscSectionGetNumFields(section, &Nf);
1265:   /* Determine which discretizations we have */
1266:   for (f = 0; f < Nf; ++f) {
1267:     PetscObject  obj;
1268:     PetscClassId id;

1270:     PetscDSGetDiscretization(prob, f, &obj);
1271:     PetscObjectGetClassId(obj, &id);
1272:     if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE;
1273:   }
1274:   /* Get local solution with boundary values */
1275:   DMGetLocalVector(dm, &locX);
1276:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1277:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1278:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1279:   /* Read DS information */
1280:   PetscDSGetTotalDimension(prob, &totDim);
1281:   PetscDSGetComponentOffsets(prob, &uOff);
1282:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1283:   ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);
1284:   PetscDSGetConstants(prob, &numConstants, &constants);
1285:   /* Read Auxiliary DS information */
1286:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1287:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1288:   if (dmAux) {
1289:     DMGetDS(dmAux, &probAux);
1290:     PetscDSGetNumFields(probAux, &NfAux);
1291:     DMGetSection(dmAux, &sectionAux);
1292:     PetscDSGetTotalDimension(probAux, &totDimAux);
1293:     PetscDSGetComponentOffsets(probAux, &aOff);
1294:   }
1295:   /* Allocate data  arrays */
1296:   PetscCalloc1(numCells*totDim, &u);
1297:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1298:   /* Read out geometry */
1299:   DMGetCoordinateField(dm,&coordField);
1300:   DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
1301:   if (maxDegree <= 1) {
1302:     DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
1303:     if (affineQuad) {
1304:       DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);
1305:     }
1306:   }
1307:   if (useFVM) {
1308:     PetscFV   fv = NULL;
1309:     Vec       grad;
1310:     PetscInt  fStart, fEnd;
1311:     PetscBool compGrad;

1313:     for (f = 0; f < Nf; ++f) {
1314:       PetscObject  obj;
1315:       PetscClassId id;

1317:       PetscDSGetDiscretization(prob, f, &obj);
1318:       PetscObjectGetClassId(obj, &id);
1319:       if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;}
1320:     }
1321:     PetscFVGetComputeGradients(fv, &compGrad);
1322:     PetscFVSetComputeGradients(fv, PETSC_TRUE);
1323:     DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1324:     DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1325:     PetscFVSetComputeGradients(fv, compGrad);
1326:     VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1327:     /* Reconstruct and limit cell gradients */
1328:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1329:     DMGetGlobalVector(dmGrad, &grad);
1330:     DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1331:     /* Communicate gradient values */
1332:     DMGetLocalVector(dmGrad, &locGrad);
1333:     DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1334:     DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1335:     DMRestoreGlobalVector(dmGrad, &grad);
1336:     /* Handle non-essential (e.g. outflow) boundary values */
1337:     DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1338:     VecGetArrayRead(locGrad, &lgrad);
1339:   }
1340:   /* Read out data from inputs */
1341:   for (c = cStart; c < cEnd; ++c) {
1342:     PetscScalar *x = NULL;
1343:     PetscInt     i;

1345:     DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1346:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1347:     DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1348:     if (dmAux) {
1349:       DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1350:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1351:       DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1352:     }
1353:   }
1354:   /* Do integration for each field */
1355:   for (f = 0; f < Nf; ++f) {
1356:     PetscObject  obj;
1357:     PetscClassId id;
1358:     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1360:     PetscDSGetDiscretization(prob, f, &obj);
1361:     PetscObjectGetClassId(obj, &id);
1362:     if (id == PETSCFE_CLASSID) {
1363:       PetscFE         fe = (PetscFE) obj;
1364:       PetscQuadrature q;
1365:       PetscFEGeom     *chunkGeom = NULL;
1366:       PetscInt        Nq, Nb;

1368:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1369:       PetscFEGetQuadrature(fe, &q);
1370:       PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1371:       PetscFEGetDimension(fe, &Nb);
1372:       blockSize = Nb*Nq;
1373:       batchSize = numBlocks * blockSize;
1374:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1375:       numChunks = numCells / (numBatches*batchSize);
1376:       Ne        = numChunks*numBatches*batchSize;
1377:       Nr        = numCells % (numBatches*batchSize);
1378:       offset    = numCells - Nr;
1379:       if (!affineQuad) {
1380:         DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);
1381:       }
1382:       PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
1383:       PetscFEIntegrate(fe, prob, f, Ne, chunkGeom, u, probAux, a, cintegral);
1384:       PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);
1385:       PetscFEIntegrate(fe, prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);
1386:       PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);
1387:       if (!affineQuad) {
1388:         PetscFEGeomDestroy(&cgeomFEM);
1389:       }
1390:     } else if (id == PETSCFV_CLASSID) {
1391:       PetscInt       foff;
1392:       PetscPointFunc obj_func;
1393:       PetscScalar    lint;

1395:       PetscDSGetObjective(prob, f, &obj_func);
1396:       PetscDSGetFieldOffset(prob, f, &foff);
1397:       if (obj_func) {
1398:         for (c = 0; c < numCells; ++c) {
1399:           PetscScalar *u_x;

1401:           DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1402:           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);
1403:           cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1404:         }
1405:       }
1406:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1407:   }
1408:   /* Cleanup data arrays */
1409:   if (useFVM) {
1410:     VecRestoreArrayRead(locGrad, &lgrad);
1411:     VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1412:     DMRestoreLocalVector(dmGrad, &locGrad);
1413:     VecDestroy(&faceGeometryFVM);
1414:     VecDestroy(&cellGeometryFVM);
1415:     DMDestroy(&dmGrad);
1416:   }
1417:   if (dmAux) {PetscFree(a);}
1418:   PetscFree(u);
1419:   /* Cleanup */
1420:   if (affineQuad) {
1421:     PetscFEGeomDestroy(&cgeomFEM);
1422:   }
1423:   PetscQuadratureDestroy(&affineQuad);
1424:   ISDestroy(&cellIS);
1425:   DMRestoreLocalVector(dm, &locX);
1426:   return(0);
1427: }

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

1432:   Input Parameters:
1433: + dm - The mesh
1434: . X  - Global input vector
1435: - user - The user context

1437:   Output Parameter:
1438: . integral - Integral for each field

1440:   Level: developer

1442: .seealso: DMPlexComputeResidualFEM()
1443: @*/
1444: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
1445: {
1446:   DM_Plex       *mesh = (DM_Plex *) dm->data;
1447:   PetscScalar   *cintegral, *lintegral;
1448:   PetscInt       Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;

1455:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1456:   DMGetNumFields(dm, &Nf);
1457:   DMPlexGetVTKCellHeight(dm, &cellHeight);
1458:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1459:   DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1460:   cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1461:   /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1462:   PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);
1463:   DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1464:   /* Sum up values */
1465:   for (cell = cStart; cell < cEnd; ++cell) {
1466:     const PetscInt c = cell - cStart;

1468:     if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1469:     for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f];
1470:   }
1471:   MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1472:   if (mesh->printFEM) {
1473:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");
1474:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));}
1475:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1476:   }
1477:   PetscFree2(lintegral, cintegral);
1478:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1479:   return(0);
1480: }

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

1485:   Input Parameters:
1486: + dm - The mesh
1487: . X  - Global input vector
1488: - user - The user context

1490:   Output Parameter:
1491: . integral - Cellwise integrals for each field

1493:   Level: developer

1495: .seealso: DMPlexComputeResidualFEM()
1496: @*/
1497: PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user)
1498: {
1499:   DM_Plex       *mesh = (DM_Plex *) dm->data;
1500:   DM             dmF;
1501:   PetscSection   sectionF;
1502:   PetscScalar   *cintegral, *af;
1503:   PetscInt       Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;

1510:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1511:   DMGetNumFields(dm, &Nf);
1512:   DMPlexGetVTKCellHeight(dm, &cellHeight);
1513:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1514:   DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1515:   cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1516:   /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1517:   PetscCalloc1((cEnd-cStart)*Nf, &cintegral);
1518:   DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1519:   /* Put values in F*/
1520:   VecGetDM(F, &dmF);
1521:   DMGetSection(dmF, &sectionF);
1522:   VecGetArray(F, &af);
1523:   for (cell = cStart; cell < cEnd; ++cell) {
1524:     const PetscInt c = cell - cStart;
1525:     PetscInt       dof, off;

1527:     if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1528:     PetscSectionGetDof(sectionF, cell, &dof);
1529:     PetscSectionGetOffset(sectionF, cell, &off);
1530:     if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf);
1531:     for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f];
1532:   }
1533:   VecRestoreArray(F, &af);
1534:   PetscFree(cintegral);
1535:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1536:   return(0);
1537: }

1539: static PetscErrorCode DMPlexComputeBdIntegral_Internal(DM dm, Vec locX, IS pointIS,
1540:                                                        void (*func)(PetscInt, PetscInt, PetscInt,
1541:                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1542:                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1543:                                                                     PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1544:                                                        PetscScalar *fintegral, void *user)
1545: {
1546:   DM                 plex = NULL, plexA = NULL;
1547:   PetscDS            prob, probAux = NULL;
1548:   PetscSection       section, sectionAux = NULL;
1549:   Vec                locA = NULL;
1550:   DMField            coordField;
1551:   PetscInt           Nf,        totDim,        *uOff, *uOff_x;
1552:   PetscInt           NfAux = 0, totDimAux = 0, *aOff = NULL;
1553:   PetscScalar       *u, *a = NULL;
1554:   const PetscScalar *constants;
1555:   PetscInt           numConstants, f;
1556:   PetscErrorCode     ierr;

1559:   DMGetCoordinateField(dm, &coordField);
1560:   DMConvert(dm, DMPLEX, &plex);
1561:   DMGetDS(dm, &prob);
1562:   DMGetDefaultSection(dm, &section);
1563:   PetscSectionGetNumFields(section, &Nf);
1564:   /* Determine which discretizations we have */
1565:   for (f = 0; f < Nf; ++f) {
1566:     PetscObject  obj;
1567:     PetscClassId id;

1569:     PetscDSGetDiscretization(prob, f, &obj);
1570:     PetscObjectGetClassId(obj, &id);
1571:     if (id == PETSCFV_CLASSID) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not supported for FVM (field %D)", f);
1572:   }
1573:   /* Read DS information */
1574:   PetscDSGetTotalDimension(prob, &totDim);
1575:   PetscDSGetComponentOffsets(prob, &uOff);
1576:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1577:   PetscDSGetConstants(prob, &numConstants, &constants);
1578:   /* Read Auxiliary DS information */
1579:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1580:   if (locA) {
1581:     DM dmAux;

1583:     VecGetDM(locA, &dmAux);
1584:     DMConvert(dmAux, DMPLEX, &plexA);
1585:     DMGetDS(dmAux, &probAux);
1586:     PetscDSGetNumFields(probAux, &NfAux);
1587:     DMGetDefaultSection(dmAux, &sectionAux);
1588:     PetscDSGetTotalDimension(probAux, &totDimAux);
1589:     PetscDSGetComponentOffsets(probAux, &aOff);
1590:   }
1591:   /* Integrate over points */
1592:   {
1593:     PetscFEGeom    *fgeom, *chunkGeom = NULL;
1594:     PetscInt        maxDegree;
1595:     PetscQuadrature qGeom = NULL;
1596:     const PetscInt *points;
1597:     PetscInt        numFaces, face, Nq, field;
1598:     PetscInt        numChunks, chunkSize, chunk, Nr, offset;

1600:     ISGetLocalSize(pointIS, &numFaces);
1601:     ISGetIndices(pointIS, &points);
1602:     PetscCalloc2(numFaces*totDim, &u, locA ? numFaces*totDimAux : 0, &a);
1603:     DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
1604:     for (field = 0; field < Nf; ++field) {
1605:       PetscFE fe;

1607:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1608:       if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, pointIS, &qGeom);}
1609:       if (!qGeom) {
1610:         PetscFEGetFaceQuadrature(fe, &qGeom);
1611:         PetscObjectReference((PetscObject) qGeom);
1612:       }
1613:       PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
1614:       DMPlexGetFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
1615:       for (face = 0; face < numFaces; ++face) {
1616:         const PetscInt point = points[face], *support, *cone;
1617:         PetscScalar    *x    = NULL;
1618:         PetscInt       i, coneSize, faceLoc;

1620:         DMPlexGetSupport(dm, point, &support);
1621:         DMPlexGetConeSize(dm, support[0], &coneSize);
1622:         DMPlexGetCone(dm, support[0], &cone);
1623:         for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break;
1624:         if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of support[0] %D", face, support[0]);
1625:         fgeom->face[face][0] = faceLoc;
1626:         DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);
1627:         for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
1628:         DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);
1629:         if (locA) {
1630:           PetscInt subp;
1631:           DMPlexGetSubpoint(plexA, support[0], &subp);
1632:           DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);
1633:           for (i = 0; i < totDimAux; ++i) a[f*totDimAux+i] = x[i];
1634:           DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);
1635:         }
1636:       }
1637:       /* Get blocking */
1638:       {
1639:         PetscQuadrature q;
1640:         PetscInt        numBatches, batchSize, numBlocks, blockSize;
1641:         PetscInt        Nq, Nb;

1643:         PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1644:         PetscFEGetQuadrature(fe, &q);
1645:         PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1646:         PetscFEGetDimension(fe, &Nb);
1647:         blockSize = Nb*Nq;
1648:         batchSize = numBlocks * blockSize;
1649:         chunkSize = numBatches*batchSize;
1650:         PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1651:         numChunks = numFaces / chunkSize;
1652:         Nr        = numFaces % chunkSize;
1653:         offset    = numFaces - Nr;
1654:       }
1655:       /* Do integration for each field */
1656:       for (chunk = 0; chunk < numChunks; ++chunk) {
1657:         PetscFEGeomGetChunk(fgeom, chunk*chunkSize, (chunk+1)*chunkSize, &chunkGeom);
1658:         PetscFEIntegrateBd(fe, prob, field, func, chunkSize, chunkGeom, u, probAux, a, fintegral);
1659:         PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);
1660:       }
1661:       PetscFEGeomGetChunk(fgeom, offset, numFaces, &chunkGeom);
1662:       PetscFEIntegrateBd(fe, prob, field, func, Nr, chunkGeom, &u[offset*totDim], probAux, a ? &a[offset*totDimAux] : NULL, &fintegral[offset*Nf]);
1663:       PetscFEGeomRestoreChunk(fgeom, offset, numFaces, &chunkGeom);
1664:       /* Cleanup data arrays */
1665:       DMPlexRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
1666:       PetscQuadratureDestroy(&qGeom);
1667:       PetscFree2(u, a);
1668:       ISRestoreIndices(pointIS, &points);
1669:     }
1670:   }
1671:   if (plex)  {DMDestroy(&plex);}
1672:   if (plexA) {DMDestroy(&plexA);}
1673:   return(0);
1674: }

1676: /*@
1677:   DMPlexComputeBdIntegral - Form the integral over the specified boundary from the global input X using pointwise functions specified by the user

1679:   Input Parameters:
1680: + dm      - The mesh
1681: . X       - Global input vector
1682: . label   - The boundary DMLabel
1683: . numVals - The number of label values to use, or PETSC_DETERMINE for all values
1684: . vals    - The label values to use, or PETSC_NULL for all values
1685: . func    = The function to integrate along the boundary
1686: - user    - The user context

1688:   Output Parameter:
1689: . integral - Integral for each field

1691:   Level: developer

1693: .seealso: DMPlexComputeIntegralFEM(), DMPlexComputeBdResidualFEM()
1694: @*/
1695: PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt numVals, const PetscInt vals[],
1696:                                        void (*func)(PetscInt, PetscInt, PetscInt,
1697:                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1698:                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1699:                                                     PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
1700:                                        PetscScalar *integral, void *user)
1701: {
1702:   Vec            locX;
1703:   PetscSection   section;
1704:   DMLabel        depthLabel;
1705:   IS             facetIS;
1706:   PetscInt       dim, Nf, f, v;

1715:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1716:   DMPlexGetDepthLabel(dm, &depthLabel);
1717:   DMGetDimension(dm, &dim);
1718:   DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);
1719:   DMGetDefaultSection(dm, &section);
1720:   PetscSectionGetNumFields(section, &Nf);
1721:   /* Get local solution with boundary values */
1722:   DMGetLocalVector(dm, &locX);
1723:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1724:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1725:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1726:   /* Loop over label values */
1727:   PetscMemzero(integral, Nf * sizeof(PetscScalar));
1728:   for (v = 0; v < numVals; ++v) {
1729:     IS           pointIS;
1730:     PetscInt     numFaces, face;
1731:     PetscScalar *fintegral;

1733:     DMLabelGetStratumIS(label, vals[v], &pointIS);
1734:     if (!pointIS) continue; /* No points with that id on this process */
1735:     {
1736:       IS isectIS;

1738:       /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */
1739:       ISIntersect_Caching_Internal(facetIS, pointIS, &isectIS);
1740:       ISDestroy(&pointIS);
1741:       pointIS = isectIS;
1742:     }
1743:     ISGetLocalSize(pointIS, &numFaces);
1744:     PetscCalloc1(numFaces*Nf, &fintegral);
1745:     DMPlexComputeBdIntegral_Internal(dm, locX, pointIS, func, fintegral, user);
1746:     /* Sum point contributions into integral */
1747:     for (f = 0; f < Nf; ++f) for (face = 0; face < numFaces; ++face) integral[f] += fintegral[face*Nf+f];
1748:     PetscFree(fintegral);
1749:     ISDestroy(&pointIS);
1750:   }
1751:   DMRestoreLocalVector(dm, &locX);
1752:   ISDestroy(&facetIS);
1753:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1754:   return(0);
1755: }

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

1760:   Input Parameters:
1761: + dmf  - The fine mesh
1762: . dmc  - The coarse mesh
1763: - user - The user context

1765:   Output Parameter:
1766: . In  - The interpolation matrix

1768:   Level: developer

1770: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1771: @*/
1772: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1773: {
1774:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1775:   const char       *name  = "Interpolator";
1776:   PetscDS           prob;
1777:   PetscFE          *feRef;
1778:   PetscFV          *fvRef;
1779:   PetscSection      fsection, fglobalSection;
1780:   PetscSection      csection, cglobalSection;
1781:   PetscScalar      *elemMat;
1782:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1783:   PetscInt          cTotDim, rTotDim = 0;
1784:   PetscErrorCode    ierr;

1787:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1788:   DMGetDimension(dmf, &dim);
1789:   DMGetSection(dmf, &fsection);
1790:   DMGetGlobalSection(dmf, &fglobalSection);
1791:   DMGetSection(dmc, &csection);
1792:   DMGetGlobalSection(dmc, &cglobalSection);
1793:   PetscSectionGetNumFields(fsection, &Nf);
1794:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1795:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1796:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1797:   DMGetDS(dmf, &prob);
1798:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1799:   for (f = 0; f < Nf; ++f) {
1800:     PetscObject  obj;
1801:     PetscClassId id;
1802:     PetscInt     rNb = 0, Nc = 0;

1804:     PetscDSGetDiscretization(prob, f, &obj);
1805:     PetscObjectGetClassId(obj, &id);
1806:     if (id == PETSCFE_CLASSID) {
1807:       PetscFE fe = (PetscFE) obj;

1809:       PetscFERefine(fe, &feRef[f]);
1810:       PetscFEGetDimension(feRef[f], &rNb);
1811:       PetscFEGetNumComponents(fe, &Nc);
1812:     } else if (id == PETSCFV_CLASSID) {
1813:       PetscFV        fv = (PetscFV) obj;
1814:       PetscDualSpace Q;

1816:       PetscFVRefine(fv, &fvRef[f]);
1817:       PetscFVGetDualSpace(fvRef[f], &Q);
1818:       PetscDualSpaceGetDimension(Q, &rNb);
1819:       PetscFVGetNumComponents(fv, &Nc);
1820:     }
1821:     rTotDim += rNb;
1822:   }
1823:   PetscDSGetTotalDimension(prob, &cTotDim);
1824:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1825:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1826:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1827:     PetscDualSpace   Qref;
1828:     PetscQuadrature  f;
1829:     const PetscReal *qpoints, *qweights;
1830:     PetscReal       *points;
1831:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1833:     /* Compose points from all dual basis functionals */
1834:     if (feRef[fieldI]) {
1835:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
1836:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
1837:     } else {
1838:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1839:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1840:     }
1841:     PetscDualSpaceGetDimension(Qref, &fpdim);
1842:     for (i = 0; i < fpdim; ++i) {
1843:       PetscDualSpaceGetFunctional(Qref, i, &f);
1844:       PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
1845:       npoints += Np;
1846:     }
1847:     PetscMalloc1(npoints*dim,&points);
1848:     for (i = 0, k = 0; i < fpdim; ++i) {
1849:       PetscDualSpaceGetFunctional(Qref, i, &f);
1850:       PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1851:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1852:     }

1854:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1855:       PetscObject  obj;
1856:       PetscClassId id;
1857:       PetscReal   *B;
1858:       PetscInt     NcJ = 0, cpdim = 0, j, qNc;

1860:       PetscDSGetDiscretization(prob, fieldJ, &obj);
1861:       PetscObjectGetClassId(obj, &id);
1862:       if (id == PETSCFE_CLASSID) {
1863:         PetscFE fe = (PetscFE) obj;

1865:         /* Evaluate basis at points */
1866:         PetscFEGetNumComponents(fe, &NcJ);
1867:         PetscFEGetDimension(fe, &cpdim);
1868:         /* For now, fields only interpolate themselves */
1869:         if (fieldI == fieldJ) {
1870:           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);
1871:           PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1872:           for (i = 0, k = 0; i < fpdim; ++i) {
1873:             PetscDualSpaceGetFunctional(Qref, i, &f);
1874:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1875:             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);
1876:             for (p = 0; p < Np; ++p, ++k) {
1877:               for (j = 0; j < cpdim; ++j) {
1878:                 /*
1879:                    cTotDim:            Total columns in element interpolation matrix, sum of number of dual basis functionals in each field
1880:                    offsetI, offsetJ:   Offsets into the larger element interpolation matrix for different fields
1881:                    fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals
1882:                    qNC, Nc, Ncj, c:    Number of components in this field
1883:                    Np, p:              Number of quad points in the fine grid functional i
1884:                    k:                  i*Np + p, overall point number for the interpolation
1885:                 */
1886:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
1887:               }
1888:             }
1889:           }
1890:           PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1891:         }
1892:       } else if (id == PETSCFV_CLASSID) {
1893:         PetscFV        fv = (PetscFV) obj;

1895:         /* Evaluate constant function at points */
1896:         PetscFVGetNumComponents(fv, &NcJ);
1897:         cpdim = 1;
1898:         /* For now, fields only interpolate themselves */
1899:         if (fieldI == fieldJ) {
1900:           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);
1901:           for (i = 0, k = 0; i < fpdim; ++i) {
1902:             PetscDualSpaceGetFunctional(Qref, i, &f);
1903:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1904:             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);
1905:             for (p = 0; p < Np; ++p, ++k) {
1906:               for (j = 0; j < cpdim; ++j) {
1907:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += 1.0*qweights[p*qNc+c];
1908:               }
1909:             }
1910:           }
1911:         }
1912:       }
1913:       offsetJ += cpdim;
1914:     }
1915:     offsetI += fpdim;
1916:     PetscFree(points);
1917:   }
1918:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1919:   /* Preallocate matrix */
1920:   {
1921:     Mat          preallocator;
1922:     PetscScalar *vals;
1923:     PetscInt    *cellCIndices, *cellFIndices;
1924:     PetscInt     locRows, locCols, cell;

1926:     MatGetLocalSize(In, &locRows, &locCols);
1927:     MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1928:     MatSetType(preallocator, MATPREALLOCATOR);
1929:     MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1930:     MatSetUp(preallocator);
1931:     PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1932:     for (cell = cStart; cell < cEnd; ++cell) {
1933:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1934:       MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1935:     }
1936:     PetscFree3(vals,cellCIndices,cellFIndices);
1937:     MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1938:     MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1939:     MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1940:     MatDestroy(&preallocator);
1941:   }
1942:   /* Fill matrix */
1943:   MatZeroEntries(In);
1944:   for (c = cStart; c < cEnd; ++c) {
1945:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1946:   }
1947:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1948:   PetscFree2(feRef,fvRef);
1949:   PetscFree(elemMat);
1950:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1951:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1952:   if (mesh->printFEM) {
1953:     PetscPrintf(PetscObjectComm((PetscObject)In), "%s:\n", name);
1954:     MatChop(In, 1.0e-10);
1955:     MatView(In, NULL);
1956:   }
1957:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1958:   return(0);
1959: }

1961: PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user)
1962: {
1963:   SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness");
1964: }

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

1969:   Input Parameters:
1970: + dmf  - The fine mesh
1971: . dmc  - The coarse mesh
1972: - user - The user context

1974:   Output Parameter:
1975: . In  - The interpolation matrix

1977:   Level: developer

1979: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1980: @*/
1981: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1982: {
1983:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
1984:   const char    *name = "Interpolator";
1985:   PetscDS        prob;
1986:   PetscSection   fsection, csection, globalFSection, globalCSection;
1987:   PetscHSetIJ    ht;
1988:   PetscLayout    rLayout;
1989:   PetscInt      *dnz, *onz;
1990:   PetscInt       locRows, rStart, rEnd;
1991:   PetscReal     *x, *v0, *J, *invJ, detJ;
1992:   PetscReal     *v0c, *Jc, *invJc, detJc;
1993:   PetscScalar   *elemMat;
1994:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

1998:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1999:   DMGetCoordinateDim(dmc, &dim);
2000:   DMGetDS(dmc, &prob);
2001:   PetscDSGetRefCoordArrays(prob, &x, NULL);
2002:   PetscDSGetNumFields(prob, &Nf);
2003:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2004:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2005:   DMGetSection(dmf, &fsection);
2006:   DMGetGlobalSection(dmf, &globalFSection);
2007:   DMGetSection(dmc, &csection);
2008:   DMGetGlobalSection(dmc, &globalCSection);
2009:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2010:   PetscDSGetTotalDimension(prob, &totDim);
2011:   PetscMalloc1(totDim, &elemMat);

2013:   MatGetLocalSize(In, &locRows, NULL);
2014:   PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
2015:   PetscLayoutSetLocalSize(rLayout, locRows);
2016:   PetscLayoutSetBlockSize(rLayout, 1);
2017:   PetscLayoutSetUp(rLayout);
2018:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2019:   PetscLayoutDestroy(&rLayout);
2020:   PetscCalloc2(locRows,&dnz,locRows,&onz);
2021:   PetscHSetIJCreate(&ht);
2022:   for (field = 0; field < Nf; ++field) {
2023:     PetscObject      obj;
2024:     PetscClassId     id;
2025:     PetscDualSpace   Q = NULL;
2026:     PetscQuadrature  f;
2027:     const PetscReal *qpoints;
2028:     PetscInt         Nc, Np, fpdim, i, d;

2030:     PetscDSGetDiscretization(prob, field, &obj);
2031:     PetscObjectGetClassId(obj, &id);
2032:     if (id == PETSCFE_CLASSID) {
2033:       PetscFE fe = (PetscFE) obj;

2035:       PetscFEGetDualSpace(fe, &Q);
2036:       PetscFEGetNumComponents(fe, &Nc);
2037:     } else if (id == PETSCFV_CLASSID) {
2038:       PetscFV fv = (PetscFV) obj;

2040:       PetscFVGetDualSpace(fv, &Q);
2041:       Nc   = 1;
2042:     }
2043:     PetscDualSpaceGetDimension(Q, &fpdim);
2044:     /* For each fine grid cell */
2045:     for (cell = cStart; cell < cEnd; ++cell) {
2046:       PetscInt *findices,   *cindices;
2047:       PetscInt  numFIndices, numCIndices;

2049:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2050:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2051:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
2052:       for (i = 0; i < fpdim; ++i) {
2053:         Vec             pointVec;
2054:         PetscScalar    *pV;
2055:         PetscSF         coarseCellSF = NULL;
2056:         const PetscSFNode *coarseCells;
2057:         PetscInt        numCoarseCells, q, c;

2059:         /* Get points from the dual basis functional quadrature */
2060:         PetscDualSpaceGetFunctional(Q, i, &f);
2061:         PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2062:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2063:         VecSetBlockSize(pointVec, dim);
2064:         VecGetArray(pointVec, &pV);
2065:         for (q = 0; q < Np; ++q) {
2066:           const PetscReal xi0[3] = {-1., -1., -1.};

2068:           /* Transform point to real space */
2069:           CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2070:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2071:         }
2072:         VecRestoreArray(pointVec, &pV);
2073:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2074:         /* OPT: Pack all quad points from fine cell */
2075:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2076:         PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2077:         /* Update preallocation info */
2078:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2079:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2080:         {
2081:           PetscHashIJKey key;
2082:           PetscBool      missing;

2084:           key.i = findices[i];
2085:           if (key.i >= 0) {
2086:             /* Get indices for coarse elements */
2087:             for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2088:               DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2089:               for (c = 0; c < numCIndices; ++c) {
2090:                 key.j = cindices[c];
2091:                 if (key.j < 0) continue;
2092:                 PetscHSetIJQueryAdd(ht, key, &missing);
2093:                 if (missing) {
2094:                   if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2095:                   else                                     ++onz[key.i-rStart];
2096:                 }
2097:               }
2098:               DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2099:             }
2100:           }
2101:         }
2102:         PetscSFDestroy(&coarseCellSF);
2103:         VecDestroy(&pointVec);
2104:       }
2105:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2106:     }
2107:   }
2108:   PetscHSetIJDestroy(&ht);
2109:   MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
2110:   MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2111:   PetscFree2(dnz,onz);
2112:   for (field = 0; field < Nf; ++field) {
2113:     PetscObject      obj;
2114:     PetscClassId     id;
2115:     PetscDualSpace   Q = NULL;
2116:     PetscQuadrature  f;
2117:     const PetscReal *qpoints, *qweights;
2118:     PetscInt         Nc, qNc, Np, fpdim, i, d;

2120:     PetscDSGetDiscretization(prob, field, &obj);
2121:     PetscObjectGetClassId(obj, &id);
2122:     if (id == PETSCFE_CLASSID) {
2123:       PetscFE fe = (PetscFE) obj;

2125:       PetscFEGetDualSpace(fe, &Q);
2126:       PetscFEGetNumComponents(fe, &Nc);
2127:     } else if (id == PETSCFV_CLASSID) {
2128:       PetscFV fv = (PetscFV) obj;

2130:       PetscFVGetDualSpace(fv, &Q);
2131:       Nc   = 1;
2132:     } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %d",field);
2133:     PetscDualSpaceGetDimension(Q, &fpdim);
2134:     /* For each fine grid cell */
2135:     for (cell = cStart; cell < cEnd; ++cell) {
2136:       PetscInt *findices,   *cindices;
2137:       PetscInt  numFIndices, numCIndices;

2139:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2140:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2141:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
2142:       for (i = 0; i < fpdim; ++i) {
2143:         Vec             pointVec;
2144:         PetscScalar    *pV;
2145:         PetscSF         coarseCellSF = NULL;
2146:         const PetscSFNode *coarseCells;
2147:         PetscInt        numCoarseCells, cpdim, q, c, j;

2149:         /* Get points from the dual basis functional quadrature */
2150:         PetscDualSpaceGetFunctional(Q, i, &f);
2151:         PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
2152:         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);
2153:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2154:         VecSetBlockSize(pointVec, dim);
2155:         VecGetArray(pointVec, &pV);
2156:         for (q = 0; q < Np; ++q) {
2157:           const PetscReal xi0[3] = {-1., -1., -1.};

2159:           /* Transform point to real space */
2160:           CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2161:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2162:         }
2163:         VecRestoreArray(pointVec, &pV);
2164:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2165:         /* OPT: Read this out from preallocation information */
2166:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2167:         /* Update preallocation info */
2168:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2169:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2170:         VecGetArray(pointVec, &pV);
2171:         for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2172:           PetscReal pVReal[3];
2173:           const PetscReal xi0[3] = {-1., -1., -1.};

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

2181:           if (id == PETSCFE_CLASSID) {
2182:             PetscFE    fe = (PetscFE) obj;
2183:             PetscReal *B;

2185:             /* Evaluate coarse basis on contained point */
2186:             PetscFEGetDimension(fe, &cpdim);
2187:             PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2188:             PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2189:             /* Get elemMat entries by multiplying by weight */
2190:             for (j = 0; j < cpdim; ++j) {
2191:               for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c];
2192:             }
2193:             PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2194:           } else {
2195:             cpdim = 1;
2196:             for (j = 0; j < cpdim; ++j) {
2197:               for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
2198:             }
2199:           }
2200:           /* Update interpolator */
2201:           if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2202:           if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2203:           MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
2204:           DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2205:         }
2206:         VecRestoreArray(pointVec, &pV);
2207:         PetscSFDestroy(&coarseCellSF);
2208:         VecDestroy(&pointVec);
2209:       }
2210:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2211:     }
2212:   }
2213:   PetscFree3(v0,J,invJ);
2214:   PetscFree3(v0c,Jc,invJc);
2215:   PetscFree(elemMat);
2216:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2217:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2218:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2219:   return(0);
2220: }

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

2225:   Input Parameters:
2226: + dmf  - The fine mesh
2227: . dmc  - The coarse mesh
2228: - user - The user context

2230:   Output Parameter:
2231: . mass  - The mass matrix

2233:   Level: developer

2235: .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2236: @*/
2237: PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user)
2238: {
2239:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
2240:   const char    *name = "Mass Matrix";
2241:   PetscDS        prob;
2242:   PetscSection   fsection, csection, globalFSection, globalCSection;
2243:   PetscHSetIJ    ht;
2244:   PetscLayout    rLayout;
2245:   PetscInt      *dnz, *onz;
2246:   PetscInt       locRows, rStart, rEnd;
2247:   PetscReal     *x, *v0, *J, *invJ, detJ;
2248:   PetscReal     *v0c, *Jc, *invJc, detJc;
2249:   PetscScalar   *elemMat;
2250:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

2254:   DMGetCoordinateDim(dmc, &dim);
2255:   DMGetDS(dmc, &prob);
2256:   PetscDSGetRefCoordArrays(prob, &x, NULL);
2257:   PetscDSGetNumFields(prob, &Nf);
2258:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2259:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2260:   DMGetSection(dmf, &fsection);
2261:   DMGetGlobalSection(dmf, &globalFSection);
2262:   DMGetSection(dmc, &csection);
2263:   DMGetGlobalSection(dmc, &globalCSection);
2264:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2265:   PetscDSGetTotalDimension(prob, &totDim);
2266:   PetscMalloc1(totDim, &elemMat);

2268:   MatGetLocalSize(mass, &locRows, NULL);
2269:   PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);
2270:   PetscLayoutSetLocalSize(rLayout, locRows);
2271:   PetscLayoutSetBlockSize(rLayout, 1);
2272:   PetscLayoutSetUp(rLayout);
2273:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2274:   PetscLayoutDestroy(&rLayout);
2275:   PetscCalloc2(locRows,&dnz,locRows,&onz);
2276:   PetscHSetIJCreate(&ht);
2277:   for (field = 0; field < Nf; ++field) {
2278:     PetscObject      obj;
2279:     PetscClassId     id;
2280:     PetscQuadrature  quad;
2281:     const PetscReal *qpoints;
2282:     PetscInt         Nq, Nc, i, d;

2284:     PetscDSGetDiscretization(prob, field, &obj);
2285:     PetscObjectGetClassId(obj, &id);
2286:     if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);}
2287:     else                       {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2288:     PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);
2289:     /* For each fine grid cell */
2290:     for (cell = cStart; cell < cEnd; ++cell) {
2291:       Vec                pointVec;
2292:       PetscScalar       *pV;
2293:       PetscSF            coarseCellSF = NULL;
2294:       const PetscSFNode *coarseCells;
2295:       PetscInt           numCoarseCells, q, c;
2296:       PetscInt          *findices,   *cindices;
2297:       PetscInt           numFIndices, numCIndices;

2299:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2300:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2301:       /* Get points from the quadrature */
2302:       VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2303:       VecSetBlockSize(pointVec, dim);
2304:       VecGetArray(pointVec, &pV);
2305:       for (q = 0; q < Nq; ++q) {
2306:         const PetscReal xi0[3] = {-1., -1., -1.};

2308:         /* Transform point to real space */
2309:         CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2310:         for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2311:       }
2312:       VecRestoreArray(pointVec, &pV);
2313:       /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2314:       DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2315:       PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2316:       /* Update preallocation info */
2317:       PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2318:       if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2319:       {
2320:         PetscHashIJKey key;
2321:         PetscBool      missing;

2323:         for (i = 0; i < numFIndices; ++i) {
2324:           key.i = findices[i];
2325:           if (key.i >= 0) {
2326:             /* Get indices for coarse elements */
2327:             for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2328:               DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2329:               for (c = 0; c < numCIndices; ++c) {
2330:                 key.j = cindices[c];
2331:                 if (key.j < 0) continue;
2332:                 PetscHSetIJQueryAdd(ht, key, &missing);
2333:                 if (missing) {
2334:                   if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2335:                   else                                     ++onz[key.i-rStart];
2336:                 }
2337:               }
2338:               DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2339:             }
2340:           }
2341:         }
2342:       }
2343:       PetscSFDestroy(&coarseCellSF);
2344:       VecDestroy(&pointVec);
2345:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2346:     }
2347:   }
2348:   PetscHSetIJDestroy(&ht);
2349:   MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
2350:   MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2351:   PetscFree2(dnz,onz);
2352:   for (field = 0; field < Nf; ++field) {
2353:     PetscObject      obj;
2354:     PetscClassId     id;
2355:     PetscQuadrature  quad;
2356:     PetscReal       *Bfine;
2357:     const PetscReal *qpoints, *qweights;
2358:     PetscInt         Nq, Nc, i, d;

2360:     PetscDSGetDiscretization(prob, field, &obj);
2361:     PetscObjectGetClassId(obj, &id);
2362:     if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);}
2363:     else                       {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2364:     PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);
2365:     /* For each fine grid cell */
2366:     for (cell = cStart; cell < cEnd; ++cell) {
2367:       Vec                pointVec;
2368:       PetscScalar       *pV;
2369:       PetscSF            coarseCellSF = NULL;
2370:       const PetscSFNode *coarseCells;
2371:       PetscInt           numCoarseCells, cpdim, q, c, j;
2372:       PetscInt          *findices,   *cindices;
2373:       PetscInt           numFIndices, numCIndices;

2375:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2376:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2377:       /* Get points from the quadrature */
2378:       VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2379:       VecSetBlockSize(pointVec, dim);
2380:       VecGetArray(pointVec, &pV);
2381:       for (q = 0; q < Nq; ++q) {
2382:         const PetscReal xi0[3] = {-1., -1., -1.};

2384:         /* Transform point to real space */
2385:         CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2386:         for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2387:       }
2388:       VecRestoreArray(pointVec, &pV);
2389:       /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2390:       DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2391:       /* Update matrix */
2392:       PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2393:       if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2394:       VecGetArray(pointVec, &pV);
2395:       for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2396:         PetscReal pVReal[3];
2397:         const PetscReal xi0[3] = {-1., -1., -1.};


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

2406:         if (id == PETSCFE_CLASSID) {
2407:           PetscFE    fe = (PetscFE) obj;
2408:           PetscReal *B;

2410:           /* Evaluate coarse basis on contained point */
2411:           PetscFEGetDimension(fe, &cpdim);
2412:           PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2413:           /* Get elemMat entries by multiplying by weight */
2414:           for (i = 0; i < numFIndices; ++i) {
2415:             PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2416:             for (j = 0; j < cpdim; ++j) {
2417:               for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ;
2418:             }
2419:             /* Update interpolator */
2420:             if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2421:             if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2422:             MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2423:           }
2424:           PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2425:         } else {
2426:           cpdim = 1;
2427:           for (i = 0; i < numFIndices; ++i) {
2428:             PetscMemzero(elemMat, cpdim * sizeof(PetscScalar));
2429:             for (j = 0; j < cpdim; ++j) {
2430:               for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ;
2431:             }
2432:             /* Update interpolator */
2433:             if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2434:             PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);
2435:             if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2436:             MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2437:           }
2438:         }
2439:         DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2440:       }
2441:       VecRestoreArray(pointVec, &pV);
2442:       PetscSFDestroy(&coarseCellSF);
2443:       VecDestroy(&pointVec);
2444:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2445:     }
2446:   }
2447:   PetscFree3(v0,J,invJ);
2448:   PetscFree3(v0c,Jc,invJc);
2449:   PetscFree(elemMat);
2450:   MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
2451:   MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
2452:   return(0);
2453: }

2455: /*@
2456:   DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns

2458:   Input Parameters:
2459: + dmc  - The coarse mesh
2460: - dmf  - The fine mesh
2461: - user - The user context

2463:   Output Parameter:
2464: . sc   - The mapping

2466:   Level: developer

2468: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2469: @*/
2470: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
2471: {
2472:   PetscDS        prob;
2473:   PetscFE       *feRef;
2474:   PetscFV       *fvRef;
2475:   Vec            fv, cv;
2476:   IS             fis, cis;
2477:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
2478:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
2479:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;
2480:   PetscBool     *needAvg;

2484:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2485:   DMGetDimension(dmf, &dim);
2486:   DMGetSection(dmf, &fsection);
2487:   DMGetGlobalSection(dmf, &fglobalSection);
2488:   DMGetSection(dmc, &csection);
2489:   DMGetGlobalSection(dmc, &cglobalSection);
2490:   PetscSectionGetNumFields(fsection, &Nf);
2491:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
2492:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
2493:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2494:   DMGetDS(dmc, &prob);
2495:   PetscCalloc3(Nf,&feRef,Nf,&fvRef,Nf,&needAvg);
2496:   for (f = 0; f < Nf; ++f) {
2497:     PetscObject  obj;
2498:     PetscClassId id;
2499:     PetscInt     fNb = 0, Nc = 0;

2501:     PetscDSGetDiscretization(prob, f, &obj);
2502:     PetscObjectGetClassId(obj, &id);
2503:     if (id == PETSCFE_CLASSID) {
2504:       PetscFE    fe = (PetscFE) obj;
2505:       PetscSpace sp;
2506:       PetscInt   maxDegree;

2508:       PetscFERefine(fe, &feRef[f]);
2509:       PetscFEGetDimension(feRef[f], &fNb);
2510:       PetscFEGetNumComponents(fe, &Nc);
2511:       PetscFEGetBasisSpace(fe, &sp);
2512:       PetscSpaceGetDegree(sp, NULL, &maxDegree);
2513:       if (!maxDegree) needAvg[f] = PETSC_TRUE;
2514:     } else if (id == PETSCFV_CLASSID) {
2515:       PetscFV        fv = (PetscFV) obj;
2516:       PetscDualSpace Q;

2518:       PetscFVRefine(fv, &fvRef[f]);
2519:       PetscFVGetDualSpace(fvRef[f], &Q);
2520:       PetscDualSpaceGetDimension(Q, &fNb);
2521:       PetscFVGetNumComponents(fv, &Nc);
2522:       needAvg[f] = PETSC_TRUE;
2523:     }
2524:     fTotDim += fNb;
2525:   }
2526:   PetscDSGetTotalDimension(prob, &cTotDim);
2527:   PetscMalloc1(cTotDim,&cmap);
2528:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
2529:     PetscFE        feC;
2530:     PetscFV        fvC;
2531:     PetscDualSpace QF, QC;
2532:     PetscInt       order = -1, NcF, NcC, fpdim, cpdim;

2534:     if (feRef[field]) {
2535:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
2536:       PetscFEGetNumComponents(feC, &NcC);
2537:       PetscFEGetNumComponents(feRef[field], &NcF);
2538:       PetscFEGetDualSpace(feRef[field], &QF);
2539:       PetscDualSpaceGetOrder(QF, &order);
2540:       PetscDualSpaceGetDimension(QF, &fpdim);
2541:       PetscFEGetDualSpace(feC, &QC);
2542:       PetscDualSpaceGetDimension(QC, &cpdim);
2543:     } else {
2544:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
2545:       PetscFVGetNumComponents(fvC, &NcC);
2546:       PetscFVGetNumComponents(fvRef[field], &NcF);
2547:       PetscFVGetDualSpace(fvRef[field], &QF);
2548:       PetscDualSpaceGetDimension(QF, &fpdim);
2549:       PetscFVGetDualSpace(fvC, &QC);
2550:       PetscDualSpaceGetDimension(QC, &cpdim);
2551:     }
2552:     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);
2553:     for (c = 0; c < cpdim; ++c) {
2554:       PetscQuadrature  cfunc;
2555:       const PetscReal *cqpoints, *cqweights;
2556:       PetscInt         NqcC, NpC;
2557:       PetscBool        found = PETSC_FALSE;

2559:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
2560:       PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);
2561:       if (NqcC != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcC, NcC);
2562:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
2563:       for (f = 0; f < fpdim; ++f) {
2564:         PetscQuadrature  ffunc;
2565:         const PetscReal *fqpoints, *fqweights;
2566:         PetscReal        sum = 0.0;
2567:         PetscInt         NqcF, NpF;

2569:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
2570:         PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);
2571:         if (NqcF != NcF) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcF, NcF);
2572:         if (NpC != NpF) continue;
2573:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
2574:         if (sum > 1.0e-9) continue;
2575:         for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]);
2576:         if (sum < 1.0e-9) continue;
2577:         cmap[offsetC+c] = offsetF+f;
2578:         found = PETSC_TRUE;
2579:         break;
2580:       }
2581:       if (!found) {
2582:         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
2583:         if (fvRef[field] || (feRef[field] && order == 0)) {
2584:           cmap[offsetC+c] = offsetF+0;
2585:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
2586:       }
2587:     }
2588:     offsetC += cpdim;
2589:     offsetF += fpdim;
2590:   }
2591:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
2592:   PetscFree3(feRef,fvRef,needAvg);

2594:   DMGetGlobalVector(dmf, &fv);
2595:   DMGetGlobalVector(dmc, &cv);
2596:   VecGetOwnershipRange(cv, &startC, &endC);
2597:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
2598:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
2599:   PetscMalloc1(m,&cindices);
2600:   PetscMalloc1(m,&findices);
2601:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
2602:   for (c = cStart; c < cEnd; ++c) {
2603:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
2604:     for (d = 0; d < cTotDim; ++d) {
2605:       if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
2606:       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]]);
2607:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
2608:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
2609:     }
2610:   }
2611:   PetscFree(cmap);
2612:   PetscFree2(cellCIndices,cellFIndices);

2614:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
2615:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
2616:   VecScatterCreate(cv, cis, fv, fis, sc);
2617:   ISDestroy(&cis);
2618:   ISDestroy(&fis);
2619:   DMRestoreGlobalVector(dmf, &fv);
2620:   DMRestoreGlobalVector(dmc, &cv);
2621:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2622:   return(0);
2623: }

2625: /*@C
2626:   DMPlexGetCellFields - Retrieve the field values values for a chunk of cells

2628:   Input Parameters:
2629: + dm     - The DM
2630: . cellIS - The cells to include
2631: . locX   - A local vector with the solution fields
2632: . locX_t - A local vector with solution field time derivatives, or NULL
2633: - locA   - A local vector with auxiliary fields, or NULL

2635:   Output Parameters:
2636: + u   - The field coefficients
2637: . u_t - The fields derivative coefficients
2638: - a   - The auxiliary field coefficients

2640:   Level: developer

2642: .seealso: DMPlexGetFaceFields()
2643: @*/
2644: PetscErrorCode DMPlexGetCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
2645: {
2646:   DM              plex, plexA = NULL;
2647:   PetscSection    section, sectionAux;
2648:   PetscDS         prob;
2649:   const PetscInt *cells;
2650:   PetscInt        cStart, cEnd, numCells, totDim, totDimAux, c;
2651:   PetscErrorCode  ierr;

2661:   DMPlexConvertPlex(dm, &plex, PETSC_FALSE);
2662:   ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
2663:   DMGetSection(dm, &section);
2664:   DMGetCellDS(dm, cStart, &prob);
2665:   PetscDSGetTotalDimension(prob, &totDim);
2666:   if (locA) {
2667:     DM      dmAux;
2668:     PetscDS probAux;

2670:     VecGetDM(locA, &dmAux);
2671:     DMPlexConvertPlex(dmAux, &plexA, PETSC_FALSE);
2672:     DMGetSection(dmAux, &sectionAux);
2673:     DMGetDS(dmAux, &probAux);
2674:     PetscDSGetTotalDimension(probAux, &totDimAux);
2675:   }
2676:   numCells = cEnd - cStart;
2677:   DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u);
2678:   if (locX_t) {DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u_t);} else {*u_t = NULL;}
2679:   if (locA)   {DMGetWorkArray(dm, numCells*totDimAux, MPIU_SCALAR, a);} else {*a = NULL;}
2680:   for (c = cStart; c < cEnd; ++c) {
2681:     const PetscInt cell = cells ? cells[c] : c;
2682:     const PetscInt cind = c - cStart;
2683:     PetscScalar   *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
2684:     PetscInt       i;

2686:     DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x);
2687:     for (i = 0; i < totDim; ++i) ul[cind*totDim+i] = x[i];
2688:     DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x);
2689:     if (locX_t) {
2690:       DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x_t);
2691:       for (i = 0; i < totDim; ++i) ul_t[cind*totDim+i] = x_t[i];
2692:       DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x_t);
2693:     }
2694:     if (locA) {
2695:       PetscInt subcell;
2696:       DMPlexGetAuxiliaryPoint(plex, plexA, cell, &subcell);
2697:       DMPlexVecGetClosure(plexA, sectionAux, locA, subcell, NULL, &x);
2698:       for (i = 0; i < totDimAux; ++i) al[cind*totDimAux+i] = x[i];
2699:       DMPlexVecRestoreClosure(plexA, sectionAux, locA, subcell, NULL, &x);
2700:     }
2701:   }
2702:   DMDestroy(&plex);
2703:   if (locA) {DMDestroy(&plexA);}
2704:   ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
2705:   return(0);
2706: }

2708: /*@C
2709:   DMPlexRestoreCellFields - Restore the field values values for a chunk of cells

2711:   Input Parameters:
2712: + dm     - The DM
2713: . cellIS - The cells to include
2714: . locX   - A local vector with the solution fields
2715: . locX_t - A local vector with solution field time derivatives, or NULL
2716: - locA   - A local vector with auxiliary fields, or NULL

2718:   Output Parameters:
2719: + u   - The field coefficients
2720: . u_t - The fields derivative coefficients
2721: - a   - The auxiliary field coefficients

2723:   Level: developer

2725: .seealso: DMPlexGetFaceFields()
2726: @*/
2727: PetscErrorCode DMPlexRestoreCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
2728: {

2732:   DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u);
2733:   if (locX_t) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u_t);}
2734:   if (locA)   {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, a);}
2735:   return(0);
2736: }

2738: /*@C
2739:   DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces

2741:   Input Parameters:
2742: + dm     - The DM
2743: . fStart - The first face to include
2744: . fEnd   - The first face to exclude
2745: . locX   - A local vector with the solution fields
2746: . locX_t - A local vector with solution field time derivatives, or NULL
2747: . faceGeometry - A local vector with face geometry
2748: . cellGeometry - A local vector with cell geometry
2749: - locaGrad - A local vector with field gradients, or NULL

2751:   Output Parameters:
2752: + Nface - The number of faces with field values
2753: . uL - The field values at the left side of the face
2754: - uR - The field values at the right side of the face

2756:   Level: developer

2758: .seealso: DMPlexGetCellFields()
2759: @*/
2760: PetscErrorCode DMPlexGetFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR)
2761: {
2762:   DM                 dmFace, dmCell, dmGrad = NULL;
2763:   PetscSection       section;
2764:   PetscDS            prob;
2765:   DMLabel            ghostLabel;
2766:   const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
2767:   PetscBool         *isFE;
2768:   PetscInt           dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
2769:   PetscErrorCode     ierr;

2780:   DMGetDimension(dm, &dim);
2781:   DMGetDS(dm, &prob);
2782:   DMGetSection(dm, &section);
2783:   PetscDSGetNumFields(prob, &Nf);
2784:   PetscDSGetTotalComponents(prob, &Nc);
2785:   PetscMalloc1(Nf, &isFE);
2786:   for (f = 0; f < Nf; ++f) {
2787:     PetscObject  obj;
2788:     PetscClassId id;

2790:     PetscDSGetDiscretization(prob, f, &obj);
2791:     PetscObjectGetClassId(obj, &id);
2792:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
2793:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
2794:     else                            {isFE[f] = PETSC_FALSE;}
2795:   }
2796:   DMGetLabel(dm, "ghost", &ghostLabel);
2797:   VecGetArrayRead(locX, &x);
2798:   VecGetDM(faceGeometry, &dmFace);
2799:   VecGetArrayRead(faceGeometry, &facegeom);
2800:   VecGetDM(cellGeometry, &dmCell);
2801:   VecGetArrayRead(cellGeometry, &cellgeom);
2802:   if (locGrad) {
2803:     VecGetDM(locGrad, &dmGrad);
2804:     VecGetArrayRead(locGrad, &lgrad);
2805:   }
2806:   DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uL);
2807:   DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uR);
2808:   /* Right now just eat the extra work for FE (could make a cell loop) */
2809:   for (face = fStart, iface = 0; face < fEnd; ++face) {
2810:     const PetscInt        *cells;
2811:     PetscFVFaceGeom       *fg;
2812:     PetscFVCellGeom       *cgL, *cgR;
2813:     PetscScalar           *xL, *xR, *gL, *gR;
2814:     PetscScalar           *uLl = *uL, *uRl = *uR;
2815:     PetscInt               ghost, nsupp, nchild;

2817:     DMLabelGetValue(ghostLabel, face, &ghost);
2818:     DMPlexGetSupportSize(dm, face, &nsupp);
2819:     DMPlexGetTreeChildren(dm, face, &nchild, NULL);
2820:     if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
2821:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
2822:     DMPlexGetSupport(dm, face, &cells);
2823:     DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
2824:     DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
2825:     for (f = 0; f < Nf; ++f) {
2826:       PetscInt off;

2828:       PetscDSGetComponentOffset(prob, f, &off);
2829:       if (isFE[f]) {
2830:         const PetscInt *cone;
2831:         PetscInt        comp, coneSizeL, coneSizeR, faceLocL, faceLocR, ldof, rdof, d;

2833:         xL = xR = NULL;
2834:         PetscSectionGetFieldComponents(section, f, &comp);
2835:         DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
2836:         DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
2837:         DMPlexGetCone(dm, cells[0], &cone);
2838:         DMPlexGetConeSize(dm, cells[0], &coneSizeL);
2839:         for (faceLocL = 0; faceLocL < coneSizeL; ++faceLocL) if (cone[faceLocL] == face) break;
2840:         DMPlexGetCone(dm, cells[1], &cone);
2841:         DMPlexGetConeSize(dm, cells[1], &coneSizeR);
2842:         for (faceLocR = 0; faceLocR < coneSizeR; ++faceLocR) if (cone[faceLocR] == face) break;
2843:         if (faceLocL == coneSizeL && faceLocR == coneSizeR) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d or cell %d", face, cells[0], cells[1]);
2844:         /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
2845:         /* TODO: this is a hack that might not be right for nonconforming */
2846:         if (faceLocL < coneSizeL) {
2847:           EvaluateFaceFields(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
2848:           if (rdof == ldof && faceLocR < coneSizeR) {EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
2849:           else              {for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
2850:         }
2851:         else {
2852:           EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);
2853:           PetscSectionGetFieldComponents(section, f, &comp);
2854:           for(d = 0; d < comp; ++d) uLl[iface*Nc+off+d] = uRl[iface*Nc+off+d];
2855:         }
2856:         DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
2857:         DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
2858:       } else {
2859:         PetscFV  fv;
2860:         PetscInt numComp, c;

2862:         PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
2863:         PetscFVGetNumComponents(fv, &numComp);
2864:         DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
2865:         DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
2866:         if (dmGrad) {
2867:           PetscReal dxL[3], dxR[3];

2869:           DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
2870:           DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
2871:           DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
2872:           DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
2873:           for (c = 0; c < numComp; ++c) {
2874:             uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
2875:             uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
2876:           }
2877:         } else {
2878:           for (c = 0; c < numComp; ++c) {
2879:             uLl[iface*Nc+off+c] = xL[c];
2880:             uRl[iface*Nc+off+c] = xR[c];
2881:           }
2882:         }
2883:       }
2884:     }
2885:     ++iface;
2886:   }
2887:   *Nface = iface;
2888:   VecRestoreArrayRead(locX, &x);
2889:   VecRestoreArrayRead(faceGeometry, &facegeom);
2890:   VecRestoreArrayRead(cellGeometry, &cellgeom);
2891:   if (locGrad) {
2892:     VecRestoreArrayRead(locGrad, &lgrad);
2893:   }
2894:   PetscFree(isFE);
2895:   return(0);
2896: }

2898: /*@C
2899:   DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces

2901:   Input Parameters:
2902: + dm     - The DM
2903: . fStart - The first face to include
2904: . fEnd   - The first face to exclude
2905: . locX   - A local vector with the solution fields
2906: . locX_t - A local vector with solution field time derivatives, or NULL
2907: . faceGeometry - A local vector with face geometry
2908: . cellGeometry - A local vector with cell geometry
2909: - locaGrad - A local vector with field gradients, or NULL

2911:   Output Parameters:
2912: + Nface - The number of faces with field values
2913: . uL - The field values at the left side of the face
2914: - uR - The field values at the right side of the face

2916:   Level: developer

2918: .seealso: DMPlexGetFaceFields()
2919: @*/
2920: PetscErrorCode DMPlexRestoreFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR)
2921: {

2925:   DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uL);
2926:   DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uR);
2927:   return(0);
2928: }

2930: /*@C
2931:   DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces

2933:   Input Parameters:
2934: + dm     - The DM
2935: . fStart - The first face to include
2936: . fEnd   - The first face to exclude
2937: . faceGeometry - A local vector with face geometry
2938: - cellGeometry - A local vector with cell geometry

2940:   Output Parameters:
2941: + Nface - The number of faces with field values
2942: . fgeom - The extract the face centroid and normal
2943: - vol   - The cell volume

2945:   Level: developer

2947: .seealso: DMPlexGetCellFields()
2948: @*/
2949: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
2950: {
2951:   DM                 dmFace, dmCell;
2952:   DMLabel            ghostLabel;
2953:   const PetscScalar *facegeom, *cellgeom;
2954:   PetscInt           dim, numFaces = fEnd - fStart, iface, face;
2955:   PetscErrorCode     ierr;

2963:   DMGetDimension(dm, &dim);
2964:   DMGetLabel(dm, "ghost", &ghostLabel);
2965:   VecGetDM(faceGeometry, &dmFace);
2966:   VecGetArrayRead(faceGeometry, &facegeom);
2967:   VecGetDM(cellGeometry, &dmCell);
2968:   VecGetArrayRead(cellGeometry, &cellgeom);
2969:   PetscMalloc1(numFaces, fgeom);
2970:   DMGetWorkArray(dm, numFaces*2, MPIU_SCALAR, vol);
2971:   for (face = fStart, iface = 0; face < fEnd; ++face) {
2972:     const PetscInt        *cells;
2973:     PetscFVFaceGeom       *fg;
2974:     PetscFVCellGeom       *cgL, *cgR;
2975:     PetscFVFaceGeom       *fgeoml = *fgeom;
2976:     PetscReal             *voll   = *vol;
2977:     PetscInt               ghost, d, nchild, nsupp;

2979:     DMLabelGetValue(ghostLabel, face, &ghost);
2980:     DMPlexGetSupportSize(dm, face, &nsupp);
2981:     DMPlexGetTreeChildren(dm, face, &nchild, NULL);
2982:     if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
2983:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
2984:     DMPlexGetSupport(dm, face, &cells);
2985:     DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
2986:     DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
2987:     for (d = 0; d < dim; ++d) {
2988:       fgeoml[iface].centroid[d] = fg->centroid[d];
2989:       fgeoml[iface].normal[d]   = fg->normal[d];
2990:     }
2991:     voll[iface*2+0] = cgL->volume;
2992:     voll[iface*2+1] = cgR->volume;
2993:     ++iface;
2994:   }
2995:   *Nface = iface;
2996:   VecRestoreArrayRead(faceGeometry, &facegeom);
2997:   VecRestoreArrayRead(cellGeometry, &cellgeom);
2998:   return(0);
2999: }

3001: /*@C
3002:   DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces

3004:   Input Parameters:
3005: + dm     - The DM
3006: . fStart - The first face to include
3007: . fEnd   - The first face to exclude
3008: . faceGeometry - A local vector with face geometry
3009: - cellGeometry - A local vector with cell geometry

3011:   Output Parameters:
3012: + Nface - The number of faces with field values
3013: . fgeom - The extract the face centroid and normal
3014: - vol   - The cell volume

3016:   Level: developer

3018: .seealso: DMPlexGetFaceFields()
3019: @*/
3020: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
3021: {

3025:   PetscFree(*fgeom);
3026:   DMRestoreWorkArray(dm, 0, MPIU_REAL, vol);
3027:   return(0);
3028: }

3030: PetscErrorCode DMSNESGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3031: {
3032:   char            composeStr[33] = {0};
3033:   PetscObjectId   id;
3034:   PetscContainer  container;
3035:   PetscErrorCode  ierr;

3038:   PetscObjectGetId((PetscObject)quad,&id);
3039:   PetscSNPrintf(composeStr, 32, "DMSNESGetFEGeom_%x\n", id);
3040:   PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
3041:   if (container) {
3042:     PetscContainerGetPointer(container, (void **) geom);
3043:   } else {
3044:     DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
3045:     PetscContainerCreate(PETSC_COMM_SELF,&container);
3046:     PetscContainerSetPointer(container, (void *) *geom);
3047:     PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
3048:     PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
3049:     PetscContainerDestroy(&container);
3050:   }
3051:   return(0);
3052: }

3054: PetscErrorCode DMSNESRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3055: {
3057:   *geom = NULL;
3058:   return(0);
3059: }

3061: PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM dm, PetscSection section, IS cellIS, PetscReal t, Vec locX, Vec locX_t, Vec locF, void *user)
3062: {
3063:   DM_Plex         *mesh       = (DM_Plex *) dm->data;
3064:   const char      *name       = "Residual";
3065:   DM               dmAux      = NULL;
3066:   DMLabel          ghostLabel = NULL;
3067:   PetscDS          prob       = NULL;
3068:   PetscDS          probAux    = NULL;
3069:   PetscBool        useFEM     = PETSC_FALSE;
3070:   PetscBool        isImplicit = (locX_t || t == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
3071:   DMField          coordField = NULL;
3072:   Vec              locA;
3073:   PetscScalar     *u = NULL, *u_t, *a, *uL = NULL, *uR = NULL;
3074:   IS               chunkIS;
3075:   const PetscInt  *cells;
3076:   PetscInt         cStart, cEnd, numCells;
3077:   PetscInt         Nf, f, totDim, totDimAux, numChunks, cellChunkSize, chunk, fStart, fEnd;
3078:   PetscInt         maxDegree = PETSC_MAX_INT;
3079:   PetscQuadrature  affineQuad = NULL, *quads = NULL;
3080:   PetscFEGeom     *affineGeom = NULL, **geoms = NULL;
3081:   PetscErrorCode   ierr;

3084:   PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
3085:   /* FEM+FVM */
3086:   /* 1: Get sizes from dm and dmAux */
3087:   DMGetLabel(dm, "ghost", &ghostLabel);
3088:   DMGetDS(dm, &prob);
3089:   PetscDSGetNumFields(prob, &Nf);
3090:   PetscDSGetTotalDimension(prob, &totDim);
3091:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
3092:   if (locA) {
3093:     VecGetDM(locA, &dmAux);
3094:     DMGetDS(dmAux, &probAux);
3095:     PetscDSGetTotalDimension(probAux, &totDimAux);
3096:   }
3097:   /* 2: Get geometric data */
3098:   for (f = 0; f < Nf; ++f) {
3099:     PetscObject  obj;
3100:     PetscClassId id;
3101:     PetscBool    fimp;

3103:     PetscDSGetImplicit(prob, f, &fimp);
3104:     if (isImplicit != fimp) continue;
3105:     PetscDSGetDiscretization(prob, f, &obj);
3106:     PetscObjectGetClassId(obj, &id);
3107:     if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
3108:     if (id == PETSCFV_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Use of FVM with PCPATCH not yet implemented");
3109:   }
3110:   if (useFEM) {
3111:     DMGetCoordinateField(dm, &coordField);
3112:     DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
3113:     if (maxDegree <= 1) {
3114:       DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
3115:       if (affineQuad) {
3116:         DMSNESGetFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3117:       }
3118:     } else {
3119:       PetscCalloc2(Nf,&quads,Nf,&geoms);
3120:       for (f = 0; f < Nf; ++f) {
3121:         PetscObject  obj;
3122:         PetscClassId id;
3123:         PetscBool    fimp;

3125:         PetscDSGetImplicit(prob, f, &fimp);
3126:         if (isImplicit != fimp) continue;
3127:         PetscDSGetDiscretization(prob, f, &obj);
3128:         PetscObjectGetClassId(obj, &id);
3129:         if (id == PETSCFE_CLASSID) {
3130:           PetscFE fe = (PetscFE) obj;

3132:           PetscFEGetQuadrature(fe, &quads[f]);
3133:           PetscObjectReference((PetscObject)quads[f]);
3134:           DMSNESGetFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3135:         }
3136:       }
3137:     }
3138:   }
3139:   /* Loop over chunks */
3140:   ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3141:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
3142:   if (useFEM) {ISCreate(PETSC_COMM_SELF, &chunkIS);}
3143:   numCells      = cEnd - cStart;
3144:   numChunks     = 1;
3145:   cellChunkSize = numCells/numChunks;
3146:   numChunks     = PetscMin(1,numCells);
3147:   for (chunk = 0; chunk < numChunks; ++chunk) {
3148:     PetscScalar     *elemVec, *fluxL = NULL, *fluxR = NULL;
3149:     PetscReal       *vol = NULL;
3150:     PetscFVFaceGeom *fgeom = NULL;
3151:     PetscInt         cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c;
3152:     PetscInt         numFaces = 0;

3154:     /* Extract field coefficients */
3155:     if (useFEM) {
3156:       ISGetPointSubrange(chunkIS, cS, cE, cells);
3157:       DMPlexGetCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3158:       DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3159:       PetscMemzero(elemVec, numCells*totDim * sizeof(PetscScalar));
3160:     }
3161:     /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
3162:     /* Loop over fields */
3163:     for (f = 0; f < Nf; ++f) {
3164:       PetscObject  obj;
3165:       PetscClassId id;
3166:       PetscBool    fimp;
3167:       PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

3169:       PetscDSGetImplicit(prob, f, &fimp);
3170:       if (isImplicit != fimp) continue;
3171:       PetscDSGetDiscretization(prob, f, &obj);
3172:       PetscObjectGetClassId(obj, &id);
3173:       if (id == PETSCFE_CLASSID) {
3174:         PetscFE         fe = (PetscFE) obj;
3175:         PetscFEGeom    *geom = affineGeom ? affineGeom : geoms[f];
3176:         PetscFEGeom    *chunkGeom = NULL;
3177:         PetscQuadrature quad = affineQuad ? affineQuad : quads[f];
3178:         PetscInt        Nq, Nb;

3180:         PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3181:         PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
3182:         PetscFEGetDimension(fe, &Nb);
3183:         blockSize = Nb;
3184:         batchSize = numBlocks * blockSize;
3185:         PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3186:         numChunks = numCells / (numBatches*batchSize);
3187:         Ne        = numChunks*numBatches*batchSize;
3188:         Nr        = numCells % (numBatches*batchSize);
3189:         offset    = numCells - Nr;
3190:         /* Integrate FE residual to get elemVec (need fields at quadrature points) */
3191:         /*   For FV, I think we use a P0 basis and the cell coefficients (for subdivided cells, we can tweak the basis tabulation to be the indicator function) */
3192:         PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);
3193:         PetscFEIntegrateResidual(fe, prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);
3194:         PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);
3195:         PetscFEIntegrateResidual(fe, prob, f, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);
3196:         PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);
3197:       } else if (id == PETSCFV_CLASSID) {
3198:         PetscFV fv = (PetscFV) obj;

3200:         Ne = numFaces;
3201:         /* Riemann solve over faces (need fields at face centroids) */
3202:         /*   We need to evaluate FE fields at those coordinates */
3203:         PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
3204:       } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
3205:     }
3206:     /* Loop over domain */
3207:     if (useFEM) {
3208:       /* Add elemVec to locX */
3209:       for (c = cS; c < cE; ++c) {
3210:         const PetscInt cell = cells ? cells[c] : c;
3211:         const PetscInt cind = c - cStart;

3213:         if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cind*totDim]);}
3214:         if (ghostLabel) {
3215:           PetscInt ghostVal;

3217:           DMLabelGetValue(ghostLabel,cell,&ghostVal);
3218:           if (ghostVal > 0) continue;
3219:         }
3220:         DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);
3221:       }
3222:     }
3223:     /* Handle time derivative */
3224:     if (locX_t) {
3225:       PetscScalar *x_t, *fa;

3227:       VecGetArray(locF, &fa);
3228:       VecGetArray(locX_t, &x_t);
3229:       for (f = 0; f < Nf; ++f) {
3230:         PetscFV      fv;
3231:         PetscObject  obj;
3232:         PetscClassId id;
3233:         PetscInt     pdim, d;

3235:         PetscDSGetDiscretization(prob, f, &obj);
3236:         PetscObjectGetClassId(obj, &id);
3237:         if (id != PETSCFV_CLASSID) continue;
3238:         fv   = (PetscFV) obj;
3239:         PetscFVGetNumComponents(fv, &pdim);
3240:         for (c = cS; c < cE; ++c) {
3241:           const PetscInt cell = cells ? cells[c] : c;
3242:           PetscScalar   *u_t, *r;

3244:           if (ghostLabel) {
3245:             PetscInt ghostVal;

3247:             DMLabelGetValue(ghostLabel, cell, &ghostVal);
3248:             if (ghostVal > 0) continue;
3249:           }
3250:           DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
3251:           DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
3252:           for (d = 0; d < pdim; ++d) r[d] += u_t[d];
3253:         }
3254:       }
3255:       VecRestoreArray(locX_t, &x_t);
3256:       VecRestoreArray(locF, &fa);
3257:     }
3258:     if (useFEM) {
3259:       DMPlexRestoreCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3260:       DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3261:     }
3262:   }
3263:   if (useFEM) {ISDestroy(&chunkIS);}
3264:   ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
3265:   /* TODO Could include boundary residual here (see DMPlexComputeResidual_Internal) */
3266:   if (useFEM) {
3267:     if (maxDegree <= 1) {
3268:       DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3269:       PetscQuadratureDestroy(&affineQuad);
3270:     } else {
3271:       for (f = 0; f < Nf; ++f) {
3272:         DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3273:         PetscQuadratureDestroy(&quads[f]);
3274:       }
3275:       PetscFree2(quads,geoms);
3276:     }
3277:   }
3278:   PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
3279:   return(0);
3280: }

3282: /*
3283:   We always assemble JacP, and if the matrix is different from Jac and two different sets of point functions are provided, we also assemble Jac

3285:   X   - The local solution vector
3286:   X_t - The local solution time derviative vector, or NULL
3287: */
3288: PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM dm, PetscSection section, PetscSection globalSection, IS cellIS,
3289:                                                     PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP, void *ctx)
3290: {
3291:   DM_Plex         *mesh  = (DM_Plex *) dm->data;
3292:   const char      *name = "Jacobian", *nameP = "JacobianPre";
3293:   DM               dmAux = NULL;
3294:   PetscDS          prob,   probAux = NULL;
3295:   PetscSection     sectionAux = NULL;
3296:   Vec              A;
3297:   DMField          coordField;
3298:   PetscFEGeom     *cgeomFEM;
3299:   PetscQuadrature  qGeom = NULL;
3300:   Mat              J = Jac, JP = JacP;
3301:   PetscScalar     *work, *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL, *elemMatP = NULL, *elemMatD = NULL;
3302:   PetscBool        hasJac, hasPrec, hasDyn, assembleJac, isMatIS, isMatISP, *isFE, hasFV = PETSC_FALSE;
3303:   const PetscInt  *cells;
3304:   PetscInt         Nf, fieldI, fieldJ, maxDegree, numCells, cStart, cEnd, numChunks, chunkSize, chunk, totDim, totDimAux = 0, sz, wsz, off = 0, offCell = 0;
3305:   PetscErrorCode   ierr;

3308:   CHKMEMQ;
3309:   ISGetLocalSize(cellIS, &numCells);
3310:   ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3311:   PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
3312:   DMGetDS(dm, &prob);
3313:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3314:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3315:   if (dmAux) {
3316:     DMGetDefaultSection(dmAux, &sectionAux);
3317:     DMGetDS(dmAux, &probAux);
3318:   }
3319:   /* Get flags */
3320:   PetscDSGetNumFields(prob, &Nf);
3321:   DMGetWorkArray(dm, Nf, MPIU_BOOL, &isFE);
3322:   for (fieldI = 0; fieldI < Nf; ++fieldI) {
3323:     PetscObject  disc;
3324:     PetscClassId id;
3325:     PetscDSGetDiscretization(prob, fieldI, &disc);
3326:     PetscObjectGetClassId(disc, &id);
3327:     if (id == PETSCFE_CLASSID)      {isFE[fieldI] = PETSC_TRUE;}
3328:     else if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; isFE[fieldI] = PETSC_FALSE;}
3329:   }
3330:   PetscDSHasJacobian(prob, &hasJac);
3331:   PetscDSHasJacobianPreconditioner(prob, &hasPrec);
3332:   PetscDSHasDynamicJacobian(prob, &hasDyn);
3333:   assembleJac = hasJac && hasPrec && (Jac != JacP) ? PETSC_TRUE : PETSC_FALSE;
3334:   hasDyn      = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
3335:   PetscObjectTypeCompare((PetscObject) Jac,  MATIS, &isMatIS);
3336:   PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
3337:   /* Setup input data and temp arrays (should be DMGetWorkArray) */
3338:   if (isMatISP || isMatISP) {DMPlexGetSubdomainSection(dm, &globalSection);}
3339:   if (isMatIS)  {MatISGetLocalMat(Jac,  &J);}
3340:   if (isMatISP) {MatISGetLocalMat(JacP, &JP);}
3341:   if (hasFV)    {MatSetOption(JP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);} /* No allocated space for FV stuff, so ignore the zero entries */
3342:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3343:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3344:   PetscDSGetTotalDimension(prob, &totDim);
3345:   if (probAux) {PetscDSGetTotalDimension(probAux, &totDimAux);}
3346:   CHKMEMQ;
3347:   /* Compute batch sizes */
3348:   if (isFE[0]) {
3349:     PetscFE         fe;
3350:     PetscQuadrature q;
3351:     PetscInt        numQuadPoints, numBatches, batchSize, numBlocks, blockSize, Nb;

3353:     PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3354:     PetscFEGetQuadrature(fe, &q);
3355:     PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
3356:     PetscFEGetDimension(fe, &Nb);
3357:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3358:     blockSize = Nb*numQuadPoints;
3359:     batchSize = numBlocks  * blockSize;
3360:     chunkSize = numBatches * batchSize;
3361:     numChunks = numCells / chunkSize + numCells % chunkSize;
3362:     PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3363:   } else {
3364:     chunkSize = numCells;
3365:     numChunks = 1;
3366:   }
3367:   /* Get work space */
3368:   wsz  = (((X?1:0) + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize;
3369:   DMGetWorkArray(dm, wsz, MPIU_SCALAR, &work);
3370:   PetscMemzero(work, wsz * sizeof(PetscScalar));
3371:   off      = 0;
3372:   u        = X       ? (sz = chunkSize*totDim,        off += sz, work+off-sz) : NULL;
3373:   u_t      = X_t     ? (sz = chunkSize*totDim,        off += sz, work+off-sz) : NULL;
3374:   a        = dmAux   ? (sz = chunkSize*totDimAux,     off += sz, work+off-sz) : NULL;
3375:   elemMat  = hasJac  ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3376:   elemMatP = hasPrec ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3377:   elemMatD = hasDyn  ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3378:   if (off != wsz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error is workspace size %D should be %D", off, wsz);
3379:   /* Setup geometry */
3380:   DMGetCoordinateField(dm, &coordField);
3381:   DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);
3382:   if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, cellIS, &qGeom);}
3383:   if (!qGeom) {
3384:     PetscFE fe;

3386:     PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3387:     PetscFEGetQuadrature(fe, &qGeom);
3388:     PetscObjectReference((PetscObject) qGeom);
3389:   }
3390:   DMSNESGetFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
3391:   /* Compute volume integrals */
3392:   if (assembleJac) {MatZeroEntries(J);}
3393:   MatZeroEntries(JP);
3394:   for (chunk = 0; chunk < numChunks; ++chunk, offCell += chunkSize) {
3395:     const PetscInt   Ncell = PetscMin(chunkSize, numCells - offCell);
3396:     PetscInt         c;

3398:     /* Extract values */
3399:     for (c = 0; c < Ncell; ++c) {
3400:       const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
3401:       PetscScalar   *x = NULL,  *x_t = NULL;
3402:       PetscInt       i;

3404:       if (X) {
3405:         DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);
3406:         for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
3407:         DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);
3408:       }
3409:       if (X_t) {
3410:         DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);
3411:         for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
3412:         DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);
3413:       }
3414:       if (dmAux) {
3415:         DMPlexVecGetClosure(dmAux, sectionAux, A, cell, NULL, &x);
3416:         for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
3417:         DMPlexVecRestoreClosure(dmAux, sectionAux, A, cell, NULL, &x);
3418:       }
3419:     }
3420:     CHKMEMQ;
3421:     for (fieldI = 0; fieldI < Nf; ++fieldI) {
3422:       PetscFE fe;
3423:       PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3424:       for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
3425:         if (hasJac)  {PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN,     fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMat);}
3426:         if (hasPrec) {PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatP);}
3427:         if (hasDyn)  {PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatD);}
3428:       }
3429:       /* For finite volume, add the identity */
3430:       if (!isFE[fieldI]) {
3431:         PetscFV  fv;
3432:         PetscInt eOffset = 0, Nc, fc, foff;

3434:         PetscDSGetFieldOffset(prob, fieldI, &foff);
3435:         PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);
3436:         PetscFVGetNumComponents(fv, &Nc);
3437:         for (c = 0; c < chunkSize; ++c, eOffset += totDim*totDim) {
3438:           for (fc = 0; fc < Nc; ++fc) {
3439:             const PetscInt i = foff + fc;
3440:             if (hasJac)  {elemMat [eOffset+i*totDim+i] = 1.0;}
3441:             if (hasPrec) {elemMatP[eOffset+i*totDim+i] = 1.0;}
3442:           }
3443:         }
3444:       }
3445:     }
3446:     CHKMEMQ;
3447:     /*   Add contribution from X_t */
3448:     if (hasDyn) {for (c = 0; c < chunkSize*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];}
3449:     /* Insert values into matrix */
3450:     for (c = 0; c < Ncell; ++c) {
3451:       const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
3452:       if (mesh->printFEM > 1) {
3453:         if (hasJac)  {DMPrintCellMatrix(cell, name,  totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
3454:         if (hasPrec) {DMPrintCellMatrix(cell, nameP, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
3455:       }
3456:       if (assembleJac) {DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);}
3457:       DMPlexMatSetClosure(dm, section, globalSection, JP, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
3458:     }
3459:     CHKMEMQ;
3460:   }
3461:   /* Cleanup */
3462:   DMSNESRestoreFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
3463:   PetscQuadratureDestroy(&qGeom);
3464:   if (hasFV) {MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);}
3465:   DMRestoreWorkArray(dm, Nf, MPIU_BOOL, &isFE);
3466:   DMRestoreWorkArray(dm, ((1 + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize, MPIU_SCALAR, &work);
3467:   /* Compute boundary integrals */
3468:   /* DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, ctx); */
3469:   /* Assemble matrix */
3470:   if (assembleJac) {MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);}
3471:   MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
3472:   PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
3473:   CHKMEMQ;
3474:   return(0);
3475: }