Actual source code: plexfem.c

petsc-master 2016-12-09
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petscsf.h>

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

  9: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
 10: {
 11:   DM_Plex *mesh = (DM_Plex*) dm->data;

 16:   *scale = mesh->scale[unit];
 17:   return(0);
 18: }

 22: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
 23: {
 24:   DM_Plex *mesh = (DM_Plex*) dm->data;

 28:   mesh->scale[unit] = scale;
 29:   return(0);
 30: }

 32: PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
 33: {
 34:   switch (i) {
 35:   case 0:
 36:     switch (j) {
 37:     case 0: return 0;
 38:     case 1:
 39:       switch (k) {
 40:       case 0: return 0;
 41:       case 1: return 0;
 42:       case 2: return 1;
 43:       }
 44:     case 2:
 45:       switch (k) {
 46:       case 0: return 0;
 47:       case 1: return -1;
 48:       case 2: return 0;
 49:       }
 50:     }
 51:   case 1:
 52:     switch (j) {
 53:     case 0:
 54:       switch (k) {
 55:       case 0: return 0;
 56:       case 1: return 0;
 57:       case 2: return -1;
 58:       }
 59:     case 1: return 0;
 60:     case 2:
 61:       switch (k) {
 62:       case 0: return 1;
 63:       case 1: return 0;
 64:       case 2: return 0;
 65:       }
 66:     }
 67:   case 2:
 68:     switch (j) {
 69:     case 0:
 70:       switch (k) {
 71:       case 0: return 0;
 72:       case 1: return 1;
 73:       case 2: return 0;
 74:       }
 75:     case 1:
 76:       switch (k) {
 77:       case 0: return -1;
 78:       case 1: return 0;
 79:       case 2: return 0;
 80:       }
 81:     case 2: return 0;
 82:     }
 83:   }
 84:   return 0;
 85: }

 89: static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx)
 90: {
 91:   PetscInt *ctxInt  = (PetscInt *) ctx;
 92:   PetscInt  dim2    = ctxInt[0];
 93:   PetscInt  d       = ctxInt[1];
 94:   PetscInt  i, j, k = dim > 2 ? d - dim : d;

 97:   if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2);
 98:   for (i = 0; i < dim; i++) mode[i] = 0.;
 99:   if (d < dim) {
100:     mode[d] = 1.;
101:   } else {
102:     for (i = 0; i < dim; i++) {
103:       for (j = 0; j < dim; j++) {
104:         mode[j] += epsilon(i, j, k)*X[i];
105:       }
106:     }
107:   }
108:   return(0);
109: }

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

116:   Collective on DM

118:   Input Arguments:
119: . dm - the DM

121:   Output Argument:
122: . sp - the null space

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

126:   Level: advanced

128: .seealso: MatNullSpaceCreate()
129: @*/
130: PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
131: {
132:   MPI_Comm       comm;
133:   Vec            mode[6];
134:   PetscSection   section, globalSection;
135:   PetscInt       dim, dimEmbed, n, m, d, i, j;

139:   PetscObjectGetComm((PetscObject)dm,&comm);
140:   DMGetDimension(dm, &dim);
141:   DMGetCoordinateDim(dm, &dimEmbed);
142:   if (dim == 1) {
143:     MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
144:     return(0);
145:   }
146:   DMGetDefaultSection(dm, &section);
147:   DMGetDefaultGlobalSection(dm, &globalSection);
148:   PetscSectionGetConstrainedStorageSize(globalSection, &n);
149:   m    = (dim*(dim+1))/2;
150:   VecCreate(comm, &mode[0]);
151:   VecSetSizes(mode[0], n, PETSC_DETERMINE);
152:   VecSetUp(mode[0]);
153:   for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
154:   for (d = 0; d < m; d++) {
155:     PetscInt         ctx[2];
156:     PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
157:     void            *voidctx = (void *) (&ctx[0]);

159:     ctx[0] = dimEmbed;
160:     ctx[1] = d;
161:     DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);
162:   }
163:   for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
164:   /* Orthonormalize system */
165:   for (i = dim; i < m; ++i) {
166:     PetscScalar dots[6];

168:     VecMDot(mode[i], i, mode, dots);
169:     for (j = 0; j < i; ++j) dots[j] *= -1.0;
170:     VecMAXPY(mode[i], i, dots, mode);
171:     VecNormalize(mode[i], NULL);
172:   }
173:   MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
174:   for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
175:   return(0);
176: }

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

189:   Input Parameters:
190: + dm - the DMPlex object
191: - height - the maximum projection height >= 0

193:   Level: advanced

195: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
196: @*/
197: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
198: {
199:   DM_Plex *plex = (DM_Plex *) dm->data;

203:   plex->maxProjectionHeight = height;
204:   return(0);
205: }

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

213:   Input Parameters:
214: . dm - the DMPlex object

216:   Output Parameters:
217: . height - the maximum projection height

219:   Level: intermediate

221: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
222: @*/
223: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
224: {
225:   DM_Plex *plex = (DM_Plex *) dm->data;

229:   *height = plex->maxProjectionHeight;
230:   return(0);
231: }

235: static PetscErrorCode DMPlexInsertBoundaryValues_FEM_Internal(DM dm, PetscReal time, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX)
236: {
237:   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
238:   void            **ctxs;
239:   PetscInt          numFields;
240:   PetscErrorCode    ierr;

243:   DMGetNumFields(dm, &numFields);
244:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
245:   funcs[field] = func;
246:   ctxs[field]  = ctx;
247:   DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);
248:   PetscFree2(funcs,ctxs);
249:   return(0);
250: }

254: static PetscErrorCode DMPlexInsertBoundaryValues_FEM_AuxField_Internal(DM dm, PetscReal time, Vec locU, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[],
255:                                                                        void (*func)(PetscInt, PetscInt, PetscInt,
256:                                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
257:                                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
258:                                                                                     PetscReal, const PetscReal[], PetscScalar[]), void *ctx, Vec locX)
259: {
260:   void (**funcs)(PetscInt, PetscInt, PetscInt,
261:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
262:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
263:                  PetscReal, const PetscReal[], PetscScalar[]);
264:   void            **ctxs;
265:   PetscInt          numFields;
266:   PetscErrorCode    ierr;

269:   DMGetNumFields(dm, &numFields);
270:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
271:   funcs[field] = func;
272:   ctxs[field]  = ctx;
273:   DMProjectFieldLabelLocal(dm, time, label, numids, ids, locU, funcs, INSERT_BC_VALUES, locX);
274:   PetscFree2(funcs,ctxs);
275:   return(0);
276: }

