Actual source code: plexfem.c

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

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

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

 10:   Not collective

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

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

 19:   Level: advanced

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

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

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

 37:   Not collective

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

 44:   Level: advanced

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

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

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

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

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

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

138:   Collective on DM

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

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

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

148:   Level: advanced

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

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

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

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

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

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

213:   Level: advanced

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

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

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

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

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

237:   Level: intermediate

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

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

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

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

264:   Output Parameter:
265: . locX   - A local vector to receives the boundary values

267:   Level: developer

269: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
270: @*/
271: PetscErrorCode DMPlexInsertBoundaryValuesEssential(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)
272: {
273:   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
274:   void            **ctxs;
275:   PetscInt          numFields;
276:   PetscErrorCode    ierr;

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

288: /*@C
289:   DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector

291:   Input Parameters:
292: + dm     - The DM, with a PetscDS that matches the problem being constrained
293: . time   - The time
294: . locU   - A local vector with the input solution values
295: . field  - The field to constrain
296: . label  - The DMLabel defining constrained points
297: . numids - The number of DMLabel ids for constrained points
298: . ids    - An array of ids for constrained points
299: . func   - A pointwise function giving boundary values
300: - ctx    - An optional user context for bcFunc

302:   Output Parameter:
303: . locX   - A local vector to receives the boundary values

305:   Level: developer

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

326:   DMGetNumFields(dm, &numFields);
327:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
328:   funcs[field] = func;
329:   ctxs[field]  = ctx;
330:   DMProjectFieldLabelLocal(dm, time, label, numids, ids, locU, funcs, INSERT_BC_VALUES, locX);
331:   PetscFree2(funcs,ctxs);
332:   return(0);
333: }

335: /*@C
336:   DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector

338:   Input Parameters:
339: + dm     - The DM, with a PetscDS that matches the problem being constrained
340: . time   - The time
341: . faceGeometry - A vector with the FVM face geometry information
342: . cellGeometry - A vector with the FVM cell geometry information
343: . Grad         - A vector with the FVM cell gradient information
344: . field  - The field to constrain
345: . label  - The DMLabel defining constrained points
346: . numids - The number of DMLabel ids for constrained points
347: . ids    - An array of ids for constrained points
348: . func   - A pointwise function giving boundary values
349: - ctx    - An optional user context for bcFunc

351:   Output Parameter:
352: . locX   - A local vector to receives the boundary values

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

356:   Level: developer

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

373:   DMGetPointSF(dm, &sf);
374:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
375:   nleaves = PetscMax(0, nleaves);
376:   DMGetDimension(dm, &dim);
377:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
378:   DMGetDS(dm, &prob);
379:   VecGetDM(faceGeometry, &dmFace);
380:   VecGetArrayRead(faceGeometry, &facegeom);
381:   if (cellGeometry) {
382:     VecGetDM(cellGeometry, &dmCell);
383:     VecGetArrayRead(cellGeometry, &cellgeom);
384:   }
385:   if (Grad) {
386:     PetscFV fv;

388:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
389:     VecGetDM(Grad, &dmGrad);
390:     VecGetArrayRead(Grad, &grad);
391:     PetscFVGetNumComponents(fv, &pdim);
392:     DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
393:   }
394:   VecGetArray(locX, &x);
395:   for (i = 0; i < numids; ++i) {
396:     IS              faceIS;
397:     const PetscInt *faces;
398:     PetscInt        numFaces, f;

400:     DMLabelGetStratumIS(label, ids[i], &faceIS);
401:     if (!faceIS) continue; /* No points with that id on this process */
402:     ISGetLocalSize(faceIS, &numFaces);
403:     ISGetIndices(faceIS, &faces);
404:     for (f = 0; f < numFaces; ++f) {
405:       const PetscInt         face = faces[f], *cells;
406:       PetscFVFaceGeom        *fg;

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

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

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

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

464:   Input Parameters:
465: + dm - The DM
466: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
467: . time - The time
468: . faceGeomFVM - Face geometry data for FV discretizations
469: . cellGeomFVM - Cell geometry data for FV discretizations
470: - gradFVM - Gradient reconstruction data for FV discretizations

472:   Output Parameters:
473: . locX - Solution updated with boundary values

475:   Level: developer

477: .seealso: DMProjectFunctionLabelLocal()
478: @*/
479: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
480: {
481:   PetscInt       numBd, b;

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

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

535: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
536: {
537:   const PetscInt   debug = 0;
538:   PetscSection     section;
539:   PetscQuadrature  quad;
540:   Vec              localX;
541:   PetscScalar     *funcVal, *interpolant;
542:   PetscReal       *coords, *detJ, *J;
543:   PetscReal        localDiff = 0.0;
544:   const PetscReal *quadWeights;
545:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
546:   PetscErrorCode   ierr;

549:   DMGetDimension(dm, &dim);
550:   DMGetCoordinateDim(dm, &coordDim);
551:   DMGetDefaultSection(dm, &section);
552:   PetscSectionGetNumFields(section, &numFields);
553:   DMGetLocalVector(dm, &localX);
554:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
555:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
556:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
557:   for (field = 0; field < numFields; ++field) {
558:     PetscObject  obj;
559:     PetscClassId id;
560:     PetscInt     Nc;

562:     DMGetField(dm, field, &obj);
563:     PetscObjectGetClassId(obj, &id);
564:     if (id == PETSCFE_CLASSID) {
565:       PetscFE fe = (PetscFE) obj;

567:       PetscFEGetQuadrature(fe, &quad);
568:       PetscFEGetNumComponents(fe, &Nc);
569:     } else if (id == PETSCFV_CLASSID) {
570:       PetscFV fv = (PetscFV) obj;

572:       PetscFVGetQuadrature(fv, &quad);
573:       PetscFVGetNumComponents(fv, &Nc);
574:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
575:     numComponents += Nc;
576:   }
577:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
578:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
579:   PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
580:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
581:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
582:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
583:   for (c = cStart; c < cEnd; ++c) {
584:     PetscScalar *x = NULL;
585:     PetscReal    elemDiff = 0.0;
586:     PetscInt     qc = 0;

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

591:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
592:       PetscObject  obj;
593:       PetscClassId id;
594:       void * const ctx = ctxs ? ctxs[field] : NULL;
595:       PetscInt     Nb, Nc, q, fc;

597:       DMGetField(dm, field, &obj);
598:       PetscObjectGetClassId(obj, &id);
599:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
600:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
601:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
602:       if (debug) {
603:         char title[1024];
604:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
605:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
606:       }
607:       for (q = 0; q < Nq; ++q) {
608:         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);
609:         (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx);
610:         if (ierr) {
611:           PetscErrorCode ierr2;
612:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
613:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
614:           ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
615: 
616:         }
617:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
618:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
619:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
620:         for (fc = 0; fc < Nc; ++fc) {
621:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
622:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
623:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
624:         }
625:       }
626:       fieldOffset += Nb;
627:       qc += Nc;
628:     }
629:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
630:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
631:     localDiff += elemDiff;
632:   }
633:   PetscFree5(funcVal,interpolant,coords,detJ,J);
634:   DMRestoreLocalVector(dm, &localX);
635:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
636:   *diff = PetscSqrtReal(*diff);
637:   return(0);
638: }

640: 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)
641: {
642:   const PetscInt   debug = 0;
643:   PetscSection     section;
644:   PetscQuadrature  quad;
645:   Vec              localX;
646:   PetscScalar     *funcVal, *interpolant;
647:   const PetscReal *quadPoints, *quadWeights;
648:   PetscReal       *coords, *realSpaceDer, *J, *invJ, *detJ;
649:   PetscReal        localDiff = 0.0;
650:   PetscInt         dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset;
651:   PetscErrorCode   ierr;

654:   DMGetDimension(dm, &dim);
655:   DMGetCoordinateDim(dm, &coordDim);
656:   DMGetDefaultSection(dm, &section);
657:   PetscSectionGetNumFields(section, &numFields);
658:   DMGetLocalVector(dm, &localX);
659:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
660:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
661:   for (field = 0; field < numFields; ++field) {
662:     PetscFE  fe;
663:     PetscInt Nc;

665:     DMGetField(dm, field, (PetscObject *) &fe);
666:     PetscFEGetQuadrature(fe, &quad);
667:     PetscFEGetNumComponents(fe, &Nc);
668:     numComponents += Nc;
669:   }
670:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
671:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
672:   /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
673:   PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,numComponents,&interpolant,Nq,&detJ);
674:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
675:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
676:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
677:   for (c = cStart; c < cEnd; ++c) {
678:     PetscScalar *x = NULL;
679:     PetscReal    elemDiff = 0.0;
680:     PetscInt     qc = 0;

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

685:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
686:       PetscFE          fe;
687:       void * const     ctx = ctxs ? ctxs[field] : NULL;
688:       PetscReal       *basisDer;
689:       PetscInt         Nb, Nc, q, fc;

691:       DMGetField(dm, field, (PetscObject *) &fe);
692:       PetscFEGetDimension(fe, &Nb);
693:       PetscFEGetNumComponents(fe, &Nc);
694:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
695:       if (debug) {
696:         char title[1024];
697:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
698:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
699:       }
700:       for (q = 0; q < Nq; ++q) {
701:         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);
702:         (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx);
703:         if (ierr) {
704:           PetscErrorCode ierr2;
705:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
706:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
707:           ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);CHKERRQ(ierr2);
708: 
709:         }
710:         PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], coordDim, invJ, n, q, interpolant);
711:         for (fc = 0; fc < Nc; ++fc) {
712:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
713:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
714:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
715:         }
716:       }
717:       fieldOffset += Nb;
718:       qc          += Nc;
719:     }
720:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
721:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
722:     localDiff += elemDiff;
723:   }
724:   PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolant,detJ);
725:   DMRestoreLocalVector(dm, &localX);
726:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
727:   *diff = PetscSqrtReal(*diff);
728:   return(0);
729: }

