Actual source code: plexfem.c

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

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

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

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

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

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

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

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

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

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

116:   Collective on DM

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

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

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

126:   Level: advanced

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

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

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

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

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

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

192:   Level: advanced

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

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

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

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

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

218:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

462:   Output Parameter:
463: . X - vector

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

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

474:   Level: developer

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

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

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

509:   DMPlexGetMaxProjectionHeight(dm,&maxHeight);
510:   if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");}
511:   DMGetDS(dm, &prob);
512:   DMGetDimension(dm, &dim);
513:   DMGetDefaultSection(dm, &section);
514:   PetscSectionGetNumFields(section, &Nf);
515:   PetscMalloc2(Nf, &sp, Nf, &Ncf);
516:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
517:   PetscDSGetTotalDimension(prob, &totDim);
518:   PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
519:   PetscDSGetRefCoordArrays(prob, &x, NULL);
520:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
521:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
522:   if (dmAux) {
523:     DMGetDS(dmAux, &probAux);
524:     DMGetDefaultSection(dmAux, &sectionAux);
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:       PetscObject  obj;
542:       PetscClassId id;

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

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

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

565:         if (funcs[f]) {
566:           PetscDualSpaceGetFunctional(sp[f], d, &quad);
567:           PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);
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:         } else {
575:           for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0;
576:         }
577:         v += Ncf[f];
578:       }
579:     }
580:     DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);
581:     if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
582:     DMPlexVecSetClosure(dm, section, localX, c, values, mode);
583:   }
584:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
585:   PetscFree3(v0,J,invJ);
586:   for (f = 0; f < Nf; f++) {PetscDualSpaceDestroy(&sp[f]);}
587:   PetscFree2(sp, Ncf);
588:   return(0);
589: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

756:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

869:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

987:   Level: developer

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

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

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

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

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

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

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

1050:       PetscReal       elemDiff = 0.0;

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

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

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

1095:   Output Parameter:
1096: . integral - Local integral for each field

1098:   Level: developer

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

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

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

1162:     PetscDSGetDiscretization(prob, f, &obj);
1163:     PetscObjectGetClassId(obj, &id);
1164:     if (id == PETSCFE_CLASSID) {
1165:       PetscFE         fe = (PetscFE) obj;
1166:       PetscQuadrature q;
1167:       PetscInt        Nq, Nb;

1169:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1170:       PetscFEGetQuadrature(fe, &q);
1171:       PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1172:       PetscFEGetDimension(fe, &Nb);
1173:       blockSize = Nb*Nq;
1174:       batchSize = numBlocks * blockSize;
1175:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1176:       numChunks = numCells / (numBatches*batchSize);
1177:       Ne        = numChunks*numBatches*batchSize;
1178:       Nr        = numCells % (numBatches*batchSize);
1179:       offset    = numCells - Nr;
1180:       PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);
1181:       PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);
1182:     } else if (id == PETSCFV_CLASSID) {
1183:       /* PetscFV  fv = (PetscFV) obj; */
1184:       PetscInt       foff;
1185:       PetscPointFunc obj_func;
1186:       PetscScalar    lint;

1188:       PetscDSGetObjective(prob, f, &obj_func);
1189:       PetscDSGetFieldOffset(prob, f, &foff);
1190:       if (obj_func) {
1191:         for (c = 0; c < numCells; ++c) {
1192:           /* TODO: Need full pointwise interpolation and get centroid for x argument */
1193:           obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint);
1194:           lintegral[f] = PetscRealPart(lint)*vol[c];
1195:         }
1196:       }
1197:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1198:   }
1199:   if (dmAux) {PetscFree(a);}
1200:   if (mesh->printFEM) {
1201:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
1202:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);}
1203:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1204:   }
1205:   DMRestoreLocalVector(dm, &localX);
1206:   MPI_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1207:   PetscFree4(lintegral,u,cgeom,vol);
1208:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1209:   return(0);
1210: }

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

1217:   Input Parameters:
1218: + dmf  - The fine mesh
1219: . dmc  - The coarse mesh
1220: - user - The user context

1222:   Output Parameter:
1223: . In  - The interpolation matrix

1225:   Note:
1226:   The first member of the user context must be an FEMContext.

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

1231:   Level: developer

