Actual source code: plexfem.c

petsc-master 2016-07-27
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
  2: #include <petscsf.h>

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

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

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

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

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

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

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

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

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

116:   Collective on DM

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

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

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

126:   Level: advanced

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

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

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

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

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

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

193:   Level: advanced

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

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

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

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

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

219:   Level: intermediate

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

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

235: PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
236: {
237:   PetscDualSpace *sp, *cellsp;
238:   PetscInt       *numComp;
239:   PetscSection    section;
240:   PetscScalar    *values;
241:   PetscBool      *fieldActive;
242:   PetscInt        numFields, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, cStart, cEnd, cEndInterior, f, d, v, i, comp, maxHeight, h;
243:   PetscErrorCode  ierr;

246:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
247:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
248:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
249:   if (cEnd <= cStart) return(0);
250:   DMGetDimension(dm, &dim);
251:   DMGetCoordinateDim(dm, &dimEmbed);
252:   DMGetDefaultSection(dm, &section);
253:   PetscSectionGetNumFields(section, &numFields);
254:   PetscMalloc2(numFields,&sp,numFields,&numComp);
255:   DMPlexGetMaxProjectionHeight(dm,&maxHeight);
256:   if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
257:   if (maxHeight > 0) {PetscMalloc1(numFields,&cellsp);}
258:   else               {cellsp = sp;}
259:   for (h = 0; h <= maxHeight; h++) {
260:     DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
261:     if (!h) {pStart = cStart; pEnd = cEnd;}
262:     if (pEnd <= pStart) continue;
263:     totDim = 0;
264:     for (f = 0; f < numFields; ++f) {
265:       PetscObject  obj;
266:       PetscClassId id;

268:       DMGetField(dm, f, &obj);
269:       PetscObjectGetClassId(obj, &id);
270:       if (id == PETSCFE_CLASSID) {
271:         PetscFE fe = (PetscFE) obj;

273:         PetscFEGetNumComponents(fe, &numComp[f]);
274:         if (!h) {
275:           PetscFEGetDualSpace(fe, &cellsp[f]);
276:           sp[f] = cellsp[f];
277:         } else {
278:           PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
279:           if (!sp[f]) continue;
280:         }
281:       } else if (id == PETSCFV_CLASSID) {
282:         PetscFV fv = (PetscFV) obj;

284:         PetscFVGetNumComponents(fv, &numComp[f]);
285:         PetscFVGetDualSpace(fv, &cellsp[f]);
286:         sp[f] = cellsp[f];
287:       } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
288:       PetscDualSpaceGetDimension(sp[f], &spDim);
289:       totDim += spDim*numComp[f];
290:     }
291:     DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);
292:     if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
293:     if (!totDim) continue;
294:     DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
295:     DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);
296:     for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
297:     for (i = 0; i < numIds; ++i) {
298:       IS              pointIS;
299:       const PetscInt *points;
300:       PetscInt        n, p;

302:       DMLabelGetStratumIS(label, ids[i], &pointIS);
303:       if (!pointIS) continue; /* No points with that id on this process */
304:       ISGetLocalSize(pointIS, &n);
305:       ISGetIndices(pointIS, &points);
306:       for (p = 0; p < n; ++p) {
307:         const PetscInt    point = points[p];
308:         PetscFECellGeom   geom;

310:         if ((point < pStart) || (point >= pEnd)) continue;
311:         DMPlexComputeCellGeometryFEM(dm, point, NULL, geom.v0, geom.J, NULL, &geom.detJ);
312:         geom.dim      = dim - h;
313:         geom.dimEmbed = dimEmbed;
314:         for (f = 0, v = 0; f < numFields; ++f) {
315:           void * const ctx = ctxs ? ctxs[f] : NULL;

317:           if (!sp[f]) continue;
318:           PetscDualSpaceGetDimension(sp[f], &spDim);
319:           for (d = 0; d < spDim; ++d) {
320:             if (funcs[f]) {
321:               PetscDualSpaceApply(sp[f], d, time, &geom, numComp[f], funcs[f], ctx, &values[v]);
322:               if (ierr) {
323:                 PetscErrorCode ierr2;
324:                 ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2);
325:                 ierr2 = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr2);
326: 
327:               }
328:             } else {
329:               for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
330:             }
331:             v += numComp[f];
332:           }
333:         }
334:         DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);
335:       }
336:       ISRestoreIndices(pointIS, &points);
337:       ISDestroy(&pointIS);
338:     }
339:     DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
340:     DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);
341:   }
342:   PetscFree2(sp, numComp);
343:   if (maxHeight > 0) {
344:     PetscFree(cellsp);
345:   }
346:   return(0);
347: }

351: PetscErrorCode DMProjectFunctionLocal_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
352: {
353:   PetscDualSpace *sp, *cellsp;
354:   PetscInt       *numComp;
355:   PetscSection    section;
356:   PetscScalar    *values;
357:   PetscInt        Nf, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight;
358:   PetscBool      *isFE, hasFE = PETSC_FALSE, hasFV = PETSC_FALSE;
359:   PetscErrorCode  ierr;

362:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
363:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
364:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
365:   DMGetDefaultSection(dm, &section);
366:   PetscSectionGetNumFields(section, &Nf);
367:   PetscMalloc3(Nf, &isFE, Nf, &sp, Nf, &numComp);
368:   DMGetDimension(dm, &dim);
369:   DMGetCoordinateDim(dm, &dimEmbed);
370:   DMPlexGetMaxProjectionHeight(dm,&maxHeight);
371:   if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
372:   if (maxHeight > 0) {
373:     PetscMalloc1(Nf,&cellsp);
374:   }
375:   else {
376:     cellsp = sp;
377:   }
378:   for (f = 0; f < Nf; ++f) {
379:     PetscObject  obj;
380:     PetscClassId id;

382:     DMGetField(dm, f, &obj);
383:     PetscObjectGetClassId(obj, &id);
384:     if (id == PETSCFE_CLASSID) {
385:       PetscFE fe = (PetscFE) obj;

387:       hasFE   = PETSC_TRUE;
388:       isFE[f] = PETSC_TRUE;
389:       PetscFEGetNumComponents(fe, &numComp[f]);
390:       PetscFEGetDualSpace(fe, &cellsp[f]);
391:     } else if (id == PETSCFV_CLASSID) {
392:       PetscFV fv = (PetscFV) obj;

394:       hasFV   = PETSC_TRUE;
395:       isFE[f] = PETSC_FALSE;
396:       PetscFVGetNumComponents(fv, &numComp[f]);
397:       PetscFVGetDualSpace(fv, &cellsp[f]);
398:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
399:   }
400:   for (h = 0; h <= maxHeight; h++) {
401:     DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
402:     if (!h) {pStart = cStart; pEnd = cEnd;}
403:     if (pEnd <= pStart) continue;
404:     totDim = 0;
405:     for (f = 0; f < Nf; ++f) {
406:       if (!h) {
407:         sp[f] = cellsp[f];
408:       }
409:       else {
410:         PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
411:         if (!sp[f]) {
412:           continue;
413:         }
414:       }
415:       PetscDualSpaceGetDimension(sp[f], &spDim);
416:       totDim += spDim*numComp[f];
417:     }
418:     DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);
419:     if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
420:     if (!totDim) continue;
421:     DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
422:     for (p = pStart; p < pEnd; ++p) {
423:       PetscFECellGeom fegeom;
424:       PetscFVCellGeom fvgeom;

426:       if (hasFE) {
427:         DMPlexComputeCellGeometryFEM(dm, p, NULL, fegeom.v0, fegeom.J, NULL, &fegeom.detJ);
428:         fegeom.dim      = dim - h;
429:         fegeom.dimEmbed = dimEmbed;
430:       }
431:       if (hasFV) {DMPlexComputeCellGeometryFVM(dm, p, &fvgeom.volume, fvgeom.centroid, NULL);}
432:       for (f = 0, v = 0; f < Nf; ++f) {
433:         void * const ctx = ctxs ? ctxs[f] : NULL;

435:         if (!sp[f]) continue;
436:         PetscDualSpaceGetDimension(sp[f], &spDim);
437:         for (d = 0; d < spDim; ++d) {
438:           if (funcs[f]) {
439:             if (isFE[f]) {PetscDualSpaceApply(sp[f], d, time, &fegeom, numComp[f], funcs[f], ctx, &values[v]);}
440:             else         {PetscDualSpaceApplyFVM(sp[f], d, time, &fvgeom, numComp[f], funcs[f], ctx, &values[v]);}
441:             if (ierr) {
442:               PetscErrorCode ierr2;
443:               ierr2 = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr2);
444: 
445:             }
446:           } else {
447:             for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
448:           }
449:           v += numComp[f];
450:         }
451:       }
452:       DMPlexVecSetClosure(dm, section, localX, p, values, mode);
453:     }
454:     DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
455:   }
456:   PetscFree3(isFE, sp, numComp);
457:   if (maxHeight > 0) {
458:     PetscFree(cellsp);
459:   }
460:   return(0);
461: }

