Actual source code: plexfem.c

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

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

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

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

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

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

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

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

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

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

114:   Collective on DM

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

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

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

124:   Level: advanced

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

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

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

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

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

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

190:   Level: advanced

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

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

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

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

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

216:   Level: intermediate

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

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

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

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

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

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

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

305:       DMLabelGetStratumIS(label, ids[i], &pointIS);
306:       ISGetLocalSize(pointIS, &n);
307:       ISGetIndices(pointIS, &points);
308:       for (p = 0; p < n; ++p) {
309:         const PetscInt    point = points[p];
310:         PetscFECellGeom   geom;

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

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

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

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

389:       DMGetField(dm, f, &obj);
390:       PetscObjectGetClassId(obj, &id);
391:       if (id == PETSCFE_CLASSID) {
392:         PetscFE fe = (PetscFE) obj;

394:         PetscFEGetNumComponents(fe, &numComp[f]);
395:         if (!h) {
396:           PetscFEGetDualSpace(fe, &cellsp[f]);
397:           sp[f] = cellsp[f];
398:           PetscObjectReference((PetscObject) sp[f]);
399:         }
400:         else {
401:           PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);
402:           if (!sp[f]) {
403:             continue;
404:           }
405:         }
406:       } else if (id == PETSCFV_CLASSID) {
407:         PetscFV         fv = (PetscFV) obj;
408:         PetscQuadrature q;

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

429:       DMPlexComputeCellGeometryFEM(dm, p, NULL, geom.v0, geom.J, NULL, &geom.detJ);
430:       geom.dim      = dim - h;
431:       geom.dimEmbed = dimEmbed;
432:       for (f = 0, v = 0; f < numFields; ++f) {
433:         void * const ctx = ctxs ? ctxs[f] : NULL;

435:         if (!sp[f]) continue;
436:         PetscDualSpaceGetDimension(sp[f], &spDim);
437:         for (d = 0; d < spDim; ++d) {
438:           if (funcs[f]) {
439:             PetscDualSpaceApply(sp[f], d, &geom, numComp[f], funcs[f], ctx, &values[v]);
440:           } else {
441:             for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
442:           }
443:           v += numComp[f];
444:         }
445:       }
446:       DMPlexVecSetClosure(dm, section, localX, p, values, mode);
447:     }
448:     DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
449:     if (h || !maxHeight) {
450:       for (f = 0; f < numFields; f++) {PetscDualSpaceDestroy(&sp[f]);}
451:     }
452:   }
453:   PetscFree2(sp, numComp);
454:   if (maxHeight > 0) {
455:     for (f = 0; f < numFields; f++) {PetscDualSpaceDestroy(&cellsp[f]);}
456:     PetscFree(cellsp);
457:   }
458:   return(0);
459: }

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

466:   Input Parameters:
467: + dm      - The DM
468: . funcs   - The coordinate functions to evaluate, one per field
469: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
470: - mode    - The insertion mode for values

472:   Output Parameter:
473: . X - vector

475:   Level: developer