280: /* This ignores numcomps/comps */
281: static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad,
282:                                                               PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
283: {
284:   PetscDS            prob;
285:   PetscSF            sf;
286:   DM                 dmFace, dmCell, dmGrad;
287:   const PetscScalar *facegeom, *cellgeom = NULL, *grad;
288:   const PetscInt    *leaves;
289:   PetscScalar       *x, *fx;
290:   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
291:   PetscErrorCode     ierr, ierru = 0;

294:   DMGetPointSF(dm, &sf);
295:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
296:   nleaves = PetscMax(0, nleaves);
297:   DMGetDimension(dm, &dim);
298:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
299:   DMGetDS(dm, &prob);
300:   VecGetDM(faceGeometry, &dmFace);
301:   VecGetArrayRead(faceGeometry, &facegeom);
302:   if (cellGeometry) {
303:     VecGetDM(cellGeometry, &dmCell);
304:     VecGetArrayRead(cellGeometry, &cellgeom);
305:   }
306:   if (Grad) {
307:     PetscFV fv;

309:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
310:     VecGetDM(Grad, &dmGrad);
311:     VecGetArrayRead(Grad, &grad);
312:     PetscFVGetNumComponents(fv, &pdim);
313:     DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
314:   }
315:   VecGetArray(locX, &x);
316:   for (i = 0; i < numids; ++i) {
317:     IS              faceIS;
318:     const PetscInt *faces;
319:     PetscInt        numFaces, f;

321:     DMLabelGetStratumIS(label, ids[i], &faceIS);
322:     if (!faceIS) continue; /* No points with that id on this process */
323:     ISGetLocalSize(faceIS, &numFaces);
324:     ISGetIndices(faceIS, &faces);
325:     for (f = 0; f < numFaces; ++f) {
326:       const PetscInt         face = faces[f], *cells;
327:       PetscFVFaceGeom        *fg;

329:       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
330:       PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
331:       if (loc >= 0) continue;
332:       DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
333:       DMPlexGetSupport(dm, face, &cells);
334:       if (Grad) {
335:         PetscFVCellGeom       *cg;
336:         PetscScalar           *cx, *cgrad;
337:         PetscScalar           *xG;
338:         PetscReal              dx[3];
339:         PetscInt               d;

341:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
342:         DMPlexPointLocalRead(dm, cells[0], x, &cx);
343:         DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
344:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
345:         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
346:         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
347:         ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
348:         if (ierru) {
349:           ISRestoreIndices(faceIS, &faces);
350:           ISDestroy(&faceIS);
351:           goto cleanup;
352:         }
353:       } else {
354:         PetscScalar       *xI;
355:         PetscScalar       *xG;

357:         DMPlexPointLocalRead(dm, cells[0], x, &xI);
358:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
359:         ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
360:         if (ierru) {
361:           ISRestoreIndices(faceIS, &faces);
362:           ISDestroy(&faceIS);
363:           goto cleanup;
364:         }
365:       }
366:     }
367:     ISRestoreIndices(faceIS, &faces);
368:     ISDestroy(&faceIS);
369:   }
370:   cleanup:
371:   VecRestoreArray(locX, &x);
372:   if (Grad) {
373:     DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);
374:     VecRestoreArrayRead(Grad, &grad);
375:   }
376:   if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
377:   VecRestoreArrayRead(faceGeometry, &facegeom);
378:   CHKERRQ(ierru);
379:   return(0);
380: }

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

387:   Input Parameters:
388: + dm - The DM
389: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
390: . time - The time
391: . faceGeomFVM - Face geometry data for FV discretizations
392: . cellGeomFVM - Cell geometry data for FV discretizations
393: - gradFVM - Gradient reconstruction data for FV discretizations

395:   Output Parameters:
396: . locX - Solution updated with boundary values

398:   Level: developer