465: PetscErrorCode DMProjectFieldLocal_Plex(DM dm, Vec localU,
466:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
467:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
468:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
469:                                                        PetscReal, const PetscReal[], PetscScalar[]),
470:                                         InsertMode mode, Vec localX)
471: {
472:   DM              dmAux;
473:   PetscDS         prob, probAux = NULL;
474:   Vec             A;
475:   PetscSection    section, sectionAux = NULL;
476:   PetscDualSpace *sp;
477:   PetscInt       *Ncf;
478:   PetscScalar    *values, *u, *u_x, *a, *a_x;
479:   PetscReal      *x, *v0, *J, *invJ, detJ;
480:   PetscInt       *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
481:   PetscInt        Nf, NfAux = 0, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight;
482:   PetscErrorCode  ierr;

485:   DMPlexGetMaxProjectionHeight(dm,&maxHeight);
486:   if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");}
487:   DMGetDS(dm, &prob);
488:   DMGetDimension(dm, &dim);
489:   DMGetDefaultSection(dm, &section);
490:   PetscSectionGetNumFields(section, &Nf);
491:   PetscMalloc2(Nf, &sp, Nf, &Ncf);
492:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
493:   PetscDSGetTotalDimension(prob, &totDim);
494:   PetscDSGetComponentOffsets(prob, &uOff);
495:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
496:   PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
497:   PetscDSGetRefCoordArrays(prob, &x, NULL);
498:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
499:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
500:   if (dmAux) {
501:     DMGetDS(dmAux, &probAux);
502:     PetscDSGetNumFields(probAux, &NfAux);
503:     DMGetDefaultSection(dmAux, &sectionAux);
504:     PetscDSGetComponentOffsets(probAux, &aOff);
505:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
506:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
507:   }
508:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localU, 0.0, NULL, NULL, NULL);
509:   DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
510:   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
511:   DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
512:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
513:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
514:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
515:   for (c = cStart; c < cEnd; ++c) {
516:     PetscScalar *coefficients = NULL, *coefficientsAux = NULL;

518:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
519:     DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);
520:     if (dmAux) {DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
521:     for (f = 0, v = 0; f < Nf; ++f) {
522:       PetscObject  obj;
523:       PetscClassId id;

525:       PetscDSGetDiscretization(prob, f, &obj);
526:       PetscObjectGetClassId(obj, &id);
527:       if (id == PETSCFE_CLASSID) {
528:         PetscFE fe = (PetscFE) obj;

530:         PetscFEGetDualSpace(fe, &sp[f]);
531:         PetscFEGetNumComponents(fe, &Ncf[f]);
532:       } else if (id == PETSCFV_CLASSID) {
533:         PetscFV fv = (PetscFV) obj;

535:         PetscFVGetNumComponents(fv, &Ncf[f]);
536:         PetscFVGetDualSpace(fv, &sp[f]);
537:       }
538:       PetscDualSpaceGetDimension(sp[f], &spDim);
539:       for (d = 0; d < spDim; ++d) {
540:         PetscQuadrature  quad;
541:         const PetscReal *points, *weights;
542:         PetscInt         numPoints, q;

544:         if (funcs[f]) {
545:           PetscDualSpaceGetFunctional(sp[f], d, &quad);
546:           PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);
547:           for (q = 0; q < numPoints; ++q) {
548:             CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x);
549:             EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, coefficients,    NULL, u, u_x, NULL);
550:             EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);
551:             (*funcs[f])(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &values[v]);
552:           }
553:         } else {
554:           for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0;
555:         }
556:         v += Ncf[f];
557:       }
558:     }
559:     DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);
560:     if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
561:     DMPlexVecSetClosure(dm, section, localX, c, values, mode);
562:   }
563:   PetscFree3(v0,J,invJ);
564:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
565:   PetscFree2(sp, Ncf);
566:   return(0);
567: }

571: 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)
572: {
573:   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
574:   void            **ctxs;
575:   PetscInt          numFields;
576:   PetscErrorCode    ierr;

579:   DMGetNumFields(dm, &numFields);
580:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
581:   funcs[field] = func;
582:   ctxs[field]  = ctx;
583:   DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);
584:   PetscFree2(funcs,ctxs);
585:   return(0);
586: }

