Actual source code: plexfem.c

petsc-master 2015-05-04
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 void DMPlexProjectRigidBody(const PetscReal X[], PetscScalar *mode, void *ctx)
 90: {
 91:   PetscInt       *ctxInt = (PetscInt *)ctx;
 92:   PetscInt       dim     = ctxInt[0];
 93:   PetscInt       d       = ctxInt[1];
 94:   PetscInt       i, j, k = dim > 2 ? d - dim : d;

 96:   for (i = 0; i < dim; i++) mode[i] = 0.;
 97:   if (d < dim) {
 98:     mode[d] = 1.;
 99:   }
100:   else {
101:     for (i = 0; i < dim; i++) {
102:       for (j = 0; j < dim; j++) {
103:         mode[j] += epsilon(i, j, k)*X[i];
104:       }
105:     }
106:   }
107: }

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

114:   Collective on DM

116:   Input Arguments:
117: . dm - the DM

119:   Output Argument:
120: . sp - the null space

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

124:   Level: advanced

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

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

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

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

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

186:   Input Parameters:
187: + dm - the DMPlex object
188: - height - the maximum projection height >= 0

190:   Level: advanced

192: .seealso: DMPlexGetMaxProjectionHeight(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
193: @*/
194: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
195: {
196:   DM_Plex *plex = (DM_Plex *) dm->data;

200:   plex->maxProjectionHeight = height;
201:   return(0);
202: }

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

210:   Input Parameters:
211: . dm - the DMPlex object

213:   Output Parameters:
214: . height - the maximum projection height

216:   Level: intermediate

218: .seealso: DMPlexSetMaxProjectionHeight(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
219: @*/
220: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
221: {
222:   DM_Plex *plex = (DM_Plex *) dm->data;

226:   *height = plex->maxProjectionHeight;
227:   return(0);
228: }

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

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

265:       DMGetField(dm, f, &obj);
266:       PetscObjectGetClassId(obj, &id);
267:       if (id == PETSCFE_CLASSID) {
268:         PetscFE fe = (PetscFE) obj;

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

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

300:       DMLabelGetStratumIS(label, ids[i], &pointIS);
301:       ISGetLocalSize(pointIS, &n);
302:       ISGetIndices(pointIS, &points);
303:       for (p = 0; p < n; ++p) {
304:         const PetscInt    point = points[p];
305:         PetscFECellGeom   geom;

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

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

346: PetscErrorCode DMPlexProjectFunctionLocal(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
347: {
348:   PetscDualSpace *sp, *cellsp;
349:   PetscInt       *numComp;
350:   PetscSection    section;
351:   PetscScalar    *values;
352:   PetscInt        Nf, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight;
353:   PetscErrorCode  ierr;

356:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
357:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
358:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
359:   if (cEnd <= cStart) return(0);
360:   DMGetDefaultSection(dm, &section);
361:   PetscSectionGetNumFields(section, &Nf);
362:   PetscMalloc2(Nf, &sp, Nf, &numComp);
363:   DMGetDimension(dm, &dim);
364:   DMGetCoordinateDim(dm, &dimEmbed);
365:   DMPlexGetMaxProjectionHeight(dm,&maxHeight);
366:   if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
367:   if (maxHeight > 0) {
368:     PetscMalloc1(Nf,&cellsp);
369:   }
370:   else {
371:     cellsp = sp;
372:   }
373:   for (h = 0; h <= maxHeight; h++) {
374:     DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
375:     if (!h) {pStart = cStart; pEnd = cEnd;}
376:     if (pEnd <= pStart) continue;
377:     totDim = 0;
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:         PetscFEGetNumComponents(fe, &numComp[f]);
388:         if (!h) {
389:           PetscFEGetDualSpace(fe, &cellsp[f]);
390:           sp[f] = cellsp[f];
391:           PetscObjectReference((PetscObject) sp[f]);
392:         }
393:         else {
394:           PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
395:           if (!sp[f]) {
396:             continue;
397:           }
398:         }
399:       } else if (id == PETSCFV_CLASSID) {
400:         PetscFV fv = (PetscFV) obj;

402:         if (h) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, "Projection height > 0 not supported for finite volume");
403:         PetscFVGetNumComponents(fv, &numComp[f]);
404:         PetscFVGetDualSpace(fv, &sp[f]);
405:         PetscObjectReference((PetscObject) sp[f]);
406:       } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
407:       PetscDualSpaceGetDimension(sp[f], &spDim);
408:       totDim += spDim*numComp[f];
409:     }
410:     DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);
411:     if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
412:     if (!totDim) continue;
413:     DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
414:     for (p = pStart; p < pEnd; ++p) {
415:       PetscFECellGeom geom;

417:       DMPlexComputeCellGeometryFEM(dm, p, NULL, geom.v0, geom.J, NULL, &geom.detJ);
418:       geom.dim      = dim - h;
419:       geom.dimEmbed = dimEmbed;
420:       for (f = 0, v = 0; f < Nf; ++f) {
421:         void * const ctx = ctxs ? ctxs[f] : NULL;

423:         if (!sp[f]) continue;
424:         PetscDualSpaceGetDimension(sp[f], &spDim);
425:         for (d = 0; d < spDim; ++d) {
426:           if (funcs[f]) {
427:             PetscDualSpaceApply(sp[f], d, &geom, numComp[f], funcs[f], ctx, &values[v]);
428:           } else {
429:             for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
430:           }
431:           v += numComp[f];
432:         }
433:       }
434:       DMPlexVecSetClosure(dm, section, localX, p, values, mode);
435:     }
436:     DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
437:     if (h || !maxHeight) {
438:       for (f = 0; f < Nf; f++) {PetscDualSpaceDestroy(&sp[f]);}
439:     }
440:   }
441:   PetscFree2(sp, numComp);
442:   if (maxHeight > 0) {
443:     for (f = 0; f < Nf; f++) {PetscDualSpaceDestroy(&cellsp[f]);}
444:     PetscFree(cellsp);
445:   }
446:   return(0);
447: }

451: /*@C
452:   DMPlexProjectFunction - This projects the given function into the function space provided.

454:   Input Parameters:
455: + dm      - The DM
456: . funcs   - The coordinate functions to evaluate, one per field
457: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
458: - mode    - The insertion mode for values

460:   Output Parameter:
461: . X - vector

463:   Level: developer

465: .seealso: DMPlexComputeL2Diff()
466: @*/
467: PetscErrorCode DMPlexProjectFunction(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
468: {
469:   Vec            localX;

474:   DMGetLocalVector(dm, &localX);
475:   DMPlexProjectFunctionLocal(dm, funcs, ctxs, mode, localX);
476:   DMLocalToGlobalBegin(dm, localX, mode, X);
477:   DMLocalToGlobalEnd(dm, localX, mode, X);
478:   DMRestoreLocalVector(dm, &localX);
479:   return(0);
480: }

484: PetscErrorCode DMPlexProjectFieldLocal(DM dm, Vec localU, void (**funcs)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal [], PetscScalar []), InsertMode mode, Vec localX)
485: {
486:   DM                dmAux;
487:   PetscDS           prob, probAux = NULL;
488:   Vec               A;
489:   PetscSection      section, sectionAux = NULL;
490:   PetscDualSpace   *sp;
491:   PetscInt         *Ncf;
492:   PetscScalar      *values, *u, *u_x, *a, *a_x;
493:   PetscReal        *x, *v0, *J, *invJ, detJ;
494:   PetscInt          Nf, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight;
495:   PetscErrorCode    ierr;

498:   DMPlexGetMaxProjectionHeight(dm,&maxHeight);
499:   if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");}
500:   DMGetDS(dm, &prob);
501:   DMGetDimension(dm, &dim);
502:   DMGetDefaultSection(dm, &section);
503:   PetscSectionGetNumFields(section, &Nf);
504:   PetscMalloc2(Nf, &sp, Nf, &Ncf);
505:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
506:   PetscDSGetTotalDimension(prob, &totDim);
507:   PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
508:   PetscDSGetRefCoordArrays(prob, &x, NULL);
509:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
510:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
511:   if (dmAux) {
512:     DMGetDS(dmAux, &probAux);
513:     DMGetDefaultSection(dmAux, &sectionAux);
514:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
515:   }
516:   DMPlexInsertBoundaryValues(dm, localU, 0.0, NULL, NULL, NULL);
517:   DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
518:   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
519:   DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
520:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
521:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
522:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
523:   for (c = cStart; c < cEnd; ++c) {
524:     PetscScalar *coefficients = NULL, *coefficientsAux = NULL;

526:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
527:     DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);
528:     if (dmAux) {DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
529:     for (f = 0, v = 0; f < Nf; ++f) {
530:       PetscObject  obj;
531:       PetscClassId id;

533:       PetscDSGetDiscretization(prob, f, &obj);
534:       PetscObjectGetClassId(obj, &id);
535:       if (id == PETSCFE_CLASSID) {
536:         PetscFE fe = (PetscFE) obj;

538:         PetscFEGetDualSpace(fe, &sp[f]);
539:         PetscObjectReference((PetscObject) sp[f]);
540:         PetscFEGetNumComponents(fe, &Ncf[f]);
541:       } else if (id == PETSCFV_CLASSID) {
542:         PetscFV fv = (PetscFV) obj;

544:         PetscFVGetNumComponents(fv, &Ncf[f]);
545:         PetscFVGetDualSpace(fv, &sp[f]);
546:         PetscObjectReference((PetscObject) sp[f]);
547:       }
548:       PetscDualSpaceGetDimension(sp[f], &spDim);
549:       for (d = 0; d < spDim; ++d) {
550:         PetscQuadrature  quad;
551:         const PetscReal *points, *weights;
552:         PetscInt         numPoints, q;

554:         if (funcs[f]) {
555:           PetscDualSpaceGetFunctional(sp[f], d, &quad);
556:           PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);
557:           for (q = 0; q < numPoints; ++q) {
558:             CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x);
559:             EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, coefficients,    NULL, u, u_x, NULL);
560:             EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);
561:             (*funcs[f])(u, NULL, u_x, a, NULL, a_x, x, &values[v]);
562:           }
563:         } else {
564:           for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0;
565:         }
566:         v += Ncf[f];
567:       }
568:     }
569:     DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);
570:     if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
571:     DMPlexVecSetClosure(dm, section, localX, c, values, mode);
572:   }
573:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
574:   PetscFree3(v0,J,invJ);
575:   for (f = 0; f < Nf; f++) {PetscDualSpaceDestroy(&sp[f]);}
576:   PetscFree2(sp, Ncf);
577:   return(0);
578: }

