Actual source code: plexfem.c

petsc-master 2020-06-03
Report Typos and Errors
  1:  #include <petsc/private/dmpleximpl.h>
  2:  #include <petscsf.h>

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

  8: static PetscErrorCode DMPlexConvertPlex(DM dm, DM *plex, PetscBool copy)
  9: {
 10:   PetscBool      isPlex;

 14:   PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
 15:   if (isPlex) {
 16:     *plex = dm;
 17:     PetscObjectReference((PetscObject) dm);
 18:   } else {
 19:     PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
 20:     if (!*plex) {
 21:       DMConvert(dm,DMPLEX,plex);
 22:       PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
 23:       if (copy) {
 24:         const char *comps[3] = {"A", "dmAux"};
 25:         PetscObject obj;
 26:         PetscInt    i;

 28:         {
 29:           /* Run the subdomain hook (this will copy the DMSNES/DMTS) */
 30:           DMSubDomainHookLink link;
 31:           for (link = dm->subdomainhook; link; link = link->next) {
 32:             if (link->ddhook) {(*link->ddhook)(dm, *plex, link->ctx);}
 33:           }
 34:         }
 35:         for (i = 0; i < 3; i++) {
 36:           PetscObjectQuery((PetscObject) dm, comps[i], &obj);
 37:           PetscObjectCompose((PetscObject) *plex, comps[i], obj);
 38:         }
 39:       }
 40:     } else {
 41:       PetscObjectReference((PetscObject) *plex);
 42:     }
 43:   }
 44:   return(0);
 45: }

 47: static PetscErrorCode PetscContainerUserDestroy_PetscFEGeom (void *ctx)
 48: {
 49:   PetscFEGeom *geom = (PetscFEGeom *) ctx;

 53:   PetscFEGeomDestroy(&geom);
 54:   return(0);
 55: }

 57: static PetscErrorCode DMPlexGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
 58: {
 59:   char            composeStr[33] = {0};
 60:   PetscObjectId   id;
 61:   PetscContainer  container;
 62:   PetscErrorCode  ierr;

 65:   PetscObjectGetId((PetscObject)quad,&id);
 66:   PetscSNPrintf(composeStr, 32, "DMPlexGetFEGeom_%x\n", id);
 67:   PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
 68:   if (container) {
 69:     PetscContainerGetPointer(container, (void **) geom);
 70:   } else {
 71:     DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
 72:     PetscContainerCreate(PETSC_COMM_SELF,&container);
 73:     PetscContainerSetPointer(container, (void *) *geom);
 74:     PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
 75:     PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
 76:     PetscContainerDestroy(&container);
 77:   }
 78:   return(0);
 79: }

 81: static PetscErrorCode DMPlexRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
 82: {
 84:   *geom = NULL;
 85:   return(0);
 86: }

 88: /*@
 89:   DMPlexGetScale - Get the scale for the specified fundamental unit

 91:   Not collective

 93:   Input Arguments:
 94: + dm   - the DM
 95: - unit - The SI unit

 97:   Output Argument:
 98: . scale - The value used to scale all quantities with this unit

100:   Level: advanced

102: .seealso: DMPlexSetScale(), PetscUnit
103: @*/
104: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
105: {
106:   DM_Plex *mesh = (DM_Plex*) dm->data;

111:   *scale = mesh->scale[unit];
112:   return(0);
113: }

115: /*@
116:   DMPlexSetScale - Set the scale for the specified fundamental unit

118:   Not collective

120:   Input Arguments:
121: + dm   - the DM
122: . unit - The SI unit
123: - scale - The value used to scale all quantities with this unit

125:   Level: advanced

127: .seealso: DMPlexGetScale(), PetscUnit
128: @*/
129: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
130: {
131:   DM_Plex *mesh = (DM_Plex*) dm->data;

135:   mesh->scale[unit] = scale;
136:   return(0);
137: }

139: static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nc, PetscScalar *mode, void *ctx)
140: {
141:   const PetscInt eps[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
142:   PetscInt *ctxInt  = (PetscInt *) ctx;
143:   PetscInt  dim2    = ctxInt[0];
144:   PetscInt  d       = ctxInt[1];
145:   PetscInt  i, j, k = dim > 2 ? d - dim : d;

148:   if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %D does not match context dimension %D", dim, dim2);
149:   for (i = 0; i < dim; i++) mode[i] = 0.;
150:   if (d < dim) {
151:     mode[d] = 1.; /* Translation along axis d */
152:   } else {
153:     for (i = 0; i < dim; i++) {
154:       for (j = 0; j < dim; j++) {
155:         mode[j] += eps[i][j][k]*X[i]; /* Rotation about axis d */
156:       }
157:     }
158:   }
159:   return(0);
160: }

162: /*@
163:   DMPlexCreateRigidBody - For the default global section, create rigid body modes by function space interpolation

165:   Collective on dm

167:   Input Arguments:
168: . dm - the DM

170:   Output Argument:
171: . sp - the null space

173:   Note: This is necessary to provide a suitable coarse space for algebraic multigrid

175:   Level: advanced

177: .seealso: MatNullSpaceCreate(), PCGAMG
178: @*/
179: PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
180: {
181:   MPI_Comm       comm;
182:   Vec            mode[6];
183:   PetscSection   section, globalSection;
184:   PetscInt       dim, dimEmbed, n, m, mmin, d, i, j;

188:   PetscObjectGetComm((PetscObject)dm,&comm);
189:   DMGetDimension(dm, &dim);
190:   DMGetCoordinateDim(dm, &dimEmbed);
191:   if (dim == 1) {
192:     MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
193:     return(0);
194:   }
195:   DMGetLocalSection(dm, &section);
196:   DMGetGlobalSection(dm, &globalSection);
197:   PetscSectionGetConstrainedStorageSize(globalSection, &n);
198:   m    = (dim*(dim+1))/2;
199:   VecCreate(comm, &mode[0]);
200:   VecSetSizes(mode[0], n, PETSC_DETERMINE);
201:   VecSetUp(mode[0]);
202:   VecGetSize(mode[0], &n);
203:   mmin = PetscMin(m, n);
204:   for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
205:   for (d = 0; d < m; d++) {
206:     PetscInt         ctx[2];
207:     PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
208:     void            *voidctx = (void *) (&ctx[0]);

210:     ctx[0] = dimEmbed;
211:     ctx[1] = d;
212:     DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);
213:   }
214:   /* Orthonormalize system */
215:   for (i = 0; i < mmin; ++i) {
216:     PetscScalar dots[6];

218:     VecNormalize(mode[i], NULL);
219:     VecMDot(mode[i], mmin-i-1, mode+i+1, dots+i+1);
220:     for (j = i+1; j < mmin; ++j) {
221:       dots[j] *= -1.0;
222:       VecAXPY(mode[j], dots[j], mode[i]);
223:     }
224:   }
225:   MatNullSpaceCreate(comm, PETSC_FALSE, mmin, mode, sp);
226:   for (i = 0; i < m; ++i) {VecDestroy(&mode[i]);}
227:   return(0);
228: }

230: /*@
231:   DMPlexCreateRigidBodies - For the default global section, create rigid body modes by function space interpolation

233:   Collective on dm

235:   Input Arguments:
236: + dm    - the DM
237: . nb    - The number of bodies
238: . label - The DMLabel marking each domain
239: . nids  - The number of ids per body
240: - ids   - An array of the label ids in sequence for each domain

242:   Output Argument:
243: . sp - the null space

245:   Note: This is necessary to provide a suitable coarse space for algebraic multigrid

247:   Level: advanced

249: .seealso: MatNullSpaceCreate()
250: @*/
251: PetscErrorCode DMPlexCreateRigidBodies(DM dm, PetscInt nb, DMLabel label, const PetscInt nids[], const PetscInt ids[], MatNullSpace *sp)
252: {
253:   MPI_Comm       comm;
254:   PetscSection   section, globalSection;
255:   Vec           *mode;
256:   PetscScalar   *dots;
257:   PetscInt       dim, dimEmbed, n, m, b, d, i, j, off;

261:   PetscObjectGetComm((PetscObject)dm,&comm);
262:   DMGetDimension(dm, &dim);
263:   DMGetCoordinateDim(dm, &dimEmbed);
264:   DMGetLocalSection(dm, &section);
265:   DMGetGlobalSection(dm, &globalSection);
266:   PetscSectionGetConstrainedStorageSize(globalSection, &n);
267:   m    = nb * (dim*(dim+1))/2;
268:   PetscMalloc2(m, &mode, m, &dots);
269:   VecCreate(comm, &mode[0]);
270:   VecSetSizes(mode[0], n, PETSC_DETERMINE);
271:   VecSetUp(mode[0]);
272:   for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
273:   for (b = 0, off = 0; b < nb; ++b) {
274:     for (d = 0; d < m/nb; ++d) {
275:       PetscInt         ctx[2];
276:       PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private;
277:       void            *voidctx = (void *) (&ctx[0]);

279:       ctx[0] = dimEmbed;
280:       ctx[1] = d;
281:       DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);
282:       off   += nids[b];
283:     }
284:   }
285:   /* Orthonormalize system */
286:   for (i = 0; i < m; ++i) {
287:     PetscScalar dots[6];

289:     VecNormalize(mode[i], NULL);
290:     VecMDot(mode[i], m-i-1, mode+i+1, dots+i+1);
291:     for (j = i+1; j < m; ++j) {
292:       dots[j] *= -1.0;
293:       VecAXPY(mode[j], dots[j], mode[i]);
294:     }
295:   }
296:   MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
297:   for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
298:   PetscFree2(mode, dots);
299:   return(0);
300: }

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

311:   Input Parameters:
312: + dm - the DMPlex object
313: - height - the maximum projection height >= 0

315:   Level: advanced

317: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
318: @*/
319: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
320: {
321:   DM_Plex *plex = (DM_Plex *) dm->data;

325:   plex->maxProjectionHeight = height;
326:   return(0);
327: }

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

333:   Input Parameters:
334: . dm - the DMPlex object

336:   Output Parameters:
337: . height - the maximum projection height

339:   Level: intermediate

341: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
342: @*/
343: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
344: {
345:   DM_Plex *plex = (DM_Plex *) dm->data;

349:   *height = plex->maxProjectionHeight;
350:   return(0);
351: }

353: typedef struct {
354:   PetscReal    alpha; /* The first Euler angle, and in 2D the only one */
355:   PetscReal    beta;  /* The second Euler angle */
356:   PetscReal    gamma; /* The third Euler angle */
357:   PetscInt     dim;   /* The dimension of R */
358:   PetscScalar *R;     /* The rotation matrix, transforming a vector in the local basis to the global basis */
359:   PetscScalar *RT;    /* The transposed rotation matrix, transforming a vector in the global basis to the local basis */
360: } RotCtx;

362: /*
363:   Note: Following https://en.wikipedia.org/wiki/Euler_angles, we will specify Euler angles by extrinsic rotations, meaning that
364:   we rotate with respect to a fixed initial coordinate system, the local basis (x-y-z). The global basis (X-Y-Z) is reached as follows:
365:   $ The XYZ system rotates about the z axis by alpha. The X axis is now at angle alpha with respect to the x axis.
366:   $ The XYZ system rotates again about the x axis by beta. The Z axis is now at angle beta with respect to the z axis.
367:   $ The XYZ system rotates a third time about the z axis by gamma.
368: */
369: static PetscErrorCode DMPlexBasisTransformSetUp_Rotation_Internal(DM dm, void *ctx)
370: {
371:   RotCtx        *rc  = (RotCtx *) ctx;
372:   PetscInt       dim = rc->dim;
373:   PetscReal      c1, s1, c2, s2, c3, s3;

377:   PetscMalloc2(PetscSqr(dim), &rc->R, PetscSqr(dim), &rc->RT);
378:   switch (dim) {
379:   case 2:
380:     c1 = PetscCosReal(rc->alpha);s1 = PetscSinReal(rc->alpha);
381:     rc->R[0] =  c1;rc->R[1] = s1;
382:     rc->R[2] = -s1;rc->R[3] = c1;
383:     PetscArraycpy(rc->RT, rc->R, PetscSqr(dim));
384:     DMPlex_Transpose2D_Internal(rc->RT);break;
385:     break;
386:   case 3:
387:     c1 = PetscCosReal(rc->alpha);s1 = PetscSinReal(rc->alpha);
388:     c2 = PetscCosReal(rc->beta); s2 = PetscSinReal(rc->beta);
389:     c3 = PetscCosReal(rc->gamma);s3 = PetscSinReal(rc->gamma);
390:     rc->R[0] =  c1*c3 - c2*s1*s3;rc->R[1] =  c3*s1    + c1*c2*s3;rc->R[2] = s2*s3;
391:     rc->R[3] = -c1*s3 - c2*c3*s1;rc->R[4] =  c1*c2*c3 - s1*s3;   rc->R[5] = c3*s2;
392:     rc->R[6] =  s1*s2;           rc->R[7] = -c1*s2;              rc->R[8] = c2;
393:     PetscArraycpy(rc->RT, rc->R, PetscSqr(dim));
394:     DMPlex_Transpose3D_Internal(rc->RT);break;
395:     break;
396:   default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not supported", dim);
397:   }
398:   return(0);
399: }

401: static PetscErrorCode DMPlexBasisTransformDestroy_Rotation_Internal(DM dm, void *ctx)
402: {
403:   RotCtx        *rc = (RotCtx *) ctx;

407:   PetscFree2(rc->R, rc->RT);
408:   PetscFree(rc);
409:   return(0);
410: }

412: static PetscErrorCode DMPlexBasisTransformGetMatrix_Rotation_Internal(DM dm, const PetscReal x[], PetscBool l2g, const PetscScalar **A, void *ctx)
413: {
414:   RotCtx *rc = (RotCtx *) ctx;

418:   if (l2g) {*A = rc->R;}
419:   else     {*A = rc->RT;}
420:   return(0);
421: }

423: PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM dm, const PetscReal x[], PetscBool l2g, PetscInt dim, const PetscReal *y, PetscReal *z, void *ctx)
424: {

428:   #if defined(PETSC_USE_COMPLEX)
429:   switch (dim) {
430:     case 2:
431:     {
432:       PetscScalar yt[2], zt[2];

434:       yt[0] = y[0]; yt[1] = y[1];
435:       DMPlexBasisTransformApply_Internal(dm, x, l2g, dim, yt, zt, ctx);
436:       z[0] = PetscRealPart(zt[0]); z[1] = PetscRealPart(zt[1]);
437:     }
438:     break;
439:     case 3:
440:     {
441:       PetscScalar yt[3], zt[3];

443:       yt[0] = y[0]; yt[1] = y[1]; yt[2] = y[2];
444:       DMPlexBasisTransformApply_Internal(dm, x, l2g, dim, yt, zt, ctx);
445:       z[0] = PetscRealPart(zt[0]); z[1] = PetscRealPart(zt[1]); z[2] = PetscRealPart(zt[2]);
446:     }
447:     break;
448:   }
449:   #else
450:   DMPlexBasisTransformApply_Internal(dm, x, l2g, dim, y, z, ctx);
451:   #endif
452:   return(0);
453: }

455: PetscErrorCode DMPlexBasisTransformApply_Internal(DM dm, const PetscReal x[], PetscBool l2g, PetscInt dim, const PetscScalar *y, PetscScalar *z, void *ctx)
456: {
457:   const PetscScalar *A;
458:   PetscErrorCode     ierr;

461:   (*dm->transformGetMatrix)(dm, x, l2g, &A, ctx);
462:   switch (dim) {
463:   case 2: DMPlex_Mult2D_Internal(A, 1, y, z);break;
464:   case 3: DMPlex_Mult3D_Internal(A, 1, y, z);break;
465:   }
466:   return(0);
467: }

469: static PetscErrorCode DMPlexBasisTransformField_Internal(DM dm, DM tdm, Vec tv, PetscInt p, PetscInt f, PetscBool l2g, PetscScalar *a)
470: {
471:   PetscSection       ts;
472:   const PetscScalar *ta, *tva;
473:   PetscInt           dof;
474:   PetscErrorCode     ierr;

477:   DMGetLocalSection(tdm, &ts);
478:   PetscSectionGetFieldDof(ts, p, f, &dof);
479:   VecGetArrayRead(tv, &ta);
480:   DMPlexPointLocalFieldRead(tdm, p, f, ta, (void *) &tva);
481:   if (l2g) {
482:     switch (dof) {
483:     case 4: DMPlex_Mult2D_Internal(tva, 1, a, a);break;
484:     case 9: DMPlex_Mult3D_Internal(tva, 1, a, a);break;
485:     }
486:   } else {
487:     switch (dof) {
488:     case 4: DMPlex_MultTranspose2D_Internal(tva, 1, a, a);break;
489:     case 9: DMPlex_MultTranspose3D_Internal(tva, 1, a, a);break;
490:     }
491:   }
492:   VecRestoreArrayRead(tv, &ta);
493:   return(0);
494: }

496: static PetscErrorCode DMPlexBasisTransformFieldTensor_Internal(DM dm, DM tdm, Vec tv, PetscInt pf, PetscInt f, PetscInt pg, PetscInt g, PetscBool l2g, PetscInt lda, PetscScalar *a)
497: {
498:   PetscSection       s, ts;
499:   const PetscScalar *ta, *tvaf, *tvag;
500:   PetscInt           fdof, gdof, fpdof, gpdof;
501:   PetscErrorCode     ierr;

504:   DMGetLocalSection(dm, &s);
505:   DMGetLocalSection(tdm, &ts);
506:   PetscSectionGetFieldDof(s, pf, f, &fpdof);
507:   PetscSectionGetFieldDof(s, pg, g, &gpdof);
508:   PetscSectionGetFieldDof(ts, pf, f, &fdof);
509:   PetscSectionGetFieldDof(ts, pg, g, &gdof);
510:   VecGetArrayRead(tv, &ta);
511:   DMPlexPointLocalFieldRead(tdm, pf, f, ta, (void *) &tvaf);
512:   DMPlexPointLocalFieldRead(tdm, pg, g, ta, (void *) &tvag);
513:   if (l2g) {
514:     switch (fdof) {
515:     case 4: DMPlex_MatMult2D_Internal(tvaf, gpdof, lda, a, a);break;
516:     case 9: DMPlex_MatMult3D_Internal(tvaf, gpdof, lda, a, a);break;
517:     }
518:     switch (gdof) {
519:     case 4: DMPlex_MatMultTransposeLeft2D_Internal(tvag, fpdof, lda, a, a);break;
520:     case 9: DMPlex_MatMultTransposeLeft3D_Internal(tvag, fpdof, lda, a, a);break;
521:     }
522:   } else {
523:     switch (fdof) {
524:     case 4: DMPlex_MatMultTranspose2D_Internal(tvaf, gpdof, lda, a, a);break;
525:     case 9: DMPlex_MatMultTranspose3D_Internal(tvaf, gpdof, lda, a, a);break;
526:     }
527:     switch (gdof) {
528:     case 4: DMPlex_MatMultLeft2D_Internal(tvag, fpdof, lda, a, a);break;
529:     case 9: DMPlex_MatMultLeft3D_Internal(tvag, fpdof, lda, a, a);break;
530:     }
531:   }
532:   VecRestoreArrayRead(tv, &ta);
533:   return(0);
534: }