590: /* This ignores numcomps/comps */
591: static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad,
592:                                                               PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
593: {
594:   PetscDS            prob;
595:   PetscSF            sf;
596:   DM                 dmFace, dmCell, dmGrad;
597:   const PetscScalar *facegeom, *cellgeom = NULL, *grad;
598:   const PetscInt    *leaves;
599:   PetscScalar       *x, *fx;
600:   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
601:   PetscErrorCode     ierr, ierru = 0;

604:   DMGetPointSF(dm, &sf);
605:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
606:   nleaves = PetscMax(0, nleaves);
607:   DMGetDimension(dm, &dim);
608:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
609:   DMGetDS(dm, &prob);
610:   VecGetDM(faceGeometry, &dmFace);
611:   VecGetArrayRead(faceGeometry, &facegeom);
612:   if (cellGeometry) {
613:     VecGetDM(cellGeometry, &dmCell);
614:     VecGetArrayRead(cellGeometry, &cellgeom);
615:   }
616:   if (Grad) {
617:     PetscFV fv;

619:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
620:     VecGetDM(Grad, &dmGrad);
621:     VecGetArrayRead(Grad, &grad);
622:     PetscFVGetNumComponents(fv, &pdim);
623:     DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
624:   }
625:   VecGetArray(locX, &x);
626:   for (i = 0; i < numids; ++i) {
627:     IS              faceIS;
628:     const PetscInt *faces;
629:     PetscInt        numFaces, f;

631:     DMLabelGetStratumIS(label, ids[i], &faceIS);
632:     if (!faceIS) continue; /* No points with that id on this process */
633:     ISGetLocalSize(faceIS, &numFaces);
634:     ISGetIndices(faceIS, &faces);
635:     for (f = 0; f < numFaces; ++f) {
636:       const PetscInt         face = faces[f], *cells;
637:       PetscFVFaceGeom        *fg;

639:       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
640:       PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
641:       if (loc >= 0) continue;
642:       DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
643:       DMPlexGetSupport(dm, face, &cells);
644:       if (Grad) {
645:         PetscFVCellGeom       *cg;
646:         PetscScalar           *cx, *cgrad;
647:         PetscScalar           *xG;
648:         PetscReal              dx[3];
649:         PetscInt               d;

651:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
652:         DMPlexPointLocalRead(dm, cells[0], x, &cx);
653:         DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
654:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
655:         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
656:         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
657:         ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
658:         if (ierru) {
659:           ISRestoreIndices(faceIS, &faces);
660:           ISDestroy(&faceIS);
661:           goto cleanup;
662:         }
663:       } else {
664:         PetscScalar       *xI;
665:         PetscScalar       *xG;

667:         DMPlexPointLocalRead(dm, cells[0], x, &xI);
668:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
669:         ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
670:         if (ierru) {
671:           ISRestoreIndices(faceIS, &faces);
672:           ISDestroy(&faceIS);
673:           goto cleanup;
674:         }
675:       }
676:     }
677:     ISRestoreIndices(faceIS, &faces);
678:     ISDestroy(&faceIS);
679:   }
680:   cleanup:
681:   VecRestoreArray(locX, &x);
682:   if (Grad) {
683:     DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);
684:     VecRestoreArrayRead(Grad, &grad);
685:   }
686:   if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
687:   VecRestoreArrayRead(faceGeometry, &facegeom);
688:   CHKERRQ(ierru);
689:   return(0);
690: }

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

697:   Input Parameters:
698: + dm - The DM
699: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
700: . time - The time
701: . faceGeomFVM - Face geometry data for FV discretizations
702: . cellGeomFVM - Cell geometry data for FV discretizations
703: - gradFVM - Gradient reconstruction data for FV discretizations

705:   Output Parameters:
706: . locX - Solution updated with boundary values

708:   Level: developer

710: .seealso: DMProjectFunctionLabelLocal()
711: @*/
712: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
713: {
714:   PetscInt       numBd, b;

723:   PetscDSGetNumBoundary(dm->prob, &numBd);
724:   for (b = 0; b < numBd; ++b) {
725:     PetscBool       isEssential;
726:     const char     *labelname;
727:     DMLabel         label;
728:     PetscInt        field;
729:     PetscObject     obj;
730:     PetscClassId    id;
731:     void          (*func)();
732:     PetscInt        numids;
733:     const PetscInt *ids;
734:     void           *ctx;

736:     DMGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);
737:     if (insertEssential != isEssential) continue;
738:     DMGetLabel(dm, labelname, &label);
739:     DMGetField(dm, field, &obj);
740:     PetscObjectGetClassId(obj, &id);
741:     if (id == PETSCFE_CLASSID) {
742:       if (!isEssential) continue; /* for FEM, there is no insertion to be done for non-essential boundary conditions */
743:       DMPlexLabelAddCells(dm,label);
744:       DMPlexInsertBoundaryValues_FEM_Internal(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
745:       DMPlexLabelClearCells(dm,label);
746:     } else if (id == PETSCFV_CLASSID) {
747:       if (!faceGeomFVM) continue;
748:       DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM,
749:                                                      field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
750:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
751:   }
752:   return(0);
753: }

757: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
758: {
759:   const PetscInt   debug = 0;
760:   PetscSection     section;
761:   PetscQuadrature  quad;
762:   Vec              localX;
763:   PetscScalar     *funcVal, *interpolant;
764:   PetscReal       *coords, *v0, *J, *invJ, detJ;
765:   PetscReal        localDiff = 0.0;
766:   const PetscReal *quadPoints, *quadWeights;
767:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
768:   PetscErrorCode   ierr;

771:   DMGetDimension(dm, &dim);
772:   DMGetDefaultSection(dm, &section);
773:   PetscSectionGetNumFields(section, &numFields);
774:   DMGetLocalVector(dm, &localX);
775:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
776:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
777:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
778:   for (field = 0; field < numFields; ++field) {
779:     PetscObject  obj;
780:     PetscClassId id;
781:     PetscInt     Nc;

783:     DMGetField(dm, field, &obj);
784:     PetscObjectGetClassId(obj, &id);
785:     if (id == PETSCFE_CLASSID) {
786:       PetscFE fe = (PetscFE) obj;

788:       PetscFEGetQuadrature(fe, &quad);
789:       PetscFEGetNumComponents(fe, &Nc);
790:     } else if (id == PETSCFV_CLASSID) {
791:       PetscFV fv = (PetscFV) obj;

793:       PetscFVGetQuadrature(fv, &quad);
794:       PetscFVGetNumComponents(fv, &Nc);
795:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
796:     numComponents += Nc;
797:   }
798:   PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
799:   PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
800:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
801:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
802:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
803:   for (c = cStart; c < cEnd; ++c) {
804:     PetscScalar *x = NULL;
805:     PetscReal    elemDiff = 0.0;

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

811:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
812:       PetscObject  obj;
813:       PetscClassId id;
814:       void * const ctx = ctxs ? ctxs[field] : NULL;
815:       PetscInt     Nb, Nc, q, fc;

817:       DMGetField(dm, field, &obj);
818:       PetscObjectGetClassId(obj, &id);
819:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
820:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
821:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
822:       if (debug) {
823:         char title[1024];
824:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
825:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
826:       }
827:       for (q = 0; q < numQuadPoints; ++q) {
828:         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
829:         (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);
830:         if (ierr) {
831:           PetscErrorCode ierr2;
832:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
833:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
834:           ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2);
835: 
836:         }
837:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
838:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
839:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
840:         for (fc = 0; fc < Nc; ++fc) {
841:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);}
842:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
843:         }
844:       }
845:       fieldOffset += Nb*Nc;
846:     }
847:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
848:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
849:     localDiff += elemDiff;
850:   }
851:   PetscFree6(funcVal,interpolant,coords,v0,J,invJ);
852:   DMRestoreLocalVector(dm, &localX);
853:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
854:   *diff = PetscSqrtReal(*diff);
855:   return(0);
856: }

