Actual source code: plexfem.c

petsc-master 2015-01-26
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: /*@C
 90:   DMPlexCreateRigidBody - create rigid body modes from coordinates

 92:   Collective on DM

 94:   Input Arguments:
 95: + dm - the DM
 96: . section - the local section associated with the rigid field, or NULL for the default section
 97: - globalSection - the global section associated with the rigid field, or NULL for the default section

 99:   Output Argument:
100: . sp - the null space

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

104:   Level: advanced

106: .seealso: MatNullSpaceCreate()
107: @*/
108: PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp)
109: {
110:   MPI_Comm       comm;
111:   Vec            coordinates, localMode, mode[6];
112:   PetscSection   coordSection;
113:   PetscScalar   *coords;
114:   PetscInt       dim, vStart, vEnd, v, n, m, d, i, j;

118:   PetscObjectGetComm((PetscObject)dm,&comm);
119:   DMGetDimension(dm, &dim);
120:   if (dim == 1) {
121:     MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
122:     return(0);
123:   }
124:   if (!section)       {DMGetDefaultSection(dm, &section);}
125:   if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
126:   PetscSectionGetConstrainedStorageSize(globalSection, &n);
127:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
128:   DMGetCoordinateSection(dm, &coordSection);
129:   DMGetCoordinatesLocal(dm, &coordinates);
130:   m    = (dim*(dim+1))/2;
131:   VecCreate(comm, &mode[0]);
132:   VecSetSizes(mode[0], n, PETSC_DETERMINE);
133:   VecSetUp(mode[0]);
134:   for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
135:   /* Assume P1 */
136:   DMGetLocalVector(dm, &localMode);
137:   for (d = 0; d < dim; ++d) {
138:     PetscScalar values[3] = {0.0, 0.0, 0.0};

140:     values[d] = 1.0;
141:     VecSet(localMode, 0.0);
142:     for (v = vStart; v < vEnd; ++v) {
143:       DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);
144:     }
145:     DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);
146:     DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);
147:   }
148:   VecGetArray(coordinates, &coords);
149:   for (d = dim; d < dim*(dim+1)/2; ++d) {
150:     PetscInt i, j, k = dim > 2 ? d - dim : d;

152:     VecSet(localMode, 0.0);
153:     for (v = vStart; v < vEnd; ++v) {
154:       PetscScalar values[3] = {0.0, 0.0, 0.0};
155:       PetscInt    off;

157:       PetscSectionGetOffset(coordSection, v, &off);
158:       for (i = 0; i < dim; ++i) {
159:         for (j = 0; j < dim; ++j) {
160:           values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]);
161:         }
162:       }
163:       DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);
164:     }
165:     DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);
166:     DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);
167:   }
168:   VecRestoreArray(coordinates, &coords);
169:   DMRestoreLocalVector(dm, &localMode);
170:   for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
171:   /* Orthonormalize system */
172:   for (i = dim; i < m; ++i) {
173:     PetscScalar dots[6];

175:     VecMDot(mode[i], i, mode, dots);
176:     for (j = 0; j < i; ++j) dots[j] *= -1.0;
177:     VecMAXPY(mode[i], i, dots, mode);
178:     VecNormalize(mode[i], NULL);
179:   }
180:   MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
181:   for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
182:   return(0);
183: }

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

196:   Input Parameters:
197: + dm - the DMPlex object
198: - height - the maximum projection height >= 0

200:   Level: advanced