536: PetscErrorCode DMPlexBasisTransformPoint_Internal(DM dm, DM tdm, Vec tv, PetscInt p, PetscBool fieldActive[], PetscBool l2g, PetscScalar *a)
537: {
538:   PetscSection    s;
539:   PetscSection    clSection;
540:   IS              clPoints;
541:   const PetscInt *clp;
542:   PetscInt       *points = NULL;
543:   PetscInt        Nf, f, Np, cp, dof, d = 0;
544:   PetscErrorCode  ierr;

547:   DMGetLocalSection(dm, &s);
548:   PetscSectionGetNumFields(s, &Nf);
549:   DMPlexGetCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);
550:   for (f = 0; f < Nf; ++f) {
551:     for (cp = 0; cp < Np*2; cp += 2) {
552:       PetscSectionGetFieldDof(s, points[cp], f, &dof);
553:       if (!dof) continue;
554:       if (fieldActive[f]) {DMPlexBasisTransformField_Internal(dm, tdm, tv, points[cp], f, l2g, &a[d]);}
555:       d += dof;
556:     }
557:   }
558:   DMPlexRestoreCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);
559:   return(0);
560: }

562: PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM dm, DM tdm, Vec tv, PetscInt p, PetscBool l2g, PetscInt lda, PetscScalar *a)
563: {
564:   PetscSection    s;
565:   PetscSection    clSection;
566:   IS              clPoints;
567:   const PetscInt *clp;
568:   PetscInt       *points = NULL;
569:   PetscInt        Nf, f, g, Np, cpf, cpg, fdof, gdof, r, c = 0;
570:   PetscErrorCode  ierr;

573:   DMGetLocalSection(dm, &s);
574:   PetscSectionGetNumFields(s, &Nf);
575:   DMPlexGetCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);
576:   for (f = 0, r = 0; f < Nf; ++f) {
577:     for (cpf = 0; cpf < Np*2; cpf += 2) {
578:       PetscSectionGetFieldDof(s, points[cpf], f, &fdof);
579:       for (g = 0, c = 0; g < Nf; ++g) {
580:         for (cpg = 0; cpg < Np*2; cpg += 2) {
581:           PetscSectionGetFieldDof(s, points[cpg], g, &gdof);
582:           DMPlexBasisTransformFieldTensor_Internal(dm, tdm, tv, points[cpf], f, points[cpg], g, l2g, lda, &a[r*lda+c]);
583:           c += gdof;
584:         }
585:       }
586:       if (c != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of columns %D should be %D", c, lda);
587:       r += fdof;
588:     }
589:   }
590:   if (r != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of rows %D should be %D", c, lda);
591:   DMPlexRestoreCompressedClosure(dm, s, p, &Np, &points, &clSection, &clPoints, &clp);
592:   return(0);
593: }

595: static PetscErrorCode DMPlexBasisTransform_Internal(DM dm, Vec lv, PetscBool l2g)
596: {
597:   DM                 tdm;
598:   Vec                tv;
599:   PetscSection       ts, s;
600:   const PetscScalar *ta;
601:   PetscScalar       *a, *va;
602:   PetscInt           pStart, pEnd, p, Nf, f;
603:   PetscErrorCode     ierr;

606:   DMGetBasisTransformDM_Internal(dm, &tdm);
607:   DMGetBasisTransformVec_Internal(dm, &tv);
608:   DMGetLocalSection(tdm, &ts);
609:   DMGetLocalSection(dm, &s);
610:   PetscSectionGetChart(s, &pStart, &pEnd);
611:   PetscSectionGetNumFields(s, &Nf);
612:   VecGetArray(lv, &a);
613:   VecGetArrayRead(tv, &ta);
614:   for (p = pStart; p < pEnd; ++p) {
615:     for (f = 0; f < Nf; ++f) {
616:       DMPlexPointLocalFieldRef(dm, p, f, a, (void *) &va);
617:       DMPlexBasisTransformField_Internal(dm, tdm, tv, p, f, l2g, va);
618:     }
619:   }
620:   VecRestoreArray(lv, &a);
621:   VecRestoreArrayRead(tv, &ta);
622:   return(0);
623: }

625: /*@
626:   DMPlexGlobalToLocalBasis - Transform the values in the given local vector from the global basis to the local basis

628:   Input Parameters:
629: + dm - The DM
630: - lv - A local vector with values in the global basis

632:   Output Parameters:
633: . lv - A local vector with values in the local basis

635:   Note: This method is only intended to be called inside DMGlobalToLocal(). It is unlikely that a user will have a local vector full of coefficients for the global basis unless they are reimplementing GlobalToLocal.

637:   Level: developer

639: .seealso: DMPlexLocalToGlobalBasis(), DMGetLocalSection(), DMPlexCreateBasisRotation()
640: @*/
641: PetscErrorCode DMPlexGlobalToLocalBasis(DM dm, Vec lv)
642: {

648:   DMPlexBasisTransform_Internal(dm, lv, PETSC_FALSE);
649:   return(0);
650: }

652: /*@
653:   DMPlexLocalToGlobalBasis - Transform the values in the given local vector from the local basis to the global basis

655:   Input Parameters:
656: + dm - The DM
657: - lv - A local vector with values in the local basis

659:   Output Parameters:
660: . lv - A local vector with values in the global basis

662:   Note: This method is only intended to be called inside DMGlobalToLocal(). It is unlikely that a user would want a local vector full of coefficients for the global basis unless they are reimplementing GlobalToLocal.

664:   Level: developer

666: .seealso: DMPlexGlobalToLocalBasis(), DMGetLocalSection(), DMPlexCreateBasisRotation()
667: @*/
668: PetscErrorCode DMPlexLocalToGlobalBasis(DM dm, Vec lv)
669: {

675:   DMPlexBasisTransform_Internal(dm, lv, PETSC_TRUE);
676:   return(0);
677: }

679: /*@
680:   DMPlexCreateBasisRotation - Create an internal transformation from the global basis, used to specify boundary conditions
681:     and global solutions, to a local basis, appropriate for discretization integrals and assembly.

683:   Input Parameters:
684: + dm    - The DM
685: . alpha - The first Euler angle, and in 2D the only one
686: . beta  - The second Euler angle
687: - gamma - The third Euler angle

689:   Note: Following https://en.wikipedia.org/wiki/Euler_angles, we will specify Euler angles by extrinsic rotations, meaning that
690:   we rotate with respect to a fixed initial coordinate system, the local basis (x-y-z). The global basis (X-Y-Z) is reached as follows:
691:   $ The XYZ system rotates about the z axis by alpha. The X axis is now at angle alpha with respect to the x axis.
692:   $ The XYZ system rotates again about the x axis by beta. The Z axis is now at angle beta with respect to the z axis.
693:   $ The XYZ system rotates a third time about the z axis by gamma.

695:   Level: developer

697: .seealso: DMPlexGlobalToLocalBasis(), DMPlexLocalToGlobalBasis()
698: @*/
699: PetscErrorCode DMPlexCreateBasisRotation(DM dm, PetscReal alpha, PetscReal beta, PetscReal gamma)
700: {
701:   RotCtx        *rc;
702:   PetscInt       cdim;

705:   DMGetCoordinateDim(dm, &cdim);
706:   PetscMalloc1(1, &rc);
707:   dm->transformCtx       = rc;
708:   dm->transformSetUp     = DMPlexBasisTransformSetUp_Rotation_Internal;
709:   dm->transformDestroy   = DMPlexBasisTransformDestroy_Rotation_Internal;
710:   dm->transformGetMatrix = DMPlexBasisTransformGetMatrix_Rotation_Internal;
711:   rc->dim   = cdim;
712:   rc->alpha = alpha;
713:   rc->beta  = beta;
714:   rc->gamma = gamma;
715:   (*dm->transformSetUp)(dm, dm->transformCtx);
716:   DMConstructBasisTransform_Internal(dm);
717:   return(0);
718: }

720: /*@C
721:   DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector using a function of the coordinates

723:   Input Parameters:
724: + dm     - The DM, with a PetscDS that matches the problem being constrained
725: . time   - The time
726: . field  - The field to constrain
727: . Nc     - The number of constrained field components, or 0 for all components
728: . comps  - An array of constrained component numbers, or NULL for all components
729: . label  - The DMLabel defining constrained points
730: . numids - The number of DMLabel ids for constrained points
731: . ids    - An array of ids for constrained points
732: . func   - A pointwise function giving boundary values
733: - ctx    - An optional user context for bcFunc

735:   Output Parameter:
736: . locX   - A local vector to receives the boundary values

738:   Level: developer

740: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMPlexInsertBoundaryValuesEssentialBdField(), DMAddBoundary()
741: @*/
742: PetscErrorCode DMPlexInsertBoundaryValuesEssential(DM dm, PetscReal time, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX)
743: {
744:   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
745:   void            **ctxs;
746:   PetscInt          numFields;
747:   PetscErrorCode    ierr;

750:   DMGetNumFields(dm, &numFields);
751:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
752:   funcs[field] = func;
753:   ctxs[field]  = ctx;
754:   DMProjectFunctionLabelLocal(dm, time, label, numids, ids, Nc, comps, funcs, ctxs, INSERT_BC_VALUES, locX);
755:   PetscFree2(funcs,ctxs);
756:   return(0);
757: }

759: /*@C
760:   DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector using a function of the coordinates and field data

762:   Input Parameters:
763: + dm     - The DM, with a PetscDS that matches the problem being constrained
764: . time   - The time
765: . locU   - A local vector with the input solution values
766: . field  - The field to constrain
767: . Nc     - The number of constrained field components, or 0 for all components
768: . comps  - An array of constrained component numbers, or NULL for all components
769: . label  - The DMLabel defining constrained points
770: . numids - The number of DMLabel ids for constrained points
771: . ids    - An array of ids for constrained points
772: . func   - A pointwise function giving boundary values
773: - ctx    - An optional user context for bcFunc

775:   Output Parameter:
776: . locX   - A local vector to receives the boundary values

778:   Level: developer

780: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialBdField(), DMAddBoundary()
781: @*/
782: PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
783:                                                         void (*func)(PetscInt, PetscInt, PetscInt,
784:                                                                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
785:                                                                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
786:                                                                      PetscReal, const PetscReal[], PetscInt, const PetscScalar[],
787:                                                                      PetscScalar[]),
788:                                                         void *ctx, Vec locX)
789: {
790:   void (**funcs)(PetscInt, PetscInt, PetscInt,
791:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
792:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
793:                  PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
794:   void            **ctxs;
795:   PetscInt          numFields;
796:   PetscErrorCode    ierr;

799:   DMGetNumFields(dm, &numFields);
800:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
801:   funcs[field] = func;
802:   ctxs[field]  = ctx;
803:   DMProjectFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
804:   PetscFree2(funcs,ctxs);
805:   return(0);
806: }

808: /*@C
809:   DMPlexInsertBoundaryValuesEssentialBdField - Insert boundary values into a local vector using a function of the coodinates and boundary field data

811:   Collective on dm

813:   Input Parameters:
814: + dm     - The DM, with a PetscDS that matches the problem being constrained
815: . time   - The time
816: . locU   - A local vector with the input solution values
817: . field  - The field to constrain
818: . Nc     - The number of constrained field components, or 0 for all components
819: . comps  - An array of constrained component numbers, or NULL for all components
820: . label  - The DMLabel defining constrained points
821: . numids - The number of DMLabel ids for constrained points
822: . ids    - An array of ids for constrained points
823: . func   - A pointwise function giving boundary values, the calling sequence is given in DMProjectBdFieldLabelLocal()
824: - ctx    - An optional user context for bcFunc

826:   Output Parameter:
827: . locX   - A local vector to receive the boundary values

829:   Level: developer

831: .seealso: DMProjectBdFieldLabelLocal(), DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
832: @*/
833: PetscErrorCode DMPlexInsertBoundaryValuesEssentialBdField(DM dm, PetscReal time, Vec locU, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
834:                                                           void (*func)(PetscInt, PetscInt, PetscInt,
835:                                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
836:                                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
837:                                                                        PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[],
838:                                                                        PetscScalar[]),
839:                                                           void *ctx, Vec locX)
840: {
841:   void (**funcs)(PetscInt, PetscInt, PetscInt,
842:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
843:                  const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
844:                  PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
845:   void            **ctxs;
846:   PetscInt          numFields;
847:   PetscErrorCode    ierr;

850:   DMGetNumFields(dm, &numFields);
851:   PetscCalloc2(numFields,&funcs,numFields,&ctxs);
852:   funcs[field] = func;
853:   ctxs[field]  = ctx;
854:   DMProjectBdFieldLabelLocal(dm, time, label, numids, ids, Nc, comps, locU, funcs, INSERT_BC_VALUES, locX);
855:   PetscFree2(funcs,ctxs);
856:   return(0);
857: }

859: /*@C
860:   DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector

862:   Input Parameters:
863: + dm     - The DM, with a PetscDS that matches the problem being constrained
864: . time   - The time
865: . faceGeometry - A vector with the FVM face geometry information
866: . cellGeometry - A vector with the FVM cell geometry information
867: . Grad         - A vector with the FVM cell gradient information
868: . field  - The field to constrain
869: . Nc     - The number of constrained field components, or 0 for all components
870: . comps  - An array of constrained component numbers, or NULL for all components
871: . label  - The DMLabel defining constrained points
872: . numids - The number of DMLabel ids for constrained points
873: . ids    - An array of ids for constrained points
874: . func   - A pointwise function giving boundary values
875: - ctx    - An optional user context for bcFunc

877:   Output Parameter:
878: . locX   - A local vector to receives the boundary values

880:   Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary()

882:   Level: developer

884: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
885: @*/
886: PetscErrorCode DMPlexInsertBoundaryValuesRiemann(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, PetscInt field, PetscInt Nc, const PetscInt comps[], DMLabel label, PetscInt numids, const PetscInt ids[],
887:                                                  PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
888: {
889:   PetscDS            prob;
890:   PetscSF            sf;
891:   DM                 dmFace, dmCell, dmGrad;
892:   const PetscScalar *facegeom, *cellgeom = NULL, *grad;
893:   const PetscInt    *leaves;
894:   PetscScalar       *x, *fx;
895:   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
896:   PetscErrorCode     ierr, ierru = 0;

899:   DMGetPointSF(dm, &sf);
900:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
901:   nleaves = PetscMax(0, nleaves);
902:   DMGetDimension(dm, &dim);
903:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
904:   DMGetDS(dm, &prob);
905:   VecGetDM(faceGeometry, &dmFace);
906:   VecGetArrayRead(faceGeometry, &facegeom);
907:   if (cellGeometry) {
908:     VecGetDM(cellGeometry, &dmCell);
909:     VecGetArrayRead(cellGeometry, &cellgeom);
910:   }
911:   if (Grad) {
912:     PetscFV fv;

914:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
915:     VecGetDM(Grad, &dmGrad);
916:     VecGetArrayRead(Grad, &grad);
917:     PetscFVGetNumComponents(fv, &pdim);
918:     DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);
919:   }
920:   VecGetArray(locX, &x);
921:   for (i = 0; i < numids; ++i) {
922:     IS              faceIS;
923:     const PetscInt *faces;
924:     PetscInt        numFaces, f;

926:     DMLabelGetStratumIS(label, ids[i], &faceIS);
927:     if (!faceIS) continue; /* No points with that id on this process */
928:     ISGetLocalSize(faceIS, &numFaces);
929:     ISGetIndices(faceIS, &faces);
930:     for (f = 0; f < numFaces; ++f) {
931:       const PetscInt         face = faces[f], *cells;
932:       PetscFVFaceGeom        *fg;

934:       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
935:       PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
936:       if (loc >= 0) continue;
937:       DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
938:       DMPlexGetSupport(dm, face, &cells);
939:       if (Grad) {
940:         PetscFVCellGeom       *cg;
941:         PetscScalar           *cx, *cgrad;
942:         PetscScalar           *xG;
943:         PetscReal              dx[3];
944:         PetscInt               d;

946:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
947:         DMPlexPointLocalRead(dm, cells[0], x, &cx);
948:         DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
949:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
950:         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
951:         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
952:         ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
953:         if (ierru) {
954:           ISRestoreIndices(faceIS, &faces);
955:           ISDestroy(&faceIS);
956:           goto cleanup;
957:         }
958:       } else {
959:         PetscScalar       *xI;
960:         PetscScalar       *xG;

962:         DMPlexPointLocalRead(dm, cells[0], x, &xI);
963:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
964:         ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
965:         if (ierru) {
966:           ISRestoreIndices(faceIS, &faces);
967:           ISDestroy(&faceIS);
968:           goto cleanup;
969:         }
970:       }
971:     }
972:     ISRestoreIndices(faceIS, &faces);
973:     ISDestroy(&faceIS);
974:   }
975:   cleanup:
976:   VecRestoreArray(locX, &x);
977:   if (Grad) {
978:     DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);
979:     VecRestoreArrayRead(Grad, &grad);
980:   }
981:   if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
982:   VecRestoreArrayRead(faceGeometry, &facegeom);
983:   CHKERRQ(ierru);
984:   return(0);
985: }

987: PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
988: {
989:   PetscDS        prob;
990:   PetscInt       numBd, b;

994:   DMGetDS(dm, &prob);
995:   PetscDSGetNumBoundary(prob, &numBd);
996:   for (b = 0; b < numBd; ++b) {
997:     DMBoundaryConditionType type;
998:     const char             *name, *labelname;
999:     DMLabel                 label;
1000:     PetscInt                field, Nc;
1001:     const PetscInt         *comps;
1002:     PetscObject             obj;
1003:     PetscClassId            id;
1004:     void                    (*func)(void);
1005:     PetscInt                numids;
1006:     const PetscInt         *ids;
1007:     void                   *ctx;

1009:     DMGetBoundary(dm, b, &type, &name, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);
1010:     if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
1011:     DMGetLabel(dm, labelname, &label);
1012:     if (!label) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "Label %s for boundary condition %s does not exist in the DM", labelname, name);
1013:     DMGetField(dm, field, NULL, &obj);
1014:     PetscObjectGetClassId(obj, &id);
1015:     if (id == PETSCFE_CLASSID) {
1016:       switch (type) {
1017:         /* for FEM, there is no insertion to be done for non-essential boundary conditions */
1018:       case DM_BC_ESSENTIAL:
1019:         DMPlexLabelAddCells(dm,label);
1020:         DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
1021:         DMPlexLabelClearCells(dm,label);
1022:         break;
1023:       case DM_BC_ESSENTIAL_FIELD:
1024:         DMPlexLabelAddCells(dm,label);
1025:         DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
1026:                                                         (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1027:                                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1028:                                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);
1029:         DMPlexLabelClearCells(dm,label);
1030:         break;
1031:       default: break;
1032:       }
1033:     } else if (id == PETSCFV_CLASSID) {
1034:       if (!faceGeomFVM) continue;
1035:       DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids,
1036:                                                (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
1037:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1038:   }
1039:   return(0);
1040: }

1042: /*@
1043:   DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector

1045:   Input Parameters:
1046: + dm - The DM
1047: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
1048: . time - The time
1049: . faceGeomFVM - Face geometry data for FV discretizations
1050: . cellGeomFVM - Cell geometry data for FV discretizations
1051: - gradFVM - Gradient reconstruction data for FV discretizations

1053:   Output Parameters:
1054: . locX - Solution updated with boundary values

1056:   Level: developer

1058: .seealso: DMProjectFunctionLabelLocal()
1059: @*/
1060: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
1061: {

1070:   PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));
1071:   return(0);
1072: }