860: 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)
861: {
862:   const PetscInt  debug = 0;
863:   PetscSection    section;
864:   PetscQuadrature quad;
865:   Vec             localX;
866:   PetscScalar    *funcVal, *interpolantVec;
867:   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
868:   PetscReal       localDiff = 0.0;
869:   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp;
870:   PetscErrorCode  ierr;

873:   DMGetDimension(dm, &dim);
874:   DMGetDefaultSection(dm, &section);
875:   PetscSectionGetNumFields(section, &numFields);
876:   DMGetLocalVector(dm, &localX);
877:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
878:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
879:   for (field = 0; field < numFields; ++field) {
880:     PetscFE  fe;
881:     PetscInt Nc;

883:     DMGetField(dm, field, (PetscObject *) &fe);
884:     PetscFEGetQuadrature(fe, &quad);
885:     PetscFEGetNumComponents(fe, &Nc);
886:     numComponents += Nc;
887:   }
888:   /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
889:   PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
890:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
891:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
892:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
893:   for (c = cStart; c < cEnd; ++c) {
894:     PetscScalar *x = NULL;
895:     PetscReal    elemDiff = 0.0;

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

901:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
902:       PetscFE          fe;
903:       void * const     ctx = ctxs ? ctxs[field] : NULL;
904:       const PetscReal *quadPoints, *quadWeights;
905:       PetscReal       *basisDer;
906:       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;

908:       DMGetField(dm, field, (PetscObject *) &fe);
909:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
910:       PetscFEGetDimension(fe, &Nb);
911:       PetscFEGetNumComponents(fe, &Ncomp);
912:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
913:       if (debug) {
914:         char title[1024];
915:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
916:         DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);
917:       }
918:       for (q = 0; q < numQuadPoints; ++q) {
919:         for (d = 0; d < dim; d++) {
920:           coords[d] = v0[d];
921:           for (e = 0; e < dim; e++) {
922:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
923:           }
924:         }
925:         (*funcs[field])(dim, time, coords, n, numFields, funcVal, ctx);
926:         if (ierr) {
927:           PetscErrorCode ierr2;
928:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
929:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
930:           ierr2 = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr2);
931: 
932:         }
933:         for (fc = 0; fc < Ncomp; ++fc) {
934:           PetscScalar interpolant = 0.0;

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

940:             for (d = 0; d < dim; ++d) {
941:               realSpaceDer[d] = 0.0;
942:               for (g = 0; g < dim; ++g) {
943:                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
944:               }
945:               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
946:             }
947:           }
948:           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
949:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
950:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
951:         }
952:       }
953:       comp        += Ncomp;
954:       fieldOffset += Nb*Ncomp;
955:     }
956:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
957:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
958:     localDiff += elemDiff;
959:   }
960:   PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
961:   DMRestoreLocalVector(dm, &localX);
962:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
963:   *diff = PetscSqrtReal(*diff);
964:   return(0);
965: }

969: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
970: {
971:   const PetscInt   debug = 0;
972:   PetscSection     section;
973:   PetscQuadrature  quad;
974:   Vec              localX;
975:   PetscScalar     *funcVal, *interpolant;
976:   PetscReal       *coords, *v0, *J, *invJ, detJ;
977:   PetscReal       *localDiff;
978:   const PetscReal *quadPoints, *quadWeights;
979:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
980:   PetscErrorCode   ierr;

983:   DMGetDimension(dm, &dim);
984:   DMGetDefaultSection(dm, &section);
985:   PetscSectionGetNumFields(section, &numFields);
986:   DMGetLocalVector(dm, &localX);
987:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
988:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
989:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
990:   for (field = 0; field < numFields; ++field) {
991:     PetscObject  obj;
992:     PetscClassId id;
993:     PetscInt     Nc;

995:     DMGetField(dm, field, &obj);
996:     PetscObjectGetClassId(obj, &id);
997:     if (id == PETSCFE_CLASSID) {
998:       PetscFE fe = (PetscFE) obj;

1000:       PetscFEGetQuadrature(fe, &quad);
1001:       PetscFEGetNumComponents(fe, &Nc);
1002:     } else if (id == PETSCFV_CLASSID) {
1003:       PetscFV fv = (PetscFV) obj;

1005:       PetscFVGetQuadrature(fv, &quad);
1006:       PetscFVGetNumComponents(fv, &Nc);
1007:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1008:     numComponents += Nc;
1009:   }
1010:   PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1011:   PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1012:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1013:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1014:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1015:   for (c = cStart; c < cEnd; ++c) {
1016:     PetscScalar *x = NULL;

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

1022:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1023:       PetscObject  obj;
1024:       PetscClassId id;
1025:       void * const ctx = ctxs ? ctxs[field] : NULL;
1026:       PetscInt     Nb, Nc, q, fc;

1028:       PetscReal       elemDiff = 0.0;

1030:       DMGetField(dm, field, &obj);
1031:       PetscObjectGetClassId(obj, &id);
1032:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1033:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1034:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1035:       if (debug) {
1036:         char title[1024];
1037:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1038:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
1039:       }
1040:       for (q = 0; q < numQuadPoints; ++q) {
1041:         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
1042:         (*funcs[field])(dim, time, coords, numFields, funcVal, ctx);
1043:         if (ierr) {
1044:           PetscErrorCode ierr2;
1045:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1046:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1047:           ierr2 = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2);
1048: 
1049:         }
1050:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
1051:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1052:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1053:         for (fc = 0; fc < Nc; ++fc) {
1054:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d point %g %g %g diff %g\n", c, field, dim > 0 ? coords[0] : 0., dim > 1 ? coords[1] : 0., dim > 2 ? coords[2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);}
1055:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
1056:         }
1057:       }
1058:       fieldOffset += Nb*Nc;
1059:       localDiff[field] += elemDiff;
1060:     }
1061:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1062:   }
1063:   DMRestoreLocalVector(dm, &localX);
1064:   MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1065:   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
1066:   PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);
1067:   return(0);
1068: }