731: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
732: {
733:   const PetscInt   debug = 0;
734:   PetscSection     section;
735:   PetscQuadrature  quad;
736:   Vec              localX;
737:   PetscScalar     *funcVal, *interpolant;
738:   PetscReal       *coords, *detJ, *J;
739:   PetscReal       *localDiff;
740:   const PetscReal *quadPoints, *quadWeights;
741:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
742:   PetscErrorCode   ierr;

745:   DMGetDimension(dm, &dim);
746:   DMGetCoordinateDim(dm, &coordDim);
747:   DMGetDefaultSection(dm, &section);
748:   PetscSectionGetNumFields(section, &numFields);
749:   DMGetLocalVector(dm, &localX);
750:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
751:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
752:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
753:   for (field = 0; field < numFields; ++field) {
754:     PetscObject  obj;
755:     PetscClassId id;
756:     PetscInt     Nc;

758:     DMGetField(dm, field, &obj);
759:     PetscObjectGetClassId(obj, &id);
760:     if (id == PETSCFE_CLASSID) {
761:       PetscFE fe = (PetscFE) obj;

763:       PetscFEGetQuadrature(fe, &quad);
764:       PetscFEGetNumComponents(fe, &Nc);
765:     } else if (id == PETSCFV_CLASSID) {
766:       PetscFV fv = (PetscFV) obj;

768:       PetscFVGetQuadrature(fv, &quad);
769:       PetscFVGetNumComponents(fv, &Nc);
770:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
771:     numComponents += Nc;
772:   }
773:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
774:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
775:   PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
776:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
777:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
778:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
779:   for (c = cStart; c < cEnd; ++c) {
780:     PetscScalar *x = NULL;
781:     PetscInt     qc = 0;

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

786:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
787:       PetscObject  obj;
788:       PetscClassId id;
789:       void * const ctx = ctxs ? ctxs[field] : NULL;
790:       PetscInt     Nb, Nc, q, fc;

792:       PetscReal       elemDiff = 0.0;

794:       DMGetField(dm, field, &obj);
795:       PetscObjectGetClassId(obj, &id);
796:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
797:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
798:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
799:       if (debug) {
800:         char title[1024];
801:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
802:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
803:       }
804:       for (q = 0; q < Nq; ++q) {
805:         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);
806:         (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx);
807:         if (ierr) {
808:           PetscErrorCode ierr2;
809:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
810:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
811:           ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2);
812: 
813:         }
814:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
815:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
816:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
817:         for (fc = 0; fc < Nc; ++fc) {
818:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
819:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d point %g %g %g diff %g\n", c, field, coordDim > 0 ? coords[0] : 0., coordDim > 1 ? coords[1] : 0., coordDim > 2 ? coords[2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q]);}
820:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
821:         }
822:       }
823:       fieldOffset += Nb;
824:       qc          += Nc;
825:       localDiff[field] += elemDiff;
826:     }
827:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
828:   }
829:   DMRestoreLocalVector(dm, &localX);
830:   MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
831:   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
832:   PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);
833:   return(0);
834: }

