Actual source code: plexfem.c

petsc-master 2015-08-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, 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, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody;

159:     ctx[0] = dimEmbed;
160:     ctx[1] = d;
161:     DMPlexProjectFunction(dm, &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(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
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(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
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 DMPlexProjectFunctionLabelLocal(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, 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:           PetscObjectReference((PetscObject) sp[f]);
278:         } else {
279:           PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
280:           if (!sp[f]) continue;
281:         }
282:       } else if (id == PETSCFV_CLASSID) {
283:         PetscFV fv = (PetscFV) obj;

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

303:       DMLabelGetStratumIS(label, ids[i], &pointIS);
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, &geom, numComp[f], funcs[f], ctx, &values[v]);
322:             } else {
323:               for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
324:             }
325:             v += numComp[f];
326:           }
327:         }
328:         DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);
329:       }
330:       ISRestoreIndices(pointIS, &points);
331:       ISDestroy(&pointIS);
332:     }
333:     if (h) {
334:       for (f = 0; f < numFields; ++f) {PetscDualSpaceDestroy(&sp[f]);}
335:     }
336:   }
337:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
338:   DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);
339:   for (f = 0; f < numFields; ++f) {PetscDualSpaceDestroy(&sp[f]);}
340:   PetscFree2(sp, numComp);
341:   if (maxHeight > 0) {
342:     PetscFree(cellsp);
343:   }
344:   return(0);
345: }

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

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

385:       DMGetField(dm, f, &obj);
386:       PetscObjectGetClassId(obj, &id);
387:       if (id == PETSCFE_CLASSID) {
388:         PetscFE fe = (PetscFE) obj;

390:         PetscFEGetNumComponents(fe, &numComp[f]);
391:         if (!h) {
392:           PetscFEGetDualSpace(fe, &cellsp[f]);
393:           sp[f] = cellsp[f];
394:           PetscObjectReference((PetscObject) sp[f]);
395:         }
396:         else {
397:           PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
398:           if (!sp[f]) {
399:             continue;
400:           }
401:         }
402:       } else if (id == PETSCFV_CLASSID) {
403:         PetscFV fv = (PetscFV) obj;

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

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

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

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

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

463:   Output Parameter:
464: . X - vector

466:    Calling sequence of func:
467: $    func(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

469: +  dim - The spatial dimension
470: .  x   - The coordinates
471: .  Nf  - The number of fields
472: .  u   - The output field values
473: -  ctx - optional user-defined function context

475:   Level: developer

477: .seealso: DMPlexComputeL2Diff()
478: @*/
479: PetscErrorCode DMPlexProjectFunction(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
480: {
481:   Vec            localX;

486:   DMGetLocalVector(dm, &localX);
487:   DMPlexProjectFunctionLocal(dm, funcs, ctxs, mode, localX);
488:   DMLocalToGlobalBegin(dm, localX, mode, X);
489:   DMLocalToGlobalEnd(dm, localX, mode, X);
490:   DMRestoreLocalVector(dm, &localX);
491:   return(0);
492: }

496: 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)
497: {
498:   DM                dmAux;
499:   PetscDS           prob, probAux = NULL;
500:   Vec               A;
501:   PetscSection      section, sectionAux = NULL;
502:   PetscDualSpace   *sp;
503:   PetscInt         *Ncf;
504:   PetscScalar      *values, *u, *u_x, *a, *a_x;
505:   PetscReal        *x, *v0, *J, *invJ, detJ;
506:   PetscInt          Nf, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight;
507:   PetscErrorCode    ierr;

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

538:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
539:     DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);
540:     if (dmAux) {DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
541:     for (f = 0, v = 0; f < Nf; ++f) {
542:       PetscObject  obj;
543:       PetscClassId id;

545:       PetscDSGetDiscretization(prob, f, &obj);
546:       PetscObjectGetClassId(obj, &id);
547:       if (id == PETSCFE_CLASSID) {
548:         PetscFE fe = (PetscFE) obj;

550:         PetscFEGetDualSpace(fe, &sp[f]);
551:         PetscObjectReference((PetscObject) sp[f]);
552:         PetscFEGetNumComponents(fe, &Ncf[f]);
553:       } else if (id == PETSCFV_CLASSID) {
554:         PetscFV fv = (PetscFV) obj;

556:         PetscFVGetNumComponents(fv, &Ncf[f]);
557:         PetscFVGetDualSpace(fv, &sp[f]);
558:         PetscObjectReference((PetscObject) sp[f]);
559:       }
560:       PetscDualSpaceGetDimension(sp[f], &spDim);
561:       for (d = 0; d < spDim; ++d) {
562:         PetscQuadrature  quad;
563:         const PetscReal *points, *weights;
564:         PetscInt         numPoints, q;

566:         if (funcs[f]) {
567:           PetscDualSpaceGetFunctional(sp[f], d, &quad);
568:           PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);
569:           for (q = 0; q < numPoints; ++q) {
570:             CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x);
571:             EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, coefficients,    NULL, u, u_x, NULL);
572:             EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);
573:             (*funcs[f])(u, NULL, u_x, a, NULL, a_x, x, &values[v]);
574:           }
575:         } else {
576:           for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0;
577:         }
578:         v += Ncf[f];
579:       }
580:     }
581:     DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);
582:     if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
583:     DMPlexVecSetClosure(dm, section, localX, c, values, mode);
584:   }
585:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
586:   PetscFree3(v0,J,invJ);
587:   for (f = 0; f < Nf; f++) {PetscDualSpaceDestroy(&sp[f]);}
588:   PetscFree2(sp, Ncf);
589:   return(0);
590: }