400: .seealso: DMProjectFunctionLabelLocal()
401: @*/
402: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
403: {
404:   PetscInt       numBd, b;

413:   PetscDSGetNumBoundary(dm->prob, &numBd);
414:   for (b = 0; b < numBd; ++b) {
415:     DMBoundaryConditionType type;
416:     const char             *labelname;
417:     DMLabel                 label;
418:     PetscInt                field;
419:     PetscObject             obj;
420:     PetscClassId            id;
421:     void                  (*func)();
422:     PetscInt                numids;
423:     const PetscInt         *ids;
424:     void                   *ctx;

426:     DMGetBoundary(dm, b, &type, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);
427:     if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
428:     DMGetLabel(dm, labelname, &label);
429:     DMGetField(dm, field, &obj);
430:     PetscObjectGetClassId(obj, &id);
431:     if (id == PETSCFE_CLASSID) {
432:       switch (type) {
433:         /* for FEM, there is no insertion to be done for non-essential boundary conditions */
434:       case DM_BC_ESSENTIAL:
435:         DMPlexLabelAddCells(dm,label);
436:         DMPlexInsertBoundaryValues_FEM_Internal(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
437:         DMPlexLabelClearCells(dm,label);
438:         break;
439:       case DM_BC_ESSENTIAL_FIELD:
440:         DMPlexLabelAddCells(dm,label);
441:         DMPlexInsertBoundaryValues_FEM_AuxField_Internal(dm, time, locX, field, label, numids, ids, (void (*)(PetscInt, PetscInt, PetscInt,
442:                                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
443:                                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
444:                                                                                     PetscReal, const PetscReal[], PetscScalar[])) func, ctx, locX);
445:         DMPlexLabelClearCells(dm,label);
446:         break;
447:       default: break;
448:       }
449:     } else if (id == PETSCFV_CLASSID) {
450:       if (!faceGeomFVM) continue;
451:       DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM,
452:                                                      field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
453:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
454:   }
455:   return(0);
456: }

460: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
461: {
462:   const PetscInt   debug = 0;
463:   PetscSection     section;
464:   PetscQuadrature  quad;
465:   Vec              localX;
466:   PetscScalar     *funcVal, *interpolant;
467:   PetscReal       *coords, *v0, *J, *invJ, detJ;
468:   PetscReal        localDiff = 0.0;
469:   const PetscReal *quadPoints, *quadWeights;
470:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
471:   PetscErrorCode   ierr;

474:   DMGetDimension(dm, &dim);
475:   DMGetDefaultSection(dm, &section);
476:   PetscSectionGetNumFields(section, &numFields);
477:   DMGetLocalVector(dm, &localX);
478:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
479:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
480:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
481:   for (field = 0; field < numFields; ++field) {
482:     PetscObject  obj;
483:     PetscClassId id;
484:     PetscInt     Nc;

486:     DMGetField(dm, field, &obj);
487:     PetscObjectGetClassId(obj, &id);
488:     if (id == PETSCFE_CLASSID) {
489:       PetscFE fe = (PetscFE) obj;

491:       PetscFEGetQuadrature(fe, &quad);
492:       PetscFEGetNumComponents(fe, &Nc);
493:     } else if (id == PETSCFV_CLASSID) {
494:       PetscFV fv = (PetscFV) obj;

496:       PetscFVGetQuadrature(fv, &quad);
497:       PetscFVGetNumComponents(fv, &Nc);
498:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
499:     numComponents += Nc;
500:   }
501:   PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
502:   PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
503:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
504:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
505:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
506:   for (c = cStart; c < cEnd; ++c) {
507:     PetscScalar *x = NULL;
508:     PetscReal    elemDiff = 0.0;

510:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
511:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
512:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

514:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
515:       PetscObject  obj;
516:       PetscClassId id;
517:       void * const ctx = ctxs ? ctxs[field] : NULL;
518:       PetscInt     Nb, Nc, q, fc;

520:       DMGetField(dm, field, &obj);
521:       PetscObjectGetClassId(obj, &id);
522:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
523:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
524:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
525:       if (debug) {
526:         char title[1024];
527:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
528:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
529:       }
530:       for (q = 0; q < numQuadPoints; ++q) {
531:         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
532:         (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);
533:         if (ierr) {
534:           PetscErrorCode ierr2;
535:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
536:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
537:           ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2);
538: 
539:         }
540:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
541:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
542:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
543:         for (fc = 0; fc < Nc; ++fc) {
544:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);}
545:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
546:         }
547:       }
548:       fieldOffset += Nb*Nc;
549:     }
550:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
551:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
552:     localDiff += elemDiff;
553:   }
554:   PetscFree6(funcVal,interpolant,coords,v0,J,invJ);
555:   DMRestoreLocalVector(dm, &localX);
556:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
557:   *diff = PetscSqrtReal(*diff);
558:   return(0);
559: }

563: 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)
564: {
565:   const PetscInt  debug = 0;
566:   PetscSection    section;
567:   PetscQuadrature quad;
568:   Vec             localX;
569:   PetscScalar    *funcVal, *interpolantVec;
570:   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
571:   PetscReal       localDiff = 0.0;
572:   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp;
573:   PetscErrorCode  ierr;

576:   DMGetDimension(dm, &dim);
577:   DMGetDefaultSection(dm, &section);
578:   PetscSectionGetNumFields(section, &numFields);
579:   DMGetLocalVector(dm, &localX);
580:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
581:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
582:   for (field = 0; field < numFields; ++field) {
583:     PetscFE  fe;
584:     PetscInt Nc;

586:     DMGetField(dm, field, (PetscObject *) &fe);
587:     PetscFEGetQuadrature(fe, &quad);
588:     PetscFEGetNumComponents(fe, &Nc);
589:     numComponents += Nc;
590:   }
591:   /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
592:   PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
593:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
594:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
595:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
596:   for (c = cStart; c < cEnd; ++c) {
597:     PetscScalar *x = NULL;
598:     PetscReal    elemDiff = 0.0;

600:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
601:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
602:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

604:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
605:       PetscFE          fe;
606:       void * const     ctx = ctxs ? ctxs[field] : NULL;
607:       const PetscReal *quadPoints, *quadWeights;
608:       PetscReal       *basisDer;
609:       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;

611:       DMGetField(dm, field, (PetscObject *) &fe);
612:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
613:       PetscFEGetDimension(fe, &Nb);
614:       PetscFEGetNumComponents(fe, &Ncomp);
615:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
616:       if (debug) {
617:         char title[1024];
618:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
619:         DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);
620:       }
621:       for (q = 0; q < numQuadPoints; ++q) {
622:         for (d = 0; d < dim; d++) {
623:           coords[d] = v0[d];
624:           for (e = 0; e < dim; e++) {
625:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
626:           }
627:         }
628:         (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);
629:         if (ierr) {
630:           PetscErrorCode ierr2;
631:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
632:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
633:           ierr2 = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr2);
634: 
635:         }
636:         for (fc = 0; fc < Ncomp; ++fc) {
637:           PetscScalar interpolant = 0.0;

639:           for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
640:           for (f = 0; f < Nb; ++f) {
641:             const PetscInt fidx = f*Ncomp+fc;

643:             for (d = 0; d < dim; ++d) {
644:               realSpaceDer[d] = 0.0;
645:               for (g = 0; g < dim; ++g) {
646:                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
647:               }
648:               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
649:             }
650:           }
651:           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
652:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
653:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
654:         }
655:       }
656:       comp        += Ncomp;
657:       fieldOffset += Nb*Ncomp;
658:     }
659:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
660:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
661:     localDiff += elemDiff;
662:   }
663:   PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
664:   DMRestoreLocalVector(dm, &localX);
665:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
666:   *diff = PetscSqrtReal(*diff);
667:   return(0);
668: }