836: /*@C
837:   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.

839:   Input Parameters:
840: + dm    - The DM
841: . time  - The time
842: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
843: . ctxs  - Optional array of contexts to pass to each function, or NULL.
844: - X     - The coefficient vector u_h

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

849:   Level: developer

851: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
852: @*/
853: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
854: {
855:   PetscSection     section;
856:   PetscQuadrature  quad;
857:   Vec              localX;
858:   PetscScalar     *funcVal, *interpolant;
859:   PetscReal       *coords, *detJ, *J;
860:   const PetscReal *quadPoints, *quadWeights;
861:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
862:   PetscErrorCode   ierr;

865:   VecSet(D, 0.0);
866:   DMGetDimension(dm, &dim);
867:   DMGetCoordinateDim(dm, &coordDim);
868:   DMGetDefaultSection(dm, &section);
869:   PetscSectionGetNumFields(section, &numFields);
870:   DMGetLocalVector(dm, &localX);
871:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
872:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
873:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
874:   for (field = 0; field < numFields; ++field) {
875:     PetscObject  obj;
876:     PetscClassId id;
877:     PetscInt     Nc;

879:     DMGetField(dm, field, &obj);
880:     PetscObjectGetClassId(obj, &id);
881:     if (id == PETSCFE_CLASSID) {
882:       PetscFE fe = (PetscFE) obj;

884:       PetscFEGetQuadrature(fe, &quad);
885:       PetscFEGetNumComponents(fe, &Nc);
886:     } else if (id == PETSCFV_CLASSID) {
887:       PetscFV fv = (PetscFV) obj;

889:       PetscFVGetQuadrature(fv, &quad);
890:       PetscFVGetNumComponents(fv, &Nc);
891:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
892:     numComponents += Nc;
893:   }
894:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
895:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
896:   PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);
897:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
898:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
899:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
900:   for (c = cStart; c < cEnd; ++c) {
901:     PetscScalar *x = NULL;
902:     PetscScalar  elemDiff = 0.0;
903:     PetscInt     qc = 0;

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

908:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
909:       PetscObject  obj;
910:       PetscClassId id;
911:       void * const ctx = ctxs ? ctxs[field] : NULL;
912:       PetscInt     Nb, Nc, q, fc;

914:       DMGetField(dm, field, &obj);
915:       PetscObjectGetClassId(obj, &id);
916:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
917:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
918:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
919:       if (funcs[field]) {
920:         for (q = 0; q < Nq; ++q) {
921:           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);
922:           (*funcs[field])(dim, time, &coords[q*coordDim], Nc, funcVal, ctx);
923:           if (ierr) {
924:             PetscErrorCode ierr2;
925:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
926:             ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);
927:             ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
928: 
929:           }
930:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
931:           else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
932:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
933:           for (fc = 0; fc < Nc; ++fc) {
934:             const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
935:             elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*detJ[q];
936:           }
937:         }
938:       }
939:       fieldOffset += Nb;
940:       qc          += Nc;
941:     }
942:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
943:     VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
944:   }
945:   PetscFree5(funcVal,interpolant,coords,detJ,J);
946:   DMRestoreLocalVector(dm, &localX);
947:   VecSqrtAbs(D);
948:   return(0);
949: }

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