202: .seealso: DMPlexGetMaxProjectionHeight(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
203: @*/
204: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
205: {
206:   DM_Plex *plex = (DM_Plex *) dm->data;

210:   plex->maxProjectionHeight = height;
211:   return(0);
212: }

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

220:   Input Parameters:
221: . dm - the DMPlex object

223:   Output Parameters:
224: . height - the maximum projection height

226:   Level: intermediate

228: .seealso: DMPlexSetMaxProjectionHeight(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
229: @*/
230: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
231: {
232:   DM_Plex *plex = (DM_Plex *) dm->data;

236:   *height = plex->maxProjectionHeight;
237:   return(0);
238: }

242: PetscErrorCode DMPlexProjectFunctionLabelLocal(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
243: {
244:   PetscFE         fe;
245:   PetscDualSpace *sp, *cellsp;
246:   PetscSection    section;
247:   PetscScalar    *values;
248:   PetscBool      *fieldActive;
249:   PetscInt        numFields, numComp, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, cStart, cEnd, cEndInterior, f, d, v, i, comp, maxHeight, h;
250:   PetscErrorCode  ierr;

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

275:       DMGetField(dm, f, &obj);
276:       PetscObjectGetClassId(obj, &id);
277:       if (id == PETSCFE_CLASSID) {
278:         PetscFE fe = (PetscFE) obj;

280:         PetscFEGetNumComponents(fe, &numComp);
281:         if (!h) {
282:           PetscFEGetDualSpace(fe, &cellsp[f]);
283:           sp[f] = cellsp[f];
284:           PetscObjectReference((PetscObject) sp[f]);
285:         } else {
286:           PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
287:           if (!sp[f]) continue;
288:         }
289:       } else if (id == PETSCFV_CLASSID) {
290:         PetscFV         fv = (PetscFV) obj;
291:         PetscQuadrature q;

293:         PetscFVGetNumComponents(fv, &numComp);
294:         PetscFVGetQuadrature(fv, &q);
295:         PetscDualSpaceCreate(PetscObjectComm((PetscObject) fv), &sp[f]);
296:         PetscDualSpaceSetDM(sp[f], dm);
297:         PetscDualSpaceSetType(sp[f], PETSCDUALSPACESIMPLE);
298:         PetscDualSpaceSimpleSetDimension(sp[f], 1);
299:         PetscDualSpaceSimpleSetFunctional(sp[f], 0, q);
300:       } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
301:       PetscDualSpaceGetDimension(sp[f], &spDim);
302:       totDim += spDim*numComp;
303:     }
304:     DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);
305:     if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
306:     if (!totDim) continue;
307:     DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
308:     DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);
309:     for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
310:     for (i = 0; i < numIds; ++i) {
311:       IS              pointIS;
312:       const PetscInt *points;
313:       PetscInt        n, p;

315:       DMLabelGetStratumIS(label, ids[i], &pointIS);
316:       ISGetLocalSize(pointIS, &n);
317:       ISGetIndices(pointIS, &points);
318:       for (p = 0; p < n; ++p) {
319:         const PetscInt    point = points[p];
320:         PetscFECellGeom   geom;

322:         if ((point < pStart) || (point >= pEnd)) continue;
323:         DMPlexComputeCellGeometryFEM(dm, point, NULL, geom.v0, geom.J, NULL, &geom.detJ);
324:         geom.dim      = dim - h;
325:         geom.dimEmbed = dimEmbed;
326:         for (f = 0, v = 0; f < numFields; ++f) {
327:           void * const ctx = ctxs ? ctxs[f] : NULL;

329:           if (!sp[f]) continue;
330:           DMGetField(dm, f, (PetscObject *) &fe);
331:           PetscFEGetNumComponents(fe, &numComp);
332:           PetscDualSpaceGetDimension(sp[f], &spDim);
333:           for (d = 0; d < spDim; ++d) {
334:             if (funcs[f]) {
335:               PetscDualSpaceApply(sp[f], d, &geom, numComp, funcs[f], ctx, &values[v]);
336:             } else {
337:               for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0;
338:             }
339:             v += numComp;
340:           }
341:         }
342:         DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);
343:       }
344:       ISRestoreIndices(pointIS, &points);
345:       ISDestroy(&pointIS);
346:     }
347:     if (h) {
348:       for (f = 0; f < numFields; ++f) {PetscDualSpaceDestroy(&sp[f]);}
349:     }
350:   }
351:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
352:   DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);
353:   for (f = 0; f < numFields; ++f) {PetscDualSpaceDestroy(&sp[f]);}
354:   PetscFree(sp);
355:   if (maxHeight > 0) {
356:     PetscFree(cellsp);
357:   }
358:   return(0);
359: }