582: static PetscErrorCode DMPlexInsertBoundaryValues_FEM_Internal(DM dm, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], void (*func)(const PetscReal[], PetscScalar *, void *), void *ctx, Vec locX)
583: {
584:   void        (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx);
585:   void         **ctxs;
586:   PetscInt       numFields;

590:   DMGetNumFields(dm, &numFields);
591:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
592:   funcs[field] = func;
593:   ctxs[field]  = ctx;
594:   DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);
595:   PetscFree2(funcs,ctxs);
596:   return(0);
597: }

601: static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad,
602:                                                               PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
603: {
604:   PetscDS            prob;
605:   PetscSF            sf;
606:   DM                 dmFace, dmCell, dmGrad;
607:   const PetscScalar *facegeom, *cellgeom, *grad;
608:   const PetscInt    *leaves;
609:   PetscScalar       *x, *fx;
610:   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
611:   PetscErrorCode     ierr;

614:   DMGetPointSF(dm, &sf);
615:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
616:   nleaves = PetscMax(0, nleaves);
617:   DMGetDimension(dm, &dim);
618:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
619:   DMGetDS(dm, &prob);
620:   VecGetDM(faceGeometry, &dmFace);
621:   VecGetArrayRead(faceGeometry, &facegeom);
622:   VecGetDM(cellGeometry, &dmCell);
623:   VecGetArrayRead(cellGeometry, &cellgeom);
624:   if (Grad) {
625:     PetscFV fv;

627:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
628:     VecGetDM(Grad, &dmGrad);
629:     VecGetArrayRead(Grad, &grad);
630:     PetscFVGetNumComponents(fv, &pdim);
631:     DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
632:   }
633:   VecGetArray(locX, &x);
634:   for (i = 0; i < numids; ++i) {
635:     IS              faceIS;
636:     const PetscInt *faces;
637:     PetscInt        numFaces, f;

639:     DMLabelGetStratumIS(label, ids[i], &faceIS);
640:     if (!faceIS) continue; /* No points with that id on this process */
641:     ISGetLocalSize(faceIS, &numFaces);
642:     ISGetIndices(faceIS, &faces);
643:     for (f = 0; f < numFaces; ++f) {
644:       const PetscInt         face = faces[f], *cells;
645:       const PetscFVFaceGeom *fg;

647:       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
648:       PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
649:       if (loc >= 0) continue;
650:       DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
651:       DMPlexGetSupport(dm, face, &cells);
652:       if (Grad) {
653:         const PetscFVCellGeom *cg;
654:         const PetscScalar     *cx, *cgrad;
655:         PetscScalar           *xG;
656:         PetscReal              dx[3];
657:         PetscInt               d;

659:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
660:         DMPlexPointLocalRead(dm, cells[0], x, &cx);
661:         DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
662:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
663:         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
664:         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
665:         (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
666:       } else {
667:         const PetscScalar *xI;
668:         PetscScalar       *xG;

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

690: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
691: {
692:   PetscInt       numBd, b;

701:   DMPlexGetNumBoundary(dm, &numBd);
702:   for (b = 0; b < numBd; ++b) {
703:     PetscBool       isEssential;
704:     const char     *labelname;
705:     DMLabel         label;
706:     PetscInt        field;
707:     PetscObject     obj;
708:     PetscClassId    id;
709:     void          (*func)();
710:     PetscInt        numids;
711:     const PetscInt *ids;
712:     void           *ctx;

714:     DMPlexGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, &func, &numids, &ids, &ctx);
715:     if (!isEssential) continue;
716:     DMPlexGetLabel(dm, labelname, &label);
717:     DMGetField(dm, field, &obj);
718:     PetscObjectGetClassId(obj, &id);
719:     if (id == PETSCFE_CLASSID) {
720:       DMPlexInsertBoundaryValues_FEM_Internal(dm, field, label, numids, ids, (void (*)(const PetscReal[], PetscScalar *, void *)) func, ctx, locX);
721:     } else if (id == PETSCFV_CLASSID) {
722:       if (!faceGeomFVM) continue;
723:       DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM,
724:                                                      field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
725:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
726:   }
727:   return(0);
728: }

732: /*@C
733:   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

735:   Input Parameters:
736: + dm    - The DM
737: . funcs - The functions to evaluate for each field component
738: . ctxs  - Optional array of contexts to pass to each function, or NULL.
739: - X     - The coefficient vector u_h

741:   Output Parameter:
742: . diff - The diff ||u - u_h||_2

744:   Level: developer

746: .seealso: DMPlexProjectFunction(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff()
747: @*/
748: PetscErrorCode DMPlexComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
749: {
750:   const PetscInt   debug = 0;
751:   PetscSection     section;
752:   PetscQuadrature  quad;
753:   Vec              localX;
754:   PetscScalar     *funcVal, *interpolant;
755:   PetscReal       *coords, *v0, *J, *invJ, detJ;
756:   PetscReal        localDiff = 0.0;
757:   const PetscReal *quadPoints, *quadWeights;
758:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
759:   PetscErrorCode   ierr;

762:   DMGetDimension(dm, &dim);
763:   DMGetDefaultSection(dm, &section);
764:   PetscSectionGetNumFields(section, &numFields);
765:   DMGetLocalVector(dm, &localX);
766:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
767:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
768:   for (field = 0; field < numFields; ++field) {
769:     PetscObject  obj;
770:     PetscClassId id;
771:     PetscInt     Nc;

773:     DMGetField(dm, field, &obj);
774:     PetscObjectGetClassId(obj, &id);
775:     if (id == PETSCFE_CLASSID) {
776:       PetscFE fe = (PetscFE) obj;

778:       PetscFEGetQuadrature(fe, &quad);
779:       PetscFEGetNumComponents(fe, &Nc);
780:     } else if (id == PETSCFV_CLASSID) {
781:       PetscFV fv = (PetscFV) obj;

783:       PetscFVGetQuadrature(fv, &quad);
784:       PetscFVGetNumComponents(fv, &Nc);
785:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
786:     numComponents += Nc;
787:   }
788:   PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
789:   DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);
790:   PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
791:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
792:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
793:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
794:   for (c = cStart; c < cEnd; ++c) {
795:     PetscScalar *x = NULL;
796:     PetscReal    elemDiff = 0.0;

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

802:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
803:       PetscObject  obj;
804:       PetscClassId id;
805:       void * const ctx = ctxs ? ctxs[field] : NULL;
806:       PetscInt     Nb, Nc, q, fc;

808:       DMGetField(dm, field, &obj);
809:       PetscObjectGetClassId(obj, &id);
810:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
811:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
812:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
813:       if (debug) {
814:         char title[1024];
815:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
816:         DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);
817:       }
818:       for (q = 0; q < numQuadPoints; ++q) {
819:         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
820:         (*funcs[field])(coords, funcVal, ctx);
821:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
822:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
823:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
824:         for (fc = 0; fc < Nc; ++fc) {
825:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ);}
826:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
827:         }
828:       }
829:       fieldOffset += Nb*Nc;
830:     }
831:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
832:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
833:     localDiff += elemDiff;
834:   }
835:   PetscFree6(funcVal,interpolant,coords,v0,J,invJ);
836:   DMRestoreLocalVector(dm, &localX);
837:   MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
838:   *diff = PetscSqrtReal(*diff);
839:   return(0);
840: }

844: /*@C
845:   DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

847:   Input Parameters:
848: + dm    - The DM
849: . funcs - The gradient functions to evaluate for each field component
850: . ctxs  - Optional array of contexts to pass to each function, or NULL.
851: . X     - The coefficient vector u_h
852: - n     - The vector to project along

854:   Output Parameter:
855: . diff - The diff ||(grad u - grad u_h) . n||_2

857:   Level: developer

859: .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
860: @*/
861: PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
862: {
863:   const PetscInt  debug = 0;
864:   PetscSection    section;
865:   PetscQuadrature quad;
866:   Vec             localX;
867:   PetscScalar    *funcVal, *interpolantVec;
868:   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
869:   PetscReal       localDiff = 0.0;
870:   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp;
871:   PetscErrorCode  ierr;

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

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

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

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

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

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

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

963: /*@C
964:   DMPlexComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

966:   Input Parameters:
967: + dm    - The DM
968: . funcs - The functions to evaluate for each field component
969: . ctxs  - Optional array of contexts to pass to each function, or NULL.
970: - X     - The coefficient vector u_h

972:   Output Parameter:
973: . diff - The array of differences, ||u^f - u^f_h||_2

975:   Level: developer

977: .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2GradientDiff()
978: @*/
979: PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
980: {
981:   const PetscInt   debug = 0;
982:   PetscSection     section;
983:   PetscQuadrature  quad;
984:   Vec              localX;
985:   PetscScalar     *funcVal, *interpolant;
986:   PetscReal       *coords, *v0, *J, *invJ, detJ;
987:   PetscReal       *localDiff;
988:   const PetscReal *quadPoints, *quadWeights;
989:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
990:   PetscErrorCode   ierr;

993:   DMGetDimension(dm, &dim);
994:   DMGetDefaultSection(dm, &section);
995:   PetscSectionGetNumFields(section, &numFields);
996:   DMGetLocalVector(dm, &localX);
997:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
998:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
999:   for (field = 0; field < numFields; ++field) {
1000:     PetscObject  obj;
1001:     PetscClassId id;
1002:     PetscInt     Nc;

1004:     DMGetField(dm, field, &obj);
1005:     PetscObjectGetClassId(obj, &id);
1006:     if (id == PETSCFE_CLASSID) {
1007:       PetscFE fe = (PetscFE) obj;

1009:       PetscFEGetQuadrature(fe, &quad);
1010:       PetscFEGetNumComponents(fe, &Nc);
1011:     } else if (id == PETSCFV_CLASSID) {
1012:       PetscFV fv = (PetscFV) obj;

1014:       PetscFVGetQuadrature(fv, &quad);
1015:       PetscFVGetNumComponents(fv, &Nc);
1016:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1017:     numComponents += Nc;
1018:   }
1019:   PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1020:   DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);
1021:   PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1022:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1023:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1024:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1025:   for (c = cStart; c < cEnd; ++c) {
1026:     PetscScalar *x = NULL;

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

1032:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1033:       PetscObject  obj;
1034:       PetscClassId id;
1035:       void * const ctx = ctxs ? ctxs[field] : NULL;
1036:       PetscInt     Nb, Nc, q, fc;

1038:       PetscReal       elemDiff = 0.0;

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

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

1078:   Input Parameters:
1079: + dm - The mesh
1080: . X  - Local input vector
1081: - user - The user context

1083:   Output Parameter:
1084: . integral - Local integral for each field

1086:   Level: developer

1088: .seealso: DMPlexComputeResidualFEM()
1089: @*/
1090: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
1091: {
1092:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1093:   DM                dmAux;
1094:   Vec               localX, A;
1095:   PetscDS           prob, probAux = NULL;
1096:   PetscSection      section, sectionAux;
1097:   PetscFECellGeom  *cgeom;
1098:   PetscScalar      *u, *a = NULL;
1099:   PetscReal        *lintegral, *vol;
1100:   PetscInt          dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c;
1101:   PetscInt          totDim, totDimAux;
1102:   PetscErrorCode    ierr;

1105:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1106:   DMGetLocalVector(dm, &localX);
1107:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1108:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1109:   DMGetDimension(dm, &dim);
1110:   DMGetDefaultSection(dm, &section);
1111:   DMGetDS(dm, &prob);
1112:   PetscDSGetTotalDimension(prob, &totDim);
1113:   PetscSectionGetNumFields(section, &Nf);
1114:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1115:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1116:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1117:   numCells = cEnd - cStart;
1118:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1119:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1120:   if (dmAux) {
1121:     DMGetDefaultSection(dmAux, &sectionAux);
1122:     DMGetDS(dmAux, &probAux);
1123:     PetscDSGetTotalDimension(probAux, &totDimAux);
1124:   }
1125:   DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);
1126:   PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);
1127:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1128:   for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
1129:   for (c = cStart; c < cEnd; ++c) {
1130:     PetscScalar *x = NULL;
1131:     PetscInt     i;

1133:     DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);
1134:     DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);
1135:     if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c);
1136:     DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
1137:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1138:     DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
1139:     if (dmAux) {
1140:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1141:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1142:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1143:     }
1144:   }
1145:   for (f = 0; f < Nf; ++f) {
1146:     PetscObject  obj;
1147:     PetscClassId id;
1148:     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1150:     PetscDSGetDiscretization(prob, f, &obj);
1151:     PetscObjectGetClassId(obj, &id);
1152:     if (id == PETSCFE_CLASSID) {
1153:       PetscFE         fe = (PetscFE) obj;
1154:       PetscQuadrature q;
1155:       PetscInt        Nq, Nb;

1157:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1158:       PetscFEGetQuadrature(fe, &q);
1159:       PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1160:       PetscFEGetDimension(fe, &Nb);
1161:       blockSize = Nb*Nq;
1162:       batchSize = numBlocks * blockSize;
1163:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1164:       numChunks = numCells / (numBatches*batchSize);
1165:       Ne        = numChunks*numBatches*batchSize;
1166:       Nr        = numCells % (numBatches*batchSize);
1167:       offset    = numCells - Nr;
1168:       PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);
1169:       PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);
1170:     } else if (id == PETSCFV_CLASSID) {
1171:       /* PetscFV  fv = (PetscFV) obj; */
1172:       PetscInt    foff;
1173:       void      (*obj_func)(const PetscScalar u[], const PetscScalar u_x[], const PetscScalar u_t[], const PetscScalar a[], const PetscScalar a_x[], const PetscScalar a_t[], const PetscReal x[], PetscScalar obj[]);
1174:       PetscScalar lint;

1176:       PetscDSGetObjective(prob, f, &obj_func);
1177:       PetscDSGetFieldOffset(prob, f, &foff);
1178:       if (obj_func) {
1179:         for (c = 0; c < numCells; ++c) {
1180:           /* TODO: Need full pointwise interpolation and get centroid for x argument */
1181:           obj_func(&u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, &lint);
1182:           lintegral[f] = PetscRealPart(lint)*vol[c];
1183:         }
1184:       }
1185:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1186:   }
1187:   if (dmAux) {PetscFree(a);}
1188:   if (mesh->printFEM) {
1189:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1190:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);}
1191:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1192:   }
1193:   DMRestoreLocalVector(dm, &localX);
1194:   MPI_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1195:   PetscFree4(lintegral,u,cgeom,vol);
1196:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1197:   return(0);
1198: }