594: static PetscErrorCode DMPlexInsertBoundaryValues_FEM_Internal(DM dm, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX)
595: {
596:   PetscErrorCode (**funcs)(PetscInt, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
597:   void            **ctxs;
598:   PetscInt          numFields;
599:   PetscErrorCode    ierr;

602:   DMGetNumFields(dm, &numFields);
603:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
604:   funcs[field] = func;
605:   ctxs[field]  = ctx;
606:   DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);
607:   PetscFree2(funcs,ctxs);
608:   return(0);
609: }

613: /* This ignores numcomps/comps */
614: static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad,
615:                                                               PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
616: {
617:   PetscDS            prob;
618:   PetscSF            sf;
619:   DM                 dmFace, dmCell, dmGrad;
620:   const PetscScalar *facegeom, *cellgeom, *grad;
621:   const PetscInt    *leaves;
622:   PetscScalar       *x, *fx;
623:   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
624:   PetscErrorCode     ierr;

627:   DMGetPointSF(dm, &sf);
628:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
629:   nleaves = PetscMax(0, nleaves);
630:   DMGetDimension(dm, &dim);
631:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
632:   DMGetDS(dm, &prob);
633:   VecGetDM(faceGeometry, &dmFace);
634:   VecGetArrayRead(faceGeometry, &facegeom);
635:   VecGetDM(cellGeometry, &dmCell);
636:   VecGetArrayRead(cellGeometry, &cellgeom);
637:   if (Grad) {
638:     PetscFV fv;

640:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
641:     VecGetDM(Grad, &dmGrad);
642:     VecGetArrayRead(Grad, &grad);
643:     PetscFVGetNumComponents(fv, &pdim);
644:     DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);
645:   }
646:   VecGetArray(locX, &x);
647:   for (i = 0; i < numids; ++i) {
648:     IS              faceIS;
649:     const PetscInt *faces;
650:     PetscInt        numFaces, f;

652:     DMLabelGetStratumIS(label, ids[i], &faceIS);
653:     if (!faceIS) continue; /* No points with that id on this process */
654:     ISGetLocalSize(faceIS, &numFaces);
655:     ISGetIndices(faceIS, &faces);
656:     for (f = 0; f < numFaces; ++f) {
657:       const PetscInt         face = faces[f], *cells;
658:       const PetscFVFaceGeom *fg;

660:       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
661:       PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
662:       if (loc >= 0) continue;
663:       DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
664:       DMPlexGetSupport(dm, face, &cells);
665:       if (Grad) {
666:         const PetscFVCellGeom *cg;
667:         const PetscScalar     *cx, *cgrad;
668:         PetscScalar           *xG;
669:         PetscReal              dx[3];
670:         PetscInt               d;

672:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
673:         DMPlexPointLocalRead(dm, cells[0], x, &cx);
674:         DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
675:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
676:         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
677:         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
678:         (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
679:       } else {
680:         const PetscScalar *xI;
681:         PetscScalar       *xG;

683:         DMPlexPointLocalRead(dm, cells[0], x, &xI);
684:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
685:         (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
686:       }
687:     }
688:     ISRestoreIndices(faceIS, &faces);
689:     ISDestroy(&faceIS);
690:   }
691:   VecRestoreArray(locX, &x);
692:   if (Grad) {
693:     DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);
694:     VecRestoreArrayRead(Grad, &grad);
695:   }
696:   VecRestoreArrayRead(cellGeometry, &cellgeom);
697:   VecRestoreArrayRead(faceGeometry, &facegeom);
698:   return(0);
699: }