363: PetscErrorCode DMPlexProjectFunctionLocal(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
364: {
365:   PetscDualSpace *sp, *cellsp;
366:   PetscInt       *numComp;
367:   PetscSection    section;
368:   PetscScalar    *values;
369:   PetscInt        numFields, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight;
370:   PetscErrorCode  ierr;

373:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
374:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
375:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
376:   if (cEnd <= cStart) return(0);
377:   DMGetDefaultSection(dm, &section);
378:   PetscSectionGetNumFields(section, &numFields);
379:   PetscMalloc2(numFields, &sp, numFields, &numComp);
380:   DMGetDimension(dm, &dim);
381:   DMGetCoordinateDim(dm, &dimEmbed);
382:   DMPlexGetMaxProjectionHeight(dm,&maxHeight);
383:   if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
384:   if (maxHeight > 0) {
385:     PetscMalloc1(numFields,&cellsp);
386:   }
387:   else {
388:     cellsp = sp;
389:   }
390:   for (h = 0; h <= maxHeight; h++) {
391:     DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);
392:     if (!h) {pStart = cStart; pEnd = cEnd;}
393:     if (pEnd <= pStart) continue;
394:     totDim = 0;
395:     for (f = 0; f < numFields; ++f) {
396:       PetscObject  obj;
397:       PetscClassId id;

399:       DMGetField(dm, f, &obj);
400:       PetscObjectGetClassId(obj, &id);
401:       if (id == PETSCFE_CLASSID) {
402:         PetscFE fe = (PetscFE) obj;

404:         PetscFEGetNumComponents(fe, &numComp[f]);
405:         if (!h) {
406:           PetscFEGetDualSpace(fe, &cellsp[f]);
407:           sp[f] = cellsp[f];
408:           PetscObjectReference((PetscObject) sp[f]);
409:         }
410:         else {
411:           PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
412:           if (!sp[f]) {
413:             continue;
414:           }
415:         }
416:       } else if (id == PETSCFV_CLASSID) {
417:         PetscFV         fv = (PetscFV) obj;
418:         PetscQuadrature q;

420:         if (h) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, "Projection height > 0 not supported for finite volume");
421:         PetscFVGetNumComponents(fv, &numComp[f]);
422:         PetscFVGetQuadrature(fv, &q);
423:         PetscDualSpaceCreate(PetscObjectComm((PetscObject) fv), &sp[f]);
424:         PetscDualSpaceSetDM(sp[f], dm);
425:         PetscDualSpaceSetType(sp[f], PETSCDUALSPACESIMPLE);
426:         PetscDualSpaceSimpleSetDimension(sp[f], 1);
427:         PetscDualSpaceSimpleSetFunctional(sp[f], 0, q);
428:       } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
429:       PetscDualSpaceGetDimension(sp[f], &spDim);
430:       totDim += spDim*numComp[f];
431:     }
432:     DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);
433:     if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
434:     if (!totDim) continue;
435:     DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
436:     for (p = pStart; p < pEnd; ++p) {
437:       PetscFECellGeom geom;

439:       DMPlexComputeCellGeometryFEM(dm, p, NULL, geom.v0, geom.J, NULL, &geom.detJ);
440:       geom.dim      = dim - h;
441:       geom.dimEmbed = dimEmbed;
442:       for (f = 0, v = 0; f < numFields; ++f) {
443:         void * const ctx = ctxs ? ctxs[f] : NULL;

445:         if (!sp[f]) continue;
446:         PetscDualSpaceGetDimension(sp[f], &spDim);
447:         for (d = 0; d < spDim; ++d) {
448:           if (funcs[f]) {
449:             PetscDualSpaceApply(sp[f], d, &geom, numComp[f], funcs[f], ctx, &values[v]);
450:           } else {
451:             for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
452:           }
453:           v += numComp[f];
454:         }
455:       }
456:       DMPlexVecSetClosure(dm, section, localX, p, values, mode);
457:     }
458:     DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
459:     if (h || !maxHeight) {
460:       for (f = 0; f < numFields; f++) {PetscDualSpaceDestroy(&sp[f]);}
461:     }
462:   }
463:   PetscFree2(sp, numComp);
464:   if (maxHeight > 0) {
465:     for (f = 0; f < numFields; f++) {PetscDualSpaceDestroy(&cellsp[f]);}
466:     PetscFree(cellsp);
467:   }
468:   return(0);
469: }

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