954:   Input Parameters:
955: + dm - The mesh
956: . X  - Local input vector
957: - user - The user context

959:   Output Parameter:
960: . integral - Local integral for each field

962:   Level: developer

964: .seealso: DMPlexComputeResidualFEM()
965: @*/
966: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
967: {
968:   DM_Plex           *mesh  = (DM_Plex *) dm->data;
969:   DM                 dmAux, dmGrad;
970:   Vec                localX, A, cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
971:   PetscDS            prob, probAux = NULL;
972:   PetscSection       section, sectionAux;
973:   PetscFV            fvm = NULL;
974:   PetscFECellGeom   *cgeomFEM;
975:   PetscFVFaceGeom   *fgeomFVM;
976:   PetscFVCellGeom   *cgeomFVM;
977:   PetscScalar       *u, *a = NULL;
978:   const PetscScalar *constants, *lgrad;
979:   PetscReal         *lintegral;
980:   PetscInt          *uOff, *uOff_x, *aOff = NULL;
981:   PetscInt           dim, numConstants, Nf, NfAux = 0, f, numCells, cStart, cEnd, cEndInterior, c;
982:   PetscInt           totDim, totDimAux;
983:   PetscBool          useFVM = PETSC_FALSE;
984:   PetscErrorCode     ierr;

987:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
988:   DMGetLocalVector(dm, &localX);
989:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);
990:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
991:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
992:   DMGetDimension(dm, &dim);
993:   DMGetDefaultSection(dm, &section);
994:   DMGetDS(dm, &prob);
995:   PetscDSGetTotalDimension(prob, &totDim);
996:   PetscDSGetComponentOffsets(prob, &uOff);
997:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
998:   PetscDSGetConstants(prob, &numConstants, &constants);
999:   PetscSectionGetNumFields(section, &Nf);
1000:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1001:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1002:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1003:   numCells = cEnd - cStart;
1004:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1005:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1006:   if (dmAux) {
1007:     DMGetDS(dmAux, &probAux);
1008:     PetscDSGetNumFields(probAux, &NfAux);
1009:     DMGetDefaultSection(dmAux, &sectionAux);
1010:     PetscDSGetTotalDimension(probAux, &totDimAux);
1011:     PetscDSGetComponentOffsets(probAux, &aOff);
1012:     PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);
1013:   }
1014:   for (f = 0; f < Nf; ++f) {
1015:     PetscObject  obj;
1016:     PetscClassId id;

1018:     PetscDSGetDiscretization(prob, f, &obj);
1019:     PetscObjectGetClassId(obj, &id);
1020:     if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1021:   }
1022:   if (useFVM) {
1023:     Vec       grad;
1024:     PetscInt  fStart, fEnd;
1025:     PetscBool compGrad;

1027:     PetscFVGetComputeGradients(fvm, &compGrad);
1028:     PetscFVSetComputeGradients(fvm, PETSC_TRUE);
1029:     DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1030:     DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1031:     PetscFVSetComputeGradients(fvm, compGrad);
1032:     VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1033:     VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1034:     /* Reconstruct and limit cell gradients */
1035:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1036:     DMGetGlobalVector(dmGrad, &grad);
1037:     DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);
1038:     /* Communicate gradient values */
1039:     DMGetLocalVector(dmGrad, &locGrad);
1040:     DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1041:     DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1042:     DMRestoreGlobalVector(dmGrad, &grad);
1043:     /* Handle non-essential (e.g. outflow) boundary values */
1044:     DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1045:     VecGetArrayRead(locGrad, &lgrad);
1046:   }
1047:   PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);
1048:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1049:   for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
1050:   for (c = cStart; c < cEnd; ++c) {
1051:     PetscScalar *x = NULL;
1052:     PetscInt     i;

1054:     DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);
1055:     if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c);
1056:     DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
1057:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1058:     DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
1059:     if (dmAux) {
1060:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1061:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1062:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1063:     }
1064:   }
1065:   for (f = 0; f < Nf; ++f) {
1066:     PetscObject  obj;
1067:     PetscClassId id;
1068:     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1070:     PetscDSGetDiscretization(prob, f, &obj);
1071:     PetscObjectGetClassId(obj, &id);
1072:     if (id == PETSCFE_CLASSID) {
1073:       PetscFE         fe = (PetscFE) obj;
1074:       PetscQuadrature q;
1075:       PetscInt        Nq, Nb;

1077:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1078:       PetscFEGetQuadrature(fe, &q);
1079:       PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1080:       PetscFEGetDimension(fe, &Nb);
1081:       blockSize = Nb*Nq;
1082:       batchSize = numBlocks * blockSize;
1083:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1084:       numChunks = numCells / (numBatches*batchSize);
1085:       Ne        = numChunks*numBatches*batchSize;
1086:       Nr        = numCells % (numBatches*batchSize);
1087:       offset    = numCells - Nr;
1088:       PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);
1089:       PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);
1090:     } else if (id == PETSCFV_CLASSID) {
1091:       /* PetscFV  fv = (PetscFV) obj; */
1092:       PetscInt       foff;
1093:       PetscPointFunc obj_func;
1094:       PetscScalar    lint;

1096:       PetscDSGetObjective(prob, f, &obj_func);
1097:       PetscDSGetFieldOffset(prob, f, &foff);
1098:       if (obj_func) {
1099:         for (c = 0; c < numCells; ++c) {
1100:           PetscScalar *u_x;

1102:           DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1103:           obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, a, NULL, NULL, 0.0, cgeomFVM[c].centroid, numConstants, constants, &lint);
1104:           lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1105:         }
1106:       }
1107:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1108:   }
1109:   if (useFVM) {
1110:     VecRestoreArrayRead(locGrad, &lgrad);
1111:     VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1112:     VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1113:     DMRestoreLocalVector(dmGrad, &locGrad);
1114:     VecDestroy(&faceGeometryFVM);
1115:     VecDestroy(&cellGeometryFVM);
1116:     DMDestroy(&dmGrad);
1117:   }
1118:   if (dmAux) {PetscFree(a);}
1119:   if (mesh->printFEM) {
1120:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1121:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);}
1122:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1123:   }
1124:   DMRestoreLocalVector(dm, &localX);
1125:   MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1126:   PetscFree3(lintegral,u,cgeomFEM);
1127:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1128:   return(0);
1129: }

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