1072: /*@C
1073:   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.

1075:   Input Parameters:
1076: + dm    - The DM
1077: . time  - The time
1078: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
1079: . ctxs  - Optional array of contexts to pass to each function, or NULL.
1080: - X     - The coefficient vector u_h

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

1085:   Level: developer

1087: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1088: @*/
1089: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
1090: {
1091:   PetscSection     section;
1092:   PetscQuadrature  quad;
1093:   Vec              localX;
1094:   PetscScalar     *funcVal, *interpolant;
1095:   PetscReal       *coords, *v0, *J, *invJ, detJ;
1096:   const PetscReal *quadPoints, *quadWeights;
1097:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1098:   PetscErrorCode   ierr;

1101:   VecSet(D, 0.0);
1102:   DMGetDimension(dm, &dim);
1103:   DMGetDefaultSection(dm, &section);
1104:   PetscSectionGetNumFields(section, &numFields);
1105:   DMGetLocalVector(dm, &localX);
1106:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1107:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1108:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1109:   for (field = 0; field < numFields; ++field) {
1110:     PetscObject  obj;
1111:     PetscClassId id;
1112:     PetscInt     Nc;

1114:     DMGetField(dm, field, &obj);
1115:     PetscObjectGetClassId(obj, &id);
1116:     if (id == PETSCFE_CLASSID) {
1117:       PetscFE fe = (PetscFE) obj;

1119:       PetscFEGetQuadrature(fe, &quad);
1120:       PetscFEGetNumComponents(fe, &Nc);
1121:     } else if (id == PETSCFV_CLASSID) {
1122:       PetscFV fv = (PetscFV) obj;

1124:       PetscFVGetQuadrature(fv, &quad);
1125:       PetscFVGetNumComponents(fv, &Nc);
1126:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1127:     numComponents += Nc;
1128:   }
1129:   PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1130:   PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1131:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1132:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1133:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1134:   for (c = cStart; c < cEnd; ++c) {
1135:     PetscScalar *x = NULL;
1136:     PetscScalar  elemDiff = 0.0;

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

1142:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1143:       PetscObject  obj;
1144:       PetscClassId id;
1145:       void * const ctx = ctxs ? ctxs[field] : NULL;
1146:       PetscInt     Nb, Nc, q, fc;

1148:       DMGetField(dm, field, &obj);
1149:       PetscObjectGetClassId(obj, &id);
1150:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1151:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1152:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1153:       if (funcs[field]) {
1154:         for (q = 0; q < numQuadPoints; ++q) {
1155:           CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
1156:           (*funcs[field])(dim, time, coords, Nc, funcVal, ctx);
1157:           if (ierr) {
1158:             PetscErrorCode ierr2;
1159:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1160:             ierr2 = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr2);
1161:             ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1162: 
1163:           }
1164:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
1165:           else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1166:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1167:           for (fc = 0; fc < Nc; ++fc) {
1168:             elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
1169:           }
1170:         }
1171:       }
1172:       fieldOffset += Nb*Nc;
1173:     }
1174:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1175:     VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1176:   }
1177:   PetscFree6(funcVal,interpolant,coords,v0,J,invJ);
1178:   DMRestoreLocalVector(dm, &localX);
1179:   VecSqrtAbs(D);
1180:   return(0);
1181: }

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

1188:   Input Parameters:
1189: + dm - The mesh
1190: . X  - Local input vector
1191: - user - The user context

1193:   Output Parameter:
1194: . integral - Local integral for each field

1196:   Level: developer

1198: .seealso: DMPlexComputeResidualFEM()
1199: @*/
1200: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
1201: {
1202:   DM_Plex           *mesh  = (DM_Plex *) dm->data;
1203:   DM                 dmAux, dmGrad;
1204:   Vec                localX, A, cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
1205:   PetscDS            prob, probAux = NULL;
1206:   PetscSection       section, sectionAux;
1207:   PetscFV            fvm = NULL;
1208:   PetscFECellGeom   *cgeomFEM;
1209:   PetscFVFaceGeom   *fgeomFVM;
1210:   PetscFVCellGeom   *cgeomFVM;
1211:   PetscScalar       *u, *a = NULL;
1212:   const PetscScalar *lgrad;
1213:   PetscReal         *lintegral;
1214:   PetscInt          *uOff, *uOff_x, *aOff = NULL;
1215:   PetscInt           dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, cEndInterior, c;
1216:   PetscInt           totDim, totDimAux;
1217:   PetscBool          useFVM = PETSC_FALSE;
1218:   PetscErrorCode     ierr;

1221:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1222:   DMGetLocalVector(dm, &localX);
1223:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);
1224:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1225:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1226:   DMGetDimension(dm, &dim);
1227:   DMGetDefaultSection(dm, &section);
1228:   DMGetDS(dm, &prob);
1229:   PetscDSGetTotalDimension(prob, &totDim);
1230:   PetscDSGetComponentOffsets(prob, &uOff);
1231:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1232:   PetscSectionGetNumFields(section, &Nf);
1233:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1234:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1235:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1236:   numCells = cEnd - cStart;
1237:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1238:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1239:   if (dmAux) {
1240:     DMGetDS(dmAux, &probAux);
1241:     PetscDSGetNumFields(probAux, &NfAux);
1242:     DMGetDefaultSection(dmAux, &sectionAux);
1243:     PetscDSGetTotalDimension(probAux, &totDimAux);
1244:     PetscDSGetComponentOffsets(probAux, &aOff);
1245:     PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);
1246:   }
1247:   for (f = 0; f < Nf; ++f) {
1248:     PetscObject  obj;
1249:     PetscClassId id;

1251:     PetscDSGetDiscretization(prob, f, &obj);
1252:     PetscObjectGetClassId(obj, &id);
1253:     if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1254:   }
1255:   if (useFVM) {
1256:     Vec       grad;
1257:     PetscInt  fStart, fEnd;
1258:     PetscBool compGrad;

1260:     PetscFVGetComputeGradients(fvm, &compGrad);
1261:     PetscFVSetComputeGradients(fvm, PETSC_TRUE);
1262:     DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1263:     DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1264:     PetscFVSetComputeGradients(fvm, compGrad);
1265:     VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1266:     VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1267:     /* Reconstruct and limit cell gradients */
1268:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1269:     DMGetGlobalVector(dmGrad, &grad);
1270:     DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);
1271:     /* Communicate gradient values */
1272:     DMGetLocalVector(dmGrad, &locGrad);
1273:     DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1274:     DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1275:     DMRestoreGlobalVector(dmGrad, &grad);
1276:     /* Handle non-essential (e.g. outflow) boundary values */
1277:     DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1278:     VecGetArrayRead(locGrad, &lgrad);
1279:   }
1280:   PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);
1281:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1282:   for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
1283:   for (c = cStart; c < cEnd; ++c) {
1284:     PetscScalar *x = NULL;
1285:     PetscInt     i;

1287:     DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);
1288:     if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c);
1289:     DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
1290:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1291:     DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
1292:     if (dmAux) {
1293:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1294:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1295:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1296:     }
1297:   }
1298:   for (f = 0; f < Nf; ++f) {
1299:     PetscObject  obj;
1300:     PetscClassId id;
1301:     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1303:     PetscDSGetDiscretization(prob, f, &obj);
1304:     PetscObjectGetClassId(obj, &id);
1305:     if (id == PETSCFE_CLASSID) {
1306:       PetscFE         fe = (PetscFE) obj;
1307:       PetscQuadrature q;
1308:       PetscInt        Nq, Nb;

1310:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1311:       PetscFEGetQuadrature(fe, &q);
1312:       PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1313:       PetscFEGetDimension(fe, &Nb);
1314:       blockSize = Nb*Nq;
1315:       batchSize = numBlocks * blockSize;
1316:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1317:       numChunks = numCells / (numBatches*batchSize);
1318:       Ne        = numChunks*numBatches*batchSize;
1319:       Nr        = numCells % (numBatches*batchSize);
1320:       offset    = numCells - Nr;
1321:       PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);
1322:       PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);
1323:     } else if (id == PETSCFV_CLASSID) {
1324:       /* PetscFV  fv = (PetscFV) obj; */
1325:       PetscInt       foff;
1326:       PetscPointFunc obj_func;
1327:       PetscScalar    lint;

1329:       PetscDSGetObjective(prob, f, &obj_func);
1330:       PetscDSGetFieldOffset(prob, f, &foff);
1331:       if (obj_func) {
1332:         for (c = 0; c < numCells; ++c) {
1333:           PetscScalar *u_x;

1335:           DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1336:           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);
1337:           lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1338:         }
1339:       }
1340:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1341:   }
1342:   if (useFVM) {
1343:     VecRestoreArrayRead(locGrad, &lgrad);
1344:     VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1345:     VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1346:     DMRestoreLocalVector(dmGrad, &locGrad);
1347:     VecDestroy(&faceGeometryFVM);
1348:     VecDestroy(&cellGeometryFVM);
1349:     DMDestroy(&dmGrad);
1350:   }
1351:   if (dmAux) {PetscFree(a);}
1352:   if (mesh->printFEM) {
1353:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1354:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);}
1355:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1356:   }
1357:   DMRestoreLocalVector(dm, &localX);
1358:   MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1359:   PetscFree3(lintegral,u,cgeomFEM);
1360:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1361:   return(0);
1362: }

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