476:   Input Parameters:
477: + dm      - The DM
478: . funcs   - The coordinate functions to evaluate, one per field
479: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
480: - mode    - The insertion mode for values

482:   Output Parameter:
483: . X - vector

485:   Level: developer

487: .seealso: DMPlexComputeL2Diff()
488: @*/
489: PetscErrorCode DMPlexProjectFunction(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
490: {
491:   Vec            localX;

496:   DMGetLocalVector(dm, &localX);
497:   DMPlexProjectFunctionLocal(dm, funcs, ctxs, mode, localX);
498:   DMLocalToGlobalBegin(dm, localX, mode, X);
499:   DMLocalToGlobalEnd(dm, localX, mode, X);
500:   DMRestoreLocalVector(dm, &localX);
501:   return(0);
502: }

506: 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)
507: {
508:   DM                dmAux;
509:   PetscDS           prob, probAux;
510:   Vec               A;
511:   PetscSection      section, sectionAux;
512:   PetscScalar      *values, *u, *u_x, *a, *a_x;
513:   PetscReal        *x, *v0, *J, *invJ, detJ, **basisField, **basisFieldDer, **basisFieldAux, **basisFieldDerAux;
514:   PetscInt          Nf, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight;
515:   PetscErrorCode    ierr;

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

547:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
548:     DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);
549:     if (dmAux) {DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
550:     for (f = 0, v = 0; f < Nf; ++f) {
551:       PetscFE        fe;
552:       PetscDualSpace sp;
553:       PetscInt       Ncf;

555:       PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
556:       PetscFEGetDualSpace(fe, &sp);
557:       PetscFEGetNumComponents(fe, &Ncf);
558:       PetscDualSpaceGetDimension(sp, &spDim);
559:       for (d = 0; d < spDim; ++d) {
560:         PetscQuadrature  quad;
561:         const PetscReal *points, *weights;
562:         PetscInt         numPoints, q;

564:         if (funcs[f]) {
565:           PetscDualSpaceGetFunctional(sp, d, &quad);
566:           PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);
567:           PetscFEGetTabulation(fe, numPoints, points, &basisField[f], &basisFieldDer[f], NULL);
568:           for (q = 0; q < numPoints; ++q) {
569:             CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x);
570:             EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, coefficients,    NULL, u, u_x, NULL);
571:             EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);
572:             (*funcs[f])(u, NULL, u_x, a, NULL, a_x, x, &values[v]);
573:           }
574:           PetscFERestoreTabulation(fe, numPoints, points, &basisField[f], &basisFieldDer[f], NULL);
575:         } else {
576:           for (comp = 0; comp < Ncf; ++comp) values[v+comp] = 0.0;
577:         }
578:         v += Ncf;
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:   return(0);
588: }