1074: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1075: {
1076:   Vec              localX;
1077:   PetscErrorCode   ierr;

1080:   DMGetLocalVector(dm, &localX);
1081:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);
1082:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1083:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1084:   DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);
1085:   DMRestoreLocalVector(dm, &localX);
1086:   return(0);
1087: }

1089: /*@C
1090:   DMComputeL2DiffLocal - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

1092:   Collective on dm

1094:   Input Parameters:
1095: + dm     - The DM
1096: . time   - The time
1097: . funcs  - The functions to evaluate for each field component
1098: . ctxs   - Optional array of contexts to pass to each function, or NULL.
1099: - localX - The coefficient vector u_h, a local vector

1101:   Output Parameter:
1102: . diff - The diff ||u - u_h||_2

1104:   Level: developer

1106: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
1107: @*/
1108: PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff)
1109: {
1110:   const PetscInt   debug = ((DM_Plex*)dm->data)->printL2;
1111:   DM               tdm;
1112:   Vec              tv;
1113:   PetscSection     section;
1114:   PetscQuadrature  quad;
1115:   PetscFEGeom      fegeom;
1116:   PetscScalar     *funcVal, *interpolant;
1117:   PetscReal       *coords, *gcoords;
1118:   PetscReal        localDiff = 0.0;
1119:   const PetscReal *quadWeights;
1120:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, c, field, fieldOffset;
1121:   PetscBool        transform;
1122:   PetscErrorCode   ierr;

1125:   DMGetDimension(dm, &dim);
1126:   DMGetCoordinateDim(dm, &coordDim);
1127:   fegeom.dimEmbed = coordDim;
1128:   DMGetLocalSection(dm, &section);
1129:   PetscSectionGetNumFields(section, &numFields);
1130:   DMGetBasisTransformDM_Internal(dm, &tdm);
1131:   DMGetBasisTransformVec_Internal(dm, &tv);
1132:   DMHasBasisTransform(dm, &transform);
1133:   for (field = 0; field < numFields; ++field) {
1134:     PetscObject  obj;
1135:     PetscClassId id;
1136:     PetscInt     Nc;

1138:     DMGetField(dm, field, NULL, &obj);
1139:     PetscObjectGetClassId(obj, &id);
1140:     if (id == PETSCFE_CLASSID) {
1141:       PetscFE fe = (PetscFE) obj;

1143:       PetscFEGetQuadrature(fe, &quad);
1144:       PetscFEGetNumComponents(fe, &Nc);
1145:     } else if (id == PETSCFV_CLASSID) {
1146:       PetscFV fv = (PetscFV) obj;

1148:       PetscFVGetQuadrature(fv, &quad);
1149:       PetscFVGetNumComponents(fv, &Nc);
1150:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1151:     numComponents += Nc;
1152:   }
1153:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
1154:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1155:   PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1156:   DMPlexGetVTKCellHeight(dm, &cellHeight);
1157:   DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd);
1158:   for (c = cStart; c < cEnd; ++c) {
1159:     PetscScalar *x = NULL;
1160:     PetscReal    elemDiff = 0.0;
1161:     PetscInt     qc = 0;

1163:     DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1164:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

1166:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1167:       PetscObject  obj;
1168:       PetscClassId id;
1169:       void * const ctx = ctxs ? ctxs[field] : NULL;
1170:       PetscInt     Nb, Nc, q, fc;

1172:       DMGetField(dm, field, NULL, &obj);
1173:       PetscObjectGetClassId(obj, &id);
1174:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1175:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1176:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1177:       if (debug) {
1178:         char title[1024];
1179:         PetscSNPrintf(title, 1023, "Solution for Field %D", field);
1180:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1181:       }
1182:       for (q = 0; q < Nq; ++q) {
1183:         PetscFEGeom qgeom;

1185:         qgeom.dimEmbed = fegeom.dimEmbed;
1186:         qgeom.J        = &fegeom.J[q*coordDim*coordDim];
1187:         qgeom.invJ     = &fegeom.invJ[q*coordDim*coordDim];
1188:         qgeom.detJ     = &fegeom.detJ[q];
1189:         if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, point %D", (double)fegeom.detJ[q], c, q);
1190:         if (transform) {
1191:           gcoords = &coords[coordDim*Nq];
1192:           DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);
1193:         } else {
1194:           gcoords = &coords[coordDim*q];
1195:         }
1196:         (*funcs[field])(coordDim, time, gcoords, Nc, funcVal, ctx);
1197:         if (ierr) {
1198:           PetscErrorCode ierr2;
1199:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1200:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1201:           ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1202:           
1203:         }
1204:         if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1205:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1206:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1207:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1208:         for (fc = 0; fc < Nc; ++fc) {
1209:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1210:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %D field %D,%D point %g %g %g diff %g\n", c, field, fc, (double)(coordDim > 0 ? coords[coordDim*q] : 0.), (double)(coordDim > 1 ? coords[coordDim*q+1] : 0.),(double)(coordDim > 2 ? coords[coordDim*q+2] : 0.), (double)(PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]));}
1211:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1212:         }
1213:       }
1214:       fieldOffset += Nb;
1215:       qc += Nc;
1216:     }
1217:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1218:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %D diff %g\n", c, (double)elemDiff);}
1219:     localDiff += elemDiff;
1220:   }
1221:   PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1222:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1223:   *diff = PetscSqrtReal(*diff);
1224:   return(0);
1225: }

1227: PetscErrorCode DMComputeL2GradientDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
1228: {
1229:   const PetscInt   debug = ((DM_Plex*)dm->data)->printL2;
1230:   DM               tdm;
1231:   PetscSection     section;
1232:   PetscQuadrature  quad;
1233:   Vec              localX, tv;
1234:   PetscScalar     *funcVal, *interpolant;
1235:   const PetscReal *quadWeights;
1236:   PetscFEGeom      fegeom;
1237:   PetscReal       *coords, *gcoords;
1238:   PetscReal        localDiff = 0.0;
1239:   PetscInt         dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset;
1240:   PetscBool        transform;
1241:   PetscErrorCode   ierr;

1244:   DMGetDimension(dm, &dim);
1245:   DMGetCoordinateDim(dm, &coordDim);
1246:   fegeom.dimEmbed = coordDim;
1247:   DMGetLocalSection(dm, &section);
1248:   PetscSectionGetNumFields(section, &numFields);
1249:   DMGetLocalVector(dm, &localX);
1250:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1251:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1252:   DMGetBasisTransformDM_Internal(dm, &tdm);
1253:   DMGetBasisTransformVec_Internal(dm, &tv);
1254:   DMHasBasisTransform(dm, &transform);
1255:   for (field = 0; field < numFields; ++field) {
1256:     PetscFE  fe;
1257:     PetscInt Nc;

1259:     DMGetField(dm, field, NULL, (PetscObject *) &fe);
1260:     PetscFEGetQuadrature(fe, &quad);
1261:     PetscFEGetNumComponents(fe, &Nc);
1262:     numComponents += Nc;
1263:   }
1264:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
1265:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1266:   /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1267:   PetscMalloc6(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ,numComponents*coordDim,&interpolant,Nq,&fegeom.detJ);
1268:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1269:   for (c = cStart; c < cEnd; ++c) {
1270:     PetscScalar *x = NULL;
1271:     PetscReal    elemDiff = 0.0;
1272:     PetscInt     qc = 0;

1274:     DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1275:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

1277:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1278:       PetscFE          fe;
1279:       void * const     ctx = ctxs ? ctxs[field] : NULL;
1280:       PetscInt         Nb, Nc, q, fc;

1282:       DMGetField(dm, field, NULL, (PetscObject *) &fe);
1283:       PetscFEGetDimension(fe, &Nb);
1284:       PetscFEGetNumComponents(fe, &Nc);
1285:       if (debug) {
1286:         char title[1024];
1287:         PetscSNPrintf(title, 1023, "Solution for Field %D", field);
1288:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1289:       }
1290:       for (q = 0; q < Nq; ++q) {
1291:         PetscFEGeom qgeom;

1293:         qgeom.dimEmbed = fegeom.dimEmbed;
1294:         qgeom.J        = &fegeom.J[q*coordDim*coordDim];
1295:         qgeom.invJ     = &fegeom.invJ[q*coordDim*coordDim];
1296:         qgeom.detJ     = &fegeom.detJ[q];
1297:         if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)fegeom.detJ[q], c, q);
1298:         if (transform) {
1299:           gcoords = &coords[coordDim*Nq];
1300:           DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);
1301:         } else {
1302:           gcoords = &coords[coordDim*q];
1303:         }
1304:         (*funcs[field])(coordDim, time, gcoords, n, Nc, funcVal, ctx);
1305:         if (ierr) {
1306:           PetscErrorCode ierr2;
1307:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1308:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1309:           ierr2 = PetscFree6(funcVal,coords,fegeom.J,fegeom.invJ,interpolant,fegeom.detJ);CHKERRQ(ierr2);
1310:           
1311:         }
1312:         if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1313:         PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], &qgeom, q, interpolant);
1314:         /* Overwrite with the dot product if the normal is given */
1315:         if (n) {
1316:           for (fc = 0; fc < Nc; ++fc) {
1317:             PetscScalar sum = 0.0;
1318:             PetscInt    d;
1319:             for (d = 0; d < dim; ++d) sum += interpolant[fc*dim+d]*n[d];
1320:             interpolant[fc] = sum;
1321:           }
1322:         }
1323:         for (fc = 0; fc < Nc; ++fc) {
1324:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1325:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %D fieldDer %D,%D diff %g\n", c, field, fc, (double)(PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]));}
1326:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1327:         }
1328:       }
1329:       fieldOffset += Nb;
1330:       qc          += Nc;
1331:     }
1332:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1333:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %D diff %g\n", c, (double)elemDiff);}
1334:     localDiff += elemDiff;
1335:   }
1336:   PetscFree6(funcVal,coords,fegeom.J,fegeom.invJ,interpolant,fegeom.detJ);
1337:   DMRestoreLocalVector(dm, &localX);
1338:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1339:   *diff = PetscSqrtReal(*diff);
1340:   return(0);
1341: }

1343: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1344: {
1345:   const PetscInt   debug = ((DM_Plex*)dm->data)->printL2;
1346:   DM               tdm;
1347:   DMLabel          depthLabel;
1348:   PetscSection     section;
1349:   Vec              localX, tv;
1350:   PetscReal       *localDiff;
1351:   PetscInt         dim, depth, dE, Nf, f, Nds, s;
1352:   PetscBool        transform;
1353:   PetscErrorCode   ierr;

1356:   DMGetDimension(dm, &dim);
1357:   DMGetCoordinateDim(dm, &dE);
1358:   DMGetLocalSection(dm, &section);
1359:   DMGetLocalVector(dm, &localX);
1360:   DMGetBasisTransformDM_Internal(dm, &tdm);
1361:   DMGetBasisTransformVec_Internal(dm, &tv);
1362:   DMHasBasisTransform(dm, &transform);
1363:   DMGetNumFields(dm, &Nf);
1364:   DMPlexGetDepthLabel(dm, &depthLabel);
1365:   DMLabelGetNumValues(depthLabel, &depth);

1367:   VecSet(localX, 0.0);
1368:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1369:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1370:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1371:   DMGetNumDS(dm, &Nds);
1372:   PetscCalloc1(Nf, &localDiff);
1373:   for (s = 0; s < Nds; ++s) {
1374:     PetscDS          ds;
1375:     DMLabel          label;
1376:     IS               fieldIS, pointIS;
1377:     const PetscInt  *fields, *points = NULL;
1378:     PetscQuadrature  quad;
1379:     const PetscReal *quadPoints, *quadWeights;
1380:     PetscFEGeom      fegeom;
1381:     PetscReal       *coords, *gcoords;
1382:     PetscScalar     *funcVal, *interpolant;
1383:     PetscInt         qNc, Nq, totNc, cStart = 0, cEnd, c, dsNf;

1385:     DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);
1386:     ISGetIndices(fieldIS, &fields);
1387:     PetscDSGetNumFields(ds, &dsNf);
1388:     PetscDSGetTotalComponents(ds, &totNc);
1389:     PetscDSGetQuadrature(ds, &quad);
1390:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1391:     if ((qNc != 1) && (qNc != totNc)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, totNc);
1392:     PetscCalloc6(totNc, &funcVal, totNc, &interpolant, dE*(Nq+1), &coords,Nq, &fegeom.detJ, dE*dE*Nq, &fegeom.J, dE*dE*Nq, &fegeom.invJ);
1393:     if (!label) {
1394:       DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1395:     } else {
1396:       DMLabelGetStratumIS(label, 1, &pointIS);
1397:       ISGetLocalSize(pointIS, &cEnd);
1398:       ISGetIndices(pointIS, &points);
1399:     }
1400:     for (c = cStart; c < cEnd; ++c) {
1401:       const PetscInt cell = points ? points[c] : c;
1402:       PetscScalar   *x    = NULL;
1403:       PetscInt       qc   = 0, fOff = 0, dep;

1405:       DMLabelGetValue(depthLabel, cell, &dep);
1406:       if (dep != depth-1) continue;
1407:       DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1408:       DMPlexVecGetClosure(dm, NULL, localX, cell, NULL, &x);
1409:       for (f = 0; f < dsNf; ++f) {
1410:         PetscObject  obj;
1411:         PetscClassId id;
1412:         void * const ctx = ctxs ? ctxs[fields[f]] : NULL;
1413:         PetscInt     Nb, Nc, q, fc;
1414:         PetscReal    elemDiff = 0.0;

1416:         PetscDSGetDiscretization(ds, f, &obj);
1417:         PetscObjectGetClassId(obj, &id);
1418:         if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1419:         else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1420:         else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", fields[f]);
1421:         if (debug) {
1422:           char title[1024];
1423:           PetscSNPrintf(title, 1023, "Solution for Field %D", fields[f]);
1424:           DMPrintCellVector(cell, title, Nb, &x[fOff]);
1425:         }
1426:         for (q = 0; q < Nq; ++q) {
1427:           PetscFEGeom qgeom;

1429:           qgeom.dimEmbed = fegeom.dimEmbed;
1430:           qgeom.J        = &fegeom.J[q*dE*dE];
1431:           qgeom.invJ     = &fegeom.invJ[q*dE*dE];
1432:           qgeom.detJ     = &fegeom.detJ[q];
1433:           if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for cell %D, quadrature point %D", (double)fegeom.detJ[q], cell, q);
1434:           if (transform) {
1435:             gcoords = &coords[dE*Nq];
1436:             DMPlexBasisTransformApplyReal_Internal(dm, &coords[dE*q], PETSC_TRUE, dE, &coords[dE*q], gcoords, dm->transformCtx);
1437:           } else {
1438:             gcoords = &coords[dE*q];
1439:           }
1440:           (*funcs[fields[f]])(dE, time, gcoords, Nc, funcVal, ctx);
1441:           if (ierr) {
1442:             PetscErrorCode ierr2;
1443:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, cell, NULL, &x);CHKERRQ(ierr2);
1444:             ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1445:             ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1446:             
1447:           }
1448:           if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[dE*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1449:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fOff], &qgeom, q, interpolant);}
1450:           else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fOff], q, interpolant);}
1451:           else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", fields[f]);
1452:           for (fc = 0; fc < Nc; ++fc) {
1453:             const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1454:             if (debug) {PetscPrintf(PETSC_COMM_SELF, "    cell %D field %D,%D point %g %g %g diff %g\n", cell, fields[f], fc, (double)(dE > 0 ? coords[dE*q] : 0.), (double)(dE > 1 ? coords[dE*q+1] : 0.),(double)(dE > 2 ? coords[dE*q+2] : 0.), (double)(PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]));}
1455:             elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1456:           }
1457:         }
1458:         fOff += Nb;
1459:         qc   += Nc;
1460:         localDiff[fields[f]] += elemDiff;
1461:         if (debug) {PetscPrintf(PETSC_COMM_SELF, "  cell %D field %D cum diff %g\n", cell, fields[f], (double)localDiff[fields[f]]);}
1462:       }
1463:       DMPlexVecRestoreClosure(dm, NULL, localX, cell, NULL, &x);
1464:     }
1465:     if (label) {
1466:       ISRestoreIndices(pointIS, &points);
1467:       ISDestroy(&pointIS);
1468:     }
1469:     ISRestoreIndices(fieldIS, &fields);
1470:     PetscFree6(funcVal, interpolant, coords, fegeom.detJ, fegeom.J, fegeom.invJ);
1471:   }
1472:   DMRestoreLocalVector(dm, &localX);
1473:   MPIU_Allreduce(localDiff, diff, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1474:   PetscFree(localDiff);
1475:   for (f = 0; f < Nf; ++f) diff[f] = PetscSqrtReal(diff[f]);
1476:   return(0);
1477: }