1369:   Input Parameters:
1370: + dmf  - The fine mesh
1371: . dmc  - The coarse mesh
1372: - user - The user context

1374:   Output Parameter:
1375: . In  - The interpolation matrix

1377:   Level: developer

1379: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1380: @*/
1381: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1382: {
1383:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1384:   const char       *name  = "Interpolator";
1385:   PetscDS           prob;
1386:   PetscFE          *feRef;
1387:   PetscFV          *fvRef;
1388:   PetscSection      fsection, fglobalSection;
1389:   PetscSection      csection, cglobalSection;
1390:   PetscScalar      *elemMat;
1391:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1392:   PetscInt          cTotDim, rTotDim = 0;
1393:   PetscErrorCode    ierr;

1396:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1397:   DMGetDimension(dmf, &dim);
1398:   DMGetDefaultSection(dmf, &fsection);
1399:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1400:   DMGetDefaultSection(dmc, &csection);
1401:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1402:   PetscSectionGetNumFields(fsection, &Nf);
1403:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1404:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1405:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1406:   DMGetDS(dmf, &prob);
1407:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1408:   for (f = 0; f < Nf; ++f) {
1409:     PetscObject  obj;
1410:     PetscClassId id;
1411:     PetscInt     rNb = 0, Nc = 0;

1413:     PetscDSGetDiscretization(prob, f, &obj);
1414:     PetscObjectGetClassId(obj, &id);
1415:     if (id == PETSCFE_CLASSID) {
1416:       PetscFE fe = (PetscFE) obj;

1418:       PetscFERefine(fe, &feRef[f]);
1419:       PetscFEGetDimension(feRef[f], &rNb);
1420:       PetscFEGetNumComponents(fe, &Nc);
1421:     } else if (id == PETSCFV_CLASSID) {
1422:       PetscFV        fv = (PetscFV) obj;
1423:       PetscDualSpace Q;

1425:       PetscFVRefine(fv, &fvRef[f]);
1426:       PetscFVGetDualSpace(fvRef[f], &Q);
1427:       PetscDualSpaceGetDimension(Q, &rNb);
1428:       PetscFVGetNumComponents(fv, &Nc);
1429:     }
1430:     rTotDim += rNb*Nc;
1431:   }
1432:   PetscDSGetTotalDimension(prob, &cTotDim);
1433:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1434:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1435:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1436:     PetscDualSpace   Qref;
1437:     PetscQuadrature  f;
1438:     const PetscReal *qpoints, *qweights;
1439:     PetscReal       *points;
1440:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1442:     /* Compose points from all dual basis functionals */
1443:     if (feRef[fieldI]) {
1444:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
1445:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
1446:     } else {
1447:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1448:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1449:     }
1450:     PetscDualSpaceGetDimension(Qref, &fpdim);
1451:     for (i = 0; i < fpdim; ++i) {
1452:       PetscDualSpaceGetFunctional(Qref, i, &f);
1453:       PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1454:       npoints += Np;
1455:     }
1456:     PetscMalloc1(npoints*dim,&points);
1457:     for (i = 0, k = 0; i < fpdim; ++i) {
1458:       PetscDualSpaceGetFunctional(Qref, i, &f);
1459:       PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1460:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1461:     }

1463:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1464:       PetscObject  obj;
1465:       PetscClassId id;
1466:       PetscReal   *B;
1467:       PetscInt     NcJ = 0, cpdim = 0, j;

1469:       PetscDSGetDiscretization(prob, fieldJ, &obj);
1470:       PetscObjectGetClassId(obj, &id);
1471:       if (id == PETSCFE_CLASSID) {
1472:         PetscFE fe = (PetscFE) obj;

1474:         /* Evaluate basis at points */
1475:         PetscFEGetNumComponents(fe, &NcJ);
1476:         PetscFEGetDimension(fe, &cpdim);
1477:         /* For now, fields only interpolate themselves */
1478:         if (fieldI == fieldJ) {
1479:           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);
1480:           PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1481:           for (i = 0, k = 0; i < fpdim; ++i) {
1482:             PetscDualSpaceGetFunctional(Qref, i, &f);
1483:             PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1484:             for (p = 0; p < Np; ++p, ++k) {
1485:               for (j = 0; j < cpdim; ++j) {
1486:                 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];
1487:               }
1488:             }
1489:           }
1490:           PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1491:         }
1492:       } else if (id == PETSCFV_CLASSID) {
1493:         PetscFV        fv = (PetscFV) obj;

1495:         /* Evaluate constant function at points */
1496:         PetscFVGetNumComponents(fv, &NcJ);
1497:         cpdim = 1;
1498:         /* For now, fields only interpolate themselves */
1499:         if (fieldI == fieldJ) {
1500:           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);
1501:           for (i = 0, k = 0; i < fpdim; ++i) {
1502:             PetscDualSpaceGetFunctional(Qref, i, &f);
1503:             PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1504:             for (p = 0; p < Np; ++p, ++k) {
1505:               for (j = 0; j < cpdim; ++j) {
1506:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p];
1507:               }
1508:             }
1509:           }
1510:         }
1511:       }
1512:       offsetJ += cpdim*NcJ;
1513:     }
1514:     offsetI += fpdim*Nc;
1515:     PetscFree(points);
1516:   }
1517:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1518:   /* Preallocate matrix */
1519:   {
1520:     Mat          preallocator;
1521:     PetscScalar *vals;
1522:     PetscInt    *cellCIndices, *cellFIndices;
1523:     PetscInt     locRows, locCols, cell;

1525:     MatGetLocalSize(In, &locRows, &locCols);
1526:     MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
1527:     MatSetType(preallocator, MATPREALLOCATOR);
1528:     MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
1529:     MatSetUp(preallocator);
1530:     PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1531:     for (cell = cStart; cell < cEnd; ++cell) {
1532:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1533:       MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
1534:     }
1535:     PetscFree3(vals,cellCIndices,cellFIndices);
1536:     MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
1537:     MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
1538:     MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
1539:     MatDestroy(&preallocator);
1540:   }
1541:   /* Fill matrix */
1542:   MatZeroEntries(In);
1543:   for (c = cStart; c < cEnd; ++c) {
1544:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1545:   }
1546:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1547:   PetscFree2(feRef,fvRef);
1548:   PetscFree(elemMat);
1549:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1550:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1551:   if (mesh->printFEM) {
1552:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1553:     MatChop(In, 1.0e-10);
1554:     MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1555:   }
1556:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1557:   return(0);
1558: }

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

