Actual source code: plexfem.c

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

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

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

 10:   Not collective

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

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

 19:   Level: advanced

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

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

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

 37:   Not collective

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

 44:   Level: advanced

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

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

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

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

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

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

138:   Collective on DM

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

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

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

148:   Level: advanced

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

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

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

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

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

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

213:   Level: advanced

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

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

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

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

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

237:   Level: intermediate

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

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

251: 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)
252: {
253:   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
254:   void            **ctxs;
255:   PetscInt          numFields;
256:   PetscErrorCode    ierr;

259:   DMGetNumFields(dm, &numFields);
260:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
261:   funcs[field] = func;
262:   ctxs[field]  = ctx;
263:   DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);
264:   PetscFree2(funcs,ctxs);
265:   return(0);
266: }

268: static PetscErrorCode DMPlexInsertBoundaryValues_FEM_AuxField_Internal(DM dm, PetscReal time, Vec locU, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[],
269:                                                                        void (*func)(PetscInt, PetscInt, PetscInt,
270:                                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
271:                                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
272:                                                                                     PetscReal, const PetscReal[], PetscScalar[]), void *ctx, Vec locX)
273: {
274:   void (**funcs)(PetscInt, PetscInt, PetscInt,
275:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
276:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
277:                  PetscReal, const PetscReal[], PetscScalar[]);
278:   void            **ctxs;
279:   PetscInt          numFields;
280:   PetscErrorCode    ierr;

283:   DMGetNumFields(dm, &numFields);
284:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
285:   funcs[field] = func;
286:   ctxs[field]  = ctx;
287:   DMProjectFieldLabelLocal(dm, time, label, numids, ids, locU, funcs, INSERT_BC_VALUES, locX);
288:   PetscFree2(funcs,ctxs);
289:   return(0);
290: }

292: /* This ignores numcomps/comps */
293: static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad,
294:                                                               PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
295: {
296:   PetscDS            prob;
297:   PetscSF            sf;
298:   DM                 dmFace, dmCell, dmGrad;
299:   const PetscScalar *facegeom, *cellgeom = NULL, *grad;
300:   const PetscInt    *leaves;
301:   PetscScalar       *x, *fx;
302:   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
303:   PetscErrorCode     ierr, ierru = 0;

306:   DMGetPointSF(dm, &sf);
307:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
308:   nleaves = PetscMax(0, nleaves);
309:   DMGetDimension(dm, &dim);
310:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
311:   DMGetDS(dm, &prob);
312:   VecGetDM(faceGeometry, &dmFace);
313:   VecGetArrayRead(faceGeometry, &facegeom);
314:   if (cellGeometry) {
315:     VecGetDM(cellGeometry, &dmCell);
316:     VecGetArrayRead(cellGeometry, &cellgeom);
317:   }
318:   if (Grad) {
319:     PetscFV fv;

321:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
322:     VecGetDM(Grad, &dmGrad);
323:     VecGetArrayRead(Grad, &grad);
324:     PetscFVGetNumComponents(fv, &pdim);
325:     DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
326:   }
327:   VecGetArray(locX, &x);
328:   for (i = 0; i < numids; ++i) {
329:     IS              faceIS;
330:     const PetscInt *faces;
331:     PetscInt        numFaces, f;

333:     DMLabelGetStratumIS(label, ids[i], &faceIS);
334:     if (!faceIS) continue; /* No points with that id on this process */
335:     ISGetLocalSize(faceIS, &numFaces);
336:     ISGetIndices(faceIS, &faces);
337:     for (f = 0; f < numFaces; ++f) {
338:       const PetscInt         face = faces[f], *cells;
339:       PetscFVFaceGeom        *fg;

341:       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
342:       PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
343:       if (loc >= 0) continue;
344:       DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
345:       DMPlexGetSupport(dm, face, &cells);
346:       if (Grad) {
347:         PetscFVCellGeom       *cg;
348:         PetscScalar           *cx, *cgrad;
349:         PetscScalar           *xG;
350:         PetscReal              dx[3];
351:         PetscInt               d;

353:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
354:         DMPlexPointLocalRead(dm, cells[0], x, &cx);
355:         DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
356:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
357:         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
358:         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
359:         ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
360:         if (ierru) {
361:           ISRestoreIndices(faceIS, &faces);
362:           ISDestroy(&faceIS);
363:           goto cleanup;
364:         }
365:       } else {
366:         PetscScalar       *xI;
367:         PetscScalar       *xG;

369:         DMPlexPointLocalRead(dm, cells[0], x, &xI);
370:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
371:         ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
372:         if (ierru) {
373:           ISRestoreIndices(faceIS, &faces);
374:           ISDestroy(&faceIS);
375:           goto cleanup;
376:         }
377:       }
378:     }
379:     ISRestoreIndices(faceIS, &faces);
380:     ISDestroy(&faceIS);
381:   }
382:   cleanup:
383:   VecRestoreArray(locX, &x);
384:   if (Grad) {
385:     DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);
386:     VecRestoreArrayRead(Grad, &grad);
387:   }
388:   if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
389:   VecRestoreArrayRead(faceGeometry, &facegeom);
390:   CHKERRQ(ierru);
391:   return(0);
392: }

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