672: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
673: {
674:   const PetscInt   debug = 0;
675:   PetscSection     section;
676:   PetscQuadrature  quad;
677:   Vec              localX;
678:   PetscScalar     *funcVal, *interpolant;
679:   PetscReal       *coords, *v0, *J, *invJ, detJ;
680:   PetscReal       *localDiff;
681:   const PetscReal *quadPoints, *quadWeights;
682:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
683:   PetscErrorCode   ierr;

686:   DMGetDimension(dm, &dim);
687:   DMGetDefaultSection(dm, &section);
688:   PetscSectionGetNumFields(section, &numFields);
689:   DMGetLocalVector(dm, &localX);
690:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
691:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
692:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
693:   for (field = 0; field < numFields; ++field) {
694:     PetscObject  obj;
695:     PetscClassId id;
696:     PetscInt     Nc;

698:     DMGetField(dm, field, &obj);
699:     PetscObjectGetClassId(obj, &id);
700:     if (id == PETSCFE_CLASSID) {
701:       PetscFE fe = (PetscFE) obj;

703:       PetscFEGetQuadrature(fe, &quad);
704:       PetscFEGetNumComponents(fe, &Nc);
705:     } else if (id == PETSCFV_CLASSID) {
706:       PetscFV fv = (PetscFV) obj;

708:       PetscFVGetQuadrature(fv, &quad);
709:       PetscFVGetNumComponents(fv, &Nc);
710:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
711:     numComponents += Nc;
712:   }
713:   PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
714:   PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
715:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
716:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
717:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
718:   for (c = cStart; c < cEnd; ++c) {
719:     PetscScalar *x = NULL;

721:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
722:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
723:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

725:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
726:       PetscObject  obj;
727:       PetscClassId id;
728:       void * const ctx = ctxs ? ctxs[field] : NULL;
729:       PetscInt     Nb, Nc, q, fc;

731:       PetscReal       elemDiff = 0.0;

733:       DMGetField(dm, field, &obj);
734:       PetscObjectGetClassId(obj, &id);
735:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
736:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
737:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
738:       if (debug) {
739:         char title[1024];
740:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
741:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
742:       }
743:       for (q = 0; q < numQuadPoints; ++q) {
744:         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
745:         (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);
746:         if (ierr) {
747:           PetscErrorCode ierr2;
748:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
749:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
750:           ierr2 = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2);
751: 
752:         }
753:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
754:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
755:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
756:         for (fc = 0; fc < Nc; ++fc) {
757:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d point %g %g %g diff %g\n", c, field, dim > 0 ? coords[0] : 0., dim > 1 ? coords[1] : 0., dim > 2 ? coords[2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);}
758:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
759:         }
760:       }
761:       fieldOffset += Nb*Nc;
762:       localDiff[field] += elemDiff;
763:     }
764:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
765:   }
766:   DMRestoreLocalVector(dm, &localX);
767:   MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
768:   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
769:   PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);
770:   return(0);
771: }

775: /*@C
776:   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.

778:   Input Parameters:
779: + dm    - The DM
780: . time  - The time
781: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
782: . ctxs  - Optional array of contexts to pass to each function, or NULL.
783: - X     - The coefficient vector u_h

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

788:   Level: developer

790: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
791: @*/
792: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
793: {
794:   PetscSection     section;
795:   PetscQuadrature  quad;
796:   Vec              localX;
797:   PetscScalar     *funcVal, *interpolant;
798:   PetscReal       *coords, *v0, *J, *invJ, detJ;
799:   const PetscReal *quadPoints, *quadWeights;
800:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
801:   PetscErrorCode   ierr;

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

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

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

827:       PetscFVGetQuadrature(fv, &quad);
828:       PetscFVGetNumComponents(fv, &Nc);
829:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
830:     numComponents += Nc;
831:   }
832:   PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
833:   PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
834:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
835:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
836:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
837:   for (c = cStart; c < cEnd; ++c) {
838:     PetscScalar *x = NULL;
839:     PetscScalar  elemDiff = 0.0;

841:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
842:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
843:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

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

851:       DMGetField(dm, field, &obj);
852:       PetscObjectGetClassId(obj, &id);
853:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
854:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
855:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
856:       if (funcs[field]) {
857:         for (q = 0; q < numQuadPoints; ++q) {
858:           CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
859:           (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);
860:           if (ierr) {
861:             PetscErrorCode ierr2;
862:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
863:             ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2);
864:             ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
865: 
866:           }
867:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
868:           else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
869:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
870:           for (fc = 0; fc < Nc; ++fc) {
871:             elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
872:           }
873:         }
874:       }
875:       fieldOffset += Nb*Nc;
876:     }
877:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
878:     VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
879:   }
880:   PetscFree6(funcVal,interpolant,coords,v0,J,invJ);
881:   DMRestoreLocalVector(dm, &localX);
882:   VecSqrtAbs(D);
883:   return(0);
884: }

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

891:   Input Parameters:
892: + dm - The mesh
893: . X  - Local input vector
894: - user - The user context

896:   Output Parameter:
897: . integral - Local integral for each field

899:   Level: developer