592: 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)
593: {
594:   void        (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx);
595:   void         **ctxs;
596:   PetscInt       numFields;

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

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

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

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

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

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

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

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

700: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
701: {
702:   PetscInt       numBd, b;

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

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

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

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

751:   Output Parameter:
752: . diff - The diff ||u - u_h||_2

754:   Level: developer

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

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

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

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

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

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

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

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

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

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

864:   Output Parameter:
865: . diff - The diff ||(grad u - grad u_h) . n||_2

867:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

985:   Level: developer

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

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

1014:     DMGetField(dm, field, &obj);
1015:     PetscObjectGetClassId(obj, &id);
1016:     if (id == PETSCFE_CLASSID) {
1017:       PetscFE fe = (PetscFE) obj;

1019:       PetscFEGetQuadrature(fe, &quad);
1020:       PetscFEGetNumComponents(fe, &Nc);
1021:     } else if (id == PETSCFV_CLASSID) {
1022:       PetscFV fv = (PetscFV) obj;

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

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

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

1048:       PetscReal       elemDiff = 0.0;

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

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

1088:   Input Parameters:
1089: + dm - The mesh
1090: . X  - Local input vector
1091: - user - The user context

1093:   Output Parameter:
1094: . integral - Local integral for each field

1096:   Level: developer

1098: .seealso: DMPlexComputeResidualFEM()
1099: @*/
1100: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
1101: {
1102:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1103:   DM                dmAux;
1104:   Vec               localX, A;
1105:   PetscDS           prob, probAux = NULL;
1106:   PetscQuadrature   q;
1107:   PetscSection      section, sectionAux;
1108:   PetscFECellGeom  *cgeom;
1109:   PetscScalar      *u, *a = NULL;
1110:   PetscInt          dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c;
1111:   PetscInt          totDim, totDimAux;
1112:   PetscErrorCode    ierr;

1115:   /*PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);*/
1116:   DMGetLocalVector(dm, &localX);
1117:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1118:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1119:   DMGetDimension(dm, &dim);
1120:   DMGetDefaultSection(dm, &section);
1121:   DMGetDS(dm, &prob);
1122:   PetscDSGetTotalDimension(prob, &totDim);
1123:   PetscSectionGetNumFields(section, &Nf);
1124:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1125:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1126:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1127:   numCells = cEnd - cStart;
1128:   for (f = 0; f < Nf; ++f) {integral[f]    = 0.0;}
1129:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1130:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1131:   if (dmAux) {
1132:     DMGetDefaultSection(dmAux, &sectionAux);
1133:     DMGetDS(dmAux, &probAux);
1134:     PetscDSGetTotalDimension(probAux, &totDimAux);
1135:   }
1136:   DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);
1137:   PetscMalloc2(numCells*totDim,&u,numCells,&cgeom);
1138:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1139:   for (c = cStart; c < cEnd; ++c) {
1140:     PetscScalar *x = NULL;
1141:     PetscInt     i;

1143:     DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);
1144:     if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c);
1145:     DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
1146:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1147:     DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
1148:     if (dmAux) {
1149:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1150:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1151:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1152:     }
1153:   }
1154:   for (f = 0; f < Nf; ++f) {
1155:     PetscFE  fe;
1156:     PetscInt numQuadPoints, Nb;
1157:     /* Conforming batches */
1158:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1159:     /* Remainder */
1160:     PetscInt Nr, offset;

1162:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1163:     PetscFEGetQuadrature(fe, &q);
1164:     PetscFEGetDimension(fe, &Nb);
1165:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1166:     PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1167:     blockSize = Nb*numQuadPoints;
1168:     batchSize = numBlocks * blockSize;
1169:      PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1170:     numChunks = numCells / (numBatches*batchSize);
1171:     Ne        = numChunks*numBatches*batchSize;
1172:     Nr        = numCells % (numBatches*batchSize);
1173:     offset    = numCells - Nr;
1174:     PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, integral);
1175:     PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], integral);
1176:   }
1177:   PetscFree2(u,cgeom);
1178:   if (dmAux) {PetscFree(a);}
1179:   if (mesh->printFEM) {
1180:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1181:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", integral[f]);}
1182:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1183:   }
1184:   DMRestoreLocalVector(dm, &localX);
1185:   /* TODO: Allreduce for integral */
1186:   /*PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);*/
1187:   return(0);
1188: }

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

1195:   Input Parameters:
1196: + dmf  - The fine mesh
1197: . dmc  - The coarse mesh
1198: - user - The user context

1200:   Output Parameter:
1201: . In  - The interpolation matrix

1203:   Note:
1204:   The first member of the user context must be an FEMContext.

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

1209:   Level: developer