1479: /*@C
1480:   DMPlexComputeL2DiffVec - This function computes the cellwise L_2 difference between a function u and an FEM interpolant solution u_h, and stores it in a Vec.

1482:   Collective on dm

1484:   Input Parameters:
1485: + dm    - The DM
1486: . time  - The time
1487: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
1488: . ctxs  - Optional array of contexts to pass to each function, or NULL.
1489: - X     - The coefficient vector u_h

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

1494:   Level: developer

1496: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1497: @*/
1498: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
1499: {
1500:   PetscSection     section;
1501:   PetscQuadrature  quad;
1502:   Vec              localX;
1503:   PetscFEGeom      fegeom;
1504:   PetscScalar     *funcVal, *interpolant;
1505:   PetscReal       *coords;
1506:   const PetscReal *quadPoints, *quadWeights;
1507:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, c, field, fieldOffset;
1508:   PetscErrorCode   ierr;

1511:   VecSet(D, 0.0);
1512:   DMGetDimension(dm, &dim);
1513:   DMGetCoordinateDim(dm, &coordDim);
1514:   DMGetLocalSection(dm, &section);
1515:   PetscSectionGetNumFields(section, &numFields);
1516:   DMGetLocalVector(dm, &localX);
1517:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1518:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1519:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1520:   for (field = 0; field < numFields; ++field) {
1521:     PetscObject  obj;
1522:     PetscClassId id;
1523:     PetscInt     Nc;

1525:     DMGetField(dm, field, NULL, &obj);
1526:     PetscObjectGetClassId(obj, &id);
1527:     if (id == PETSCFE_CLASSID) {
1528:       PetscFE fe = (PetscFE) obj;

1530:       PetscFEGetQuadrature(fe, &quad);
1531:       PetscFEGetNumComponents(fe, &Nc);
1532:     } else if (id == PETSCFV_CLASSID) {
1533:       PetscFV fv = (PetscFV) obj;

1535:       PetscFVGetQuadrature(fv, &quad);
1536:       PetscFVGetNumComponents(fv, &Nc);
1537:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1538:     numComponents += Nc;
1539:   }
1540:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1541:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1542:   PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1543:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1544:   for (c = cStart; c < cEnd; ++c) {
1545:     PetscScalar *x = NULL;
1546:     PetscScalar  elemDiff = 0.0;
1547:     PetscInt     qc = 0;

1549:     DMPlexComputeCellGeometryFEM(dm, c, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1550:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

1552:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1553:       PetscObject  obj;
1554:       PetscClassId id;
1555:       void * const ctx = ctxs ? ctxs[field] : NULL;
1556:       PetscInt     Nb, Nc, q, fc;

1558:       DMGetField(dm, field, NULL, &obj);
1559:       PetscObjectGetClassId(obj, &id);
1560:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1561:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1562:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1563:       if (funcs[field]) {
1564:         for (q = 0; q < Nq; ++q) {
1565:           PetscFEGeom qgeom;

1567:           qgeom.dimEmbed = fegeom.dimEmbed;
1568:           qgeom.J        = &fegeom.J[q*coordDim*coordDim];
1569:           qgeom.invJ     = &fegeom.invJ[q*coordDim*coordDim];
1570:           qgeom.detJ     = &fegeom.detJ[q];
1571:           if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)fegeom.detJ[q], c, q);
1572:           (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx);
1573:           if (ierr) {
1574:             PetscErrorCode ierr2;
1575:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1576:             ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1577:             ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1578:             
1579:           }
1580:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1581:           else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1582:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1583:           for (fc = 0; fc < Nc; ++fc) {
1584:             const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1585:             elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1586:           }
1587:         }
1588:       }
1589:       fieldOffset += Nb;
1590:       qc          += Nc;
1591:     }
1592:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1593:     VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1594:   }
1595:   PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1596:   DMRestoreLocalVector(dm, &localX);
1597:   VecSqrtAbs(D);
1598:   return(0);
1599: }

1601: /*@C
1602:   DMPlexComputeGradientClementInterpolant - This function computes the L2 projection of the cellwise gradient of a function u onto P1, and stores it in a Vec.

1604:   Collective on dm

1606:   Input Parameters:
1607: + dm - The DM
1608: - LocX  - The coefficient vector u_h

1610:   Output Parameter:
1611: . locC - A Vec which holds the Clement interpolant of the gradient

1613:   Notes:
1614:     Add citation to (Clement, 1975) and definition of the interpolant
1615:   \nabla u_h(v_i) = \sum_{T_i \in support(v_i)} |T_i| \nabla u_h(T_i) / \sum_{T_i \in support(v_i)} |T_i| where |T_i| is the cell volume

1617:   Level: developer

1619: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1620: @*/
1621: PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC)
1622: {
1623:   DM_Plex         *mesh  = (DM_Plex *) dm->data;
1624:   PetscInt         debug = mesh->printFEM;
1625:   DM               dmC;
1626:   PetscSection     section;
1627:   PetscQuadrature  quad;
1628:   PetscScalar     *interpolant, *gradsum;
1629:   PetscFEGeom      fegeom;
1630:   PetscReal       *coords;
1631:   const PetscReal *quadPoints, *quadWeights;
1632:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, vStart, vEnd, v, field, fieldOffset;
1633:   PetscErrorCode   ierr;

1636:   VecGetDM(locC, &dmC);
1637:   VecSet(locC, 0.0);
1638:   DMGetDimension(dm, &dim);
1639:   DMGetCoordinateDim(dm, &coordDim);
1640:   fegeom.dimEmbed = coordDim;
1641:   DMGetLocalSection(dm, &section);
1642:   PetscSectionGetNumFields(section, &numFields);
1643:   for (field = 0; field < numFields; ++field) {
1644:     PetscObject  obj;
1645:     PetscClassId id;
1646:     PetscInt     Nc;

1648:     DMGetField(dm, field, NULL, &obj);
1649:     PetscObjectGetClassId(obj, &id);
1650:     if (id == PETSCFE_CLASSID) {
1651:       PetscFE fe = (PetscFE) obj;

1653:       PetscFEGetQuadrature(fe, &quad);
1654:       PetscFEGetNumComponents(fe, &Nc);
1655:     } else if (id == PETSCFV_CLASSID) {
1656:       PetscFV fv = (PetscFV) obj;

1658:       PetscFVGetQuadrature(fv, &quad);
1659:       PetscFVGetNumComponents(fv, &Nc);
1660:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1661:     numComponents += Nc;
1662:   }
1663:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1664:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1665:   PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1666:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1667:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1668:   for (v = vStart; v < vEnd; ++v) {
1669:     PetscScalar volsum = 0.0;
1670:     PetscInt   *star = NULL;
1671:     PetscInt    starSize, st, d, fc;

1673:     PetscArrayzero(gradsum, coordDim*numComponents);
1674:     DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1675:     for (st = 0; st < starSize*2; st += 2) {
1676:       const PetscInt cell = star[st];
1677:       PetscScalar   *grad = &gradsum[coordDim*numComponents];
1678:       PetscScalar   *x    = NULL;
1679:       PetscReal      vol  = 0.0;

1681:       if ((cell < cStart) || (cell >= cEnd)) continue;
1682:       DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1683:       DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
1684:       for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1685:         PetscObject  obj;
1686:         PetscClassId id;
1687:         PetscInt     Nb, Nc, q, qc = 0;

1689:         PetscArrayzero(grad, coordDim*numComponents);
1690:         DMGetField(dm, field, NULL, &obj);
1691:         PetscObjectGetClassId(obj, &id);
1692:         if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1693:         else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1694:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1695:         for (q = 0; q < Nq; ++q) {
1696:           PetscFEGeom qgeom;

1698:           qgeom.dimEmbed = fegeom.dimEmbed;
1699:           qgeom.J        = &fegeom.J[q*coordDim*coordDim];
1700:           qgeom.invJ     = &fegeom.invJ[q*coordDim*coordDim];
1701:           qgeom.detJ     = &fegeom.detJ[q];
1702:           if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", (double)fegeom.detJ[q], cell, q);
1703:           if (ierr) {
1704:             PetscErrorCode ierr2;
1705:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
1706:             ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
1707:             ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1708:             
1709:           }
1710:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1711:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", field);
1712:           for (fc = 0; fc < Nc; ++fc) {
1713:             const PetscReal wt = quadWeights[q*qNc+qc+fc];

1715:             for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q];
1716:           }
1717:           vol += quadWeights[q*qNc]*fegeom.detJ[q];
1718:         }
1719:         fieldOffset += Nb;
1720:         qc          += Nc;
1721:       }
1722:       DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
1723:       for (fc = 0; fc < numComponents; ++fc) {
1724:         for (d = 0; d < coordDim; ++d) {
1725:           gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
1726:         }
1727:       }
1728:       volsum += vol;
1729:       if (debug) {
1730:         PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
1731:         for (fc = 0; fc < numComponents; ++fc) {
1732:           for (d = 0; d < coordDim; ++d) {
1733:             if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
1734:             PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
1735:           }
1736:         }
1737:         PetscPrintf(PETSC_COMM_SELF, "]\n");
1738:       }
1739:     }
1740:     for (fc = 0; fc < numComponents; ++fc) {
1741:       for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum;
1742:     }
1743:     DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1744:     DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
1745:   }
1746:   PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1747:   return(0);
1748: }

1750: static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user)
1751: {
1752:   DM                 dmAux = NULL;
1753:   PetscDS            prob,    probAux = NULL;
1754:   PetscSection       section, sectionAux;
1755:   Vec                locX,    locA;
1756:   PetscInt           dim, numCells = cEnd - cStart, c, f;
1757:   PetscBool          useFVM = PETSC_FALSE;
1758:   /* DS */
1759:   PetscInt           Nf,    totDim,    *uOff, *uOff_x, numConstants;
1760:   PetscInt           NfAux, totDimAux, *aOff;
1761:   PetscScalar       *u, *a;
1762:   const PetscScalar *constants;
1763:   /* Geometry */
1764:   PetscFEGeom       *cgeomFEM;
1765:   DM                 dmGrad;
1766:   PetscQuadrature    affineQuad = NULL;
1767:   Vec                cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
1768:   PetscFVCellGeom   *cgeomFVM;
1769:   const PetscScalar *lgrad;
1770:   PetscInt           maxDegree;
1771:   DMField            coordField;
1772:   IS                 cellIS;
1773:   PetscErrorCode     ierr;

1776:   DMGetDS(dm, &prob);
1777:   DMGetDimension(dm, &dim);
1778:   DMGetLocalSection(dm, &section);
1779:   PetscSectionGetNumFields(section, &Nf);
1780:   /* Determine which discretizations we have */
1781:   for (f = 0; f < Nf; ++f) {
1782:     PetscObject  obj;
1783:     PetscClassId id;

1785:     PetscDSGetDiscretization(prob, f, &obj);
1786:     PetscObjectGetClassId(obj, &id);
1787:     if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE;
1788:   }
1789:   /* Get local solution with boundary values */
1790:   DMGetLocalVector(dm, &locX);
1791:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1792:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1793:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1794:   /* Read DS information */
1795:   PetscDSGetTotalDimension(prob, &totDim);
1796:   PetscDSGetComponentOffsets(prob, &uOff);
1797:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1798:   ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);
1799:   PetscDSGetConstants(prob, &numConstants, &constants);
1800:   /* Read Auxiliary DS information */
1801:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1802:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1803:   if (dmAux) {
1804:     DMGetDS(dmAux, &probAux);
1805:     PetscDSGetNumFields(probAux, &NfAux);
1806:     DMGetLocalSection(dmAux, &sectionAux);
1807:     PetscDSGetTotalDimension(probAux, &totDimAux);
1808:     PetscDSGetComponentOffsets(probAux, &aOff);
1809:   }
1810:   /* Allocate data  arrays */
1811:   PetscCalloc1(numCells*totDim, &u);
1812:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1813:   /* Read out geometry */
1814:   DMGetCoordinateField(dm,&coordField);
1815:   DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
1816:   if (maxDegree <= 1) {
1817:     DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
1818:     if (affineQuad) {
1819:       DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);
1820:     }
1821:   }
1822:   if (useFVM) {
1823:     PetscFV   fv = NULL;
1824:     Vec       grad;
1825:     PetscInt  fStart, fEnd;
1826:     PetscBool compGrad;

1828:     for (f = 0; f < Nf; ++f) {
1829:       PetscObject  obj;
1830:       PetscClassId id;

1832:       PetscDSGetDiscretization(prob, f, &obj);
1833:       PetscObjectGetClassId(obj, &id);
1834:       if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;}
1835:     }
1836:     PetscFVGetComputeGradients(fv, &compGrad);
1837:     PetscFVSetComputeGradients(fv, PETSC_TRUE);
1838:     DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1839:     DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1840:     PetscFVSetComputeGradients(fv, compGrad);
1841:     VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1842:     /* Reconstruct and limit cell gradients */
1843:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1844:     DMGetGlobalVector(dmGrad, &grad);
1845:     DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1846:     /* Communicate gradient values */
1847:     DMGetLocalVector(dmGrad, &locGrad);
1848:     DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1849:     DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1850:     DMRestoreGlobalVector(dmGrad, &grad);
1851:     /* Handle non-essential (e.g. outflow) boundary values */
1852:     DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1853:     VecGetArrayRead(locGrad, &lgrad);
1854:   }
1855:   /* Read out data from inputs */
1856:   for (c = cStart; c < cEnd; ++c) {
1857:     PetscScalar *x = NULL;
1858:     PetscInt     i;

1860:     DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1861:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1862:     DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1863:     if (dmAux) {
1864:       DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1865:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1866:       DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1867:     }
1868:   }
1869:   /* Do integration for each field */
1870:   for (f = 0; f < Nf; ++f) {
1871:     PetscObject  obj;
1872:     PetscClassId id;
1873:     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1875:     PetscDSGetDiscretization(prob, f, &obj);
1876:     PetscObjectGetClassId(obj, &id);
1877:     if (id == PETSCFE_CLASSID) {
1878:       PetscFE         fe = (PetscFE) obj;
1879:       PetscQuadrature q;
1880:       PetscFEGeom     *chunkGeom = NULL;
1881:       PetscInt        Nq, Nb;

1883:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1884:       PetscFEGetQuadrature(fe, &q);
1885:       PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1886:       PetscFEGetDimension(fe, &Nb);
1887:       blockSize = Nb*Nq;
1888:       batchSize = numBlocks * blockSize;
1889:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1890:       numChunks = numCells / (numBatches*batchSize);
1891:       Ne        = numChunks*numBatches*batchSize;
1892:       Nr        = numCells % (numBatches*batchSize);
1893:       offset    = numCells - Nr;
1894:       if (!affineQuad) {
1895:         DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);
1896:       }
1897:       PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
1898:       PetscFEIntegrate(prob, f, Ne, chunkGeom, u, probAux, a, cintegral);
1899:       PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);
1900:       PetscFEIntegrate(prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);
1901:       PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);
1902:       if (!affineQuad) {
1903:         PetscFEGeomDestroy(&cgeomFEM);
1904:       }
1905:     } else if (id == PETSCFV_CLASSID) {
1906:       PetscInt       foff;
1907:       PetscPointFunc obj_func;
1908:       PetscScalar    lint;

1910:       PetscDSGetObjective(prob, f, &obj_func);
1911:       PetscDSGetFieldOffset(prob, f, &foff);
1912:       if (obj_func) {
1913:         for (c = 0; c < numCells; ++c) {
1914:           PetscScalar *u_x;

1916:           DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1917:           obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, &a[totDimAux*c], NULL, NULL, 0.0, cgeomFVM[c].centroid, numConstants, constants, &lint);
1918:           cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1919:         }
1920:       }
1921:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
1922:   }
1923:   /* Cleanup data arrays */
1924:   if (useFVM) {
1925:     VecRestoreArrayRead(locGrad, &lgrad);
1926:     VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1927:     DMRestoreLocalVector(dmGrad, &locGrad);
1928:     VecDestroy(&faceGeometryFVM);
1929:     VecDestroy(&cellGeometryFVM);
1930:     DMDestroy(&dmGrad);
1931:   }
1932:   if (dmAux) {PetscFree(a);}
1933:   PetscFree(u);
1934:   /* Cleanup */
1935:   if (affineQuad) {
1936:     PetscFEGeomDestroy(&cgeomFEM);
1937:   }
1938:   PetscQuadratureDestroy(&affineQuad);
1939:   ISDestroy(&cellIS);
1940:   DMRestoreLocalVector(dm, &locX);
1941:   return(0);
1942: }

1944: /*@
1945:   DMPlexComputeIntegralFEM - Form the integral over the domain from the global input X using pointwise functions specified by the user

1947:   Input Parameters:
1948: + dm - The mesh
1949: . X  - Global input vector
1950: - user - The user context

1952:   Output Parameter:
1953: . integral - Integral for each field

1955:   Level: developer

1957: .seealso: DMPlexComputeResidualFEM()
1958: @*/
1959: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
1960: {
1961:   DM_Plex       *mesh = (DM_Plex *) dm->data;
1962:   PetscScalar   *cintegral, *lintegral;
1963:   PetscInt       Nf, f, cellHeight, cStart, cEnd, cell;

1970:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1971:   DMGetNumFields(dm, &Nf);
1972:   DMPlexGetVTKCellHeight(dm, &cellHeight);
1973:   DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd);
1974:   /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1975:   PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);
1976:   DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1977:   /* Sum up values */
1978:   for (cell = cStart; cell < cEnd; ++cell) {
1979:     const PetscInt c = cell - cStart;

1981:     if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1982:     for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f];
1983:   }
1984:   MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1985:   if (mesh->printFEM) {
1986:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");
1987:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));}
1988:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1989:   }
1990:   PetscFree2(lintegral, cintegral);
1991:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1992:   return(0);
1993: }

1995: /*@
1996:   DMPlexComputeCellwiseIntegralFEM - Form the vector of cellwise integrals F from the global input X using pointwise functions specified by the user

1998:   Input Parameters:
1999: + dm - The mesh
2000: . X  - Global input vector
2001: - user - The user context

2003:   Output Parameter:
2004: . integral - Cellwise integrals for each field

2006:   Level: developer

2008: .seealso: DMPlexComputeResidualFEM()
2009: @*/
2010: PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user)
2011: {
2012:   DM_Plex       *mesh = (DM_Plex *) dm->data;
2013:   DM             dmF;
2014:   PetscSection   sectionF;
2015:   PetscScalar   *cintegral, *af;
2016:   PetscInt       Nf, f, cellHeight, cStart, cEnd, cell;

2023:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
2024:   DMGetNumFields(dm, &Nf);
2025:   DMPlexGetVTKCellHeight(dm, &cellHeight);
2026:   DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd);
2027:   /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
2028:   PetscCalloc1((cEnd-cStart)*Nf, &cintegral);
2029:   DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
2030:   /* Put values in F*/
2031:   VecGetDM(F, &dmF);
2032:   DMGetLocalSection(dmF, &sectionF);
2033:   VecGetArray(F, &af);
2034:   for (cell = cStart; cell < cEnd; ++cell) {
2035:     const PetscInt c = cell - cStart;
2036:     PetscInt       dof, off;

2038:     if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
2039:     PetscSectionGetDof(sectionF, cell, &dof);
2040:     PetscSectionGetOffset(sectionF, cell, &off);
2041:     if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf);
2042:     for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f];
2043:   }
2044:   VecRestoreArray(F, &af);
2045:   PetscFree(cintegral);
2046:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
2047:   return(0);
2048: }