901: .seealso: DMPlexComputeResidualFEM()
902: @*/
903: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
904: {
905:   DM_Plex           *mesh  = (DM_Plex *) dm->data;
906:   DM                 dmAux, dmGrad;
907:   Vec                localX, A, cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
908:   PetscDS            prob, probAux = NULL;
909:   PetscSection       section, sectionAux;
910:   PetscFV            fvm = NULL;
911:   PetscFECellGeom   *cgeomFEM;
912:   PetscFVFaceGeom   *fgeomFVM;
913:   PetscFVCellGeom   *cgeomFVM;
914:   PetscScalar       *u, *a = NULL;
915:   const PetscScalar *lgrad;
916:   PetscReal         *lintegral;
917:   PetscInt          *uOff, *uOff_x, *aOff = NULL;
918:   PetscInt           dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, cEndInterior, c;
919:   PetscInt           totDim, totDimAux;
920:   PetscBool          useFVM = PETSC_FALSE;
921:   PetscErrorCode     ierr;

924:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
925:   DMGetLocalVector(dm, &localX);
926:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);
927:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
928:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
929:   DMGetDimension(dm, &dim);
930:   DMGetDefaultSection(dm, &section);
931:   DMGetDS(dm, &prob);
932:   PetscDSGetTotalDimension(prob, &totDim);
933:   PetscDSGetComponentOffsets(prob, &uOff);
934:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
935:   PetscSectionGetNumFields(section, &Nf);
936:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
937:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
938:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
939:   numCells = cEnd - cStart;
940:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
941:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
942:   if (dmAux) {
943:     DMGetDS(dmAux, &probAux);
944:     PetscDSGetNumFields(probAux, &NfAux);
945:     DMGetDefaultSection(dmAux, &sectionAux);
946:     PetscDSGetTotalDimension(probAux, &totDimAux);
947:     PetscDSGetComponentOffsets(probAux, &aOff);
948:     PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);
949:   }
950:   for (f = 0; f < Nf; ++f) {
951:     PetscObject  obj;
952:     PetscClassId id;

954:     PetscDSGetDiscretization(prob, f, &obj);
955:     PetscObjectGetClassId(obj, &id);
956:     if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
957:   }
958:   if (useFVM) {
959:     Vec       grad;
960:     PetscInt  fStart, fEnd;
961:     PetscBool compGrad;

963:     PetscFVGetComputeGradients(fvm, &compGrad);
964:     PetscFVSetComputeGradients(fvm, PETSC_TRUE);
965:     DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
966:     DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);
967:     PetscFVSetComputeGradients(fvm, compGrad);
968:     VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
969:     VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
970:     /* Reconstruct and limit cell gradients */
971:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
972:     DMGetGlobalVector(dmGrad, &grad);
973:     DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);
974:     /* Communicate gradient values */
975:     DMGetLocalVector(dmGrad, &locGrad);
976:     DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
977:     DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
978:     DMRestoreGlobalVector(dmGrad, &grad);
979:     /* Handle non-essential (e.g. outflow) boundary values */
980:     DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
981:     VecGetArrayRead(locGrad, &lgrad);
982:   }
983:   PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);
984:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
985:   for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
986:   for (c = cStart; c < cEnd; ++c) {
987:     PetscScalar *x = NULL;
988:     PetscInt     i;

990:     DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);
991:     if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c);
992:     DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
993:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
994:     DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
995:     if (dmAux) {
996:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
997:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
998:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
999:     }
1000:   }
1001:   for (f = 0; f < Nf; ++f) {
1002:     PetscObject  obj;
1003:     PetscClassId id;
1004:     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1006:     PetscDSGetDiscretization(prob, f, &obj);
1007:     PetscObjectGetClassId(obj, &id);
1008:     if (id == PETSCFE_CLASSID) {
1009:       PetscFE         fe = (PetscFE) obj;
1010:       PetscQuadrature q;
1011:       PetscInt        Nq, Nb;

1013:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1014:       PetscFEGetQuadrature(fe, &q);
1015:       PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1016:       PetscFEGetDimension(fe, &Nb);
1017:       blockSize = Nb*Nq;
1018:       batchSize = numBlocks * blockSize;
1019:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1020:       numChunks = numCells / (numBatches*batchSize);
1021:       Ne        = numChunks*numBatches*batchSize;
1022:       Nr        = numCells % (numBatches*batchSize);
1023:       offset    = numCells - Nr;
1024:       PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);
1025:       PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);
1026:     } else if (id == PETSCFV_CLASSID) {
1027:       /* PetscFV  fv = (PetscFV) obj; */
1028:       PetscInt       foff;
1029:       PetscPointFunc obj_func;
1030:       PetscScalar    lint;

1032:       PetscDSGetObjective(prob, f, &obj_func);
1033:       PetscDSGetFieldOffset(prob, f, &foff);
1034:       if (obj_func) {
1035:         for (c = 0; c < numCells; ++c) {
1036:           PetscScalar *u_x;

1038:           DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1039:           obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, a, NULL, NULL, 0.0, cgeomFVM[c].centroid, &lint);
1040:           lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1041:         }
1042:       }
1043:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1044:   }
1045:   if (useFVM) {
1046:     VecRestoreArrayRead(locGrad, &lgrad);
1047:     VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1048:     VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1049:     DMRestoreLocalVector(dmGrad, &locGrad);
1050:     VecDestroy(&faceGeometryFVM);
1051:     VecDestroy(&cellGeometryFVM);
1052:     DMDestroy(&dmGrad);
1053:   }
1054:   if (dmAux) {PetscFree(a);}
1055:   if (mesh->printFEM) {
1056:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1057:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);}
1058:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1059:   }
1060:   DMRestoreLocalVector(dm, &localX);
1061:   MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1062:   PetscFree3(lintegral,u,cgeomFEM);
1063:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1064:   return(0);
1065: }

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

1072:   Input Parameters:
1073: + dmf  - The fine mesh
1074: . dmc  - The coarse mesh
1075: - user - The user context

1077:   Output Parameter:
1078: . In  - The interpolation matrix

1080:   Level: developer