1134:   Input Parameters:
1135: + dmf  - The fine mesh
1136: . dmc  - The coarse mesh
1137: - user - The user context

1139:   Output Parameter:
1140: . In  - The interpolation matrix

1142:   Level: developer

1144: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1145: @*/
1146: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1147: {
1148:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1149:   const char       *name  = "Interpolator";
1150:   PetscDS           prob;
1151:   PetscFE          *feRef;
1152:   PetscFV          *fvRef;
1153:   PetscSection      fsection, fglobalSection;
1154:   PetscSection      csection, cglobalSection;
1155:   PetscScalar      *elemMat;
1156:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1157:   PetscInt          cTotDim, rTotDim = 0;
1158:   PetscErrorCode    ierr;

1161:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1162:   DMGetDimension(dmf, &dim);
1163:   DMGetDefaultSection(dmf, &fsection);
1164:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1165:   DMGetDefaultSection(dmc, &csection);
1166:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1167:   PetscSectionGetNumFields(fsection, &Nf);
1168:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1169:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1170:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1171:   DMGetDS(dmf, &prob);
1172:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1173:   for (f = 0; f < Nf; ++f) {
1174:     PetscObject  obj;
1175:     PetscClassId id;
1176:     PetscInt     rNb = 0, Nc = 0;

1178:     PetscDSGetDiscretization(prob, f, &obj);
1179:     PetscObjectGetClassId(obj, &id);
1180:     if (id == PETSCFE_CLASSID) {
1181:       PetscFE fe = (PetscFE) obj;

1183:       PetscFERefine(fe, &feRef[f]);
1184:       PetscFEGetDimension(feRef[f], &rNb);
1185:       PetscFEGetNumComponents(fe, &Nc);
1186:     } else if (id == PETSCFV_CLASSID) {
1187:       PetscFV        fv = (PetscFV) obj;
1188:       PetscDualSpace Q;

1190:       PetscFVRefine(fv, &fvRef[f]);
1191:       PetscFVGetDualSpace(fvRef[f], &Q);
1192:       PetscDualSpaceGetDimension(Q, &rNb);
1193:       PetscFVGetNumComponents(fv, &Nc);
1194:     }
1195:     rTotDim += rNb;
1196:   }
1197:   PetscDSGetTotalDimension(prob, &cTotDim);
1198:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1199:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1200:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1201:     PetscDualSpace   Qref;
1202:     PetscQuadrature  f;
1203:     const PetscReal *qpoints, *qweights;
1204:     PetscReal       *points;
1205:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1207:     /* Compose points from all dual basis functionals */
1208:     if (feRef[fieldI]) {
1209:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
1210:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
1211:     } else {
1212:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1213:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1214:     }
1215:     PetscDualSpaceGetDimension(Qref, &fpdim);
1216:     for (i = 0; i < fpdim; ++i) {
1217:       PetscDualSpaceGetFunctional(Qref, i, &f);
1218:       PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
1219:       npoints += Np;
1220:     }
1221:     PetscMalloc1(npoints*dim,&points);
1222:     for (i = 0, k = 0; i < fpdim; ++i) {
1223:       PetscDualSpaceGetFunctional(Qref, i, &f);
1224:       PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1225:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1226:     }

1228:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1229:       PetscObject  obj;
1230:       PetscClassId id;
1231:       PetscReal   *B;
1232:       PetscInt     NcJ = 0, cpdim = 0, j, qNc;

1234:       PetscDSGetDiscretization(prob, fieldJ, &obj);
1235:       PetscObjectGetClassId(obj, &id);
1236:       if (id == PETSCFE_CLASSID) {
1237:         PetscFE fe = (PetscFE) obj;

1239:         /* Evaluate basis at points */
1240:         PetscFEGetNumComponents(fe, &NcJ);
1241:         PetscFEGetDimension(fe, &cpdim);
1242:         /* For now, fields only interpolate themselves */
1243:         if (fieldI == fieldJ) {
1244:           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);
1245:           PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1246:           for (i = 0, k = 0; i < fpdim; ++i) {
1247:             PetscDualSpaceGetFunctional(Qref, i, &f);
1248:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1249:             if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ);
1250:             for (p = 0; p < Np; ++p, ++k) {
1251:               for (j = 0; j < cpdim; ++j) {
1252:                 /* NOTE: This is not quite right, unless fpdim == number of fine grid functional quad points */
1253:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
1254:               }
1255:             }
1256:           }
1257:           PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1258:         }
1259:       } else if (id == PETSCFV_CLASSID) {
1260:         PetscFV        fv = (PetscFV) obj;

1262:         /* Evaluate constant function at points */
1263:         PetscFVGetNumComponents(fv, &NcJ);
1264:         cpdim = 1;
1265:         /* For now, fields only interpolate themselves */
1266:         if (fieldI == fieldJ) {
1267:           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);
1268:           for (i = 0, k = 0; i < fpdim; ++i) {
1269:             PetscDualSpaceGetFunctional(Qref, i, &f);
1270:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
1271:             if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ);
1272:             for (p = 0; p < Np; ++p, ++k) {
1273:               for (j = 0; j < cpdim; ++j) {
1274:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p*qNc+c];
1275:               }
1276:             }
1277:           }
1278:         }
1279:       }
1280:       offsetJ += cpdim*NcJ;
1281:     }
1282:     offsetI += fpdim*Nc;
1283:     PetscFree(points);
1284:   }
1285:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1286:   /* Preallocate matrix */
1287:   {
1288:     Mat          preallocator;
1289:     PetscScalar *vals;
1290:     PetscInt    *cellCIndices, *cellFIndices;
1291:     PetscInt     locRows, locCols, cell;

1293:     MatGetLocalSize(In, &locRows, &locCols);
1294:     MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1295:     MatSetType(preallocator, MATPREALLOCATOR);
1296:     MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1297:     MatSetUp(preallocator);
1298:     PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1299:     for (cell = cStart; cell < cEnd; ++cell) {
1300:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1301:       MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1302:     }
1303:     PetscFree3(vals,cellCIndices,cellFIndices);
1304:     MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1305:     MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1306:     MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1307:     MatDestroy(&preallocator);
1308:   }
1309:   /* Fill matrix */
1310:   MatZeroEntries(In);
1311:   for (c = cStart; c < cEnd; ++c) {
1312:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1313:   }
1314:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1315:   PetscFree2(feRef,fvRef);
1316:   PetscFree(elemMat);
1317:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1318:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1319:   if (mesh->printFEM) {
1320:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1321:     MatChop(In, 1.0e-10);
1322:     MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1323:   }
1324:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1325:   return(0);
1326: }

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