477: .seealso: DMPlexComputeL2Diff()
478: @*/
479: PetscErrorCode DMPlexProjectFunction(DM dm, void (**funcs)(const PetscReal [], 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;
500:   Vec               A;
501:   PetscSection      section, sectionAux;
502:   PetscScalar      *values, *u, *u_x, *a, *a_x;
503:   PetscReal        *x, *v0, *J, *invJ, detJ, **basisField, **basisFieldDer, **basisFieldAux, **basisFieldDerAux;
504:   PetscInt          Nf, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight;
505:   PetscErrorCode    ierr;

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

537:     DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
538:     DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);
539:     if (dmAux) {DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
540:     for (f = 0, v = 0; f < Nf; ++f) {
541:       PetscFE        fe;
542:       PetscDualSpace sp;
543:       PetscInt       Ncf;

545:       PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
546:       PetscFEGetDualSpace(fe, &sp);
547:       PetscFEGetNumComponents(fe, &Ncf);
548:       PetscDualSpaceGetDimension(sp, &spDim);
549:       for (d = 0; d < spDim; ++d) {
550:         PetscQuadrature  quad;
551:         const PetscReal *points, *weights;
552:         PetscInt         numPoints, q;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

744:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

857:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

975:   Level: developer

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

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

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

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

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

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

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

1038:       PetscReal       elemDiff = 0.0;

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

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

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

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

1086:   Level: developer

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

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

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

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

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

1185:   Input Parameters:
1186: + dmf  - The fine mesh
1187: . dmc  - The coarse mesh
1188: - user - The user context

1190:   Output Parameter:
1191: . In  - The interpolation matrix

1193:   Note:
1194:   The first member of the user context must be an FEMContext.

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

1199:   Level: developer

1201: .seealso: DMPlexComputeJacobianFEM()
1202: @*/
1203: PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1204: {
1205:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1206:   const char       *name  = "Interpolator";
1207:   PetscDS           prob;
1208:   PetscFE          *feRef;
1209:   PetscSection      fsection, fglobalSection;
1210:   PetscSection      csection, cglobalSection;
1211:   PetscScalar      *elemMat;
1212:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1213:   PetscInt          cTotDim, rTotDim = 0;
1214:   PetscErrorCode    ierr;

1217:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1218:   DMGetDimension(dmf, &dim);
1219:   DMGetDefaultSection(dmf, &fsection);
1220:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1221:   DMGetDefaultSection(dmc, &csection);
1222:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1223:   PetscSectionGetNumFields(fsection, &Nf);
1224:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1225:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1226:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1227:   DMGetDS(dmf, &prob);
1228:   PetscMalloc1(Nf,&feRef);
1229:   for (f = 0; f < Nf; ++f) {
1230:     PetscFE  fe;
1231:     PetscInt rNb, Nc;

1233:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1234:     PetscFERefine(fe, &feRef[f]);
1235:     PetscFEGetDimension(feRef[f], &rNb);
1236:     PetscFEGetNumComponents(fe, &Nc);
1237:     rTotDim += rNb*Nc;
1238:   }
1239:   PetscDSGetTotalDimension(prob, &cTotDim);
1240:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1241:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1242:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1243:     PetscDualSpace   Qref;
1244:     PetscQuadrature  f;
1245:     const PetscReal *qpoints, *qweights;
1246:     PetscReal       *points;
1247:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1249:     /* Compose points from all dual basis functionals */
1250:     PetscFEGetDualSpace(feRef[fieldI], &Qref);
1251:     PetscFEGetNumComponents(feRef[fieldI], &Nc);
1252:     PetscDualSpaceGetDimension(Qref, &fpdim);
1253:     for (i = 0; i < fpdim; ++i) {
1254:       PetscDualSpaceGetFunctional(Qref, i, &f);
1255:       PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1256:       npoints += Np;
1257:     }
1258:     PetscMalloc1(npoints*dim,&points);
1259:     for (i = 0, k = 0; i < fpdim; ++i) {
1260:       PetscDualSpaceGetFunctional(Qref, i, &f);
1261:       PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1262:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1263:     }

1265:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1266:       PetscFE    fe;
1267:       PetscReal *B;
1268:       PetscInt   NcJ, cpdim, j;

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

1302:     MatGetLocalSize(In, &locRows, NULL);
1303:     PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1304:     PetscLayoutSetLocalSize(rLayout, locRows);
1305:     PetscLayoutSetBlockSize(rLayout, 1);
1306:     PetscLayoutSetUp(rLayout);
1307:     PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1308:     PetscLayoutDestroy(&rLayout);
1309:     PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1310:     PetscHashJKCreate(&ht);
1311:     for (cell = cStart; cell < cEnd; ++cell) {
1312:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1313:       for (r = 0; r < rTotDim; ++r) {
1314:         PetscHashJKKey  key;
1315:         PetscHashJKIter missing, iter;

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

1357: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1358: {
1359:   PetscDS        prob;
1360:   PetscFE       *feRef;
1361:   Vec            fv, cv;
1362:   IS             fis, cis;
1363:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1364:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1365:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m;

1369:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1370:   DMGetDimension(dmf, &dim);
1371:   DMGetDefaultSection(dmf, &fsection);
1372:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1373:   DMGetDefaultSection(dmc, &csection);
1374:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1375:   PetscSectionGetNumFields(fsection, &Nf);
1376:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1377:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1378:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1379:   DMGetDS(dmc, &prob);
1380:   PetscMalloc1(Nf,&feRef);
1381:   for (f = 0; f < Nf; ++f) {
1382:     PetscFE  fe;
1383:     PetscInt fNb, Nc;

1385:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1386:     PetscFERefine(fe, &feRef[f]);
1387:     PetscFEGetDimension(feRef[f], &fNb);
1388:     PetscFEGetNumComponents(fe, &Nc);
1389:     fTotDim += fNb*Nc;
1390:   }
1391:   PetscDSGetTotalDimension(prob, &cTotDim);
1392:   PetscMalloc1(cTotDim,&cmap);
1393:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1394:     PetscFE        feC;
1395:     PetscDualSpace QF, QC;
1396:     PetscInt       NcF, NcC, fpdim, cpdim;

1398:     PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1399:     PetscFEGetNumComponents(feC, &NcC);
1400:     PetscFEGetNumComponents(feRef[field], &NcF);
1401:     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);
1402:     PetscFEGetDualSpace(feRef[field], &QF);
1403:     PetscDualSpaceGetDimension(QF, &fpdim);
1404:     PetscFEGetDualSpace(feC, &QC);
1405:     PetscDualSpaceGetDimension(QC, &cpdim);
1406:     for (c = 0; c < cpdim; ++c) {
1407:       PetscQuadrature  cfunc;
1408:       const PetscReal *cqpoints;
1409:       PetscInt         NpC;

1411:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
1412:       PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);
1413:       if (NpC != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1414:       for (f = 0; f < fpdim; ++f) {
1415:         PetscQuadrature  ffunc;
1416:         const PetscReal *fqpoints;
1417:         PetscReal        sum = 0.0;
1418:         PetscInt         NpF, comp;

1420:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
1421:         PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);
1422:         if (NpC != NpF) continue;
1423:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1424:         if (sum > 1.0e-9) continue;
1425:         for (comp = 0; comp < NcC; ++comp) {
1426:           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1427:         }
1428:         break;
1429:       }
1430:     }
1431:     offsetC += cpdim*NcC;
1432:     offsetF += fpdim*NcF;
1433:   }
1434:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1435:   PetscFree(feRef);

1437:   DMGetGlobalVector(dmf, &fv);
1438:   DMGetGlobalVector(dmc, &cv);
1439:   VecGetOwnershipRange(cv, &startC, NULL);
1440:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1441:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1442:   PetscMalloc1(m,&cindices);
1443:   PetscMalloc1(m,&findices);
1444:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1445:   for (c = cStart; c < cEnd; ++c) {
1446:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1447:     for (d = 0; d < cTotDim; ++d) {
1448:       if (cellCIndices[d] < 0) continue;
1449:       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]]);
1450:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1451:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1452:     }
1453:   }
1454:   PetscFree(cmap);
1455:   PetscFree2(cellCIndices,cellFIndices);

1457:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1458:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1459:   VecScatterCreate(cv, cis, fv, fis, sc);
1460:   ISDestroy(&cis);
1461:   ISDestroy(&fis);
1462:   DMRestoreGlobalVector(dmf, &fv);
1463:   DMRestoreGlobalVector(dmc, &cv);
1464:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1465:   return(0);
1466: }