1082: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1083: @*/
1084: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1085: {
1086:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1087:   const char       *name  = "Interpolator";
1088:   PetscDS           prob;
1089:   PetscFE          *feRef;
1090:   PetscFV          *fvRef;
1091:   PetscSection      fsection, fglobalSection;
1092:   PetscSection      csection, cglobalSection;
1093:   PetscScalar      *elemMat;
1094:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1095:   PetscInt          cTotDim, rTotDim = 0;
1096:   PetscErrorCode    ierr;

1099:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1100:   DMGetDimension(dmf, &dim);
1101:   DMGetDefaultSection(dmf, &fsection);
1102:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1103:   DMGetDefaultSection(dmc, &csection);
1104:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1105:   PetscSectionGetNumFields(fsection, &Nf);
1106:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1107:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1108:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1109:   DMGetDS(dmf, &prob);
1110:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1111:   for (f = 0; f < Nf; ++f) {
1112:     PetscObject  obj;
1113:     PetscClassId id;
1114:     PetscInt     rNb = 0, Nc = 0;

1116:     PetscDSGetDiscretization(prob, f, &obj);
1117:     PetscObjectGetClassId(obj, &id);
1118:     if (id == PETSCFE_CLASSID) {
1119:       PetscFE fe = (PetscFE) obj;

1121:       PetscFERefine(fe, &feRef[f]);
1122:       PetscFEGetDimension(feRef[f], &rNb);
1123:       PetscFEGetNumComponents(fe, &Nc);
1124:     } else if (id == PETSCFV_CLASSID) {
1125:       PetscFV        fv = (PetscFV) obj;
1126:       PetscDualSpace Q;

1128:       PetscFVRefine(fv, &fvRef[f]);
1129:       PetscFVGetDualSpace(fvRef[f], &Q);
1130:       PetscDualSpaceGetDimension(Q, &rNb);
1131:       PetscFVGetNumComponents(fv, &Nc);
1132:     }
1133:     rTotDim += rNb*Nc;
1134:   }
1135:   PetscDSGetTotalDimension(prob, &cTotDim);
1136:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1137:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1138:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1139:     PetscDualSpace   Qref;
1140:     PetscQuadrature  f;
1141:     const PetscReal *qpoints, *qweights;
1142:     PetscReal       *points;
1143:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1145:     /* Compose points from all dual basis functionals */
1146:     if (feRef[fieldI]) {
1147:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
1148:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
1149:     } else {
1150:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1151:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1152:     }
1153:     PetscDualSpaceGetDimension(Qref, &fpdim);
1154:     for (i = 0; i < fpdim; ++i) {
1155:       PetscDualSpaceGetFunctional(Qref, i, &f);
1156:       PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1157:       npoints += Np;
1158:     }
1159:     PetscMalloc1(npoints*dim,&points);
1160:     for (i = 0, k = 0; i < fpdim; ++i) {
1161:       PetscDualSpaceGetFunctional(Qref, i, &f);
1162:       PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1163:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1164:     }

1166:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1167:       PetscObject  obj;
1168:       PetscClassId id;
1169:       PetscReal   *B;
1170:       PetscInt     NcJ = 0, cpdim = 0, j;

1172:       PetscDSGetDiscretization(prob, fieldJ, &obj);
1173:       PetscObjectGetClassId(obj, &id);
1174:       if (id == PETSCFE_CLASSID) {
1175:         PetscFE fe = (PetscFE) obj;

1177:         /* Evaluate basis at points */
1178:         PetscFEGetNumComponents(fe, &NcJ);
1179:         PetscFEGetDimension(fe, &cpdim);
1180:         /* For now, fields only interpolate themselves */
1181:         if (fieldI == fieldJ) {
1182:           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);
1183:           PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1184:           for (i = 0, k = 0; i < fpdim; ++i) {
1185:             PetscDualSpaceGetFunctional(Qref, i, &f);
1186:             PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1187:             for (p = 0; p < Np; ++p, ++k) {
1188:               for (j = 0; j < cpdim; ++j) {
1189:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p];
1190:               }
1191:             }
1192:           }
1193:           PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1194:         }
1195:       } else if (id == PETSCFV_CLASSID) {
1196:         PetscFV        fv = (PetscFV) obj;

1198:         /* Evaluate constant function at points */
1199:         PetscFVGetNumComponents(fv, &NcJ);
1200:         cpdim = 1;
1201:         /* For now, fields only interpolate themselves */
1202:         if (fieldI == fieldJ) {
1203:           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);
1204:           for (i = 0, k = 0; i < fpdim; ++i) {
1205:             PetscDualSpaceGetFunctional(Qref, i, &f);
1206:             PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1207:             for (p = 0; p < Np; ++p, ++k) {
1208:               for (j = 0; j < cpdim; ++j) {
1209:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p];
1210:               }
1211:             }
1212:           }
1213:         }
1214:       }
1215:       offsetJ += cpdim*NcJ;
1216:     }
1217:     offsetI += fpdim*Nc;
1218:     PetscFree(points);
1219:   }
1220:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1221:   /* Preallocate matrix */
1222:   {
1223:     Mat          preallocator;
1224:     PetscScalar *vals;
1225:     PetscInt    *cellCIndices, *cellFIndices;
1226:     PetscInt     locRows, locCols, cell;

1228:     MatGetLocalSize(In, &locRows, &locCols);
1229:     MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1230:     MatSetType(preallocator, MATPREALLOCATOR);
1231:     MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1232:     MatSetUp(preallocator);
1233:     PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1234:     for (cell = cStart; cell < cEnd; ++cell) {
1235:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1236:       MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1237:     }
1238:     PetscFree3(vals,cellCIndices,cellFIndices);
1239:     MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1240:     MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1241:     MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1242:     MatDestroy(&preallocator);
1243:   }
1244:   /* Fill matrix */
1245:   MatZeroEntries(In);
1246:   for (c = cStart; c < cEnd; ++c) {
1247:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1248:   }
1249:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1250:   PetscFree2(feRef,fvRef);
1251:   PetscFree(elemMat);
1252:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1253:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1254:   if (mesh->printFEM) {
1255:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1256:     MatChop(In, 1.0e-10);
1257:     MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1258:   }
1259:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1260:   return(0);
1261: }

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

1268:   Input Parameters:
1269: + dmf  - The fine mesh
1270: . dmc  - The coarse mesh
1271: - user - The user context

1273:   Output Parameter:
1274: . In  - The interpolation matrix

1276:   Level: developer