1331:   Input Parameters:
1332: + dmf  - The fine mesh
1333: . dmc  - The coarse mesh
1334: - user - The user context

1336:   Output Parameter:
1337: . In  - The interpolation matrix

1339:   Level: developer

1341: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1342: @*/
1343: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1344: {
1345:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
1346:   const char    *name = "Interpolator";
1347:   PetscDS        prob;
1348:   PetscSection   fsection, csection, globalFSection, globalCSection;
1349:   PetscHashJK    ht;
1350:   PetscLayout    rLayout;
1351:   PetscInt      *dnz, *onz;
1352:   PetscInt       locRows, rStart, rEnd;
1353:   PetscReal     *x, *v0, *J, *invJ, detJ;
1354:   PetscReal     *v0c, *Jc, *invJc, detJc;
1355:   PetscScalar   *elemMat;
1356:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

1360:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1361:   DMGetCoordinateDim(dmc, &dim);
1362:   DMGetDS(dmc, &prob);
1363:   PetscDSGetRefCoordArrays(prob, &x, NULL);
1364:   PetscDSGetNumFields(prob, &Nf);
1365:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1366:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1367:   DMGetDefaultSection(dmf, &fsection);
1368:   DMGetDefaultGlobalSection(dmf, &globalFSection);
1369:   DMGetDefaultSection(dmc, &csection);
1370:   DMGetDefaultGlobalSection(dmc, &globalCSection);
1371:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1372:   PetscDSGetTotalDimension(prob, &totDim);
1373:   PetscMalloc1(totDim, &elemMat);

1375:   MatGetLocalSize(In, &locRows, NULL);
1376:   PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1377:   PetscLayoutSetLocalSize(rLayout, locRows);
1378:   PetscLayoutSetBlockSize(rLayout, 1);
1379:   PetscLayoutSetUp(rLayout);
1380:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1381:   PetscLayoutDestroy(&rLayout);
1382:   PetscCalloc2(locRows,&dnz,locRows,&onz);
1383:   PetscHashJKCreate(&ht);
1384:   for (field = 0; field < Nf; ++field) {
1385:     PetscObject      obj;
1386:     PetscClassId     id;
1387:     PetscDualSpace   Q = NULL;
1388:     PetscQuadrature  f;
1389:     const PetscReal *qpoints;
1390:     PetscInt         Nc, Np, fpdim, i, d;

1392:     PetscDSGetDiscretization(prob, field, &obj);
1393:     PetscObjectGetClassId(obj, &id);
1394:     if (id == PETSCFE_CLASSID) {
1395:       PetscFE fe = (PetscFE) obj;

1397:       PetscFEGetDualSpace(fe, &Q);
1398:       PetscFEGetNumComponents(fe, &Nc);
1399:     } else if (id == PETSCFV_CLASSID) {
1400:       PetscFV fv = (PetscFV) obj;

1402:       PetscFVGetDualSpace(fv, &Q);
1403:       Nc   = 1;
1404:     }
1405:     PetscDualSpaceGetDimension(Q, &fpdim);
1406:     /* For each fine grid cell */
1407:     for (cell = cStart; cell < cEnd; ++cell) {
1408:       PetscInt *findices,   *cindices;
1409:       PetscInt  numFIndices, numCIndices;

1411:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1412:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1413:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
1414:       for (i = 0; i < fpdim; ++i) {
1415:         Vec             pointVec;
1416:         PetscScalar    *pV;
1417:         PetscSF         coarseCellSF = NULL;
1418:         const PetscSFNode *coarseCells;
1419:         PetscInt        numCoarseCells, q, c;

1421:         /* Get points from the dual basis functional quadrature */
1422:         PetscDualSpaceGetFunctional(Q, i, &f);
1423:         PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
1424:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1425:         VecSetBlockSize(pointVec, dim);
1426:         VecGetArray(pointVec, &pV);
1427:         for (q = 0; q < Np; ++q) {
1428:           /* Transform point to real space */
1429:           CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1430:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1431:         }
1432:         VecRestoreArray(pointVec, &pV);
1433:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1434:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1435:         PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
1436:         /* Update preallocation info */
1437:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1438:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1439:         {
1440:           PetscHashJKKey  key;
1441:           PetscHashJKIter missing, iter;

1443:           key.j = findices[i];
1444:           if (key.j >= 0) {
1445:             /* Get indices for coarse elements */
1446:             for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1447:               DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1448:               for (c = 0; c < numCIndices; ++c) {
1449:                 key.k = cindices[c];
1450:                 if (key.k < 0) continue;
1451:                 PetscHashJKPut(ht, key, &missing, &iter);
1452:                 if (missing) {
1453:                   PetscHashJKSet(ht, iter, 1);
1454:                   if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1455:                   else                                     ++onz[key.j-rStart];
1456:                 }
1457:               }
1458:               DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1459:             }
1460:           }
1461:         }
1462:         PetscSFDestroy(&coarseCellSF);
1463:         VecDestroy(&pointVec);
1464:       }
1465:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1466:     }
1467:   }
1468:   PetscHashJKDestroy(&ht);
1469:   MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1470:   MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1471:   PetscFree2(dnz,onz);
1472:   for (field = 0; field < Nf; ++field) {
1473:     PetscObject      obj;
1474:     PetscClassId     id;
1475:     PetscDualSpace   Q = NULL;
1476:     PetscQuadrature  f;
1477:     const PetscReal *qpoints, *qweights;
1478:     PetscInt         Nc, qNc, Np, fpdim, i, d;

1480:     PetscDSGetDiscretization(prob, field, &obj);
1481:     PetscObjectGetClassId(obj, &id);
1482:     if (id == PETSCFE_CLASSID) {
1483:       PetscFE fe = (PetscFE) obj;

1485:       PetscFEGetDualSpace(fe, &Q);
1486:       PetscFEGetNumComponents(fe, &Nc);
1487:     } else if (id == PETSCFV_CLASSID) {
1488:       PetscFV fv = (PetscFV) obj;

1490:       PetscFVGetDualSpace(fv, &Q);
1491:       Nc   = 1;
1492:     }
1493:     PetscDualSpaceGetDimension(Q, &fpdim);
1494:     /* For each fine grid cell */
1495:     for (cell = cStart; cell < cEnd; ++cell) {
1496:       PetscInt *findices,   *cindices;
1497:       PetscInt  numFIndices, numCIndices;

1499:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1500:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1501:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
1502:       for (i = 0; i < fpdim; ++i) {
1503:         Vec             pointVec;
1504:         PetscScalar    *pV;
1505:         PetscSF         coarseCellSF = NULL;
1506:         const PetscSFNode *coarseCells;
1507:         PetscInt        numCoarseCells, cpdim, q, c, j;

1509:         /* Get points from the dual basis functional quadrature */
1510:         PetscDualSpaceGetFunctional(Q, i, &f);
1511:         PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
1512:         if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, Nc);
1513:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1514:         VecSetBlockSize(pointVec, dim);
1515:         VecGetArray(pointVec, &pV);
1516:         for (q = 0; q < Np; ++q) {
1517:           /* Transform point to real space */
1518:           CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1519:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1520:         }
1521:         VecRestoreArray(pointVec, &pV);
1522:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1523:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1524:         /* Update preallocation info */
1525:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1526:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1527:         VecGetArray(pointVec, &pV);
1528:         for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1529:           PetscReal pVReal[3];

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

1537:           if (id == PETSCFE_CLASSID) {
1538:             PetscFE    fe = (PetscFE) obj;
1539:             PetscReal *B;

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

1578: /*@
1579:   DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns

1581:   Input Parameters:
1582: + dmc  - The coarse mesh
1583: - dmf  - The fine mesh
1584: - user - The user context

1586:   Output Parameter:
1587: . sc   - The mapping

1589:   Level: developer

1591: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1592: @*/
1593: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1594: {
1595:   PetscDS        prob;
1596:   PetscFE       *feRef;
1597:   PetscFV       *fvRef;
1598:   Vec            fv, cv;
1599:   IS             fis, cis;
1600:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1601:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1602:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;

1606:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1607:   DMGetDimension(dmf, &dim);
1608:   DMGetDefaultSection(dmf, &fsection);
1609:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1610:   DMGetDefaultSection(dmc, &csection);
1611:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1612:   PetscSectionGetNumFields(fsection, &Nf);
1613:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1614:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1615:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1616:   DMGetDS(dmc, &prob);
1617:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1618:   for (f = 0; f < Nf; ++f) {
1619:     PetscObject  obj;
1620:     PetscClassId id;
1621:     PetscInt     fNb = 0, Nc = 0;

1623:     PetscDSGetDiscretization(prob, f, &obj);
1624:     PetscObjectGetClassId(obj, &id);
1625:     if (id == PETSCFE_CLASSID) {
1626:       PetscFE fe = (PetscFE) obj;

1628:       PetscFERefine(fe, &feRef[f]);
1629:       PetscFEGetDimension(feRef[f], &fNb);
1630:       PetscFEGetNumComponents(fe, &Nc);
1631:     } else if (id == PETSCFV_CLASSID) {
1632:       PetscFV        fv = (PetscFV) obj;
1633:       PetscDualSpace Q;

1635:       PetscFVRefine(fv, &fvRef[f]);
1636:       PetscFVGetDualSpace(fvRef[f], &Q);
1637:       PetscDualSpaceGetDimension(Q, &fNb);
1638:       PetscFVGetNumComponents(fv, &Nc);
1639:     }
1640:     fTotDim += fNb*Nc;
1641:   }
1642:   PetscDSGetTotalDimension(prob, &cTotDim);
1643:   PetscMalloc1(cTotDim,&cmap);
1644:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1645:     PetscFE        feC;
1646:     PetscFV        fvC;
1647:     PetscDualSpace QF, QC;
1648:     PetscInt       NcF, NcC, fpdim, cpdim;

1650:     if (feRef[field]) {
1651:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1652:       PetscFEGetNumComponents(feC, &NcC);
1653:       PetscFEGetNumComponents(feRef[field], &NcF);
1654:       PetscFEGetDualSpace(feRef[field], &QF);
1655:       PetscDualSpaceGetDimension(QF, &fpdim);
1656:       PetscFEGetDualSpace(feC, &QC);
1657:       PetscDualSpaceGetDimension(QC, &cpdim);
1658:     } else {
1659:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
1660:       PetscFVGetNumComponents(fvC, &NcC);
1661:       PetscFVGetNumComponents(fvRef[field], &NcF);
1662:       PetscFVGetDualSpace(fvRef[field], &QF);
1663:       PetscDualSpaceGetDimension(QF, &fpdim);
1664:       PetscFVGetDualSpace(fvC, &QC);
1665:       PetscDualSpaceGetDimension(QC, &cpdim);
1666:     }
1667:     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);
1668:     for (c = 0; c < cpdim; ++c) {
1669:       PetscQuadrature  cfunc;
1670:       const PetscReal *cqpoints;
1671:       PetscInt         NpC;
1672:       PetscBool        found = PETSC_FALSE;

1674:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
1675:       PetscQuadratureGetData(cfunc, NULL, NULL, &NpC, &cqpoints, NULL);
1676:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1677:       for (f = 0; f < fpdim; ++f) {
1678:         PetscQuadrature  ffunc;
1679:         const PetscReal *fqpoints;
1680:         PetscReal        sum = 0.0;
1681:         PetscInt         NpF, comp;

1683:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
1684:         PetscQuadratureGetData(ffunc, NULL, NULL, &NpF, &fqpoints, NULL);
1685:         if (NpC != NpF) continue;
1686:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1687:         if (sum > 1.0e-9) continue;
1688:         for (comp = 0; comp < NcC; ++comp) {
1689:           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1690:         }
1691:         found = PETSC_TRUE;
1692:         break;
1693:       }
1694:       if (!found) {
1695:         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1696:         if (fvRef[field]) {
1697:           PetscInt comp;
1698:           for (comp = 0; comp < NcC; ++comp) {
1699:             cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1700:           }
1701:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1702:       }
1703:     }
1704:     offsetC += cpdim*NcC;
1705:     offsetF += fpdim*NcF;
1706:   }
1707:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
1708:   PetscFree2(feRef,fvRef);

1710:   DMGetGlobalVector(dmf, &fv);
1711:   DMGetGlobalVector(dmc, &cv);
1712:   VecGetOwnershipRange(cv, &startC, &endC);
1713:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1714:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1715:   PetscMalloc1(m,&cindices);
1716:   PetscMalloc1(m,&findices);
1717:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1718:   for (c = cStart; c < cEnd; ++c) {
1719:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1720:     for (d = 0; d < cTotDim; ++d) {
1721:       if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
1722:       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]]);
1723:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1724:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1725:     }
1726:   }
1727:   PetscFree(cmap);
1728:   PetscFree2(cellCIndices,cellFIndices);

1730:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1731:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1732:   VecScatterCreate(cv, cis, fv, fis, sc);
1733:   ISDestroy(&cis);
1734:   ISDestroy(&fis);
1735:   DMRestoreGlobalVector(dmf, &fv);
1736:   DMRestoreGlobalVector(dmc, &cv);
1737:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1738:   return(0);
1739: }