703: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
704: {
705:   PetscInt       numBd, b;

714:   DMPlexGetNumBoundary(dm, &numBd);
715:   for (b = 0; b < numBd; ++b) {
716:     PetscBool       isEssential;
717:     const char     *labelname;
718:     DMLabel         label;
719:     PetscInt        field;
720:     PetscObject     obj;
721:     PetscClassId    id;
722:     void          (*func)();
723:     PetscInt        numids;
724:     const PetscInt *ids;
725:     void           *ctx;

727:     DMPlexGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);
728:     if (!isEssential) continue;
729:     DMPlexGetLabel(dm, labelname, &label);
730:     DMGetField(dm, field, &obj);
731:     PetscObjectGetClassId(obj, &id);
732:     if (id == PETSCFE_CLASSID) {
733:       DMPlexInsertBoundaryValues_FEM_Internal(dm, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
734:     } else if (id == PETSCFV_CLASSID) {
735:       if (!faceGeomFVM) continue;
736:       DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM,
737:                                                      field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
738:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
739:   }
740:   return(0);
741: }

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

748:   Input Parameters:
749: + dm    - The DM
750: . funcs - The functions to evaluate for each field component
751: . ctxs  - Optional array of contexts to pass to each function, or NULL.
752: - X     - The coefficient vector u_h

754:   Output Parameter:
755: . diff - The diff ||u - u_h||_2

757:   Level: developer

759: .seealso: DMPlexProjectFunction(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff()
760: @*/
761: PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
762: {
763:   const PetscInt   debug = 0;
764:   PetscSection     section;
765:   PetscQuadrature  quad;
766:   Vec              localX;
767:   PetscScalar     *funcVal, *interpolant;
768:   PetscReal       *coords, *v0, *J, *invJ, detJ;
769:   PetscReal        localDiff = 0.0;
770:   const PetscReal *quadPoints, *quadWeights;
771:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
772:   PetscErrorCode   ierr;

775:   DMGetDimension(dm, &dim);
776:   DMGetDefaultSection(dm, &section);
777:   PetscSectionGetNumFields(section, &numFields);
778:   DMGetLocalVector(dm, &localX);
779:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
780:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
781:   for (field = 0; field < numFields; ++field) {
782:     PetscObject  obj;
783:     PetscClassId id;
784:     PetscInt     Nc;

786:     DMGetField(dm, field, &obj);
787:     PetscObjectGetClassId(obj, &id);
788:     if (id == PETSCFE_CLASSID) {
789:       PetscFE fe = (PetscFE) obj;

791:       PetscFEGetQuadrature(fe, &quad);
792:       PetscFEGetNumComponents(fe, &Nc);
793:     } else if (id == PETSCFV_CLASSID) {
794:       PetscFV fv = (PetscFV) obj;

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

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

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

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

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

860:   Input Parameters:
861: + dm    - The DM
862: . funcs - The gradient functions to evaluate for each field component
863: . ctxs  - Optional array of contexts to pass to each function, or NULL.
864: . X     - The coefficient vector u_h
865: - n     - The vector to project along

867:   Output Parameter:
868: . diff - The diff ||(grad u - grad u_h) . n||_2

870:   Level: developer

872: .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
873: @*/
874: PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
875: {
876:   const PetscInt  debug = 0;
877:   PetscSection    section;
878:   PetscQuadrature quad;
879:   Vec             localX;
880:   PetscScalar    *funcVal, *interpolantVec;
881:   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
882:   PetscReal       localDiff = 0.0;
883:   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp;
884:   PetscErrorCode  ierr;

887:   DMGetDimension(dm, &dim);
888:   DMGetDefaultSection(dm, &section);
889:   PetscSectionGetNumFields(section, &numFields);
890:   DMGetLocalVector(dm, &localX);
891:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
892:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
893:   for (field = 0; field < numFields; ++field) {
894:     PetscFE  fe;
895:     PetscInt Nc;

897:     DMGetField(dm, field, (PetscObject *) &fe);
898:     PetscFEGetQuadrature(fe, &quad);
899:     PetscFEGetNumComponents(fe, &Nc);
900:     numComponents += Nc;
901:   }
902:   /* DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
903:   PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
904:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
905:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
906:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
907:   for (c = cStart; c < cEnd; ++c) {
908:     PetscScalar *x = NULL;
909:     PetscReal    elemDiff = 0.0;

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

915:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
916:       PetscFE          fe;
917:       void * const     ctx = ctxs ? ctxs[field] : NULL;
918:       const PetscReal *quadPoints, *quadWeights;
919:       PetscReal       *basisDer;
920:       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;

922:       DMGetField(dm, field, (PetscObject *) &fe);
923:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
924:       PetscFEGetDimension(fe, &Nb);
925:       PetscFEGetNumComponents(fe, &Ncomp);
926:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
927:       if (debug) {
928:         char title[1024];
929:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
930:         DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);
931:       }
932:       for (q = 0; q < numQuadPoints; ++q) {
933:         for (d = 0; d < dim; d++) {
934:           coords[d] = v0[d];
935:           for (e = 0; e < dim; e++) {
936:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
937:           }
938:         }
939:         (*funcs[field])(dim, coords, n, numFields, funcVal, ctx);
940:         for (fc = 0; fc < Ncomp; ++fc) {
941:           PetscScalar interpolant = 0.0;

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

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

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

979:   Input Parameters:
980: + dm    - The DM
981: . funcs - The functions to evaluate for each field component
982: . ctxs  - Optional array of contexts to pass to each function, or NULL.
983: - X     - The coefficient vector u_h

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

988:   Level: developer

990: .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2GradientDiff()
991: @*/
992: PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
993: {
994:   const PetscInt   debug = 0;
995:   PetscSection     section;
996:   PetscQuadrature  quad;
997:   Vec              localX;
998:   PetscScalar     *funcVal, *interpolant;
999:   PetscReal       *coords, *v0, *J, *invJ, detJ;
1000:   PetscReal       *localDiff;
1001:   const PetscReal *quadPoints, *quadWeights;
1002:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1003:   PetscErrorCode   ierr;

1006:   DMGetDimension(dm, &dim);
1007:   DMGetDefaultSection(dm, &section);
1008:   PetscSectionGetNumFields(section, &numFields);
1009:   DMGetLocalVector(dm, &localX);
1010:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1011:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1012:   for (field = 0; field < numFields; ++field) {
1013:     PetscObject  obj;
1014:     PetscClassId id;
1015:     PetscInt     Nc;

1017:     DMGetField(dm, field, &obj);
1018:     PetscObjectGetClassId(obj, &id);
1019:     if (id == PETSCFE_CLASSID) {
1020:       PetscFE fe = (PetscFE) obj;

1022:       PetscFEGetQuadrature(fe, &quad);
1023:       PetscFEGetNumComponents(fe, &Nc);
1024:     } else if (id == PETSCFV_CLASSID) {
1025:       PetscFV fv = (PetscFV) obj;

1027:       PetscFVGetQuadrature(fv, &quad);
1028:       PetscFVGetNumComponents(fv, &Nc);
1029:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1030:     numComponents += Nc;
1031:   }
1032:   PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1033:   DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);
1034:   PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1035:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1036:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1037:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1038:   for (c = cStart; c < cEnd; ++c) {
1039:     PetscScalar *x = NULL;

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

1045:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1046:       PetscObject  obj;
1047:       PetscClassId id;
1048:       void * const ctx = ctxs ? ctxs[field] : NULL;
1049:       PetscInt     Nb, Nc, q, fc;

1051:       PetscReal       elemDiff = 0.0;

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

1088: /*@C
1089:   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.

1091:   Input Parameters:
1092: + dm    - The DM
1093: . funcs - The functions to evaluate for each field component
1094: . ctxs  - Optional array of contexts to pass to each function, or NULL.
1095: - X     - The coefficient vector u_h

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

1100:   Level: developer

1102: .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff()
1103: @*/
1104: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscErrorCode (**funcs)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
1105: {
1106:   PetscSection     section;
1107:   PetscQuadrature  quad;
1108:   Vec              localX;
1109:   PetscScalar     *funcVal, *interpolant;
1110:   PetscReal       *coords, *v0, *J, *invJ, detJ;
1111:   const PetscReal *quadPoints, *quadWeights;
1112:   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1113:   PetscErrorCode   ierr;

1116:   VecSet(D, 0.0);
1117:   DMGetDimension(dm, &dim);
1118:   DMGetDefaultSection(dm, &section);
1119:   PetscSectionGetNumFields(section, &numFields);
1120:   DMGetLocalVector(dm, &localX);
1121:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1122:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1123:   for (field = 0; field < numFields; ++field) {
1124:     PetscObject  obj;
1125:     PetscClassId id;
1126:     PetscInt     Nc;

1128:     DMGetField(dm, field, &obj);
1129:     PetscObjectGetClassId(obj, &id);
1130:     if (id == PETSCFE_CLASSID) {
1131:       PetscFE fe = (PetscFE) obj;

1133:       PetscFEGetQuadrature(fe, &quad);
1134:       PetscFEGetNumComponents(fe, &Nc);
1135:     } else if (id == PETSCFV_CLASSID) {
1136:       PetscFV fv = (PetscFV) obj;

1138:       PetscFVGetQuadrature(fv, &quad);
1139:       PetscFVGetNumComponents(fv, &Nc);
1140:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1141:     numComponents += Nc;
1142:   }
1143:   PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
1144:   DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);
1145:   PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1146:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1147:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1148:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1149:   for (c = cStart; c < cEnd; ++c) {
1150:     PetscScalar *x = NULL;
1151:     PetscScalar  elemDiff = 0.0;

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

1157:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1158:       PetscObject  obj;
1159:       PetscClassId id;
1160:       void * const ctx = ctxs ? ctxs[field] : NULL;
1161:       PetscInt     Nb, Nc, q, fc;

1163:       DMGetField(dm, field, &obj);
1164:       PetscObjectGetClassId(obj, &id);
1165:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1166:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1167:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1168:       for (q = 0; q < numQuadPoints; ++q) {
1169:         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
1170:         (*funcs[field])(dim, coords, Nc, funcVal, ctx);
1171:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);}
1172:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1173:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1174:         for (fc = 0; fc < Nc; ++fc) {
1175:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
1176:         }
1177:       }
1178:       fieldOffset += Nb*Nc;
1179:     }
1180:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1181:     VecSetValuesSection(D, section, c, &elemDiff, ADD_VALUES);
1182:   }
1183:   PetscFree6(funcVal,interpolant,coords,v0,J,invJ);
1184:   DMRestoreLocalVector(dm, &localX);
1185:   VecSqrtAbs(D);
1186:   return(0);
1187: }

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

1194:   Input Parameters:
1195: + dm - The mesh
1196: . X  - Local input vector
1197: - user - The user context

1199:   Output Parameter:
1200: . integral - Local integral for each field

1202:   Level: developer

1204: .seealso: DMPlexComputeResidualFEM()
1205: @*/
1206: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
1207: {
1208:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1209:   DM                dmAux;
1210:   Vec               localX, A;
1211:   PetscDS           prob, probAux = NULL;
1212:   PetscSection      section, sectionAux;
1213:   PetscFECellGeom  *cgeom;
1214:   PetscScalar      *u, *a = NULL;
1215:   PetscReal        *lintegral, *vol;
1216:   PetscInt          dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c;
1217:   PetscInt          totDim, totDimAux;
1218:   PetscErrorCode    ierr;

1221:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1222:   DMGetLocalVector(dm, &localX);
1223:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1224:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1225:   DMGetDimension(dm, &dim);
1226:   DMGetDefaultSection(dm, &section);
1227:   DMGetDS(dm, &prob);
1228:   PetscDSGetTotalDimension(prob, &totDim);
1229:   PetscSectionGetNumFields(section, &Nf);
1230:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1231:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1232:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1233:   numCells = cEnd - cStart;
1234:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1235:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1236:   if (dmAux) {
1237:     DMGetDefaultSection(dmAux, &sectionAux);
1238:     DMGetDS(dmAux, &probAux);
1239:     PetscDSGetTotalDimension(probAux, &totDimAux);
1240:   }
1241:   DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);
1242:   PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);
1243:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1244:   for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
1245:   for (c = cStart; c < cEnd; ++c) {
1246:     PetscScalar *x = NULL;
1247:     PetscInt     i;

1249:     DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);
1250:     DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);
1251:     if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c);
1252:     DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
1253:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1254:     DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
1255:     if (dmAux) {
1256:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1257:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1258:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1259:     }
1260:   }
1261:   for (f = 0; f < Nf; ++f) {
1262:     PetscObject  obj;
1263:     PetscClassId id;
1264:     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1266:     PetscDSGetDiscretization(prob, f, &obj);
1267:     PetscObjectGetClassId(obj, &id);
1268:     if (id == PETSCFE_CLASSID) {
1269:       PetscFE         fe = (PetscFE) obj;
1270:       PetscQuadrature q;
1271:       PetscInt        Nq, Nb;

1273:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1274:       PetscFEGetQuadrature(fe, &q);
1275:       PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1276:       PetscFEGetDimension(fe, &Nb);
1277:       blockSize = Nb*Nq;
1278:       batchSize = numBlocks * blockSize;
1279:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1280:       numChunks = numCells / (numBatches*batchSize);
1281:       Ne        = numChunks*numBatches*batchSize;
1282:       Nr        = numCells % (numBatches*batchSize);
1283:       offset    = numCells - Nr;
1284:       PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);
1285:       PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);
1286:     } else if (id == PETSCFV_CLASSID) {
1287:       /* PetscFV  fv = (PetscFV) obj; */
1288:       PetscInt       foff;
1289:       PetscPointFunc obj_func;
1290:       PetscScalar    lint;

1292:       PetscDSGetObjective(prob, f, &obj_func);
1293:       PetscDSGetFieldOffset(prob, f, &foff);
1294:       if (obj_func) {
1295:         for (c = 0; c < numCells; ++c) {
1296:           /* TODO: Need full pointwise interpolation and get centroid for x argument */
1297:           obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint);
1298:           lintegral[f] = PetscRealPart(lint)*vol[c];
1299:         }
1300:       }
1301:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1302:   }
1303:   if (dmAux) {PetscFree(a);}
1304:   if (mesh->printFEM) {
1305:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1306:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);}
1307:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1308:   }
1309:   DMRestoreLocalVector(dm, &localX);
1310:   MPI_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1311:   PetscFree4(lintegral,u,cgeom,vol);
1312:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1313:   return(0);
1314: }

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