1233: .seealso: DMPlexComputeJacobianFEM()
1234: @*/
1235: PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1236: {
1237:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1238:   const char       *name  = "Interpolator";
1239:   PetscDS           prob;
1240:   PetscFE          *feRef;
1241:   PetscFV          *fvRef;
1242:   PetscSection      fsection, fglobalSection;
1243:   PetscSection      csection, cglobalSection;
1244:   PetscScalar      *elemMat;
1245:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1246:   PetscInt          cTotDim, rTotDim = 0;
1247:   PetscErrorCode    ierr;

1250:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1251:   DMGetDimension(dmf, &dim);
1252:   DMGetDefaultSection(dmf, &fsection);
1253:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1254:   DMGetDefaultSection(dmc, &csection);
1255:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1256:   PetscSectionGetNumFields(fsection, &Nf);
1257:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1258:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1259:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1260:   DMGetDS(dmf, &prob);
1261:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1262:   for (f = 0; f < Nf; ++f) {
1263:     PetscObject  obj;
1264:     PetscClassId id;
1265:     PetscInt     rNb = 0, Nc = 0;

1267:     PetscDSGetDiscretization(prob, f, &obj);
1268:     PetscObjectGetClassId(obj, &id);
1269:     if (id == PETSCFE_CLASSID) {
1270:       PetscFE fe = (PetscFE) obj;

1272:       PetscFERefine(fe, &feRef[f]);
1273:       PetscFEGetDimension(feRef[f], &rNb);
1274:       PetscFEGetNumComponents(fe, &Nc);
1275:     } else if (id == PETSCFV_CLASSID) {
1276:       PetscFV        fv = (PetscFV) obj;
1277:       PetscDualSpace Q;

1279:       PetscFVRefine(fv, &fvRef[f]);
1280:       PetscFVGetDualSpace(fvRef[f], &Q);
1281:       PetscDualSpaceGetDimension(Q, &rNb);
1282:       PetscFVGetNumComponents(fv, &Nc);
1283:     }
1284:     rTotDim += rNb*Nc;
1285:   }
1286:   PetscDSGetTotalDimension(prob, &cTotDim);
1287:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1288:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1289:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1290:     PetscDualSpace   Qref;
1291:     PetscQuadrature  f;
1292:     const PetscReal *qpoints, *qweights;
1293:     PetscReal       *points;
1294:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1296:     /* Compose points from all dual basis functionals */
1297:     if (feRef[fieldI]) {
1298:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
1299:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
1300:     } else {
1301:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
1302:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
1303:     }
1304:     PetscDualSpaceGetDimension(Qref, &fpdim);
1305:     for (i = 0; i < fpdim; ++i) {
1306:       PetscDualSpaceGetFunctional(Qref, i, &f);
1307:       PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1308:       npoints += Np;
1309:     }
1310:     PetscMalloc1(npoints*dim,&points);
1311:     for (i = 0, k = 0; i < fpdim; ++i) {
1312:       PetscDualSpaceGetFunctional(Qref, i, &f);
1313:       PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1314:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1315:     }

1317:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1318:       PetscObject  obj;
1319:       PetscClassId id;
1320:       PetscReal   *B;
1321:       PetscInt     NcJ = 0, cpdim = 0, j;

1323:       PetscDSGetDiscretization(prob, fieldJ, &obj);
1324:       PetscObjectGetClassId(obj, &id);
1325:       if (id == PETSCFE_CLASSID) {
1326:         PetscFE fe = (PetscFE) obj;

1328:         /* Evaluate basis at points */
1329:         PetscFEGetNumComponents(fe, &NcJ);
1330:         PetscFEGetDimension(fe, &cpdim);
1331:         /* For now, fields only interpolate themselves */
1332:         if (fieldI == fieldJ) {
1333:           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);
1334:           PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1335:           for (i = 0, k = 0; i < fpdim; ++i) {
1336:             PetscDualSpaceGetFunctional(Qref, i, &f);
1337:             PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1338:             for (p = 0; p < Np; ++p, ++k) {
1339:               for (j = 0; j < cpdim; ++j) {
1340:                 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];
1341:               }
1342:             }
1343:           }
1344:           PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1345:         }
1346:       } else if (id == PETSCFV_CLASSID) {
1347:         PetscFV        fv = (PetscFV) obj;

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

1379:     MatGetLocalSize(In, &locRows, NULL);
1380:     PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1381:     PetscLayoutSetLocalSize(rLayout, locRows);
1382:     PetscLayoutSetBlockSize(rLayout, 1);
1383:     PetscLayoutSetUp(rLayout);
1384:     PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1385:     PetscLayoutDestroy(&rLayout);
1386:     PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1387:     PetscHashJKCreate(&ht);
1388:     for (cell = cStart; cell < cEnd; ++cell) {
1389:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1390:       for (r = 0; r < rTotDim; ++r) {
1391:         PetscHashJKKey  key;
1392:         PetscHashJKIter missing, iter;

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

1434: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1435: {
1436:   PetscDS        prob;
1437:   PetscFE       *feRef;
1438:   PetscFV       *fvRef;
1439:   Vec            fv, cv;
1440:   IS             fis, cis;
1441:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1442:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1443:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m;

1447:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1448:   DMGetDimension(dmf, &dim);
1449:   DMGetDefaultSection(dmf, &fsection);
1450:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1451:   DMGetDefaultSection(dmc, &csection);
1452:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1453:   PetscSectionGetNumFields(fsection, &Nf);
1454:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1455:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
1456:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1457:   DMGetDS(dmc, &prob);
1458:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
1459:   for (f = 0; f < Nf; ++f) {
1460:     PetscObject  obj;
1461:     PetscClassId id;
1462:     PetscInt     fNb = 0, Nc = 0;

1464:     PetscDSGetDiscretization(prob, f, &obj);
1465:     PetscObjectGetClassId(obj, &id);
1466:     if (id == PETSCFE_CLASSID) {
1467:       PetscFE fe = (PetscFE) obj;

1469:       PetscFERefine(fe, &feRef[f]);
1470:       PetscFEGetDimension(feRef[f], &fNb);
1471:       PetscFEGetNumComponents(fe, &Nc);
1472:     } else if (id == PETSCFV_CLASSID) {
1473:       PetscFV        fv = (PetscFV) obj;
1474:       PetscDualSpace Q;

1476:       PetscFVRefine(fv, &fvRef[f]);
1477:       PetscFVGetDualSpace(fvRef[f], &Q);
1478:       PetscDualSpaceGetDimension(Q, &fNb);
1479:       PetscFVGetNumComponents(fv, &Nc);
1480:     }
1481:     fTotDim += fNb*Nc;
1482:   }
1483:   PetscDSGetTotalDimension(prob, &cTotDim);
1484:   PetscMalloc1(cTotDim,&cmap);
1485:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1486:     PetscFE        feC;
1487:     PetscFV        fvC;
1488:     PetscDualSpace QF, QC;
1489:     PetscInt       NcF, NcC, fpdim, cpdim;

1491:     if (feRef[field]) {
1492:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1493:       PetscFEGetNumComponents(feC, &NcC);
1494:       PetscFEGetNumComponents(feRef[field], &NcF);
1495:       PetscFEGetDualSpace(feRef[field], &QF);
1496:       PetscDualSpaceGetDimension(QF, &fpdim);
1497:       PetscFEGetDualSpace(feC, &QC);
1498:       PetscDualSpaceGetDimension(QC, &cpdim);
1499:     } else {
1500:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
1501:       PetscFVGetNumComponents(fvC, &NcC);
1502:       PetscFVGetNumComponents(fvRef[field], &NcF);
1503:       PetscFVGetDualSpace(fvRef[field], &QF);
1504:       PetscDualSpaceGetDimension(QF, &fpdim);
1505:       PetscFVGetDualSpace(fvC, &QC);
1506:       PetscDualSpaceGetDimension(QC, &cpdim);
1507:     }
1508:     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);
1509:     for (c = 0; c < cpdim; ++c) {
1510:       PetscQuadrature  cfunc;
1511:       const PetscReal *cqpoints;
1512:       PetscInt         NpC;
1513:       PetscBool        found = PETSC_FALSE;

1515:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
1516:       PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);
1517:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1518:       for (f = 0; f < fpdim; ++f) {
1519:         PetscQuadrature  ffunc;
1520:         const PetscReal *fqpoints;
1521:         PetscReal        sum = 0.0;
1522:         PetscInt         NpF, comp;

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

1551:   DMGetGlobalVector(dmf, &fv);
1552:   DMGetGlobalVector(dmc, &cv);
1553:   VecGetOwnershipRange(cv, &startC, NULL);
1554:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1555:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1556:   PetscMalloc1(m,&cindices);
1557:   PetscMalloc1(m,&findices);
1558:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1559:   for (c = cStart; c < cEnd; ++c) {
1560:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1561:     for (d = 0; d < cTotDim; ++d) {
1562:       if (cellCIndices[d] < 0) continue;
1563:       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]]);
1564:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1565:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1566:     }
1567:   }
1568:   PetscFree(cmap);
1569:   PetscFree2(cellCIndices,cellFIndices);

1571:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1572:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1573:   VecScatterCreate(cv, cis, fv, fis, sc);
1574:   ISDestroy(&cis);
1575:   ISDestroy(&fis);
1576:   DMRestoreGlobalVector(dmf, &fv);
1577:   DMRestoreGlobalVector(dmc, &cv);
1578:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1579:   return(0);
1580: }