1211: .seealso: DMPlexComputeJacobianFEM()
1212: @*/
1213: PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1214: {
1215:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1216:   const char       *name  = "Interpolator";
1217:   PetscDS           prob;
1218:   PetscFE          *feRef;
1219:   PetscSection      fsection, fglobalSection;
1220:   PetscSection      csection, cglobalSection;
1221:   PetscScalar      *elemMat;
1222:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1223:   PetscInt          cTotDim, rTotDim = 0;
1224:   PetscErrorCode    ierr;

1227:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1228:   DMGetDimension(dmf, &dim);
1229:   DMGetDefaultSection(dmf, &fsection);
1230:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1231:   DMGetDefaultSection(dmc, &csection);
1232:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1233:   PetscSectionGetNumFields(fsection, &Nf);
1234:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1235:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1236:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1237:   DMGetDS(dmf, &prob);
1238:   PetscMalloc1(Nf,&feRef);
1239:   for (f = 0; f < Nf; ++f) {
1240:     PetscFE  fe;
1241:     PetscInt rNb, Nc;

1243:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1244:     PetscFERefine(fe, &feRef[f]);
1245:     PetscFEGetDimension(feRef[f], &rNb);
1246:     PetscFEGetNumComponents(fe, &Nc);
1247:     rTotDim += rNb*Nc;
1248:   }
1249:   PetscDSGetTotalDimension(prob, &cTotDim);
1250:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1251:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1252:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1253:     PetscDualSpace   Qref;
1254:     PetscQuadrature  f;
1255:     const PetscReal *qpoints, *qweights;
1256:     PetscReal       *points;
1257:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1259:     /* Compose points from all dual basis functionals */
1260:     PetscFEGetDualSpace(feRef[fieldI], &Qref);
1261:     PetscFEGetNumComponents(feRef[fieldI], &Nc);
1262:     PetscDualSpaceGetDimension(Qref, &fpdim);
1263:     for (i = 0; i < fpdim; ++i) {
1264:       PetscDualSpaceGetFunctional(Qref, i, &f);
1265:       PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1266:       npoints += Np;
1267:     }
1268:     PetscMalloc1(npoints*dim,&points);
1269:     for (i = 0, k = 0; i < fpdim; ++i) {
1270:       PetscDualSpaceGetFunctional(Qref, i, &f);
1271:       PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1272:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1273:     }

1275:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1276:       PetscFE    fe;
1277:       PetscReal *B;
1278:       PetscInt   NcJ, cpdim, j;

1280:       /* Evaluate basis at points */
1281:       PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
1282:       PetscFEGetNumComponents(fe, &NcJ);
1283:       PetscFEGetDimension(fe, &cpdim);
1284:       /* For now, fields only interpolate themselves */
1285:       if (fieldI == fieldJ) {
1286:         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);
1287:         PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1288:         for (i = 0, k = 0; i < fpdim; ++i) {
1289:           PetscDualSpaceGetFunctional(Qref, i, &f);
1290:           PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1291:           for (p = 0; p < Np; ++p, ++k) {
1292:             for (j = 0; j < cpdim; ++j) {
1293:               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];
1294:             }
1295:           }
1296:         }
1297:         PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1298:       }
1299:       offsetJ += cpdim*NcJ;
1300:     }
1301:     offsetI += fpdim*Nc;
1302:     PetscFree(points);
1303:   }
1304:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1305:   /* Preallocate matrix */
1306:   {
1307:     PetscHashJK ht;
1308:     PetscLayout rLayout;
1309:     PetscInt   *dnz, *onz, *cellCIndices, *cellFIndices;
1310:     PetscInt    locRows, rStart, rEnd, cell, r;

1312:     MatGetLocalSize(In, &locRows, NULL);
1313:     PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1314:     PetscLayoutSetLocalSize(rLayout, locRows);
1315:     PetscLayoutSetBlockSize(rLayout, 1);
1316:     PetscLayoutSetUp(rLayout);
1317:     PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1318:     PetscLayoutDestroy(&rLayout);
1319:     PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1320:     PetscHashJKCreate(&ht);
1321:     for (cell = cStart; cell < cEnd; ++cell) {
1322:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1323:       for (r = 0; r < rTotDim; ++r) {
1324:         PetscHashJKKey  key;
1325:         PetscHashJKIter missing, iter;

1327:         key.j = cellFIndices[r];
1328:         if (key.j < 0) continue;
1329:         for (c = 0; c < cTotDim; ++c) {
1330:           key.k = cellCIndices[c];
1331:           if (key.k < 0) continue;
1332:           PetscHashJKPut(ht, key, &missing, &iter);
1333:           if (missing) {
1334:             PetscHashJKSet(ht, iter, 1);
1335:             if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1336:             else                                     ++onz[key.j-rStart];
1337:           }
1338:         }
1339:       }
1340:     }
1341:     PetscHashJKDestroy(&ht);
1342:     MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1343:     MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1344:     PetscFree4(dnz,onz,cellCIndices,cellFIndices);
1345:   }
1346:   /* Fill matrix */
1347:   MatZeroEntries(In);
1348:   for (c = cStart; c < cEnd; ++c) {
1349:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1350:   }
1351:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1352:   PetscFree(feRef);
1353:   PetscFree(elemMat);
1354:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1355:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1356:   if (mesh->printFEM) {
1357:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1358:     MatChop(In, 1.0e-10);
1359:     MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1360:   }
1361:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1362:   return(0);
1363: }