397:   Input Parameters:
398: + dm - The DM
399: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
400: . time - The time
401: . faceGeomFVM - Face geometry data for FV discretizations
402: . cellGeomFVM - Cell geometry data for FV discretizations
403: - gradFVM - Gradient reconstruction data for FV discretizations

405:   Output Parameters:
406: . locX - Solution updated with boundary values

408:   Level: developer

410: .seealso: DMProjectFunctionLabelLocal()
411: @*/
412: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
413: {
414:   PetscInt       numBd, b;

423:   PetscDSGetNumBoundary(dm->prob, &numBd);
424:   for (b = 0; b < numBd; ++b) {
425:     DMBoundaryConditionType type;
426:     const char             *labelname;
427:     DMLabel                 label;
428:     PetscInt                field;
429:     PetscObject             obj;
430:     PetscClassId            id;
431:     void                    (*func)(void);
432:     PetscInt                numids;
433:     const PetscInt         *ids;
434:     void                   *ctx;

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

468: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
469: {
470:   const PetscInt   debug = 0;
471:   PetscSection     section;
472:   PetscQuadrature  quad;
473:   Vec              localX;
474:   PetscScalar     *funcVal, *interpolant;
475:   PetscReal       *coords, *detJ, *J;
476:   PetscReal        localDiff = 0.0;
477:   const PetscReal *quadWeights;
478:   PetscInt         dim, coordDim, numFields, numComponents = 0, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
479:   PetscErrorCode   ierr;

482:   DMGetDimension(dm, &dim);
483:   DMGetCoordinateDim(dm, &coordDim);
484:   DMGetDefaultSection(dm, &section);
485:   PetscSectionGetNumFields(section, &numFields);
486:   DMGetLocalVector(dm, &localX);
487:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
488:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
489:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
490:   for (field = 0; field < numFields; ++field) {
491:     PetscObject  obj;
492:     PetscClassId id;
493:     PetscInt     Nc;

495:     DMGetField(dm, field, &obj);
496:     PetscObjectGetClassId(obj, &id);
497:     if (id == PETSCFE_CLASSID) {
498:       PetscFE fe = (PetscFE) obj;

500:       PetscFEGetQuadrature(fe, &quad);
501:       PetscFEGetNumComponents(fe, &Nc);
502:     } else if (id == PETSCFV_CLASSID) {
503:       PetscFV fv = (PetscFV) obj;

505:       PetscFVGetQuadrature(fv, &quad);
506:       PetscFVGetNumComponents(fv, &Nc);
507:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
508:     numComponents += Nc;
509:   }
510:   PetscQuadratureGetData(quad, NULL, &Nq, NULL, &quadWeights);
511:   PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
512:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
513:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
514:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
515:   for (c = cStart; c < cEnd; ++c) {
516:     PetscScalar *x = NULL;
517:     PetscReal    elemDiff = 0.0;

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

522:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
523:       PetscObject  obj;
524:       PetscClassId id;
525:       void * const ctx = ctxs ? ctxs[field] : NULL;
526:       PetscInt     Nb, Nc, q, fc;

528:       DMGetField(dm, field, &obj);
529:       PetscObjectGetClassId(obj, &id);
530:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
531:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
532:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
533:       if (debug) {
534:         char title[1024];
535:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
536:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
537:       }
538:       for (q = 0; q < Nq; ++q) {
539:         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);
540:         (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx);
541:         if (ierr) {
542:           PetscErrorCode ierr2;
543:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
544:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
545:           ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
546: 
547:         }
548:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
549:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
550:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
551:         for (fc = 0; fc < Nc; ++fc) {
552:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]);}
553:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q];
554:         }
555:       }
556:       fieldOffset += Nb*Nc;
557:     }
558:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
559:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
560:     localDiff += elemDiff;
561:   }
562:   PetscFree5(funcVal,interpolant,coords,detJ,J);
563:   DMRestoreLocalVector(dm, &localX);
564:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
565:   *diff = PetscSqrtReal(*diff);
566:   return(0);
567: }