1565:   Input Parameters:
1566: + dmf  - The fine mesh
1567: . dmc  - The coarse mesh
1568: - user - The user context

1570:   Output Parameter:
1571: . In  - The interpolation matrix

1573:   Level: developer

1575: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1576: @*/
1577: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1578: {
1579:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
1580:   const char    *name = "Interpolator";
1581:   PetscDS        prob;
1582:   PetscSection   fsection, csection, globalFSection, globalCSection;
1583:   PetscHashJK    ht;
1584:   PetscLayout    rLayout;
1585:   PetscInt      *dnz, *onz;
1586:   PetscInt       locRows, rStart, rEnd;
1587:   PetscReal     *x, *v0, *J, *invJ, detJ;
1588:   PetscReal     *v0c, *Jc, *invJc, detJc;
1589:   PetscScalar   *elemMat;
1590:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

1594:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1595:   DMGetCoordinateDim(dmc, &dim);
1596:   DMGetDS(dmc, &prob);
1597:   PetscDSGetRefCoordArrays(prob, &x, NULL);
1598:   PetscDSGetNumFields(prob, &Nf);
1599:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1600:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1601:   DMGetDefaultSection(dmf, &fsection);
1602:   DMGetDefaultGlobalSection(dmf, &globalFSection);
1603:   DMGetDefaultSection(dmc, &csection);
1604:   DMGetDefaultGlobalSection(dmc, &globalCSection);
1605:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1606:   PetscDSGetTotalDimension(prob, &totDim);
1607:   PetscMalloc1(totDim*totDim, &elemMat);

1609:   MatGetLocalSize(In, &locRows, NULL);
1610:   PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1611:   PetscLayoutSetLocalSize(rLayout, locRows);
1612:   PetscLayoutSetBlockSize(rLayout, 1);
1613:   PetscLayoutSetUp(rLayout);
1614:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1615:   PetscLayoutDestroy(&rLayout);
1616:   PetscCalloc2(locRows,&dnz,locRows,&onz);
1617:   PetscHashJKCreate(&ht);
1618:   for (field = 0; field < Nf; ++field) {
1619:     PetscObject      obj;
1620:     PetscClassId     id;
1621:     PetscDualSpace   Q = NULL;
1622:     PetscQuadrature  f;
1623:     const PetscReal *qpoints;
1624:     PetscInt         Nc, Np, fpdim, i, d;

1626:     PetscDSGetDiscretization(prob, field, &obj);
1627:     PetscObjectGetClassId(obj, &id);
1628:     if (id == PETSCFE_CLASSID) {
1629:       PetscFE fe = (PetscFE) obj;

1631:       PetscFEGetDualSpace(fe, &Q);
1632:       PetscFEGetNumComponents(fe, &Nc);
1633:     } else if (id == PETSCFV_CLASSID) {
1634:       PetscFV fv = (PetscFV) obj;

1636:       PetscFVGetDualSpace(fv, &Q);
1637:       Nc   = 1;
1638:     }
1639:     PetscDualSpaceGetDimension(Q, &fpdim);
1640:     /* For each fine grid cell */
1641:     for (cell = cStart; cell < cEnd; ++cell) {
1642:       PetscInt *findices,   *cindices;
1643:       PetscInt  numFIndices, numCIndices;

1645:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1646:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1647:       if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc);
1648:       for (i = 0; i < fpdim; ++i) {
1649:         Vec             pointVec;
1650:         PetscScalar    *pV;
1651:         PetscSF         coarseCellSF = NULL;
1652:         const PetscSFNode *coarseCells;
1653:         PetscInt        numCoarseCells, q, r, c;

1655:         /* Get points from the dual basis functional quadrature */
1656:         PetscDualSpaceGetFunctional(Q, i, &f);
1657:         PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1658:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1659:         VecSetBlockSize(pointVec, dim);
1660:         VecGetArray(pointVec, &pV);
1661:         for (q = 0; q < Np; ++q) {
1662:           /* Transform point to real space */
1663:           CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1664:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1665:         }
1666:         VecRestoreArray(pointVec, &pV);
1667:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1668:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1669:         PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
1670:         /* Update preallocation info */
1671:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1672:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1673:         for (r = 0; r < Nc; ++r) {
1674:           PetscHashJKKey  key;
1675:           PetscHashJKIter missing, iter;

1677:           key.j = findices[i*Nc+r];
1678:           if (key.j < 0) continue;
1679:           /* Get indices for coarse elements */
1680:           for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1681:             DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1682:             for (c = 0; c < numCIndices; ++c) {
1683:               key.k = cindices[c];
1684:               if (key.k < 0) continue;
1685:               PetscHashJKPut(ht, key, &missing, &iter);
1686:               if (missing) {
1687:                 PetscHashJKSet(ht, iter, 1);
1688:                 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1689:                 else                                     ++onz[key.j-rStart];
1690:               }
1691:             }
1692:             DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1693:           }
1694:         }
1695:         PetscSFDestroy(&coarseCellSF);
1696:         VecDestroy(&pointVec);
1697:       }
1698:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1699:     }
1700:   }
1701:   PetscHashJKDestroy(&ht);
1702:   MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1703:   MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1704:   PetscFree2(dnz,onz);
1705:   for (field = 0; field < Nf; ++field) {
1706:     PetscObject      obj;
1707:     PetscClassId     id;
1708:     PetscDualSpace   Q = NULL;
1709:     PetscQuadrature  f;
1710:     const PetscReal *qpoints, *qweights;
1711:     PetscInt         Nc, Np, fpdim, i, d;

1713:     PetscDSGetDiscretization(prob, field, &obj);
1714:     PetscObjectGetClassId(obj, &id);
1715:     if (id == PETSCFE_CLASSID) {
1716:       PetscFE fe = (PetscFE) obj;

1718:       PetscFEGetDualSpace(fe, &Q);
1719:       PetscFEGetNumComponents(fe, &Nc);
1720:     } else if (id == PETSCFV_CLASSID) {
1721:       PetscFV fv = (PetscFV) obj;

1723:       PetscFVGetDualSpace(fv, &Q);
1724:       Nc   = 1;
1725:     }
1726:     PetscDualSpaceGetDimension(Q, &fpdim);
1727:     /* For each fine grid cell */
1728:     for (cell = cStart; cell < cEnd; ++cell) {
1729:       PetscInt *findices,   *cindices;
1730:       PetscInt  numFIndices, numCIndices;

1732:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1733:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1734:       if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc);
1735:       for (i = 0; i < fpdim; ++i) {
1736:         Vec             pointVec;
1737:         PetscScalar    *pV;
1738:         PetscSF         coarseCellSF = NULL;
1739:         const PetscSFNode *coarseCells;
1740:         PetscInt        numCoarseCells, cpdim, q, c, j;

1742:         /* Get points from the dual basis functional quadrature */
1743:         PetscDualSpaceGetFunctional(Q, i, &f);
1744:         PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);
1745:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1746:         VecSetBlockSize(pointVec, dim);
1747:         VecGetArray(pointVec, &pV);
1748:         for (q = 0; q < Np; ++q) {
1749:           /* Transform point to real space */
1750:           CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1751:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1752:         }
1753:         VecRestoreArray(pointVec, &pV);
1754:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1755:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
1756:         /* Update preallocation info */
1757:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
1758:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
1759:         VecGetArray(pointVec, &pV);
1760:         for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1761:           PetscReal pVReal[3];

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

1769:           if (id == PETSCFE_CLASSID) {
1770:             PetscFE    fe = (PetscFE) obj;
1771:             PetscReal *B;

1773:             /* Evaluate coarse basis on contained point */
1774:             PetscFEGetDimension(fe, &cpdim);
1775:             PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
1776:             PetscMemzero(elemMat, cpdim*Nc*Nc * sizeof(PetscScalar));
1777:             /* Get elemMat entries by multiplying by weight */
1778:             for (j = 0; j < cpdim; ++j) {
1779:               for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell];
1780:             }
1781:             PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
1782:           } else {
1783:             cpdim = 1;
1784:             for (j = 0; j < cpdim; ++j) {
1785:               for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell];
1786:             }
1787:           }
1788:           /* Update interpolator */
1789:           if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);}
1790:           if (numCIndices != cpdim*Nc) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D*%D", numCIndices, cpdim, Nc);
1791:           MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);
1792:           DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
1793:         }
1794:         VecRestoreArray(pointVec, &pV);
1795:         PetscSFDestroy(&coarseCellSF);
1796:         VecDestroy(&pointVec);
1797:       }
1798:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
1799:     }
1800:   }
1801:   PetscFree3(v0,J,invJ);
1802:   PetscFree3(v0c,Jc,invJc);
1803:   PetscFree(elemMat);
1804:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1805:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1806:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1807:   return(0);
1808: }