2050: static PetscErrorCode DMPlexComputeBdIntegral_Internal(DM dm, Vec locX, IS pointIS,
2051:                                                        void (*func)(PetscInt, PetscInt, PetscInt,
2052:                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2053:                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2054:                                                                     PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
2055:                                                        PetscScalar *fintegral, void *user)
2056: {
2057:   DM                 plex = NULL, plexA = NULL;
2058:   DMEnclosureType    encAux;
2059:   PetscDS            prob, probAux = NULL;
2060:   PetscSection       section, sectionAux = NULL;
2061:   Vec                locA = NULL;
2062:   DMField            coordField;
2063:   PetscInt           Nf,        totDim,        *uOff, *uOff_x;
2064:   PetscInt           NfAux = 0, totDimAux = 0, *aOff = NULL;
2065:   PetscScalar       *u, *a = NULL;
2066:   const PetscScalar *constants;
2067:   PetscInt           numConstants, f;
2068:   PetscErrorCode     ierr;

2071:   DMGetCoordinateField(dm, &coordField);
2072:   DMConvert(dm, DMPLEX, &plex);
2073:   DMGetDS(dm, &prob);
2074:   DMGetLocalSection(dm, &section);
2075:   PetscSectionGetNumFields(section, &Nf);
2076:   /* Determine which discretizations we have */
2077:   for (f = 0; f < Nf; ++f) {
2078:     PetscObject  obj;
2079:     PetscClassId id;

2081:     PetscDSGetDiscretization(prob, f, &obj);
2082:     PetscObjectGetClassId(obj, &id);
2083:     if (id == PETSCFV_CLASSID) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not supported for FVM (field %D)", f);
2084:   }
2085:   /* Read DS information */
2086:   PetscDSGetTotalDimension(prob, &totDim);
2087:   PetscDSGetComponentOffsets(prob, &uOff);
2088:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
2089:   PetscDSGetConstants(prob, &numConstants, &constants);
2090:   /* Read Auxiliary DS information */
2091:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
2092:   if (locA) {
2093:     DM dmAux;

2095:     VecGetDM(locA, &dmAux);
2096:     DMGetEnclosureRelation(dmAux, dm, &encAux);
2097:     DMConvert(dmAux, DMPLEX, &plexA);
2098:     DMGetDS(dmAux, &probAux);
2099:     PetscDSGetNumFields(probAux, &NfAux);
2100:     DMGetLocalSection(dmAux, &sectionAux);
2101:     PetscDSGetTotalDimension(probAux, &totDimAux);
2102:     PetscDSGetComponentOffsets(probAux, &aOff);
2103:   }
2104:   /* Integrate over points */
2105:   {
2106:     PetscFEGeom    *fgeom, *chunkGeom = NULL;
2107:     PetscInt        maxDegree;
2108:     PetscQuadrature qGeom = NULL;
2109:     const PetscInt *points;
2110:     PetscInt        numFaces, face, Nq, field;
2111:     PetscInt        numChunks, chunkSize, chunk, Nr, offset;

2113:     ISGetLocalSize(pointIS, &numFaces);
2114:     ISGetIndices(pointIS, &points);
2115:     PetscCalloc2(numFaces*totDim, &u, locA ? numFaces*totDimAux : 0, &a);
2116:     DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
2117:     for (field = 0; field < Nf; ++field) {
2118:       PetscFE fe;

2120:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
2121:       if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, pointIS, &qGeom);}
2122:       if (!qGeom) {
2123:         PetscFEGetFaceQuadrature(fe, &qGeom);
2124:         PetscObjectReference((PetscObject) qGeom);
2125:       }
2126:       PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
2127:       DMPlexGetFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
2128:       for (face = 0; face < numFaces; ++face) {
2129:         const PetscInt point = points[face], *support, *cone;
2130:         PetscScalar    *x    = NULL;
2131:         PetscInt       i, coneSize, faceLoc;

2133:         DMPlexGetSupport(dm, point, &support);
2134:         DMPlexGetConeSize(dm, support[0], &coneSize);
2135:         DMPlexGetCone(dm, support[0], &cone);
2136:         for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break;
2137:         if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of support[0] %D", face, support[0]);
2138:         fgeom->face[face][0] = faceLoc;
2139:         DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);
2140:         for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
2141:         DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);
2142:         if (locA) {
2143:           PetscInt subp;
2144:           DMGetEnclosurePoint(plexA, dm, encAux, support[0], &subp);
2145:           DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);
2146:           for (i = 0; i < totDimAux; ++i) a[f*totDimAux+i] = x[i];
2147:           DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);
2148:         }
2149:       }
2150:       /* Get blocking */
2151:       {
2152:         PetscQuadrature q;
2153:         PetscInt        numBatches, batchSize, numBlocks, blockSize;
2154:         PetscInt        Nq, Nb;

2156:         PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2157:         PetscFEGetQuadrature(fe, &q);
2158:         PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
2159:         PetscFEGetDimension(fe, &Nb);
2160:         blockSize = Nb*Nq;
2161:         batchSize = numBlocks * blockSize;
2162:         chunkSize = numBatches*batchSize;
2163:         PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2164:         numChunks = numFaces / chunkSize;
2165:         Nr        = numFaces % chunkSize;
2166:         offset    = numFaces - Nr;
2167:       }
2168:       /* Do integration for each field */
2169:       for (chunk = 0; chunk < numChunks; ++chunk) {
2170:         PetscFEGeomGetChunk(fgeom, chunk*chunkSize, (chunk+1)*chunkSize, &chunkGeom);
2171:         PetscFEIntegrateBd(prob, field, func, chunkSize, chunkGeom, u, probAux, a, fintegral);
2172:         PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);
2173:       }
2174:       PetscFEGeomGetChunk(fgeom, offset, numFaces, &chunkGeom);
2175:       PetscFEIntegrateBd(prob, field, func, Nr, chunkGeom, &u[offset*totDim], probAux, a ? &a[offset*totDimAux] : NULL, &fintegral[offset*Nf]);
2176:       PetscFEGeomRestoreChunk(fgeom, offset, numFaces, &chunkGeom);
2177:       /* Cleanup data arrays */
2178:       DMPlexRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
2179:       PetscQuadratureDestroy(&qGeom);
2180:       PetscFree2(u, a);
2181:       ISRestoreIndices(pointIS, &points);
2182:     }
2183:   }
2184:   if (plex)  {DMDestroy(&plex);}
2185:   if (plexA) {DMDestroy(&plexA);}
2186:   return(0);
2187: }

2189: /*@
2190:   DMPlexComputeBdIntegral - Form the integral over the specified boundary from the global input X using pointwise functions specified by the user

2192:   Input Parameters:
2193: + dm      - The mesh
2194: . X       - Global input vector
2195: . label   - The boundary DMLabel
2196: . numVals - The number of label values to use, or PETSC_DETERMINE for all values
2197: . vals    - The label values to use, or PETSC_NULL for all values
2198: . func    = The function to integrate along the boundary
2199: - user    - The user context

2201:   Output Parameter:
2202: . integral - Integral for each field

2204:   Level: developer

2206: .seealso: DMPlexComputeIntegralFEM(), DMPlexComputeBdResidualFEM()
2207: @*/
2208: PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt numVals, const PetscInt vals[],
2209:                                        void (*func)(PetscInt, PetscInt, PetscInt,
2210:                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2211:                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2212:                                                     PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
2213:                                        PetscScalar *integral, void *user)
2214: {
2215:   Vec            locX;
2216:   PetscSection   section;
2217:   DMLabel        depthLabel;
2218:   IS             facetIS;
2219:   PetscInt       dim, Nf, f, v;

2228:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
2229:   DMPlexGetDepthLabel(dm, &depthLabel);
2230:   DMGetDimension(dm, &dim);
2231:   DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);
2232:   DMGetLocalSection(dm, &section);
2233:   PetscSectionGetNumFields(section, &Nf);
2234:   /* Get local solution with boundary values */
2235:   DMGetLocalVector(dm, &locX);
2236:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
2237:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
2238:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
2239:   /* Loop over label values */
2240:   PetscArrayzero(integral, Nf);
2241:   for (v = 0; v < numVals; ++v) {
2242:     IS           pointIS;
2243:     PetscInt     numFaces, face;
2244:     PetscScalar *fintegral;

2246:     DMLabelGetStratumIS(label, vals[v], &pointIS);
2247:     if (!pointIS) continue; /* No points with that id on this process */
2248:     {
2249:       IS isectIS;

2251:       /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */
2252:       ISIntersect_Caching_Internal(facetIS, pointIS, &isectIS);
2253:       ISDestroy(&pointIS);
2254:       pointIS = isectIS;
2255:     }
2256:     ISGetLocalSize(pointIS, &numFaces);
2257:     PetscCalloc1(numFaces*Nf, &fintegral);
2258:     DMPlexComputeBdIntegral_Internal(dm, locX, pointIS, func, fintegral, user);
2259:     /* Sum point contributions into integral */
2260:     for (f = 0; f < Nf; ++f) for (face = 0; face < numFaces; ++face) integral[f] += fintegral[face*Nf+f];
2261:     PetscFree(fintegral);
2262:     ISDestroy(&pointIS);
2263:   }
2264:   DMRestoreLocalVector(dm, &locX);
2265:   ISDestroy(&facetIS);
2266:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
2267:   return(0);
2268: }

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

2273:   Input Parameters:
2274: + dmc  - The coarse mesh
2275: . dmf  - The fine mesh
2276: . isRefined - Flag indicating regular refinement, rather than the same topology
2277: - user - The user context

2279:   Output Parameter:
2280: . In  - The interpolation matrix

2282:   Level: developer

2284: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2285: @*/
2286: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, PetscBool isRefined, Mat In, void *user)
2287: {
2288:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
2289:   const char       *name  = "Interpolator";
2290:   PetscDS           cds, rds;
2291:   PetscFE          *feRef;
2292:   PetscFV          *fvRef;
2293:   PetscSection      fsection, fglobalSection;
2294:   PetscSection      csection, cglobalSection;
2295:   PetscScalar      *elemMat;
2296:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c;
2297:   PetscInt          cTotDim, rTotDim = 0;
2298:   PetscErrorCode    ierr;

2301:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2302:   DMGetDimension(dmf, &dim);
2303:   DMGetLocalSection(dmf, &fsection);
2304:   DMGetGlobalSection(dmf, &fglobalSection);
2305:   DMGetLocalSection(dmc, &csection);
2306:   DMGetGlobalSection(dmc, &cglobalSection);
2307:   PetscSectionGetNumFields(fsection, &Nf);
2308:   DMPlexGetSimplexOrBoxCells(dmc, 0, &cStart, &cEnd);
2309:   DMGetDS(dmc, &cds);
2310:   DMGetDS(dmf, &rds);
2311:   PetscCalloc2(Nf, &feRef, Nf, &fvRef);
2312:   for (f = 0; f < Nf; ++f) {
2313:     PetscObject  obj;
2314:     PetscClassId id;
2315:     PetscInt     rNb = 0, Nc = 0;

2317:     PetscDSGetDiscretization(rds, f, &obj);
2318:     PetscObjectGetClassId(obj, &id);
2319:     if (id == PETSCFE_CLASSID) {
2320:       PetscFE fe = (PetscFE) obj;

2322:       if (isRefined) {
2323:         PetscFERefine(fe, &feRef[f]);
2324:       } else {
2325:         PetscObjectReference((PetscObject) fe);
2326:         feRef[f] = fe;
2327:       }
2328:       PetscFEGetDimension(feRef[f], &rNb);
2329:       PetscFEGetNumComponents(fe, &Nc);
2330:     } else if (id == PETSCFV_CLASSID) {
2331:       PetscFV        fv = (PetscFV) obj;
2332:       PetscDualSpace Q;

2334:       if (isRefined) {
2335:         PetscFVRefine(fv, &fvRef[f]);
2336:       } else {
2337:         PetscObjectReference((PetscObject) fv);
2338:         fvRef[f] = fv;
2339:       }
2340:       PetscFVGetDualSpace(fvRef[f], &Q);
2341:       PetscDualSpaceGetDimension(Q, &rNb);
2342:       PetscFVGetNumComponents(fv, &Nc);
2343:     }
2344:     rTotDim += rNb;
2345:   }
2346:   PetscDSGetTotalDimension(cds, &cTotDim);
2347:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
2348:   PetscArrayzero(elemMat, rTotDim*cTotDim);
2349:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
2350:     PetscDualSpace   Qref;
2351:     PetscQuadrature  f;
2352:     const PetscReal *qpoints, *qweights;
2353:     PetscReal       *points;
2354:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

2356:     /* Compose points from all dual basis functionals */
2357:     if (feRef[fieldI]) {
2358:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
2359:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
2360:     } else {
2361:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
2362:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
2363:     }
2364:     PetscDualSpaceGetDimension(Qref, &fpdim);
2365:     for (i = 0; i < fpdim; ++i) {
2366:       PetscDualSpaceGetFunctional(Qref, i, &f);
2367:       PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
2368:       npoints += Np;
2369:     }
2370:     PetscMalloc1(npoints*dim,&points);
2371:     for (i = 0, k = 0; i < fpdim; ++i) {
2372:       PetscDualSpaceGetFunctional(Qref, i, &f);
2373:       PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2374:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
2375:     }

2377:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
2378:       PetscObject  obj;
2379:       PetscClassId id;
2380:       PetscInt     NcJ = 0, cpdim = 0, j, qNc;

2382:       PetscDSGetDiscretization(cds, fieldJ, &obj);
2383:       PetscObjectGetClassId(obj, &id);
2384:       if (id == PETSCFE_CLASSID) {
2385:         PetscFE           fe = (PetscFE) obj;
2386:         PetscTabulation T  = NULL;

2388:         /* Evaluate basis at points */
2389:         PetscFEGetNumComponents(fe, &NcJ);
2390:         PetscFEGetDimension(fe, &cpdim);
2391:         /* For now, fields only interpolate themselves */
2392:         if (fieldI == fieldJ) {
2393:           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);
2394:           PetscFECreateTabulation(fe, 1, npoints, points, 0, &T);
2395:           for (i = 0, k = 0; i < fpdim; ++i) {
2396:             PetscDualSpaceGetFunctional(Qref, i, &f);
2397:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
2398:             if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ);
2399:             for (p = 0; p < Np; ++p, ++k) {
2400:               for (j = 0; j < cpdim; ++j) {
2401:                 /*
2402:                    cTotDim:            Total columns in element interpolation matrix, sum of number of dual basis functionals in each field
2403:                    offsetI, offsetJ:   Offsets into the larger element interpolation matrix for different fields
2404:                    fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals
2405:                    qNC, Nc, Ncj, c:    Number of components in this field
2406:                    Np, p:              Number of quad points in the fine grid functional i
2407:                    k:                  i*Np + p, overall point number for the interpolation
2408:                 */
2409:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += T->T[0][k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
2410:               }
2411:             }
2412:           }
2413:           PetscTabulationDestroy(&T);
2414:         }
2415:       } else if (id == PETSCFV_CLASSID) {
2416:         PetscFV        fv = (PetscFV) obj;

2418:         /* Evaluate constant function at points */
2419:         PetscFVGetNumComponents(fv, &NcJ);
2420:         cpdim = 1;
2421:         /* For now, fields only interpolate themselves */
2422:         if (fieldI == fieldJ) {
2423:           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);
2424:           for (i = 0, k = 0; i < fpdim; ++i) {
2425:             PetscDualSpaceGetFunctional(Qref, i, &f);
2426:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
2427:             if (qNc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, NcJ);
2428:             for (p = 0; p < Np; ++p, ++k) {
2429:               for (j = 0; j < cpdim; ++j) {
2430:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += 1.0*qweights[p*qNc+c];
2431:               }
2432:             }
2433:           }
2434:         }
2435:       }
2436:       offsetJ += cpdim;
2437:     }
2438:     offsetI += fpdim;
2439:     PetscFree(points);
2440:   }
2441:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
2442:   /* Preallocate matrix */
2443:   {
2444:     Mat          preallocator;
2445:     PetscScalar *vals;
2446:     PetscInt    *cellCIndices, *cellFIndices;
2447:     PetscInt     locRows, locCols, cell;

2449:     MatGetLocalSize(In, &locRows, &locCols);
2450:     MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
2451:     MatSetType(preallocator, MATPREALLOCATOR);
2452:     MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
2453:     MatSetUp(preallocator);
2454:     PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
2455:     for (cell = cStart; cell < cEnd; ++cell) {
2456:       if (isRefined) {
2457:         DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
2458:         MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
2459:       } else {
2460:         DMPlexMatSetClosureGeneral(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, preallocator, cell, vals, INSERT_VALUES);
2461:       }
2462:     }
2463:     PetscFree3(vals,cellCIndices,cellFIndices);
2464:     MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
2465:     MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
2466:     MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
2467:     MatDestroy(&preallocator);
2468:   }
2469:   /* Fill matrix */
2470:   MatZeroEntries(In);
2471:   for (c = cStart; c < cEnd; ++c) {
2472:     if (isRefined) {
2473:       DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
2474:     } else {
2475:       DMPlexMatSetClosureGeneral(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
2476:     }
2477:   }
2478:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
2479:   PetscFree2(feRef,fvRef);
2480:   PetscFree(elemMat);
2481:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2482:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2483:   if (mesh->printFEM) {
2484:     PetscPrintf(PetscObjectComm((PetscObject)In), "%s:\n", name);
2485:     MatChop(In, 1.0e-10);
2486:     MatView(In, NULL);
2487:   }
2488:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2489:   return(0);
2490: }

2492: PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user)
2493: {
2494:   SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness");
2495: }

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

2500:   Input Parameters:
2501: + dmf  - The fine mesh
2502: . dmc  - The coarse mesh
2503: - user - The user context

2505:   Output Parameter:
2506: . In  - The interpolation matrix

2508:   Level: developer

2510: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2511: @*/
2512: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
2513: {
2514:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
2515:   const char    *name = "Interpolator";
2516:   PetscDS        prob;
2517:   PetscSection   fsection, csection, globalFSection, globalCSection;
2518:   PetscHSetIJ    ht;
2519:   PetscLayout    rLayout;
2520:   PetscInt      *dnz, *onz;
2521:   PetscInt       locRows, rStart, rEnd;
2522:   PetscReal     *x, *v0, *J, *invJ, detJ;
2523:   PetscReal     *v0c, *Jc, *invJc, detJc;
2524:   PetscScalar   *elemMat;
2525:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

2529:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2530:   DMGetCoordinateDim(dmc, &dim);
2531:   DMGetDS(dmc, &prob);
2532:   PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);
2533:   PetscDSGetNumFields(prob, &Nf);
2534:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2535:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2536:   DMGetLocalSection(dmf, &fsection);
2537:   DMGetGlobalSection(dmf, &globalFSection);
2538:   DMGetLocalSection(dmc, &csection);
2539:   DMGetGlobalSection(dmc, &globalCSection);
2540:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2541:   PetscDSGetTotalDimension(prob, &totDim);
2542:   PetscMalloc1(totDim, &elemMat);