1321:   Input Parameters:
1322: + dmf  - The fine mesh
1323: . dmc  - The coarse mesh
1324: - user - The user context

1326:   Output Parameter:
1327: . In  - The interpolation matrix

1329:   Level: developer

1331: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
1332: @*/
1333: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
1334: {
1335:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1336:   const char       *name  = "Interpolator";
1337:   PetscDS           prob;
1338:   PetscFE          *feRef;
1339:   PetscFV          *fvRef;
1340:   PetscSection      fsection, fglobalSection;
1341:   PetscSection      csection, cglobalSection;
1342:   PetscScalar      *elemMat;
1343:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1344:   PetscInt          cTotDim, rTotDim = 0;
1345:   PetscErrorCode    ierr;

1348:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1349:   DMGetDimension(dmf, &dim);
1350:   DMGetDefaultSection(dmf, &fsection);
1351:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1352:   DMGetDefaultSection(dmc, &csection);
1353:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1354:   PetscSectionGetNumFields(fsection, &Nf);
1355:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1356:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1357:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1358:   DMGetDS(dmf, &prob);
1359:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1360:   for (f = 0; f < Nf; ++f) {
1361:     PetscObject  obj;
1362:     PetscClassId id;
1363:     PetscInt     rNb = 0, Nc = 0;

1365:     PetscDSGetDiscretization(prob, f, &obj);
1366:     PetscObjectGetClassId(obj, &id);
1367:     if (id == PETSCFE_CLASSID) {
1368:       PetscFE fe = (PetscFE) obj;

1370:       PetscFERefine(fe, &feRef[f]);
1371:       PetscFEGetDimension(feRef[f], &rNb);
1372:       PetscFEGetNumComponents(fe, &Nc);
1373:     } else if (id == PETSCFV_CLASSID) {
1374:       PetscFV        fv = (PetscFV) obj;
1375:       PetscDualSpace Q;

1377:       PetscFVRefine(fv, &fvRef[f]);
1378:       PetscFVGetDualSpace(fvRef[f], &Q);
1379:       PetscDualSpaceGetDimension(Q, &rNb);
1380:       PetscFVGetNumComponents(fv, &Nc);
1381:     }
1382:     rTotDim += rNb*Nc;
1383:   }
1384:   PetscDSGetTotalDimension(prob, &cTotDim);
1385:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1386:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1387:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1388:     PetscDualSpace   Qref;
1389:     PetscQuadrature  f;
1390:     const PetscReal *qpoints, *qweights;
1391:     PetscReal       *points;
1392:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1394:     /* Compose points from all dual basis functionals */
1395:     if (feRef[fieldI]) {
1396:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
1397:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
1398:     } else {
1399:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1400:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1401:     }
1402:     PetscDualSpaceGetDimension(Qref, &fpdim);
1403:     for (i = 0; i < fpdim; ++i) {
1404:       PetscDualSpaceGetFunctional(Qref, i, &f);
1405:       PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1406:       npoints += Np;
1407:     }
1408:     PetscMalloc1(npoints*dim,&points);
1409:     for (i = 0, k = 0; i < fpdim; ++i) {
1410:       PetscDualSpaceGetFunctional(Qref, i, &f);
1411:       PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1412:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1413:     }

1415:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1416:       PetscObject  obj;
1417:       PetscClassId id;
1418:       PetscReal   *B;
1419:       PetscInt     NcJ = 0, cpdim = 0, j;

1421:       PetscDSGetDiscretization(prob, fieldJ, &obj);
1422:       PetscObjectGetClassId(obj, &id);
1423:       if (id == PETSCFE_CLASSID) {
1424:         PetscFE fe = (PetscFE) obj;

1426:         /* Evaluate basis at points */
1427:         PetscFEGetNumComponents(fe, &NcJ);
1428:         PetscFEGetDimension(fe, &cpdim);
1429:         /* For now, fields only interpolate themselves */
1430:         if (fieldI == fieldJ) {
1431:           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);
1432:           PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1433:           for (i = 0, k = 0; i < fpdim; ++i) {
1434:             PetscDualSpaceGetFunctional(Qref, i, &f);
1435:             PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1436:             for (p = 0; p < Np; ++p, ++k) {
1437:               for (j = 0; j < cpdim; ++j) {
1438:                 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];
1439:               }
1440:             }
1441:           }
1442:           PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1443:         }
1444:       } else if (id == PETSCFV_CLASSID) {
1445:         PetscFV        fv = (PetscFV) obj;

1447:         /* Evaluate constant function at points */
1448:         PetscFVGetNumComponents(fv, &NcJ);
1449:         cpdim = 1;
1450:         /* For now, fields only interpolate themselves */
1451:         if (fieldI == fieldJ) {
1452:           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);
1453:           for (i = 0, k = 0; i < fpdim; ++i) {
1454:             PetscDualSpaceGetFunctional(Qref, i, &f);
1455:             PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1456:             for (p = 0; p < Np; ++p, ++k) {
1457:               for (j = 0; j < cpdim; ++j) {
1458:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p];
1459:               }
1460:             }
1461:           }
1462:         }
1463:       }
1464:       offsetJ += cpdim*NcJ;
1465:     }
1466:     offsetI += fpdim*Nc;
1467:     PetscFree(points);
1468:   }
1469:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1470:   /* Preallocate matrix */
1471:   {
1472:     PetscHashJK ht;
1473:     PetscLayout rLayout;
1474:     PetscInt   *dnz, *onz, *cellCIndices, *cellFIndices;
1475:     PetscInt    locRows, rStart, rEnd, cell, r;

1477:     MatGetLocalSize(In, &locRows, NULL);
1478:     PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1479:     PetscLayoutSetLocalSize(rLayout, locRows);
1480:     PetscLayoutSetBlockSize(rLayout, 1);
1481:     PetscLayoutSetUp(rLayout);
1482:     PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1483:     PetscLayoutDestroy(&rLayout);
1484:     PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1485:     PetscHashJKCreate(&ht);
1486:     for (cell = cStart; cell < cEnd; ++cell) {
1487:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1488:       for (r = 0; r < rTotDim; ++r) {
1489:         PetscHashJKKey  key;
1490:         PetscHashJKIter missing, iter;

1492:         key.j = cellFIndices[r];
1493:         if (key.j < 0) continue;
1494:         for (c = 0; c < cTotDim; ++c) {
1495:           key.k = cellCIndices[c];
1496:           if (key.k < 0) continue;
1497:           PetscHashJKPut(ht, key, &missing, &iter);
1498:           if (missing) {
1499:             PetscHashJKSet(ht, iter, 1);
1500:             if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1501:             else                                     ++onz[key.j-rStart];
1502:           }
1503:         }
1504:       }
1505:     }
1506:     PetscHashJKDestroy(&ht);
1507:     MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1508:     MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1509:     PetscFree4(dnz,onz,cellCIndices,cellFIndices);
1510:   }
1511:   /* Fill matrix */
1512:   MatZeroEntries(In);
1513:   for (c = cStart; c < cEnd; ++c) {
1514:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1515:   }
1516:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1517:   PetscFree2(feRef,fvRef);
1518:   PetscFree(elemMat);
1519:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1520:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1521:   if (mesh->printFEM) {
1522:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1523:     MatChop(In, 1.0e-10);
1524:     MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1525:   }
1526:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1527:   return(0);
1528: }

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