1278: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1279: @*/
1280: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1281: {
1282:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
1283:   const char    *name = "Interpolator";
1284:   PetscDS        prob;
1285:   PetscSection   fsection, csection, globalFSection, globalCSection;
1286:   PetscHashJK    ht;
1287:   PetscLayout    rLayout;
1288:   PetscInt      *dnz, *onz;
1289:   PetscInt       locRows, rStart, rEnd;
1290:   PetscReal     *x, *v0, *J, *invJ, detJ;
1291:   PetscReal     *v0c, *Jc, *invJc, detJc;
1292:   PetscScalar   *elemMat;
1293:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

1297:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1298:   DMGetCoordinateDim(dmc, &dim);
1299:   DMGetDS(dmc, &prob);
1300:   PetscDSGetRefCoordArrays(prob, &x, NULL);
1301:   PetscDSGetNumFields(prob, &Nf);
1302:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1303:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1304:   DMGetDefaultSection(dmf, &fsection);
1305:   DMGetDefaultGlobalSection(dmf, &globalFSection);
1306:   DMGetDefaultSection(dmc, &csection);
1307:   DMGetDefaultGlobalSection(dmc, &globalCSection);
1308:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1309:   PetscDSGetTotalDimension(prob, &totDim);
1310:   PetscMalloc1(totDim*totDim, &elemMat);

1312:   MatGetLocalSize(In, &locRows, NULL);
1313:   PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1314:   PetscLayoutSetLocalSize(rLayout, locRows);
1315:   PetscLayoutSetBlockSize(rLayout, 1);
1316:   PetscLayoutSetUp(rLayout);
1317:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1318:   PetscLayoutDestroy(&rLayout);
1319:   PetscCalloc2(locRows,&dnz,locRows,&onz);
1320:   PetscHashJKCreate(&ht);
1321:   for (field = 0; field < Nf; ++field) {
1322:     PetscObject      obj;
1323:     PetscClassId     id;
1324:     PetscDualSpace   Q = NULL;
1325:     PetscQuadrature  f;
1326:     const PetscReal *qpoints;
1327:     PetscInt         Nc, Np, fpdim, i, d;

1329:     PetscDSGetDiscretization(prob, field, &obj);
1330:     PetscObjectGetClassId(obj, &id);
1331:     if (id == PETSCFE_CLASSID) {
1332:       PetscFE fe = (PetscFE) obj;

1334:       PetscFEGetDualSpace(fe, &Q);
1335:       PetscFEGetNumComponents(fe, &Nc);
1336:     } else if (id == PETSCFV_CLASSID) {
1337:       PetscFV fv = (PetscFV) obj;

1339:       PetscFVGetDualSpace(fv, &Q);
1340:       Nc   = 1;
1341:     }
1342:     PetscDualSpaceGetDimension(Q, &fpdim);
1343:     /* For each fine grid cell */
1344:     for (cell = cStart; cell < cEnd; ++cell) {
1345:       PetscInt *findices,   *cindices;
1346:       PetscInt  numFIndices, numCIndices;

1348:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1349:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1350:       if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc);
1351:       for (i = 0; i < fpdim; ++i) {
1352:         Vec             pointVec;
1353:         PetscScalar    *pV;
1354:         PetscSF         coarseCellSF = NULL;
1355:         const PetscSFNode *coarseCells;
1356:         PetscInt        numCoarseCells, q, r, c;

1358:         /* Get points from the dual basis functional quadrature */
1359:         PetscDualSpaceGetFunctional(Q, i, &f);
1360:         PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1361:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1362:         VecSetBlockSize(pointVec, dim);
1363:         VecGetArray(pointVec, &pV);
1364:         for (q = 0; q < Np; ++q) {
1365:           /* Transform point to real space */
1366:           CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1367:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1368:         }
1369:         VecRestoreArray(pointVec, &pV);
1370:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1371:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1372:         PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
1373:         /* Update preallocation info */
1374:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1375:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1376:         for (r = 0; r < Nc; ++r) {
1377:           PetscHashJKKey  key;
1378:           PetscHashJKIter missing, iter;

1380:           key.j = findices[i*Nc+r];
1381:           if (key.j < 0) continue;
1382:           /* Get indices for coarse elements */
1383:           for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1384:             DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1385:             for (c = 0; c < numCIndices; ++c) {
1386:               key.k = cindices[c];
1387:               if (key.k < 0) continue;
1388:               PetscHashJKPut(ht, key, &missing, &iter);
1389:               if (missing) {
1390:                 PetscHashJKSet(ht, iter, 1);
1391:                 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1392:                 else                                     ++onz[key.j-rStart];
1393:               }
1394:             }
1395:             DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1396:           }
1397:         }
1398:         PetscSFDestroy(&coarseCellSF);
1399:         VecDestroy(&pointVec);
1400:       }
1401:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1402:     }
1403:   }
1404:   PetscHashJKDestroy(&ht);
1405:   MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1406:   MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1407:   PetscFree2(dnz,onz);
1408:   for (field = 0; field < Nf; ++field) {
1409:     PetscObject      obj;
1410:     PetscClassId     id;
1411:     PetscDualSpace   Q = NULL;
1412:     PetscQuadrature  f;
1413:     const PetscReal *qpoints, *qweights;
1414:     PetscInt         Nc, Np, fpdim, i, d;

1416:     PetscDSGetDiscretization(prob, field, &obj);
1417:     PetscObjectGetClassId(obj, &id);
1418:     if (id == PETSCFE_CLASSID) {
1419:       PetscFE fe = (PetscFE) obj;

1421:       PetscFEGetDualSpace(fe, &Q);
1422:       PetscFEGetNumComponents(fe, &Nc);
1423:     } else if (id == PETSCFV_CLASSID) {
1424:       PetscFV fv = (PetscFV) obj;

1426:       PetscFVGetDualSpace(fv, &Q);
1427:       Nc   = 1;
1428:     }
1429:     PetscDualSpaceGetDimension(Q, &fpdim);
1430:     /* For each fine grid cell */
1431:     for (cell = cStart; cell < cEnd; ++cell) {
1432:       PetscInt *findices,   *cindices;
1433:       PetscInt  numFIndices, numCIndices;

1435:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1436:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1437:       if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc);
1438:       for (i = 0; i < fpdim; ++i) {
1439:         Vec             pointVec;
1440:         PetscScalar    *pV;
1441:         PetscSF         coarseCellSF = NULL;
1442:         const PetscSFNode *coarseCells;
1443:         PetscInt        numCoarseCells, cpdim, q, c, j;

1445:         /* Get points from the dual basis functional quadrature */
1446:         PetscDualSpaceGetFunctional(Q, i, &f);
1447:         PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);
1448:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1449:         VecSetBlockSize(pointVec, dim);
1450:         VecGetArray(pointVec, &pV);
1451:         for (q = 0; q < Np; ++q) {
1452:           /* Transform point to real space */
1453:           CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1454:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1455:         }
1456:         VecRestoreArray(pointVec, &pV);
1457:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1458:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1459:         /* Update preallocation info */
1460:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1461:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1462:         VecGetArray(pointVec, &pV);
1463:         for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1464:           PetscReal pVReal[3];

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

1472:           if (id == PETSCFE_CLASSID) {
1473:             PetscFE    fe = (PetscFE) obj;
1474:             PetscReal *B;

1476:             /* Evaluate coarse basis on contained point */
1477:             PetscFEGetDimension(fe, &cpdim);
1478:             PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
1479:             PetscMemzero(elemMat, cpdim*Nc*Nc * sizeof(PetscScalar));
1480:             /* Get elemMat entries by multiplying by weight */
1481:             for (j = 0; j < cpdim; ++j) {
1482:               for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell];
1483:             }
1484:             PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
1485:           } else {
1486:             cpdim = 1;
1487:             for (j = 0; j < cpdim; ++j) {
1488:               for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell];
1489:             }
1490:           }
1491:           /* Update interpolator */
1492:           if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);}
1493:           if (numCIndices != cpdim*Nc) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D*%D", numCIndices, cpdim, Nc);
1494:           MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);
1495:           DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1496:         }
1497:         VecRestoreArray(pointVec, &pV);
1498:         PetscSFDestroy(&coarseCellSF);
1499:         VecDestroy(&pointVec);
1500:       }
1501:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1502:     }
1503:   }
1504:   PetscFree3(v0,J,invJ);
1505:   PetscFree3(v0c,Jc,invJc);
1506:   PetscFree(elemMat);
1507:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1508:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1509:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1510:   return(0);
1511: }