2544:   MatGetLocalSize(In, &locRows, NULL);
2545:   PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
2546:   PetscLayoutSetLocalSize(rLayout, locRows);
2547:   PetscLayoutSetBlockSize(rLayout, 1);
2548:   PetscLayoutSetUp(rLayout);
2549:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2550:   PetscLayoutDestroy(&rLayout);
2551:   PetscCalloc2(locRows,&dnz,locRows,&onz);
2552:   PetscHSetIJCreate(&ht);
2553:   for (field = 0; field < Nf; ++field) {
2554:     PetscObject      obj;
2555:     PetscClassId     id;
2556:     PetscDualSpace   Q = NULL;
2557:     PetscQuadrature  f;
2558:     const PetscReal *qpoints;
2559:     PetscInt         Nc, Np, fpdim, i, d;

2561:     PetscDSGetDiscretization(prob, field, &obj);
2562:     PetscObjectGetClassId(obj, &id);
2563:     if (id == PETSCFE_CLASSID) {
2564:       PetscFE fe = (PetscFE) obj;

2566:       PetscFEGetDualSpace(fe, &Q);
2567:       PetscFEGetNumComponents(fe, &Nc);
2568:     } else if (id == PETSCFV_CLASSID) {
2569:       PetscFV fv = (PetscFV) obj;

2571:       PetscFVGetDualSpace(fv, &Q);
2572:       Nc   = 1;
2573:     }
2574:     PetscDualSpaceGetDimension(Q, &fpdim);
2575:     /* For each fine grid cell */
2576:     for (cell = cStart; cell < cEnd; ++cell) {
2577:       PetscInt *findices,   *cindices;
2578:       PetscInt  numFIndices, numCIndices;

2580:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2581:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2582:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %D != %D dual basis vecs", numFIndices, fpdim);
2583:       for (i = 0; i < fpdim; ++i) {
2584:         Vec             pointVec;
2585:         PetscScalar    *pV;
2586:         PetscSF         coarseCellSF = NULL;
2587:         const PetscSFNode *coarseCells;
2588:         PetscInt        numCoarseCells, q, c;

2590:         /* Get points from the dual basis functional quadrature */
2591:         PetscDualSpaceGetFunctional(Q, i, &f);
2592:         PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2593:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2594:         VecSetBlockSize(pointVec, dim);
2595:         VecGetArray(pointVec, &pV);
2596:         for (q = 0; q < Np; ++q) {
2597:           const PetscReal xi0[3] = {-1., -1., -1.};

2599:           /* Transform point to real space */
2600:           CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2601:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2602:         }
2603:         VecRestoreArray(pointVec, &pV);
2604:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2605:         /* OPT: Pack all quad points from fine cell */
2606:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2607:         PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2608:         /* Update preallocation info */
2609:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2610:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2611:         {
2612:           PetscHashIJKey key;
2613:           PetscBool      missing;

2615:           key.i = findices[i];
2616:           if (key.i >= 0) {
2617:             /* Get indices for coarse elements */
2618:             for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2619:               DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2620:               for (c = 0; c < numCIndices; ++c) {
2621:                 key.j = cindices[c];
2622:                 if (key.j < 0) continue;
2623:                 PetscHSetIJQueryAdd(ht, key, &missing);
2624:                 if (missing) {
2625:                   if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2626:                   else                                     ++onz[key.i-rStart];
2627:                 }
2628:               }
2629:               DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2630:             }
2631:           }
2632:         }
2633:         PetscSFDestroy(&coarseCellSF);
2634:         VecDestroy(&pointVec);
2635:       }
2636:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2637:     }
2638:   }
2639:   PetscHSetIJDestroy(&ht);
2640:   MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
2641:   MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2642:   PetscFree2(dnz,onz);
2643:   for (field = 0; field < Nf; ++field) {
2644:     PetscObject       obj;
2645:     PetscClassId      id;
2646:     PetscDualSpace    Q = NULL;
2647:     PetscTabulation T = NULL;
2648:     PetscQuadrature   f;
2649:     const PetscReal  *qpoints, *qweights;
2650:     PetscInt          Nc, qNc, Np, fpdim, i, d;

2652:     PetscDSGetDiscretization(prob, field, &obj);
2653:     PetscObjectGetClassId(obj, &id);
2654:     if (id == PETSCFE_CLASSID) {
2655:       PetscFE fe = (PetscFE) obj;

2657:       PetscFEGetDualSpace(fe, &Q);
2658:       PetscFEGetNumComponents(fe, &Nc);
2659:       PetscFECreateTabulation(fe, 1, 1, x, 0, &T);
2660:     } else if (id == PETSCFV_CLASSID) {
2661:       PetscFV fv = (PetscFV) obj;

2663:       PetscFVGetDualSpace(fv, &Q);
2664:       Nc   = 1;
2665:     } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %D",field);
2666:     PetscDualSpaceGetDimension(Q, &fpdim);
2667:     /* For each fine grid cell */
2668:     for (cell = cStart; cell < cEnd; ++cell) {
2669:       PetscInt *findices,   *cindices;
2670:       PetscInt  numFIndices, numCIndices;

2672:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2673:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2674:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %D != %D dual basis vecs", numFIndices, fpdim);
2675:       for (i = 0; i < fpdim; ++i) {
2676:         Vec             pointVec;
2677:         PetscScalar    *pV;
2678:         PetscSF         coarseCellSF = NULL;
2679:         const PetscSFNode *coarseCells;
2680:         PetscInt        numCoarseCells, cpdim, q, c, j;

2682:         /* Get points from the dual basis functional quadrature */
2683:         PetscDualSpaceGetFunctional(Q, i, &f);
2684:         PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
2685:         if (qNc != Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in quadrature %D does not match coarse field %D", qNc, Nc);
2686:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2687:         VecSetBlockSize(pointVec, dim);
2688:         VecGetArray(pointVec, &pV);
2689:         for (q = 0; q < Np; ++q) {
2690:           const PetscReal xi0[3] = {-1., -1., -1.};

2692:           /* Transform point to real space */
2693:           CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2694:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2695:         }
2696:         VecRestoreArray(pointVec, &pV);
2697:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2698:         /* OPT: Read this out from preallocation information */
2699:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2700:         /* Update preallocation info */
2701:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2702:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2703:         VecGetArray(pointVec, &pV);
2704:         for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2705:           PetscReal pVReal[3];
2706:           const PetscReal xi0[3] = {-1., -1., -1.};

2708:           DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2709:           /* Transform points from real space to coarse reference space */
2710:           DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
2711:           for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
2712:           CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x);

2714:           if (id == PETSCFE_CLASSID) {
2715:             PetscFE fe = (PetscFE) obj;

2717:             /* Evaluate coarse basis on contained point */
2718:             PetscFEGetDimension(fe, &cpdim);
2719:             PetscFEComputeTabulation(fe, 1, x, 0, T);
2720:             PetscArrayzero(elemMat, cpdim);
2721:             /* Get elemMat entries by multiplying by weight */
2722:             for (j = 0; j < cpdim; ++j) {
2723:               for (c = 0; c < Nc; ++c) elemMat[j] += T->T[0][j*Nc + c]*qweights[ccell*qNc + c];
2724:             }
2725:           } else {
2726:             cpdim = 1;
2727:             for (j = 0; j < cpdim; ++j) {
2728:               for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
2729:             }
2730:           }
2731:           /* Update interpolator */
2732:           if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2733:           if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2734:           MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
2735:           DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2736:         }
2737:         VecRestoreArray(pointVec, &pV);
2738:         PetscSFDestroy(&coarseCellSF);
2739:         VecDestroy(&pointVec);
2740:       }
2741:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2742:     }
2743:     if (id == PETSCFE_CLASSID) {PetscTabulationDestroy(&T);}
2744:   }
2745:   PetscFree3(v0,J,invJ);
2746:   PetscFree3(v0c,Jc,invJc);
2747:   PetscFree(elemMat);
2748:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2749:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2750:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2751:   return(0);
2752: }

2754: /*@
2755:   DMPlexComputeMassMatrixGeneral - Form the local portion of the mass matrix M from the coarse DM to a non-nested fine DM.

2757:   Input Parameters:
2758: + dmf  - The fine mesh
2759: . dmc  - The coarse mesh
2760: - user - The user context

2762:   Output Parameter:
2763: . mass  - The mass matrix

2765:   Level: developer

2767: .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2768: @*/
2769: PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user)
2770: {
2771:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
2772:   const char    *name = "Mass Matrix";
2773:   PetscDS        prob;
2774:   PetscSection   fsection, csection, globalFSection, globalCSection;
2775:   PetscHSetIJ    ht;
2776:   PetscLayout    rLayout;
2777:   PetscInt      *dnz, *onz;
2778:   PetscInt       locRows, rStart, rEnd;
2779:   PetscReal     *x, *v0, *J, *invJ, detJ;
2780:   PetscReal     *v0c, *Jc, *invJc, detJc;
2781:   PetscScalar   *elemMat;
2782:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

2786:   DMGetCoordinateDim(dmc, &dim);
2787:   DMGetDS(dmc, &prob);
2788:   PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);
2789:   PetscDSGetNumFields(prob, &Nf);
2790:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2791:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2792:   DMGetLocalSection(dmf, &fsection);
2793:   DMGetGlobalSection(dmf, &globalFSection);
2794:   DMGetLocalSection(dmc, &csection);
2795:   DMGetGlobalSection(dmc, &globalCSection);
2796:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2797:   PetscDSGetTotalDimension(prob, &totDim);
2798:   PetscMalloc1(totDim, &elemMat);

2800:   MatGetLocalSize(mass, &locRows, NULL);
2801:   PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);
2802:   PetscLayoutSetLocalSize(rLayout, locRows);
2803:   PetscLayoutSetBlockSize(rLayout, 1);
2804:   PetscLayoutSetUp(rLayout);
2805:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2806:   PetscLayoutDestroy(&rLayout);
2807:   PetscCalloc2(locRows,&dnz,locRows,&onz);
2808:   PetscHSetIJCreate(&ht);
2809:   for (field = 0; field < Nf; ++field) {
2810:     PetscObject      obj;
2811:     PetscClassId     id;
2812:     PetscQuadrature  quad;
2813:     const PetscReal *qpoints;
2814:     PetscInt         Nq, Nc, i, d;

2816:     PetscDSGetDiscretization(prob, field, &obj);
2817:     PetscObjectGetClassId(obj, &id);
2818:     if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);}
2819:     else                       {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2820:     PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);
2821:     /* For each fine grid cell */
2822:     for (cell = cStart; cell < cEnd; ++cell) {
2823:       Vec                pointVec;
2824:       PetscScalar       *pV;
2825:       PetscSF            coarseCellSF = NULL;
2826:       const PetscSFNode *coarseCells;
2827:       PetscInt           numCoarseCells, q, c;
2828:       PetscInt          *findices,   *cindices;
2829:       PetscInt           numFIndices, numCIndices;

2831:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2832:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2833:       /* Get points from the quadrature */
2834:       VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2835:       VecSetBlockSize(pointVec, dim);
2836:       VecGetArray(pointVec, &pV);
2837:       for (q = 0; q < Nq; ++q) {
2838:         const PetscReal xi0[3] = {-1., -1., -1.};

2840:         /* Transform point to real space */
2841:         CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2842:         for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2843:       }
2844:       VecRestoreArray(pointVec, &pV);
2845:       /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2846:       DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2847:       PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2848:       /* Update preallocation info */
2849:       PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2850:       if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2851:       {
2852:         PetscHashIJKey key;
2853:         PetscBool      missing;

2855:         for (i = 0; i < numFIndices; ++i) {
2856:           key.i = findices[i];
2857:           if (key.i >= 0) {
2858:             /* Get indices for coarse elements */
2859:             for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2860:               DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2861:               for (c = 0; c < numCIndices; ++c) {
2862:                 key.j = cindices[c];
2863:                 if (key.j < 0) continue;
2864:                 PetscHSetIJQueryAdd(ht, key, &missing);
2865:                 if (missing) {
2866:                   if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2867:                   else                                     ++onz[key.i-rStart];
2868:                 }
2869:               }
2870:               DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2871:             }
2872:           }
2873:         }
2874:       }
2875:       PetscSFDestroy(&coarseCellSF);
2876:       VecDestroy(&pointVec);
2877:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2878:     }
2879:   }
2880:   PetscHSetIJDestroy(&ht);
2881:   MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
2882:   MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2883:   PetscFree2(dnz,onz);
2884:   for (field = 0; field < Nf; ++field) {
2885:     PetscObject       obj;
2886:     PetscClassId      id;
2887:     PetscTabulation T, Tfine;
2888:     PetscQuadrature   quad;
2889:     const PetscReal  *qpoints, *qweights;
2890:     PetscInt          Nq, Nc, i, d;

2892:     PetscDSGetDiscretization(prob, field, &obj);
2893:     PetscObjectGetClassId(obj, &id);
2894:     if (id == PETSCFE_CLASSID) {
2895:       PetscFEGetQuadrature((PetscFE) obj, &quad);
2896:       PetscFEGetCellTabulation((PetscFE) obj, &Tfine);
2897:       PetscFECreateTabulation((PetscFE) obj, 1, 1, x, 0, &T);
2898:     } else {
2899:       PetscFVGetQuadrature((PetscFV) obj, &quad);
2900:     }
2901:     PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, &qweights);
2902:     /* For each fine grid cell */
2903:     for (cell = cStart; cell < cEnd; ++cell) {
2904:       Vec                pointVec;
2905:       PetscScalar       *pV;
2906:       PetscSF            coarseCellSF = NULL;
2907:       const PetscSFNode *coarseCells;
2908:       PetscInt           numCoarseCells, cpdim, q, c, j;
2909:       PetscInt          *findices,   *cindices;
2910:       PetscInt           numFIndices, numCIndices;

2912:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2913:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2914:       /* Get points from the quadrature */
2915:       VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2916:       VecSetBlockSize(pointVec, dim);
2917:       VecGetArray(pointVec, &pV);
2918:       for (q = 0; q < Nq; ++q) {
2919:         const PetscReal xi0[3] = {-1., -1., -1.};

2921:         /* Transform point to real space */
2922:         CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2923:         for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2924:       }
2925:       VecRestoreArray(pointVec, &pV);
2926:       /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2927:       DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2928:       /* Update matrix */
2929:       PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2930:       if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2931:       VecGetArray(pointVec, &pV);
2932:       for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2933:         PetscReal pVReal[3];
2934:         const PetscReal xi0[3] = {-1., -1., -1.};


2937:         DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2938:         /* Transform points from real space to coarse reference space */
2939:         DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);
2940:         for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]);
2941:         CoordinatesRealToRef(dim, dim, xi0, v0c, invJc, pVReal, x);

2943:         if (id == PETSCFE_CLASSID) {
2944:           PetscFE fe = (PetscFE) obj;

2946:           /* Evaluate coarse basis on contained point */
2947:           PetscFEGetDimension(fe, &cpdim);
2948:           PetscFEComputeTabulation(fe, 1, x, 0, T);
2949:           /* Get elemMat entries by multiplying by weight */
2950:           for (i = 0; i < numFIndices; ++i) {
2951:             PetscArrayzero(elemMat, cpdim);
2952:             for (j = 0; j < cpdim; ++j) {
2953:               for (c = 0; c < Nc; ++c) elemMat[j] += T->T[0][j*Nc + c]*Tfine->T[0][(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ;
2954:             }
2955:             /* Update interpolator */
2956:             if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2957:             if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2958:             MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2959:           }
2960:         } else {
2961:           cpdim = 1;
2962:           for (i = 0; i < numFIndices; ++i) {
2963:             PetscArrayzero(elemMat, cpdim);
2964:             for (j = 0; j < cpdim; ++j) {
2965:               for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ;
2966:             }
2967:             /* Update interpolator */
2968:             if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2969:             PetscPrintf(PETSC_COMM_SELF, "Nq: %D %D Nf: %D %D Nc: %D %D\n", ccell, Nq, i, numFIndices, j, numCIndices);
2970:             if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2971:             MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2972:           }
2973:         }
2974:         DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, PETSC_FALSE, &numCIndices, &cindices, NULL, NULL);
2975:       }
2976:       VecRestoreArray(pointVec, &pV);
2977:       PetscSFDestroy(&coarseCellSF);
2978:       VecDestroy(&pointVec);
2979:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
2980:     }
2981:     if (id == PETSCFE_CLASSID) {PetscTabulationDestroy(&T);}
2982:   }
2983:   PetscFree3(v0,J,invJ);
2984:   PetscFree3(v0c,Jc,invJc);
2985:   PetscFree(elemMat);
2986:   MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
2987:   MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
2988:   return(0);
2989: }