569: 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)
570: {
571:   const PetscInt  debug = 0;
572:   PetscSection    section;
573:   PetscQuadrature quad;
574:   Vec             localX;
575:   PetscScalar    *funcVal, *interpolantVec;
576:   PetscReal      *coords, *realSpaceDer, *J, *invJ, *detJ;
577:   PetscReal       localDiff = 0.0;
578:   PetscInt        dim, coordDim, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp;
579:   PetscErrorCode  ierr;

582:   DMGetDimension(dm, &dim);
583:   DMGetCoordinateDim(dm, &coordDim);
584:   DMGetDefaultSection(dm, &section);
585:   PetscSectionGetNumFields(section, &numFields);
586:   DMGetLocalVector(dm, &localX);
587:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
588:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
589:   for (field = 0; field < numFields; ++field) {
590:     PetscFE  fe;
591:     PetscInt Nc;

593:     DMGetField(dm, field, (PetscObject *) &fe);
594:     PetscFEGetQuadrature(fe, &quad);
595:     PetscQuadratureGetData(quad, NULL, &Nq, NULL, NULL);
596:     PetscFEGetNumComponents(fe, &Nc);
597:     numComponents += Nc;
598:   }
599:   /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
600:   PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,coordDim,&interpolantVec,Nq,&detJ);
601:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
602:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
603:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
604:   for (c = cStart; c < cEnd; ++c) {
605:     PetscScalar *x = NULL;
606:     PetscReal    elemDiff = 0.0;

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

611:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
612:       PetscFE          fe;
613:       void * const     ctx = ctxs ? ctxs[field] : NULL;
614:       const PetscReal *quadPoints, *quadWeights;
615:       PetscReal       *basisDer;
616:       PetscInt         numQuadPoints, Nb, Ncomp, q, d, fc, f, g;

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

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

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

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, *detJ, *J;
680:   PetscReal       *localDiff;
681:   const PetscReal *quadPoints, *quadWeights;
682:   PetscInt         dim, coordDim, numFields, numComponents = 0, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
683:   PetscErrorCode   ierr;

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

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

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

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

722:     DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);
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 < Nq; ++q) {
744:         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);
745:         (*funcs[field])(coordDim, time, &coords[coordDim*q], 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 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);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, coordDim > 0 ? coords[0] : 0., coordDim > 1 ? coords[1] : 0., coordDim > 2 ? coords[2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]);}
758:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q];
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:   PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);
770:   return(0);
771: }