1535:   Input Parameters:
1536: + dmf  - The fine mesh
1537: . dmc  - The coarse mesh
1538: - user - The user context

1540:   Output Parameter:
1541: . In  - The interpolation matrix

1543:   Level: developer

1545: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
1546: @*/
1547: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
1548: {
1549:   PetscDS        prob;
1550:   PetscSection   fsection, csection, globalFSection, globalCSection;
1551:   PetscHashJK    ht;
1552:   PetscLayout    rLayout;
1553:   PetscInt      *dnz, *onz;
1554:   PetscInt       locRows, rStart, rEnd;
1555:   PetscReal     *x, *v0, *J, *invJ, detJ;
1556:   PetscReal     *v0c, *Jc, *invJc, detJc;
1557:   PetscScalar   *elemMat;
1558:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

1562:   DMGetCoordinateDim(dmc, &dim);
1563:   DMGetDS(dmc, &prob);
1564:   PetscDSGetRefCoordArrays(prob, &x, NULL);
1565:   PetscDSGetNumFields(prob, &Nf);
1566:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
1567:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
1568:   DMGetDefaultSection(dmf, &fsection);
1569:   DMGetDefaultGlobalSection(dmf, &globalFSection);
1570:   DMGetDefaultSection(dmc, &csection);
1571:   DMGetDefaultGlobalSection(dmc, &globalCSection);
1572:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
1573:   PetscDSGetTotalDimension(prob, &totDim);
1574:   PetscMalloc1(totDim*totDim, &elemMat);

1576:   MatGetLocalSize(In, &locRows, NULL);
1577:   PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1578:   PetscLayoutSetLocalSize(rLayout, locRows);
1579:   PetscLayoutSetBlockSize(rLayout, 1);
1580:   PetscLayoutSetUp(rLayout);
1581:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1582:   PetscLayoutDestroy(&rLayout);
1583:   PetscCalloc2(locRows,&dnz,locRows,&onz);
1584:   PetscHashJKCreate(&ht);
1585:   for (field = 0; field < Nf; ++field) {
1586:     PetscObject      obj;
1587:     PetscClassId     id;
1588:     PetscDualSpace   Q = NULL;
1589:     PetscQuadrature  f;
1590:     const PetscReal *qpoints;
1591:     PetscInt         Nc, Np, fpdim, i, d;

1593:     PetscDSGetDiscretization(prob, field, &obj);
1594:     PetscObjectGetClassId(obj, &id);
1595:     if (id == PETSCFE_CLASSID) {
1596:       PetscFE fe = (PetscFE) obj;

1598:       PetscFEGetDualSpace(fe, &Q);
1599:       PetscFEGetNumComponents(fe, &Nc);
1600:     } else if (id == PETSCFV_CLASSID) {
1601:       PetscFV fv = (PetscFV) obj;

1603:       PetscFVGetDualSpace(fv, &Q);
1604:       Nc   = 1;
1605:     }
1606:     PetscDualSpaceGetDimension(Q, &fpdim);
1607:     /* For each fine grid cell */
1608:     for (cell = cStart; cell < cEnd; ++cell) {
1609:       PetscInt *findices,   *cindices;
1610:       PetscInt  numFIndices, numCIndices;

1612:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);
1613:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1614:       if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc);
1615:       for (i = 0; i < fpdim; ++i) {
1616:         Vec             pointVec;
1617:         PetscScalar    *pV;
1618:         IS              coarseCellIS;
1619:         const PetscInt *coarseCells;
1620:         PetscInt        numCoarseCells, q, r, c;

1622:         /* Get points from the dual basis functional quadrature */
1623:         PetscDualSpaceGetFunctional(Q, i, &f);
1624:         PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1625:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1626:         VecSetBlockSize(pointVec, dim);
1627:         VecGetArray(pointVec, &pV);
1628:         for (q = 0; q < Np; ++q) {
1629:           /* Transform point to real space */
1630:           CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1631:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1632:         }
1633:         VecRestoreArray(pointVec, &pV);
1634:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1635:         DMLocatePoints(dmc, pointVec, &coarseCellIS);
1636:         ISViewFromOptions(coarseCellIS, NULL, "-interp_is_view");
1637:         /* Update preallocation info */
1638:         ISGetLocalSize(coarseCellIS, &numCoarseCells);
1639:         ISGetIndices(coarseCellIS, &coarseCells);
1640:         for (r = 0; r < Nc; ++r) {
1641:           PetscHashJKKey  key;
1642:           PetscHashJKIter missing, iter;

1644:           key.j = findices[i*Nc+r];
1645:           if (key.j < 0) continue;
1646:           /* Get indices for coarse elements */
1647:           for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1648:             DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);
1649:             for (c = 0; c < numCIndices; ++c) {
1650:               key.k = cindices[c];
1651:               if (key.k < 0) continue;
1652:               PetscHashJKPut(ht, key, &missing, &iter);
1653:               if (missing) {
1654:                 PetscHashJKSet(ht, iter, 1);
1655:                 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1656:                 else                                     ++onz[key.j-rStart];
1657:               }
1658:             }
1659:             DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);
1660:           }
1661:         }
1662:         ISRestoreIndices(coarseCellIS, &coarseCells);
1663:         ISDestroy(&coarseCellIS);
1664:         VecDestroy(&pointVec);
1665:       }
1666:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);
1667:     }
1668:   }
1669:   PetscHashJKDestroy(&ht);
1670:   MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1671:   MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1672:   PetscFree2(dnz,onz);
1673:   for (field = 0; field < Nf; ++field) {
1674:     PetscObject      obj;
1675:     PetscClassId     id;
1676:     PetscDualSpace   Q = NULL;
1677:     PetscQuadrature  f;
1678:     const PetscReal *qpoints, *qweights;
1679:     PetscInt         Nc, Np, fpdim, i, d;

1681:     PetscDSGetDiscretization(prob, field, &obj);
1682:     PetscObjectGetClassId(obj, &id);
1683:     if (id == PETSCFE_CLASSID) {
1684:       PetscFE fe = (PetscFE) obj;

1686:       PetscFEGetDualSpace(fe, &Q);
1687:       PetscFEGetNumComponents(fe, &Nc);
1688:     } else if (id == PETSCFV_CLASSID) {
1689:       PetscFV fv = (PetscFV) obj;

1691:       PetscFVGetDualSpace(fv, &Q);
1692:       Nc   = 1;
1693:     }
1694:     PetscDualSpaceGetDimension(Q, &fpdim);
1695:     /* For each fine grid cell */
1696:     for (cell = cStart; cell < cEnd; ++cell) {
1697:       PetscInt *findices,   *cindices;
1698:       PetscInt  numFIndices, numCIndices;

1700:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);
1701:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
1702:       if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc);
1703:       for (i = 0; i < fpdim; ++i) {
1704:         Vec             pointVec;
1705:         PetscScalar    *pV;
1706:         IS              coarseCellIS;
1707:         const PetscInt *coarseCells;
1708:         PetscInt        numCoarseCells, cpdim, q, c, j;

1710:         /* Get points from the dual basis functional quadrature */
1711:         PetscDualSpaceGetFunctional(Q, i, &f);
1712:         PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);
1713:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
1714:         VecSetBlockSize(pointVec, dim);
1715:         VecGetArray(pointVec, &pV);
1716:         for (q = 0; q < Np; ++q) {
1717:           /* Transform point to real space */
1718:           CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x);
1719:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
1720:         }
1721:         VecRestoreArray(pointVec, &pV);
1722:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
1723:         DMLocatePoints(dmc, pointVec, &coarseCellIS);
1724:         /* Update preallocation info */
1725:         ISGetLocalSize(coarseCellIS, &numCoarseCells);
1726:         ISGetIndices(coarseCellIS, &coarseCells);
1727:         VecGetArray(pointVec, &pV);
1728:         for (ccell = 0; ccell < numCoarseCells; ++ccell) {
1729:           PetscReal pVReal[3];

1731:           DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);
1732:           /* Transform points from real space to coarse reference space */
1733:           DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell], NULL, v0c, Jc, invJc, &detJc);
1734:           CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x);
1735:           for (d = 0; d < dim; ++d) pV[ccell*dim+d] = pVReal[d];