2991: /*@
2992:   DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns

2994:   Input Parameters:
2995: + dmc  - The coarse mesh
2996: - dmf  - The fine mesh
2997: - user - The user context

2999:   Output Parameter:
3000: . sc   - The mapping

3002:   Level: developer

3004: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
3005: @*/
3006: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
3007: {
3008:   PetscDS        prob;
3009:   PetscFE       *feRef;
3010:   PetscFV       *fvRef;
3011:   Vec            fv, cv;
3012:   IS             fis, cis;
3013:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
3014:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
3015:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, c, dim, d, startC, endC, offsetC, offsetF, m;
3016:   PetscBool     *needAvg;

3020:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
3021:   DMGetDimension(dmf, &dim);
3022:   DMGetLocalSection(dmf, &fsection);
3023:   DMGetGlobalSection(dmf, &fglobalSection);
3024:   DMGetLocalSection(dmc, &csection);
3025:   DMGetGlobalSection(dmc, &cglobalSection);
3026:   PetscSectionGetNumFields(fsection, &Nf);
3027:   DMPlexGetSimplexOrBoxCells(dmc, 0, &cStart, &cEnd);
3028:   DMGetDS(dmc, &prob);
3029:   PetscCalloc3(Nf,&feRef,Nf,&fvRef,Nf,&needAvg);
3030:   for (f = 0; f < Nf; ++f) {
3031:     PetscObject  obj;
3032:     PetscClassId id;
3033:     PetscInt     fNb = 0, Nc = 0;

3035:     PetscDSGetDiscretization(prob, f, &obj);
3036:     PetscObjectGetClassId(obj, &id);
3037:     if (id == PETSCFE_CLASSID) {
3038:       PetscFE    fe = (PetscFE) obj;
3039:       PetscSpace sp;
3040:       PetscInt   maxDegree;

3042:       PetscFERefine(fe, &feRef[f]);
3043:       PetscFEGetDimension(feRef[f], &fNb);
3044:       PetscFEGetNumComponents(fe, &Nc);
3045:       PetscFEGetBasisSpace(fe, &sp);
3046:       PetscSpaceGetDegree(sp, NULL, &maxDegree);
3047:       if (!maxDegree) needAvg[f] = PETSC_TRUE;
3048:     } else if (id == PETSCFV_CLASSID) {
3049:       PetscFV        fv = (PetscFV) obj;
3050:       PetscDualSpace Q;

3052:       PetscFVRefine(fv, &fvRef[f]);
3053:       PetscFVGetDualSpace(fvRef[f], &Q);
3054:       PetscDualSpaceGetDimension(Q, &fNb);
3055:       PetscFVGetNumComponents(fv, &Nc);
3056:       needAvg[f] = PETSC_TRUE;
3057:     }
3058:     fTotDim += fNb;
3059:   }
3060:   PetscDSGetTotalDimension(prob, &cTotDim);
3061:   PetscMalloc1(cTotDim,&cmap);
3062:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
3063:     PetscFE        feC;
3064:     PetscFV        fvC;
3065:     PetscDualSpace QF, QC;
3066:     PetscInt       order = -1, NcF, NcC, fpdim, cpdim;

3068:     if (feRef[field]) {
3069:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
3070:       PetscFEGetNumComponents(feC, &NcC);
3071:       PetscFEGetNumComponents(feRef[field], &NcF);
3072:       PetscFEGetDualSpace(feRef[field], &QF);
3073:       PetscDualSpaceGetOrder(QF, &order);
3074:       PetscDualSpaceGetDimension(QF, &fpdim);
3075:       PetscFEGetDualSpace(feC, &QC);
3076:       PetscDualSpaceGetDimension(QC, &cpdim);
3077:     } else {
3078:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
3079:       PetscFVGetNumComponents(fvC, &NcC);
3080:       PetscFVGetNumComponents(fvRef[field], &NcF);
3081:       PetscFVGetDualSpace(fvRef[field], &QF);
3082:       PetscDualSpaceGetDimension(QF, &fpdim);
3083:       PetscFVGetDualSpace(fvC, &QC);
3084:       PetscDualSpaceGetDimension(QC, &cpdim);
3085:     }
3086:     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);
3087:     for (c = 0; c < cpdim; ++c) {
3088:       PetscQuadrature  cfunc;
3089:       const PetscReal *cqpoints, *cqweights;
3090:       PetscInt         NqcC, NpC;
3091:       PetscBool        found = PETSC_FALSE;

3093:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
3094:       PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);
3095:       if (NqcC != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components %D", NqcC, NcC);
3096:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
3097:       for (f = 0; f < fpdim; ++f) {
3098:         PetscQuadrature  ffunc;
3099:         const PetscReal *fqpoints, *fqweights;
3100:         PetscReal        sum = 0.0;
3101:         PetscInt         NqcF, NpF;

3103:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
3104:         PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);
3105:         if (NqcF != NcF) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components %D", NqcF, NcF);
3106:         if (NpC != NpF) continue;
3107:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
3108:         if (sum > 1.0e-9) continue;
3109:         for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]);
3110:         if (sum < 1.0e-9) continue;
3111:         cmap[offsetC+c] = offsetF+f;
3112:         found = PETSC_TRUE;
3113:         break;
3114:       }
3115:       if (!found) {
3116:         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
3117:         if (fvRef[field] || (feRef[field] && order == 0)) {
3118:           cmap[offsetC+c] = offsetF+0;
3119:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
3120:       }
3121:     }
3122:     offsetC += cpdim;
3123:     offsetF += fpdim;
3124:   }
3125:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
3126:   PetscFree3(feRef,fvRef,needAvg);

3128:   DMGetGlobalVector(dmf, &fv);
3129:   DMGetGlobalVector(dmc, &cv);
3130:   VecGetOwnershipRange(cv, &startC, &endC);
3131:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
3132:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
3133:   PetscMalloc1(m,&cindices);
3134:   PetscMalloc1(m,&findices);
3135:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
3136:   for (c = cStart; c < cEnd; ++c) {
3137:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
3138:     for (d = 0; d < cTotDim; ++d) {
3139:       if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
3140:       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]]);
3141:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
3142:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
3143:     }
3144:   }
3145:   PetscFree(cmap);
3146:   PetscFree2(cellCIndices,cellFIndices);

3148:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
3149:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
3150:   VecScatterCreate(cv, cis, fv, fis, sc);
3151:   ISDestroy(&cis);
3152:   ISDestroy(&fis);
3153:   DMRestoreGlobalVector(dmf, &fv);
3154:   DMRestoreGlobalVector(dmc, &cv);
3155:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
3156:   return(0);
3157: }

3159: /*@C
3160:   DMPlexGetCellFields - Retrieve the field values values for a chunk of cells

3162:   Input Parameters:
3163: + dm     - The DM
3164: . cellIS - The cells to include
3165: . locX   - A local vector with the solution fields
3166: . locX_t - A local vector with solution field time derivatives, or NULL
3167: - locA   - A local vector with auxiliary fields, or NULL

3169:   Output Parameters:
3170: + u   - The field coefficients
3171: . u_t - The fields derivative coefficients
3172: - a   - The auxiliary field coefficients

3174:   Level: developer

3176: .seealso: DMPlexGetFaceFields()
3177: @*/
3178: PetscErrorCode DMPlexGetCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
3179: {
3180:   DM              plex, plexA = NULL;
3181:   DMEnclosureType encAux;
3182:   PetscSection    section, sectionAux;
3183:   PetscDS         prob;
3184:   const PetscInt *cells;
3185:   PetscInt        cStart, cEnd, numCells, totDim, totDimAux, c;
3186:   PetscErrorCode  ierr;

3196:   DMPlexConvertPlex(dm, &plex, PETSC_FALSE);
3197:   ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3198:   DMGetLocalSection(dm, &section);
3199:   DMGetCellDS(dm, cells ? cells[cStart] : cStart, &prob);
3200:   PetscDSGetTotalDimension(prob, &totDim);
3201:   if (locA) {
3202:     DM      dmAux;
3203:     PetscDS probAux;

3205:     VecGetDM(locA, &dmAux);
3206:     DMGetEnclosureRelation(dmAux, dm, &encAux);
3207:     DMPlexConvertPlex(dmAux, &plexA, PETSC_FALSE);
3208:     DMGetLocalSection(dmAux, &sectionAux);
3209:     DMGetDS(dmAux, &probAux);
3210:     PetscDSGetTotalDimension(probAux, &totDimAux);
3211:   }
3212:   numCells = cEnd - cStart;
3213:   DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u);
3214:   if (locX_t) {DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u_t);} else {*u_t = NULL;}
3215:   if (locA)   {DMGetWorkArray(dm, numCells*totDimAux, MPIU_SCALAR, a);} else {*a = NULL;}
3216:   for (c = cStart; c < cEnd; ++c) {
3217:     const PetscInt cell = cells ? cells[c] : c;
3218:     const PetscInt cind = c - cStart;
3219:     PetscScalar   *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
3220:     PetscInt       i;

3222:     DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x);
3223:     for (i = 0; i < totDim; ++i) ul[cind*totDim+i] = x[i];
3224:     DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x);
3225:     if (locX_t) {
3226:       DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x_t);
3227:       for (i = 0; i < totDim; ++i) ul_t[cind*totDim+i] = x_t[i];
3228:       DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x_t);
3229:     }
3230:     if (locA) {
3231:       PetscInt subcell;
3232:       DMGetEnclosurePoint(plexA, dm, encAux, cell, &subcell);
3233:       DMPlexVecGetClosure(plexA, sectionAux, locA, subcell, NULL, &x);
3234:       for (i = 0; i < totDimAux; ++i) al[cind*totDimAux+i] = x[i];
3235:       DMPlexVecRestoreClosure(plexA, sectionAux, locA, subcell, NULL, &x);
3236:     }
3237:   }
3238:   DMDestroy(&plex);
3239:   if (locA) {DMDestroy(&plexA);}
3240:   ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
3241:   return(0);
3242: }

3244: /*@C
3245:   DMPlexRestoreCellFields - Restore the field values values for a chunk of cells

3247:   Input Parameters:
3248: + dm     - The DM
3249: . cellIS - The cells to include
3250: . locX   - A local vector with the solution fields
3251: . locX_t - A local vector with solution field time derivatives, or NULL
3252: - locA   - A local vector with auxiliary fields, or NULL

3254:   Output Parameters:
3255: + u   - The field coefficients
3256: . u_t - The fields derivative coefficients
3257: - a   - The auxiliary field coefficients

3259:   Level: developer

3261: .seealso: DMPlexGetFaceFields()
3262: @*/
3263: PetscErrorCode DMPlexRestoreCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
3264: {

3268:   DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u);
3269:   if (locX_t) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u_t);}
3270:   if (locA)   {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, a);}
3271:   return(0);
3272: }

3274: /*@C
3275:   DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces

3277:   Input Parameters:
3278: + dm     - The DM
3279: . fStart - The first face to include
3280: . fEnd   - The first face to exclude
3281: . locX   - A local vector with the solution fields
3282: . locX_t - A local vector with solution field time derivatives, or NULL
3283: . faceGeometry - A local vector with face geometry
3284: . cellGeometry - A local vector with cell geometry
3285: - locaGrad - A local vector with field gradients, or NULL

3287:   Output Parameters:
3288: + Nface - The number of faces with field values
3289: . uL - The field values at the left side of the face
3290: - uR - The field values at the right side of the face

3292:   Level: developer

3294: .seealso: DMPlexGetCellFields()
3295: @*/
3296: PetscErrorCode DMPlexGetFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR)
3297: {
3298:   DM                 dmFace, dmCell, dmGrad = NULL;
3299:   PetscSection       section;
3300:   PetscDS            prob;
3301:   DMLabel            ghostLabel;
3302:   const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
3303:   PetscBool         *isFE;
3304:   PetscInt           dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
3305:   PetscErrorCode     ierr;

3316:   DMGetDimension(dm, &dim);
3317:   DMGetDS(dm, &prob);
3318:   DMGetLocalSection(dm, &section);
3319:   PetscDSGetNumFields(prob, &Nf);
3320:   PetscDSGetTotalComponents(prob, &Nc);
3321:   PetscMalloc1(Nf, &isFE);
3322:   for (f = 0; f < Nf; ++f) {
3323:     PetscObject  obj;
3324:     PetscClassId id;

3326:     PetscDSGetDiscretization(prob, f, &obj);
3327:     PetscObjectGetClassId(obj, &id);
3328:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
3329:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
3330:     else                            {isFE[f] = PETSC_FALSE;}
3331:   }
3332:   DMGetLabel(dm, "ghost", &ghostLabel);
3333:   VecGetArrayRead(locX, &x);
3334:   VecGetDM(faceGeometry, &dmFace);
3335:   VecGetArrayRead(faceGeometry, &facegeom);
3336:   VecGetDM(cellGeometry, &dmCell);
3337:   VecGetArrayRead(cellGeometry, &cellgeom);
3338:   if (locGrad) {
3339:     VecGetDM(locGrad, &dmGrad);
3340:     VecGetArrayRead(locGrad, &lgrad);
3341:   }
3342:   DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uL);
3343:   DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uR);
3344:   /* Right now just eat the extra work for FE (could make a cell loop) */
3345:   for (face = fStart, iface = 0; face < fEnd; ++face) {
3346:     const PetscInt        *cells;
3347:     PetscFVFaceGeom       *fg;
3348:     PetscFVCellGeom       *cgL, *cgR;
3349:     PetscScalar           *xL, *xR, *gL, *gR;
3350:     PetscScalar           *uLl = *uL, *uRl = *uR;
3351:     PetscInt               ghost, nsupp, nchild;

3353:     DMLabelGetValue(ghostLabel, face, &ghost);
3354:     DMPlexGetSupportSize(dm, face, &nsupp);
3355:     DMPlexGetTreeChildren(dm, face, &nchild, NULL);
3356:     if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
3357:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
3358:     DMPlexGetSupport(dm, face, &cells);
3359:     DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
3360:     DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
3361:     for (f = 0; f < Nf; ++f) {
3362:       PetscInt off;

3364:       PetscDSGetComponentOffset(prob, f, &off);
3365:       if (isFE[f]) {
3366:         const PetscInt *cone;
3367:         PetscInt        comp, coneSizeL, coneSizeR, faceLocL, faceLocR, ldof, rdof, d;

3369:         xL = xR = NULL;
3370:         PetscSectionGetFieldComponents(section, f, &comp);
3371:         DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
3372:         DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
3373:         DMPlexGetCone(dm, cells[0], &cone);
3374:         DMPlexGetConeSize(dm, cells[0], &coneSizeL);
3375:         for (faceLocL = 0; faceLocL < coneSizeL; ++faceLocL) if (cone[faceLocL] == face) break;
3376:         DMPlexGetCone(dm, cells[1], &cone);
3377:         DMPlexGetConeSize(dm, cells[1], &coneSizeR);
3378:         for (faceLocR = 0; faceLocR < coneSizeR; ++faceLocR) if (cone[faceLocR] == face) break;
3379:         if (faceLocL == coneSizeL && faceLocR == coneSizeR) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of cell %D or cell %D", face, cells[0], cells[1]);
3380:         /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
3381:         /* TODO: this is a hack that might not be right for nonconforming */
3382:         if (faceLocL < coneSizeL) {
3383:           PetscFEEvaluateFaceFields_Internal(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
3384:           if (rdof == ldof && faceLocR < coneSizeR) {PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
3385:           else              {for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
3386:         }
3387:         else {
3388:           PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);
3389:           PetscSectionGetFieldComponents(section, f, &comp);
3390:           for(d = 0; d < comp; ++d) uLl[iface*Nc+off+d] = uRl[iface*Nc+off+d];
3391:         }
3392:         DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
3393:         DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
3394:       } else {
3395:         PetscFV  fv;
3396:         PetscInt numComp, c;

3398:         PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
3399:         PetscFVGetNumComponents(fv, &numComp);
3400:         DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
3401:         DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
3402:         if (dmGrad) {
3403:           PetscReal dxL[3], dxR[3];

3405:           DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
3406:           DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
3407:           DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
3408:           DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
3409:           for (c = 0; c < numComp; ++c) {
3410:             uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
3411:             uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
3412:           }
3413:         } else {
3414:           for (c = 0; c < numComp; ++c) {
3415:             uLl[iface*Nc+off+c] = xL[c];
3416:             uRl[iface*Nc+off+c] = xR[c];
3417:           }
3418:         }
3419:       }
3420:     }
3421:     ++iface;
3422:   }
3423:   *Nface = iface;
3424:   VecRestoreArrayRead(locX, &x);
3425:   VecRestoreArrayRead(faceGeometry, &facegeom);
3426:   VecRestoreArrayRead(cellGeometry, &cellgeom);
3427:   if (locGrad) {
3428:     VecRestoreArrayRead(locGrad, &lgrad);
3429:   }
3430:   PetscFree(isFE);
3431:   return(0);
3432: }

3434: /*@C
3435:   DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces

3437:   Input Parameters:
3438: + dm     - The DM
3439: . fStart - The first face to include
3440: . fEnd   - The first face to exclude
3441: . locX   - A local vector with the solution fields
3442: . locX_t - A local vector with solution field time derivatives, or NULL
3443: . faceGeometry - A local vector with face geometry
3444: . cellGeometry - A local vector with cell geometry
3445: - locaGrad - A local vector with field gradients, or NULL

3447:   Output Parameters:
3448: + Nface - The number of faces with field values
3449: . uL - The field values at the left side of the face
3450: - uR - The field values at the right side of the face

3452:   Level: developer

3454: .seealso: DMPlexGetFaceFields()
3455: @*/
3456: PetscErrorCode DMPlexRestoreFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR)
3457: {

3461:   DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uL);
3462:   DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uR);
3463:   return(0);
3464: }

3466: /*@C
3467:   DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces

3469:   Input Parameters:
3470: + dm     - The DM
3471: . fStart - The first face to include
3472: . fEnd   - The first face to exclude
3473: . faceGeometry - A local vector with face geometry
3474: - cellGeometry - A local vector with cell geometry

3476:   Output Parameters:
3477: + Nface - The number of faces with field values
3478: . fgeom - The extract the face centroid and normal
3479: - vol   - The cell volume

3481:   Level: developer

3483: .seealso: DMPlexGetCellFields()
3484: @*/
3485: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
3486: {
3487:   DM                 dmFace, dmCell;
3488:   DMLabel            ghostLabel;
3489:   const PetscScalar *facegeom, *cellgeom;
3490:   PetscInt           dim, numFaces = fEnd - fStart, iface, face;
3491:   PetscErrorCode     ierr;

3499:   DMGetDimension(dm, &dim);
3500:   DMGetLabel(dm, "ghost", &ghostLabel);
3501:   VecGetDM(faceGeometry, &dmFace);
3502:   VecGetArrayRead(faceGeometry, &facegeom);
3503:   VecGetDM(cellGeometry, &dmCell);
3504:   VecGetArrayRead(cellGeometry, &cellgeom);
3505:   PetscMalloc1(numFaces, fgeom);
3506:   DMGetWorkArray(dm, numFaces*2, MPIU_SCALAR, vol);
3507:   for (face = fStart, iface = 0; face < fEnd; ++face) {
3508:     const PetscInt        *cells;
3509:     PetscFVFaceGeom       *fg;
3510:     PetscFVCellGeom       *cgL, *cgR;
3511:     PetscFVFaceGeom       *fgeoml = *fgeom;
3512:     PetscReal             *voll   = *vol;
3513:     PetscInt               ghost, d, nchild, nsupp;

3515:     DMLabelGetValue(ghostLabel, face, &ghost);
3516:     DMPlexGetSupportSize(dm, face, &nsupp);
3517:     DMPlexGetTreeChildren(dm, face, &nchild, NULL);
3518:     if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
3519:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
3520:     DMPlexGetSupport(dm, face, &cells);
3521:     DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
3522:     DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
3523:     for (d = 0; d < dim; ++d) {
3524:       fgeoml[iface].centroid[d] = fg->centroid[d];
3525:       fgeoml[iface].normal[d]   = fg->normal[d];
3526:     }
3527:     voll[iface*2+0] = cgL->volume;
3528:     voll[iface*2+1] = cgR->volume;
3529:     ++iface;
3530:   }
3531:   *Nface = iface;
3532:   VecRestoreArrayRead(faceGeometry, &facegeom);
3533:   VecRestoreArrayRead(cellGeometry, &cellgeom);
3534:   return(0);
3535: }

3537: /*@C
3538:   DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces

3540:   Input Parameters:
3541: + dm     - The DM
3542: . fStart - The first face to include
3543: . fEnd   - The first face to exclude
3544: . faceGeometry - A local vector with face geometry
3545: - cellGeometry - A local vector with cell geometry

3547:   Output Parameters:
3548: + Nface - The number of faces with field values
3549: . fgeom - The extract the face centroid and normal
3550: - vol   - The cell volume

3552:   Level: developer

3554: .seealso: DMPlexGetFaceFields()
3555: @*/
3556: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
3557: {

3561:   PetscFree(*fgeom);
3562:   DMRestoreWorkArray(dm, 0, MPIU_REAL, vol);
3563:   return(0);
3564: }