1367: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1368: {
1369:   PetscDS        prob;
1370:   PetscFE       *feRef;
1371:   Vec            fv, cv;
1372:   IS             fis, cis;
1373:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1374:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1375:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m;

1379:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1380:   DMGetDimension(dmf, &dim);
1381:   DMGetDefaultSection(dmf, &fsection);
1382:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1383:   DMGetDefaultSection(dmc, &csection);
1384:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1385:   PetscSectionGetNumFields(fsection, &Nf);
1386:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1387:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1388:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1389:   DMGetDS(dmc, &prob);
1390:   PetscMalloc1(Nf,&feRef);
1391:   for (f = 0; f < Nf; ++f) {
1392:     PetscFE  fe;
1393:     PetscInt fNb, Nc;

1395:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1396:     PetscFERefine(fe, &feRef[f]);
1397:     PetscFEGetDimension(feRef[f], &fNb);
1398:     PetscFEGetNumComponents(fe, &Nc);
1399:     fTotDim += fNb*Nc;
1400:   }
1401:   PetscDSGetTotalDimension(prob, &cTotDim);
1402:   PetscMalloc1(cTotDim,&cmap);
1403:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1404:     PetscFE        feC;
1405:     PetscDualSpace QF, QC;
1406:     PetscInt       NcF, NcC, fpdim, cpdim;

1408:     PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1409:     PetscFEGetNumComponents(feC, &NcC);
1410:     PetscFEGetNumComponents(feRef[field], &NcF);
1411:     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);
1412:     PetscFEGetDualSpace(feRef[field], &QF);
1413:     PetscDualSpaceGetDimension(QF, &fpdim);
1414:     PetscFEGetDualSpace(feC, &QC);
1415:     PetscDualSpaceGetDimension(QC, &cpdim);
1416:     for (c = 0; c < cpdim; ++c) {
1417:       PetscQuadrature  cfunc;
1418:       const PetscReal *cqpoints;
1419:       PetscInt         NpC;

1421:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
1422:       PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);
1423:       if (NpC != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1424:       for (f = 0; f < fpdim; ++f) {
1425:         PetscQuadrature  ffunc;
1426:         const PetscReal *fqpoints;
1427:         PetscReal        sum = 0.0;
1428:         PetscInt         NpF, comp;

1430:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
1431:         PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);
1432:         if (NpC != NpF) continue;
1433:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1434:         if (sum > 1.0e-9) continue;
1435:         for (comp = 0; comp < NcC; ++comp) {
1436:           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1437:         }
1438:         break;
1439:       }
1440:     }
1441:     offsetC += cpdim*NcC;
1442:     offsetF += fpdim*NcF;
1443:   }
1444:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1445:   PetscFree(feRef);

1447:   DMGetGlobalVector(dmf, &fv);
1448:   DMGetGlobalVector(dmc, &cv);
1449:   VecGetOwnershipRange(cv, &startC, NULL);
1450:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1451:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1452:   PetscMalloc1(m,&cindices);
1453:   PetscMalloc1(m,&findices);
1454:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1455:   for (c = cStart; c < cEnd; ++c) {
1456:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1457:     for (d = 0; d < cTotDim; ++d) {
1458:       if (cellCIndices[d] < 0) continue;
1459:       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]]);
1460:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1461:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1462:     }
1463:   }
1464:   PetscFree(cmap);
1465:   PetscFree2(cellCIndices,cellFIndices);

1467:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1468:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1469:   VecScatterCreate(cv, cis, fv, fis, sc);
1470:   ISDestroy(&cis);
1471:   ISDestroy(&fis);
1472:   DMRestoreGlobalVector(dmf, &fv);
1473:   DMRestoreGlobalVector(dmc, &cv);
1474:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1475:   return(0);
1476: }