773: /*@C
774:   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.

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

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

786:   Level: developer

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

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

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

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

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

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

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

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

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

887:   Input Parameters:
888: + dm - The mesh
889: . X  - Local input vector
890: - user - The user context

892:   Output Parameter:
893: . integral - Local integral for each field

895:   Level: developer

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

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

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

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

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

1002:     PetscDSGetDiscretization(prob, f, &obj);
1003:     PetscObjectGetClassId(obj, &id);
1004:     if (id == PETSCFE_CLASSID) {
1005:       PetscFE         fe = (PetscFE) obj;
1006:       PetscQuadrature q;
1007:       PetscInt        Nq, Nb;

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

1028:       PetscDSGetObjective(prob, f, &obj_func);
1029:       PetscDSGetFieldOffset(prob, f, &foff);
1030:       if (obj_func) {
1031:         for (c = 0; c < numCells; ++c) {
1032:           PetscScalar *u_x;

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

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

1066:   Input Parameters:
1067: + dmf  - The fine mesh
1068: . dmc  - The coarse mesh
1069: - user - The user context

1071:   Output Parameter:
1072: . In  - The interpolation matrix

1074:   Level: developer

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

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

1110:     PetscDSGetDiscretization(prob, f, &obj);
1111:     PetscObjectGetClassId(obj, &id);
1112:     if (id == PETSCFE_CLASSID) {
1113:       PetscFE fe = (PetscFE) obj;

1115:       PetscFERefine(fe, &feRef[f]);
1116:       PetscFEGetDimension(feRef[f], &rNb);
1117:       PetscFEGetNumComponents(fe, &Nc);
1118:     } else if (id == PETSCFV_CLASSID) {
1119:       PetscFV        fv = (PetscFV) obj;
1120:       PetscDualSpace Q;

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

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

1160:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1161:       PetscObject  obj;
1162:       PetscClassId id;
1163:       PetscReal   *B;
1164:       PetscInt     NcJ = 0, cpdim = 0, j;

1166:       PetscDSGetDiscretization(prob, fieldJ, &obj);
1167:       PetscObjectGetClassId(obj, &id);
1168:       if (id == PETSCFE_CLASSID) {
1169:         PetscFE fe = (PetscFE) obj;

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

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

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

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

1260:   Input Parameters:
1261: + dmf  - The fine mesh
1262: . dmc  - The coarse mesh
1263: - user - The user context

1265:   Output Parameter:
1266: . In  - The interpolation matrix

1268:   Level: developer

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

1289:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1290:   DMGetCoordinateDim(dmc, &dim);
1291:   DMGetDS(dmc, &prob);
1292:   PetscDSGetRefCoordArrays(prob, &x, NULL);
1293:   PetscDSGetNumFields(prob, &Nf);
1294:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1295:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1296:   DMGetDefaultSection(dmf, &fsection);
1297:   DMGetDefaultGlobalSection(dmf, &globalFSection);
1298:   DMGetDefaultSection(dmc, &csection);
1299:   DMGetDefaultGlobalSection(dmc, &globalCSection);
1300:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1301:   PetscDSGetTotalDimension(prob, &totDim);
1302:   PetscMalloc1(totDim*totDim, &elemMat);

1304:   MatGetLocalSize(In, &locRows, NULL);
1305:   PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1306:   PetscLayoutSetLocalSize(rLayout, locRows);
1307:   PetscLayoutSetBlockSize(rLayout, 1);
1308:   PetscLayoutSetUp(rLayout);
1309:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1310:   PetscLayoutDestroy(&rLayout);
1311:   PetscCalloc2(locRows,&dnz,locRows,&onz);
1312:   PetscHashJKCreate(&ht);
1313:   for (field = 0; field < Nf; ++field) {
1314:     PetscObject      obj;
1315:     PetscClassId     id;
1316:     PetscDualSpace   Q = NULL;
1317:     PetscQuadrature  f;
1318:     const PetscReal *qpoints;
1319:     PetscInt         Nc, Np, fpdim, i, d;

1321:     PetscDSGetDiscretization(prob, field, &obj);
1322:     PetscObjectGetClassId(obj, &id);
1323:     if (id == PETSCFE_CLASSID) {
1324:       PetscFE fe = (PetscFE) obj;

1326:       PetscFEGetDualSpace(fe, &Q);
1327:       PetscFEGetNumComponents(fe, &Nc);
1328:     } else if (id == PETSCFV_CLASSID) {
1329:       PetscFV fv = (PetscFV) obj;

1331:       PetscFVGetDualSpace(fv, &Q);
1332:       Nc   = 1;
1333:     }
1334:     PetscDualSpaceGetDimension(Q, &fpdim);
1335:     /* For each fine grid cell */
1336:     for (cell = cStart; cell < cEnd; ++cell) {
1337:       PetscInt *findices,   *cindices;
1338:       PetscInt  numFIndices, numCIndices;

1340:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1341:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1342:       if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc);
1343:       for (i = 0; i < fpdim; ++i) {
1344:         Vec             pointVec;
1345:         PetscScalar    *pV;
1346:         PetscSF         coarseCellSF = NULL;
1347:         const PetscSFNode *coarseCells;
1348:         PetscInt        numCoarseCells, q, r, c;

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

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

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

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

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

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

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

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

1464:           if (id == PETSCFE_CLASSID) {
1465:             PetscFE    fe = (PetscFE) obj;
1466:             PetscReal *B;

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

1505: /*@
1506:   DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns

1508:   Input Parameters:
1509: + dmc  - The coarse mesh
1510: - dmf  - The fine mesh
1511: - user - The user context

1513:   Output Parameter:
1514: . sc   - The mapping

1516:   Level: developer

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

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

1550:     PetscDSGetDiscretization(prob, f, &obj);
1551:     PetscObjectGetClassId(obj, &id);
1552:     if (id == PETSCFE_CLASSID) {
1553:       PetscFE fe = (PetscFE) obj;

1555:       PetscFERefine(fe, &feRef[f]);
1556:       PetscFEGetDimension(feRef[f], &fNb);
1557:       PetscFEGetNumComponents(fe, &Nc);
1558:     } else if (id == PETSCFV_CLASSID) {
1559:       PetscFV        fv = (PetscFV) obj;
1560:       PetscDualSpace Q;

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

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

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

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

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

1657:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1658:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1659:   VecScatterCreate(cv, cis, fv, fis, sc);
1660:   ISDestroy(&cis);
1661:   ISDestroy(&fis);
1662:   DMRestoreGlobalVector(dmf, &fv);
1663:   DMRestoreGlobalVector(dmc, &cv);
1664:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1665:   return(0);
1666: }