3566: PetscErrorCode DMSNESGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3567: {
3568:   char            composeStr[33] = {0};
3569:   PetscObjectId   id;
3570:   PetscContainer  container;
3571:   PetscErrorCode  ierr;

3574:   PetscObjectGetId((PetscObject)quad,&id);
3575:   PetscSNPrintf(composeStr, 32, "DMSNESGetFEGeom_%x\n", id);
3576:   PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
3577:   if (container) {
3578:     PetscContainerGetPointer(container, (void **) geom);
3579:   } else {
3580:     DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
3581:     PetscContainerCreate(PETSC_COMM_SELF,&container);
3582:     PetscContainerSetPointer(container, (void *) *geom);
3583:     PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
3584:     PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
3585:     PetscContainerDestroy(&container);
3586:   }
3587:   return(0);
3588: }

3590: PetscErrorCode DMSNESRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3591: {
3593:   *geom = NULL;
3594:   return(0);
3595: }

3597: PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM dm, PetscSection section, IS cellIS, PetscReal t, Vec locX, Vec locX_t, Vec locF, void *user)
3598: {
3599:   DM_Plex         *mesh       = (DM_Plex *) dm->data;
3600:   const char      *name       = "Residual";
3601:   DM               dmAux      = NULL;
3602:   DMLabel          ghostLabel = NULL;
3603:   PetscDS          prob       = NULL;
3604:   PetscDS          probAux    = NULL;
3605:   PetscBool        useFEM     = PETSC_FALSE;
3606:   PetscBool        isImplicit = (locX_t || t == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
3607:   DMField          coordField = NULL;
3608:   Vec              locA;
3609:   PetscScalar     *u = NULL, *u_t, *a, *uL = NULL, *uR = NULL;
3610:   IS               chunkIS;
3611:   const PetscInt  *cells;
3612:   PetscInt         cStart, cEnd, numCells;
3613:   PetscInt         Nf, f, totDim, totDimAux, numChunks, cellChunkSize, chunk, fStart, fEnd;
3614:   PetscInt         maxDegree = PETSC_MAX_INT;
3615:   PetscQuadrature  affineQuad = NULL, *quads = NULL;
3616:   PetscFEGeom     *affineGeom = NULL, **geoms = NULL;
3617:   PetscErrorCode   ierr;

3620:   PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
3621:   /* FEM+FVM */
3622:   /* 1: Get sizes from dm and dmAux */
3623:   DMGetLabel(dm, "ghost", &ghostLabel);
3624:   DMGetDS(dm, &prob);
3625:   PetscDSGetNumFields(prob, &Nf);
3626:   PetscDSGetTotalDimension(prob, &totDim);
3627:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
3628:   if (locA) {
3629:     VecGetDM(locA, &dmAux);
3630:     DMGetDS(dmAux, &probAux);
3631:     PetscDSGetTotalDimension(probAux, &totDimAux);
3632:   }
3633:   /* 2: Get geometric data */
3634:   for (f = 0; f < Nf; ++f) {
3635:     PetscObject  obj;
3636:     PetscClassId id;
3637:     PetscBool    fimp;

3639:     PetscDSGetImplicit(prob, f, &fimp);
3640:     if (isImplicit != fimp) continue;
3641:     PetscDSGetDiscretization(prob, f, &obj);
3642:     PetscObjectGetClassId(obj, &id);
3643:     if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
3644:     if (id == PETSCFV_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Use of FVM with PCPATCH not yet implemented");
3645:   }
3646:   if (useFEM) {
3647:     DMGetCoordinateField(dm, &coordField);
3648:     DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
3649:     if (maxDegree <= 1) {
3650:       DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
3651:       if (affineQuad) {
3652:         DMSNESGetFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3653:       }
3654:     } else {
3655:       PetscCalloc2(Nf,&quads,Nf,&geoms);
3656:       for (f = 0; f < Nf; ++f) {
3657:         PetscObject  obj;
3658:         PetscClassId id;
3659:         PetscBool    fimp;

3661:         PetscDSGetImplicit(prob, f, &fimp);
3662:         if (isImplicit != fimp) continue;
3663:         PetscDSGetDiscretization(prob, f, &obj);
3664:         PetscObjectGetClassId(obj, &id);
3665:         if (id == PETSCFE_CLASSID) {
3666:           PetscFE fe = (PetscFE) obj;

3668:           PetscFEGetQuadrature(fe, &quads[f]);
3669:           PetscObjectReference((PetscObject)quads[f]);
3670:           DMSNESGetFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3671:         }
3672:       }
3673:     }
3674:   }
3675:   /* Loop over chunks */
3676:   ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3677:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
3678:   if (useFEM) {ISCreate(PETSC_COMM_SELF, &chunkIS);}
3679:   numCells      = cEnd - cStart;
3680:   numChunks     = 1;
3681:   cellChunkSize = numCells/numChunks;
3682:   numChunks     = PetscMin(1,numCells);
3683:   for (chunk = 0; chunk < numChunks; ++chunk) {
3684:     PetscScalar     *elemVec, *fluxL = NULL, *fluxR = NULL;
3685:     PetscReal       *vol = NULL;
3686:     PetscFVFaceGeom *fgeom = NULL;
3687:     PetscInt         cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c;
3688:     PetscInt         numFaces = 0;

3690:     /* Extract field coefficients */
3691:     if (useFEM) {
3692:       ISGetPointSubrange(chunkIS, cS, cE, cells);
3693:       DMPlexGetCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3694:       DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3695:       PetscArrayzero(elemVec, numCells*totDim);
3696:     }
3697:     /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
3698:     /* Loop over fields */
3699:     for (f = 0; f < Nf; ++f) {
3700:       PetscObject  obj;
3701:       PetscClassId id;
3702:       PetscBool    fimp;
3703:       PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

3705:       PetscDSGetImplicit(prob, f, &fimp);
3706:       if (isImplicit != fimp) continue;
3707:       PetscDSGetDiscretization(prob, f, &obj);
3708:       PetscObjectGetClassId(obj, &id);
3709:       if (id == PETSCFE_CLASSID) {
3710:         PetscFE         fe = (PetscFE) obj;
3711:         PetscFEGeom    *geom = affineGeom ? affineGeom : geoms[f];
3712:         PetscFEGeom    *chunkGeom = NULL;
3713:         PetscQuadrature quad = affineQuad ? affineQuad : quads[f];
3714:         PetscInt        Nq, Nb;

3716:         PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3717:         PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
3718:         PetscFEGetDimension(fe, &Nb);
3719:         blockSize = Nb;
3720:         batchSize = numBlocks * blockSize;
3721:         PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3722:         numChunks = numCells / (numBatches*batchSize);
3723:         Ne        = numChunks*numBatches*batchSize;
3724:         Nr        = numCells % (numBatches*batchSize);
3725:         offset    = numCells - Nr;
3726:         /* Integrate FE residual to get elemVec (need fields at quadrature points) */
3727:         /*   For FV, I think we use a P0 basis and the cell coefficients (for subdivided cells, we can tweak the basis tabulation to be the indicator function) */
3728:         PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);
3729:         PetscFEIntegrateResidual(prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);
3730:         PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);
3731:         PetscFEIntegrateResidual(prob, f, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);
3732:         PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);
3733:       } else if (id == PETSCFV_CLASSID) {
3734:         PetscFV fv = (PetscFV) obj;

3736:         Ne = numFaces;
3737:         /* Riemann solve over faces (need fields at face centroids) */
3738:         /*   We need to evaluate FE fields at those coordinates */
3739:         PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
3740:       } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
3741:     }
3742:     /* Loop over domain */
3743:     if (useFEM) {
3744:       /* Add elemVec to locX */
3745:       for (c = cS; c < cE; ++c) {
3746:         const PetscInt cell = cells ? cells[c] : c;
3747:         const PetscInt cind = c - cStart;

3749:         if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cind*totDim]);}
3750:         if (ghostLabel) {
3751:           PetscInt ghostVal;

3753:           DMLabelGetValue(ghostLabel,cell,&ghostVal);
3754:           if (ghostVal > 0) continue;
3755:         }
3756:         DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);
3757:       }
3758:     }
3759:     /* Handle time derivative */
3760:     if (locX_t) {
3761:       PetscScalar *x_t, *fa;

3763:       VecGetArray(locF, &fa);
3764:       VecGetArray(locX_t, &x_t);
3765:       for (f = 0; f < Nf; ++f) {
3766:         PetscFV      fv;
3767:         PetscObject  obj;
3768:         PetscClassId id;
3769:         PetscInt     pdim, d;

3771:         PetscDSGetDiscretization(prob, f, &obj);
3772:         PetscObjectGetClassId(obj, &id);
3773:         if (id != PETSCFV_CLASSID) continue;
3774:         fv   = (PetscFV) obj;
3775:         PetscFVGetNumComponents(fv, &pdim);
3776:         for (c = cS; c < cE; ++c) {
3777:           const PetscInt cell = cells ? cells[c] : c;
3778:           PetscScalar   *u_t, *r;

3780:           if (ghostLabel) {
3781:             PetscInt ghostVal;

3783:             DMLabelGetValue(ghostLabel, cell, &ghostVal);
3784:             if (ghostVal > 0) continue;
3785:           }
3786:           DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
3787:           DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
3788:           for (d = 0; d < pdim; ++d) r[d] += u_t[d];
3789:         }
3790:       }
3791:       VecRestoreArray(locX_t, &x_t);
3792:       VecRestoreArray(locF, &fa);
3793:     }
3794:     if (useFEM) {
3795:       DMPlexRestoreCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3796:       DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3797:     }
3798:   }
3799:   if (useFEM) {ISDestroy(&chunkIS);}
3800:   ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
3801:   /* TODO Could include boundary residual here (see DMPlexComputeResidual_Internal) */
3802:   if (useFEM) {
3803:     if (maxDegree <= 1) {
3804:       DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3805:       PetscQuadratureDestroy(&affineQuad);
3806:     } else {
3807:       for (f = 0; f < Nf; ++f) {
3808:         DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3809:         PetscQuadratureDestroy(&quads[f]);
3810:       }
3811:       PetscFree2(quads,geoms);
3812:     }
3813:   }
3814:   PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
3815:   return(0);
3816: }

3818: /*
3819:   We always assemble JacP, and if the matrix is different from Jac and two different sets of point functions are provided, we also assemble Jac

3821:   X   - The local solution vector
3822:   X_t - The local solution time derviative vector, or NULL
3823: */
3824: PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM dm, PetscSection section, PetscSection globalSection, IS cellIS,
3825:                                                     PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP, void *ctx)
3826: {
3827:   DM_Plex         *mesh  = (DM_Plex *) dm->data;
3828:   const char      *name = "Jacobian", *nameP = "JacobianPre";
3829:   DM               dmAux = NULL;
3830:   PetscDS          prob,   probAux = NULL;
3831:   PetscSection     sectionAux = NULL;
3832:   Vec              A;
3833:   DMField          coordField;
3834:   PetscFEGeom     *cgeomFEM;
3835:   PetscQuadrature  qGeom = NULL;
3836:   Mat              J = Jac, JP = JacP;
3837:   PetscScalar     *work, *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL, *elemMatP = NULL, *elemMatD = NULL;
3838:   PetscBool        hasJac, hasPrec, hasDyn, assembleJac, isMatIS, isMatISP, *isFE, hasFV = PETSC_FALSE;
3839:   const PetscInt  *cells;
3840:   PetscInt         Nf, fieldI, fieldJ, maxDegree, numCells, cStart, cEnd, numChunks, chunkSize, chunk, totDim, totDimAux = 0, sz, wsz, off = 0, offCell = 0;
3841:   PetscErrorCode   ierr;

3844:   CHKMEMQ;
3845:   ISGetLocalSize(cellIS, &numCells);
3846:   ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3847:   PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
3848:   DMGetDS(dm, &prob);
3849:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3850:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3851:   if (dmAux) {
3852:     DMGetLocalSection(dmAux, &sectionAux);
3853:     DMGetDS(dmAux, &probAux);
3854:   }
3855:   /* Get flags */
3856:   PetscDSGetNumFields(prob, &Nf);
3857:   DMGetWorkArray(dm, Nf, MPIU_BOOL, &isFE);
3858:   for (fieldI = 0; fieldI < Nf; ++fieldI) {
3859:     PetscObject  disc;
3860:     PetscClassId id;
3861:     PetscDSGetDiscretization(prob, fieldI, &disc);
3862:     PetscObjectGetClassId(disc, &id);
3863:     if (id == PETSCFE_CLASSID)      {isFE[fieldI] = PETSC_TRUE;}
3864:     else if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; isFE[fieldI] = PETSC_FALSE;}
3865:   }
3866:   PetscDSHasJacobian(prob, &hasJac);
3867:   PetscDSHasJacobianPreconditioner(prob, &hasPrec);
3868:   PetscDSHasDynamicJacobian(prob, &hasDyn);
3869:   assembleJac = hasJac && hasPrec && (Jac != JacP) ? PETSC_TRUE : PETSC_FALSE;
3870:   hasDyn      = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
3871:   PetscObjectTypeCompare((PetscObject) Jac,  MATIS, &isMatIS);
3872:   PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
3873:   /* Setup input data and temp arrays (should be DMGetWorkArray) */
3874:   if (isMatISP || isMatISP) {DMPlexGetSubdomainSection(dm, &globalSection);}
3875:   if (isMatIS)  {MatISGetLocalMat(Jac,  &J);}
3876:   if (isMatISP) {MatISGetLocalMat(JacP, &JP);}
3877:   if (hasFV)    {MatSetOption(JP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);} /* No allocated space for FV stuff, so ignore the zero entries */
3878:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3879:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3880:   PetscDSGetTotalDimension(prob, &totDim);
3881:   if (probAux) {PetscDSGetTotalDimension(probAux, &totDimAux);}
3882:   CHKMEMQ;
3883:   /* Compute batch sizes */
3884:   if (isFE[0]) {
3885:     PetscFE         fe;
3886:     PetscQuadrature q;
3887:     PetscInt        numQuadPoints, numBatches, batchSize, numBlocks, blockSize, Nb;

3889:     PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3890:     PetscFEGetQuadrature(fe, &q);
3891:     PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
3892:     PetscFEGetDimension(fe, &Nb);
3893:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3894:     blockSize = Nb*numQuadPoints;
3895:     batchSize = numBlocks  * blockSize;
3896:     chunkSize = numBatches * batchSize;
3897:     numChunks = numCells / chunkSize + numCells % chunkSize;
3898:     PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3899:   } else {
3900:     chunkSize = numCells;
3901:     numChunks = 1;
3902:   }
3903:   /* Get work space */
3904:   wsz  = (((X?1:0) + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize;
3905:   DMGetWorkArray(dm, wsz, MPIU_SCALAR, &work);
3906:   PetscArrayzero(work, wsz);
3907:   off      = 0;
3908:   u        = X       ? (sz = chunkSize*totDim,        off += sz, work+off-sz) : NULL;
3909:   u_t      = X_t     ? (sz = chunkSize*totDim,        off += sz, work+off-sz) : NULL;
3910:   a        = dmAux   ? (sz = chunkSize*totDimAux,     off += sz, work+off-sz) : NULL;
3911:   elemMat  = hasJac  ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3912:   elemMatP = hasPrec ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3913:   elemMatD = hasDyn  ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3914:   if (off != wsz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error is workspace size %D should be %D", off, wsz);
3915:   /* Setup geometry */
3916:   DMGetCoordinateField(dm, &coordField);
3917:   DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);
3918:   if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, cellIS, &qGeom);}
3919:   if (!qGeom) {
3920:     PetscFE fe;

3922:     PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3923:     PetscFEGetQuadrature(fe, &qGeom);
3924:     PetscObjectReference((PetscObject) qGeom);
3925:   }
3926:   DMSNESGetFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
3927:   /* Compute volume integrals */
3928:   if (assembleJac) {MatZeroEntries(J);}
3929:   MatZeroEntries(JP);
3930:   for (chunk = 0; chunk < numChunks; ++chunk, offCell += chunkSize) {
3931:     const PetscInt   Ncell = PetscMin(chunkSize, numCells - offCell);
3932:     PetscInt         c;

3934:     /* Extract values */
3935:     for (c = 0; c < Ncell; ++c) {
3936:       const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
3937:       PetscScalar   *x = NULL,  *x_t = NULL;
3938:       PetscInt       i;

3940:       if (X) {
3941:         DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);
3942:         for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
3943:         DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);
3944:       }
3945:       if (X_t) {
3946:         DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);
3947:         for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
3948:         DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);
3949:       }
3950:       if (dmAux) {
3951:         DMPlexVecGetClosure(dmAux, sectionAux, A, cell, NULL, &x);
3952:         for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
3953:         DMPlexVecRestoreClosure(dmAux, sectionAux, A, cell, NULL, &x);
3954:       }
3955:     }
3956:     CHKMEMQ;
3957:     for (fieldI = 0; fieldI < Nf; ++fieldI) {
3958:       PetscFE fe;
3959:       PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3960:       for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
3961:         if (hasJac)  {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN,     fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMat);}
3962:         if (hasPrec) {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatP);}
3963:         if (hasDyn)  {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatD);}
3964:       }
3965:       /* For finite volume, add the identity */
3966:       if (!isFE[fieldI]) {
3967:         PetscFV  fv;
3968:         PetscInt eOffset = 0, Nc, fc, foff;

3970:         PetscDSGetFieldOffset(prob, fieldI, &foff);
3971:         PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);
3972:         PetscFVGetNumComponents(fv, &Nc);
3973:         for (c = 0; c < chunkSize; ++c, eOffset += totDim*totDim) {
3974:           for (fc = 0; fc < Nc; ++fc) {
3975:             const PetscInt i = foff + fc;
3976:             if (hasJac)  {elemMat [eOffset+i*totDim+i] = 1.0;}
3977:             if (hasPrec) {elemMatP[eOffset+i*totDim+i] = 1.0;}
3978:           }
3979:         }
3980:       }
3981:     }
3982:     CHKMEMQ;
3983:     /*   Add contribution from X_t */
3984:     if (hasDyn) {for (c = 0; c < chunkSize*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];}
3985:     /* Insert values into matrix */
3986:     for (c = 0; c < Ncell; ++c) {
3987:       const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
3988:       if (mesh->printFEM > 1) {
3989:         if (hasJac)  {DMPrintCellMatrix(cell, name,  totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
3990:         if (hasPrec) {DMPrintCellMatrix(cell, nameP, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
3991:       }
3992:       if (assembleJac) {DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);}
3993:       DMPlexMatSetClosure(dm, section, globalSection, JP, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
3994:     }
3995:     CHKMEMQ;
3996:   }
3997:   /* Cleanup */
3998:   DMSNESRestoreFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
3999:   PetscQuadratureDestroy(&qGeom);
4000:   if (hasFV) {MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);}
4001:   DMRestoreWorkArray(dm, Nf, MPIU_BOOL, &isFE);
4002:   DMRestoreWorkArray(dm, ((1 + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize, MPIU_SCALAR, &work);
4003:   /* Compute boundary integrals */
4004:   /* DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, ctx); */
4005:   /* Assemble matrix */
4006:   if (assembleJac) {MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);}
4007:   MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
4008:   PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
4009:   CHKMEMQ;
4010:   return(0);
4011: }