1202: /*@
1203:   DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.

1205:   Input Parameters:
1206: + dmf  - The fine mesh
1207: . dmc  - The coarse mesh
1208: - user - The user context

1210:   Output Parameter:
1211: . In  - The interpolation matrix

1213:   Note:
1214:   The first member of the user context must be an FEMContext.

1216:   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1217:   like a GPU, or vectorize on a multicore machine.

1219:   Level: developer

1221: .seealso: DMPlexComputeJacobianFEM()
1222: @*/
1223: PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1224: {
1225:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1226:   const char       *name  = "Interpolator";
1227:   PetscDS           prob;
1228:   PetscFE          *feRef;
1229:   PetscFV          *fvRef;
1230:   PetscSection      fsection, fglobalSection;
1231:   PetscSection      csection, cglobalSection;
1232:   PetscScalar      *elemMat;
1233:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1234:   PetscInt          cTotDim, rTotDim = 0;
1235:   PetscErrorCode    ierr;

1238:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1239:   DMGetDimension(dmf, &dim);
1240:   DMGetDefaultSection(dmf, &fsection);
1241:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1242:   DMGetDefaultSection(dmc, &csection);
1243:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1244:   PetscSectionGetNumFields(fsection, &Nf);
1245:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1246:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1247:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1248:   DMGetDS(dmf, &prob);
1249:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1250:   for (f = 0; f < Nf; ++f) {
1251:     PetscObject  obj;
1252:     PetscClassId id;
1253:     PetscInt     rNb = 0, Nc = 0;

1255:     PetscDSGetDiscretization(prob, f, &obj);
1256:     PetscObjectGetClassId(obj, &id);
1257:     if (id == PETSCFE_CLASSID) {
1258:       PetscFE fe = (PetscFE) obj;

1260:       PetscFERefine(fe, &feRef[f]);
1261:       PetscFEGetDimension(feRef[f], &rNb);
1262:       PetscFEGetNumComponents(fe, &Nc);
1263:     } else if (id == PETSCFV_CLASSID) {
1264:       PetscFV        fv = (PetscFV) obj;
1265:       PetscDualSpace Q;

1267:       PetscFVRefine(fv, &fvRef[f]);
1268:       PetscFVGetDualSpace(fvRef[f], &Q);
1269:       PetscDualSpaceGetDimension(Q, &rNb);
1270:       PetscFVGetNumComponents(fv, &Nc);
1271:     }
1272:     rTotDim += rNb*Nc;
1273:   }
1274:   PetscDSGetTotalDimension(prob, &cTotDim);
1275:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1276:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1277:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1278:     PetscDualSpace   Qref;
1279:     PetscQuadrature  f;
1280:     const PetscReal *qpoints, *qweights;
1281:     PetscReal       *points;
1282:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1284:     /* Compose points from all dual basis functionals */
1285:     if (feRef[fieldI]) {
1286:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
1287:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
1288:     } else {
1289:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1290:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1291:     }
1292:     PetscDualSpaceGetDimension(Qref, &fpdim);
1293:     for (i = 0; i < fpdim; ++i) {
1294:       PetscDualSpaceGetFunctional(Qref, i, &f);
1295:       PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1296:       npoints += Np;
1297:     }
1298:     PetscMalloc1(npoints*dim,&points);
1299:     for (i = 0, k = 0; i < fpdim; ++i) {
1300:       PetscDualSpaceGetFunctional(Qref, i, &f);
1301:       PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1302:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1303:     }

1305:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1306:       PetscObject  obj;
1307:       PetscClassId id;
1308:       PetscReal   *B;
1309:       PetscInt     NcJ = 0, cpdim = 0, j;

1311:       PetscDSGetDiscretization(prob, fieldJ, &obj);
1312:       PetscObjectGetClassId(obj, &id);
1313:       if (id == PETSCFE_CLASSID) {
1314:         PetscFE fe = (PetscFE) obj;

1316:         /* Evaluate basis at points */
1317:         PetscFEGetNumComponents(fe, &NcJ);
1318:         PetscFEGetDimension(fe, &cpdim);
1319:         /* For now, fields only interpolate themselves */
1320:         if (fieldI == fieldJ) {
1321:           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);
1322:           PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1323:           for (i = 0, k = 0; i < fpdim; ++i) {
1324:             PetscDualSpaceGetFunctional(Qref, i, &f);
1325:             PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1326:             for (p = 0; p < Np; ++p, ++k) {
1327:               for (j = 0; j < cpdim; ++j) {
1328:                 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];
1329:               }
1330:             }
1331:           }
1332:           PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1333:         }
1334:       } else if (id == PETSCFV_CLASSID) {
1335:         PetscFV        fv = (PetscFV) obj;

1337:         /* Evaluate constant function at points */
1338:         PetscFVGetNumComponents(fv, &NcJ);
1339:         cpdim = 1;
1340:         /* For now, fields only interpolate themselves */
1341:         if (fieldI == fieldJ) {
1342:           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);
1343:           for (i = 0, k = 0; i < fpdim; ++i) {
1344:             PetscDualSpaceGetFunctional(Qref, i, &f);
1345:             PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1346:             for (p = 0; p < Np; ++p, ++k) {
1347:               for (j = 0; j < cpdim; ++j) {
1348:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p];
1349:               }
1350:             }
1351:           }
1352:         }
1353:       }
1354:       offsetJ += cpdim*NcJ;
1355:     }
1356:     offsetI += fpdim*Nc;
1357:     PetscFree(points);
1358:   }
1359:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1360:   /* Preallocate matrix */
1361:   {
1362:     PetscHashJK ht;
1363:     PetscLayout rLayout;
1364:     PetscInt   *dnz, *onz, *cellCIndices, *cellFIndices;
1365:     PetscInt    locRows, rStart, rEnd, cell, r;

1367:     MatGetLocalSize(In, &locRows, NULL);
1368:     PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1369:     PetscLayoutSetLocalSize(rLayout, locRows);
1370:     PetscLayoutSetBlockSize(rLayout, 1);
1371:     PetscLayoutSetUp(rLayout);
1372:     PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1373:     PetscLayoutDestroy(&rLayout);
1374:     PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1375:     PetscHashJKCreate(&ht);
1376:     for (cell = cStart; cell < cEnd; ++cell) {
1377:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1378:       for (r = 0; r < rTotDim; ++r) {
1379:         PetscHashJKKey  key;
1380:         PetscHashJKIter missing, iter;

1382:         key.j = cellFIndices[r];
1383:         if (key.j < 0) continue;
1384:         for (c = 0; c < cTotDim; ++c) {
1385:           key.k = cellCIndices[c];
1386:           if (key.k < 0) continue;
1387:           PetscHashJKPut(ht, key, &missing, &iter);
1388:           if (missing) {
1389:             PetscHashJKSet(ht, iter, 1);
1390:             if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1391:             else                                     ++onz[key.j-rStart];
1392:           }
1393:         }
1394:       }
1395:     }
1396:     PetscHashJKDestroy(&ht);
1397:     MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1398:     MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1399:     PetscFree4(dnz,onz,cellCIndices,cellFIndices);
1400:   }
1401:   /* Fill matrix */
1402:   MatZeroEntries(In);
1403:   for (c = cStart; c < cEnd; ++c) {
1404:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1405:   }
1406:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1407:   PetscFree2(feRef,fvRef);
1408:   PetscFree(elemMat);
1409:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1410:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1411:   if (mesh->printFEM) {
1412:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1413:     MatChop(In, 1.0e-10);
1414:     MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1415:   }
1416:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1417:   return(0);
1418: }