1812: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1813: {
1814:   PetscDS        prob;
1815:   PetscFE       *feRef;
1816:   PetscFV       *fvRef;
1817:   Vec            fv, cv;
1818:   IS             fis, cis;
1819:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1820:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1821:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;

1825:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1826:   DMGetDimension(dmf, &dim);
1827:   DMGetDefaultSection(dmf, &fsection);
1828:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1829:   DMGetDefaultSection(dmc, &csection);
1830:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1831:   PetscSectionGetNumFields(fsection, &Nf);
1832:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1833:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1834:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1835:   DMGetDS(dmc, &prob);
1836:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1837:   for (f = 0; f < Nf; ++f) {
1838:     PetscObject  obj;
1839:     PetscClassId id;
1840:     PetscInt     fNb = 0, Nc = 0;

1842:     PetscDSGetDiscretization(prob, f, &obj);
1843:     PetscObjectGetClassId(obj, &id);
1844:     if (id == PETSCFE_CLASSID) {
1845:       PetscFE fe = (PetscFE) obj;

1847:       PetscFERefine(fe, &feRef[f]);
1848:       PetscFEGetDimension(feRef[f], &fNb);
1849:       PetscFEGetNumComponents(fe, &Nc);
1850:     } else if (id == PETSCFV_CLASSID) {
1851:       PetscFV        fv = (PetscFV) obj;
1852:       PetscDualSpace Q;

1854:       PetscFVRefine(fv, &fvRef[f]);
1855:       PetscFVGetDualSpace(fvRef[f], &Q);
1856:       PetscDualSpaceGetDimension(Q, &fNb);
1857:       PetscFVGetNumComponents(fv, &Nc);
1858:     }
1859:     fTotDim += fNb*Nc;
1860:   }
1861:   PetscDSGetTotalDimension(prob, &cTotDim);
1862:   PetscMalloc1(cTotDim,&cmap);
1863:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1864:     PetscFE        feC;
1865:     PetscFV        fvC;
1866:     PetscDualSpace QF, QC;
1867:     PetscInt       NcF, NcC, fpdim, cpdim;

1869:     if (feRef[field]) {
1870:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1871:       PetscFEGetNumComponents(feC, &NcC);
1872:       PetscFEGetNumComponents(feRef[field], &NcF);
1873:       PetscFEGetDualSpace(feRef[field], &QF);
1874:       PetscDualSpaceGetDimension(QF, &fpdim);
1875:       PetscFEGetDualSpace(feC, &QC);
1876:       PetscDualSpaceGetDimension(QC, &cpdim);
1877:     } else {
1878:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
1879:       PetscFVGetNumComponents(fvC, &NcC);
1880:       PetscFVGetNumComponents(fvRef[field], &NcF);
1881:       PetscFVGetDualSpace(fvRef[field], &QF);
1882:       PetscDualSpaceGetDimension(QF, &fpdim);
1883:       PetscFVGetDualSpace(fvC, &QC);
1884:       PetscDualSpaceGetDimension(QC, &cpdim);
1885:     }
1886:     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);
1887:     for (c = 0; c < cpdim; ++c) {
1888:       PetscQuadrature  cfunc;
1889:       const PetscReal *cqpoints;
1890:       PetscInt         NpC;
1891:       PetscBool        found = PETSC_FALSE;

1893:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
1894:       PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);
1895:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1896:       for (f = 0; f < fpdim; ++f) {
1897:         PetscQuadrature  ffunc;
1898:         const PetscReal *fqpoints;
1899:         PetscReal        sum = 0.0;
1900:         PetscInt         NpF, comp;

1902:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
1903:         PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);
1904:         if (NpC != NpF) continue;
1905:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1906:         if (sum > 1.0e-9) continue;
1907:         for (comp = 0; comp < NcC; ++comp) {
1908:           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1909:         }
1910:         found = PETSC_TRUE;
1911:         break;
1912:       }
1913:       if (!found) {
1914:         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1915:         if (fvRef[field]) {
1916:           PetscInt comp;
1917:           for (comp = 0; comp < NcC; ++comp) {
1918:             cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1919:           }
1920:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1921:       }
1922:     }
1923:     offsetC += cpdim*NcC;
1924:     offsetF += fpdim*NcF;
1925:   }
1926:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
1927:   PetscFree2(feRef,fvRef);

1929:   DMGetGlobalVector(dmf, &fv);
1930:   DMGetGlobalVector(dmc, &cv);
1931:   VecGetOwnershipRange(cv, &startC, &endC);
1932:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1933:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1934:   PetscMalloc1(m,&cindices);
1935:   PetscMalloc1(m,&findices);
1936:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1937:   for (c = cStart; c < cEnd; ++c) {
1938:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1939:     for (d = 0; d < cTotDim; ++d) {
1940:       if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
1941:       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]]);
1942:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1943:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1944:     }
1945:   }
1946:   PetscFree(cmap);
1947:   PetscFree2(cellCIndices,cellFIndices);

1949:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1950:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1951:   VecScatterCreate(cv, cis, fv, fis, sc);
1952:   ISDestroy(&cis);
1953:   ISDestroy(&fis);
1954:   DMRestoreGlobalVector(dmf, &fv);
1955:   DMRestoreGlobalVector(dmc, &cv);
1956:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1957:   return(0);
1958: }