1737:           if (id == PETSCFE_CLASSID) {
1738:             PetscFE    fe = (PetscFE) obj;
1739:             PetscReal *B;

1741:             /* Evaluate coarse basis on contained point */
1742:             PetscFEGetDimension(fe, &cpdim);
1743:             PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
1744:             /* Get elemMat entries by multiplying by weight */
1745:             for (j = 0; j < cpdim; ++j) {
1746:               for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell];
1747:             }
1748:             PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
1749:           } else {
1750:             cpdim = 1;
1751:             for (j = 0; j < cpdim; ++j) {
1752:               for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell];
1753:             }
1754:           }
1755:           /* Update interpolator */
1756:           MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);
1757:           DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell], &numCIndices, &cindices);
1758:         }
1759:         VecRestoreArray(pointVec, &pV);
1760:         ISRestoreIndices(coarseCellIS, &coarseCells);
1761:         ISDestroy(&coarseCellIS);
1762:         VecDestroy(&pointVec);
1763:       }
1764:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices);
1765:     }
1766:   }
1767:   PetscFree3(v0,J,invJ);
1768:   PetscFree3(v0c,Jc,invJc);
1769:   PetscFree(elemMat);
1770:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1771:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1772:   return(0);
1773: }

1777: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1778: {
1779:   PetscDS        prob;
1780:   PetscFE       *feRef;
1781:   PetscFV       *fvRef;
1782:   Vec            fv, cv;
1783:   IS             fis, cis;
1784:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1785:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1786:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m;

1790:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1791:   DMGetDimension(dmf, &dim);
1792:   DMGetDefaultSection(dmf, &fsection);
1793:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1794:   DMGetDefaultSection(dmc, &csection);
1795:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1796:   PetscSectionGetNumFields(fsection, &Nf);
1797:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1798:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1799:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1800:   DMGetDS(dmc, &prob);
1801:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1802:   for (f = 0; f < Nf; ++f) {
1803:     PetscObject  obj;
1804:     PetscClassId id;
1805:     PetscInt     fNb = 0, Nc = 0;

1807:     PetscDSGetDiscretization(prob, f, &obj);
1808:     PetscObjectGetClassId(obj, &id);
1809:     if (id == PETSCFE_CLASSID) {
1810:       PetscFE fe = (PetscFE) obj;

1812:       PetscFERefine(fe, &feRef[f]);
1813:       PetscFEGetDimension(feRef[f], &fNb);
1814:       PetscFEGetNumComponents(fe, &Nc);
1815:     } else if (id == PETSCFV_CLASSID) {
1816:       PetscFV        fv = (PetscFV) obj;
1817:       PetscDualSpace Q;

1819:       PetscFVRefine(fv, &fvRef[f]);
1820:       PetscFVGetDualSpace(fvRef[f], &Q);
1821:       PetscDualSpaceGetDimension(Q, &fNb);
1822:       PetscFVGetNumComponents(fv, &Nc);
1823:     }
1824:     fTotDim += fNb*Nc;
1825:   }
1826:   PetscDSGetTotalDimension(prob, &cTotDim);
1827:   PetscMalloc1(cTotDim,&cmap);
1828:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1829:     PetscFE        feC;
1830:     PetscFV        fvC;
1831:     PetscDualSpace QF, QC;
1832:     PetscInt       NcF, NcC, fpdim, cpdim;

1834:     if (feRef[field]) {
1835:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1836:       PetscFEGetNumComponents(feC, &NcC);
1837:       PetscFEGetNumComponents(feRef[field], &NcF);
1838:       PetscFEGetDualSpace(feRef[field], &QF);
1839:       PetscDualSpaceGetDimension(QF, &fpdim);
1840:       PetscFEGetDualSpace(feC, &QC);
1841:       PetscDualSpaceGetDimension(QC, &cpdim);
1842:     } else {
1843:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
1844:       PetscFVGetNumComponents(fvC, &NcC);
1845:       PetscFVGetNumComponents(fvRef[field], &NcF);
1846:       PetscFVGetDualSpace(fvRef[field], &QF);
1847:       PetscDualSpaceGetDimension(QF, &fpdim);
1848:       PetscFVGetDualSpace(fvC, &QC);
1849:       PetscDualSpaceGetDimension(QC, &cpdim);
1850:     }
1851:     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);
1852:     for (c = 0; c < cpdim; ++c) {
1853:       PetscQuadrature  cfunc;
1854:       const PetscReal *cqpoints;
1855:       PetscInt         NpC;
1856:       PetscBool        found = PETSC_FALSE;

1858:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
1859:       PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);
1860:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1861:       for (f = 0; f < fpdim; ++f) {
1862:         PetscQuadrature  ffunc;
1863:         const PetscReal *fqpoints;
1864:         PetscReal        sum = 0.0;
1865:         PetscInt         NpF, comp;

1867:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
1868:         PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);
1869:         if (NpC != NpF) continue;
1870:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1871:         if (sum > 1.0e-9) continue;
1872:         for (comp = 0; comp < NcC; ++comp) {
1873:           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1874:         }
1875:         found = PETSC_TRUE;
1876:         break;
1877:       }
1878:       if (!found) {
1879:         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1880:         if (fvRef[field]) {
1881:           PetscInt comp;
1882:           for (comp = 0; comp < NcC; ++comp) {
1883:             cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1884:           }
1885:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1886:       }
1887:     }
1888:     offsetC += cpdim*NcC;
1889:     offsetF += fpdim*NcF;
1890:   }
1891:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
1892:   PetscFree2(feRef,fvRef);

1894:   DMGetGlobalVector(dmf, &fv);
1895:   DMGetGlobalVector(dmc, &cv);
1896:   VecGetOwnershipRange(cv, &startC, NULL);
1897:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1898:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1899:   PetscMalloc1(m,&cindices);
1900:   PetscMalloc1(m,&findices);
1901:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1902:   for (c = cStart; c < cEnd; ++c) {
1903:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1904:     for (d = 0; d < cTotDim; ++d) {
1905:       if (cellCIndices[d] < 0) continue;
1906:       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]]);
1907:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1908:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1909:     }
1910:   }
1911:   PetscFree(cmap);
1912:   PetscFree2(cellCIndices,cellFIndices);

1914:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1915:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1916:   VecScatterCreate(cv, cis, fv, fis, sc);
1917:   ISDestroy(&cis);
1918:   ISDestroy(&fis);
1919:   DMRestoreGlobalVector(dmf, &fv);
1920:   DMRestoreGlobalVector(dmc, &cv);
1921:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1922:   return(0);
1923: }