1422: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1423: {
1424:   PetscDS        prob;
1425:   PetscFE       *feRef;
1426:   PetscFV       *fvRef;
1427:   Vec            fv, cv;
1428:   IS             fis, cis;
1429:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1430:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1431:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m;

1435:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1436:   DMGetDimension(dmf, &dim);
1437:   DMGetDefaultSection(dmf, &fsection);
1438:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1439:   DMGetDefaultSection(dmc, &csection);
1440:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1441:   PetscSectionGetNumFields(fsection, &Nf);
1442:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1443:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1444:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1445:   DMGetDS(dmc, &prob);
1446:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1447:   for (f = 0; f < Nf; ++f) {
1448:     PetscObject  obj;
1449:     PetscClassId id;
1450:     PetscInt     fNb = 0, Nc = 0;

1452:     PetscDSGetDiscretization(prob, f, &obj);
1453:     PetscObjectGetClassId(obj, &id);
1454:     if (id == PETSCFE_CLASSID) {
1455:       PetscFE fe = (PetscFE) obj;

1457:       PetscFERefine(fe, &feRef[f]);
1458:       PetscFEGetDimension(feRef[f], &fNb);
1459:       PetscFEGetNumComponents(fe, &Nc);
1460:     } else if (id == PETSCFV_CLASSID) {
1461:       PetscFV        fv = (PetscFV) obj;
1462:       PetscDualSpace Q;

1464:       PetscFVRefine(fv, &fvRef[f]);
1465:       PetscFVGetDualSpace(fvRef[f], &Q);
1466:       PetscDualSpaceGetDimension(Q, &fNb);
1467:       PetscFVGetNumComponents(fv, &Nc);
1468:     }
1469:     fTotDim += fNb*Nc;
1470:   }
1471:   PetscDSGetTotalDimension(prob, &cTotDim);
1472:   PetscMalloc1(cTotDim,&cmap);
1473:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1474:     PetscFE        feC;
1475:     PetscFV        fvC;
1476:     PetscDualSpace QF, QC;
1477:     PetscInt       NcF, NcC, fpdim, cpdim;

1479:     if (feRef[field]) {
1480:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1481:       PetscFEGetNumComponents(feC, &NcC);
1482:       PetscFEGetNumComponents(feRef[field], &NcF);
1483:       PetscFEGetDualSpace(feRef[field], &QF);
1484:       PetscDualSpaceGetDimension(QF, &fpdim);
1485:       PetscFEGetDualSpace(feC, &QC);
1486:       PetscDualSpaceGetDimension(QC, &cpdim);
1487:     } else {
1488:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
1489:       PetscFVGetNumComponents(fvC, &NcC);
1490:       PetscFVGetNumComponents(fvRef[field], &NcF);
1491:       PetscFVGetDualSpace(fvRef[field], &QF);
1492:       PetscDualSpaceGetDimension(QF, &fpdim);
1493:       PetscFVGetDualSpace(fvC, &QC);
1494:       PetscDualSpaceGetDimension(QC, &cpdim);
1495:     }
1496:     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);
1497:     for (c = 0; c < cpdim; ++c) {
1498:       PetscQuadrature  cfunc;
1499:       const PetscReal *cqpoints;
1500:       PetscInt         NpC;
1501:       PetscBool        found = PETSC_FALSE;

1503:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
1504:       PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);
1505:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1506:       for (f = 0; f < fpdim; ++f) {
1507:         PetscQuadrature  ffunc;
1508:         const PetscReal *fqpoints;
1509:         PetscReal        sum = 0.0;
1510:         PetscInt         NpF, comp;

1512:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
1513:         PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);
1514:         if (NpC != NpF) continue;
1515:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1516:         if (sum > 1.0e-9) continue;
1517:         for (comp = 0; comp < NcC; ++comp) {
1518:           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1519:         }
1520:         found = PETSC_TRUE;
1521:         break;
1522:       }
1523:       if (!found) {
1524:         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1525:         if (fvRef[field]) {
1526:           PetscInt comp;
1527:           for (comp = 0; comp < NcC; ++comp) {
1528:             cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1529:           }
1530:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1531:       }
1532:     }
1533:     offsetC += cpdim*NcC;
1534:     offsetF += fpdim*NcF;
1535:   }
1536:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
1537:   PetscFree2(feRef,fvRef);

1539:   DMGetGlobalVector(dmf, &fv);
1540:   DMGetGlobalVector(dmc, &cv);
1541:   VecGetOwnershipRange(cv, &startC, NULL);
1542:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1543:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1544:   PetscMalloc1(m,&cindices);
1545:   PetscMalloc1(m,&findices);
1546:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1547:   for (c = cStart; c < cEnd; ++c) {
1548:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1549:     for (d = 0; d < cTotDim; ++d) {
1550:       if (cellCIndices[d] < 0) continue;
1551:       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]]);
1552:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1553:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1554:     }
1555:   }
1556:   PetscFree(cmap);
1557:   PetscFree2(cellCIndices,cellFIndices);

1559:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1560:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1561:   VecScatterCreate(cv, cis, fv, fis, sc);
1562:   ISDestroy(&cis);
1563:   ISDestroy(&fis);
1564:   DMRestoreGlobalVector(dmf, &fv);
1565:   DMRestoreGlobalVector(dmc, &cv);
1566:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1567:   return(0);
1568: }