1515: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1516: {
1517:   PetscDS        prob;
1518:   PetscFE       *feRef;
1519:   PetscFV       *fvRef;
1520:   Vec            fv, cv;
1521:   IS             fis, cis;
1522:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1523:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1524:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;

1528:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1529:   DMGetDimension(dmf, &dim);
1530:   DMGetDefaultSection(dmf, &fsection);
1531:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1532:   DMGetDefaultSection(dmc, &csection);
1533:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1534:   PetscSectionGetNumFields(fsection, &Nf);
1535:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1536:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1537:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1538:   DMGetDS(dmc, &prob);
1539:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1540:   for (f = 0; f < Nf; ++f) {
1541:     PetscObject  obj;
1542:     PetscClassId id;
1543:     PetscInt     fNb = 0, Nc = 0;

1545:     PetscDSGetDiscretization(prob, f, &obj);
1546:     PetscObjectGetClassId(obj, &id);
1547:     if (id == PETSCFE_CLASSID) {
1548:       PetscFE fe = (PetscFE) obj;

1550:       PetscFERefine(fe, &feRef[f]);
1551:       PetscFEGetDimension(feRef[f], &fNb);
1552:       PetscFEGetNumComponents(fe, &Nc);
1553:     } else if (id == PETSCFV_CLASSID) {
1554:       PetscFV        fv = (PetscFV) obj;
1555:       PetscDualSpace Q;

1557:       PetscFVRefine(fv, &fvRef[f]);
1558:       PetscFVGetDualSpace(fvRef[f], &Q);
1559:       PetscDualSpaceGetDimension(Q, &fNb);
1560:       PetscFVGetNumComponents(fv, &Nc);
1561:     }
1562:     fTotDim += fNb*Nc;
1563:   }
1564:   PetscDSGetTotalDimension(prob, &cTotDim);
1565:   PetscMalloc1(cTotDim,&cmap);
1566:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1567:     PetscFE        feC;
1568:     PetscFV        fvC;
1569:     PetscDualSpace QF, QC;
1570:     PetscInt       NcF, NcC, fpdim, cpdim;

1572:     if (feRef[field]) {
1573:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1574:       PetscFEGetNumComponents(feC, &NcC);
1575:       PetscFEGetNumComponents(feRef[field], &NcF);
1576:       PetscFEGetDualSpace(feRef[field], &QF);
1577:       PetscDualSpaceGetDimension(QF, &fpdim);
1578:       PetscFEGetDualSpace(feC, &QC);
1579:       PetscDualSpaceGetDimension(QC, &cpdim);
1580:     } else {
1581:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
1582:       PetscFVGetNumComponents(fvC, &NcC);
1583:       PetscFVGetNumComponents(fvRef[field], &NcF);
1584:       PetscFVGetDualSpace(fvRef[field], &QF);
1585:       PetscDualSpaceGetDimension(QF, &fpdim);
1586:       PetscFVGetDualSpace(fvC, &QC);
1587:       PetscDualSpaceGetDimension(QC, &cpdim);
1588:     }
1589:     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);
1590:     for (c = 0; c < cpdim; ++c) {
1591:       PetscQuadrature  cfunc;
1592:       const PetscReal *cqpoints;
1593:       PetscInt         NpC;
1594:       PetscBool        found = PETSC_FALSE;

1596:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
1597:       PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);
1598:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1599:       for (f = 0; f < fpdim; ++f) {
1600:         PetscQuadrature  ffunc;
1601:         const PetscReal *fqpoints;
1602:         PetscReal        sum = 0.0;
1603:         PetscInt         NpF, comp;

1605:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
1606:         PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);
1607:         if (NpC != NpF) continue;
1608:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1609:         if (sum > 1.0e-9) continue;
1610:         for (comp = 0; comp < NcC; ++comp) {
1611:           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1612:         }
1613:         found = PETSC_TRUE;
1614:         break;
1615:       }
1616:       if (!found) {
1617:         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1618:         if (fvRef[field]) {
1619:           PetscInt comp;
1620:           for (comp = 0; comp < NcC; ++comp) {
1621:             cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1622:           }
1623:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1624:       }
1625:     }
1626:     offsetC += cpdim*NcC;
1627:     offsetF += fpdim*NcF;
1628:   }
1629:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
1630:   PetscFree2(feRef,fvRef);

1632:   DMGetGlobalVector(dmf, &fv);
1633:   DMGetGlobalVector(dmc, &cv);
1634:   VecGetOwnershipRange(cv, &startC, &endC);
1635:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1636:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1637:   PetscMalloc1(m,&cindices);
1638:   PetscMalloc1(m,&findices);
1639:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1640:   for (c = cStart; c < cEnd; ++c) {
1641:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1642:     for (d = 0; d < cTotDim; ++d) {
1643:       if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
1644:       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]]);
1645:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1646:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1647:     }
1648:   }
1649:   PetscFree(cmap);
1650:   PetscFree2(cellCIndices,cellFIndices);

1652:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1653:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1654:   VecScatterCreate(cv, cis, fv, fis, sc);
1655:   ISDestroy(&cis);
1656:   ISDestroy(&fis);
1657:   DMRestoreGlobalVector(dmf, &fv);
1658:   DMRestoreGlobalVector(dmc, &cv);
1659:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1660:   return(0);
1661: }