Actual source code: plexfem.c

petsc-master 2019-08-18
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 Nf, 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:   for (i = 0; i < PetscMin(dim, mmin); ++i) {VecNormalize(mode[i], NULL);}
215:   /* Orthonormalize system */
216:   for (i = dim; i < mmin; ++i) {
217:     PetscScalar dots[6];

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

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

232:   Collective on dm

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

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

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

246:   Level: advanced

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

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

278:       ctx[0] = dimEmbed;
279:       ctx[1] = d;
280:       DMProjectFunctionLabel(dm, 0.0, label, nids[b], &ids[off], 0, NULL, &func, &voidctx, INSERT_VALUES, mode[d]);
281:       off   += nids[b];
282:     }
283:   }
284:   for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
285:   /* Orthonormalize system */
286:   for (i = 0; i < m; ++i) {
287:     VecMDot(mode[i], i, mode, dots);
288:     for (j = 0; j < i; ++j) dots[j] *= -1.0;
289:     VecMAXPY(mode[i], i, dots, mode);
290:     VecNormalize(mode[i], NULL);
291:   }
292:   MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
293:   for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
294:   PetscFree2(mode, dots);
295:   return(0);
296: }

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

307:   Input Parameters:
308: + dm - the DMPlex object
309: - height - the maximum projection height >= 0

311:   Level: advanced

313: .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
314: @*/
315: PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
316: {
317:   DM_Plex *plex = (DM_Plex *) dm->data;

321:   plex->maxProjectionHeight = height;
322:   return(0);
323: }

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

329:   Input Parameters:
330: . dm - the DMPlex object

332:   Output Parameters:
333: . height - the maximum projection height

335:   Level: intermediate

337: .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal()
338: @*/
339: PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
340: {
341:   DM_Plex *plex = (DM_Plex *) dm->data;

345:   *height = plex->maxProjectionHeight;
346:   return(0);
347: }

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

358: /*
359:   Note: Following https://en.wikipedia.org/wiki/Euler_angles, we will specify Euler angles by extrinsic rotations, meaning that
360:   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:
361:   $ The XYZ system rotates about the z axis by alpha. The X axis is now at angle alpha with respect to the x axis.
362:   $ 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.
363:   $ The XYZ system rotates a third time about the z axis by gamma.
364: */
365: static PetscErrorCode DMPlexBasisTransformSetUp_Rotation_Internal(DM dm, void *ctx)
366: {
367:   RotCtx        *rc  = (RotCtx *) ctx;
368:   PetscInt       dim = rc->dim;
369:   PetscReal      c1, s1, c2, s2, c3, s3;

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

397: static PetscErrorCode DMPlexBasisTransformDestroy_Rotation_Internal(DM dm, void *ctx)
398: {
399:   RotCtx        *rc = (RotCtx *) ctx;

403:   PetscFree2(rc->R, rc->RT);
404:   PetscFree(rc);
405:   return(0);
406: }

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

414:   if (l2g) {*A = rc->R;}
415:   else     {*A = rc->RT;}
416:   return(0);
417: }

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

424:   #if defined(PETSC_USE_COMPLEX)
425:   switch (dim) {
426:     case 2:
427:     {
428:       PetscScalar yt[2], zt[2];

430:       yt[0] = y[0]; yt[1] = y[1];
431:       DMPlexBasisTransformApply_Internal(dm, x, l2g, dim, yt, zt, ctx);
432:       z[0] = PetscRealPart(zt[0]); z[1] = PetscRealPart(zt[1]);
433:     }
434:     break;
435:     case 3:
436:     {
437:       PetscScalar yt[3], zt[3];

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

451: PetscErrorCode DMPlexBasisTransformApply_Internal(DM dm, const PetscReal x[], PetscBool l2g, PetscInt dim, const PetscScalar *y, PetscScalar *z, void *ctx)
452: {
453:   const PetscScalar *A;
454:   PetscErrorCode     ierr;

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

465: static PetscErrorCode DMPlexBasisTransformField_Internal(DM dm, DM tdm, Vec tv, PetscInt p, PetscInt f, PetscBool l2g, PetscScalar *a)
466: {
467:   PetscSection       ts;
468:   const PetscScalar *ta, *tva;
469:   PetscInt           dof;
470:   PetscErrorCode     ierr;

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

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

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

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

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

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

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

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

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

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

624:   Input Parameters:
625: + dm - The DM
626: - lv - A local vector with values in the global basis

628:   Output Parameters:
629: . lv - A local vector with values in the local basis

631:   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.

633:   Level: developer

635: .seealso: DMPlexLocalToGlobalBasis(), DMGetLocalSection()
636: @*/
637: PetscErrorCode DMPlexGlobalToLocalBasis(DM dm, Vec lv)
638: {

644:   DMPlexBasisTransform_Internal(dm, lv, PETSC_FALSE);
645:   return(0);
646: }

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

651:   Input Parameters:
652: + dm - The DM
653: - lv - A local vector with values in the local basis

655:   Output Parameters:
656: . lv - A local vector with values in the global basis

658:   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.

660:   Level: developer

662: .seealso: DMPlexGlobalToLocalBasis(), DMGetLocalSection()
663: @*/
664: PetscErrorCode DMPlexLocalToGlobalBasis(DM dm, Vec lv)
665: {

671:   DMPlexBasisTransform_Internal(dm, lv, PETSC_TRUE);
672:   return(0);
673: }

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

679:   Input Parameters:
680: + dm    - The DM
681: . alpha - The first Euler angle, and in 2D the only one
682: . beta  - The second Euler angle
683: . gamma - The third Euler angle

685:   Note: Following https://en.wikipedia.org/wiki/Euler_angles, we will specify Euler angles by extrinsic rotations, meaning that
686:   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:
687:   $ The XYZ system rotates about the z axis by alpha. The X axis is now at angle alpha with respect to the x axis.
688:   $ 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.
689:   $ The XYZ system rotates a third time about the z axis by gamma.

691:   Level: developer

693: .seealso: DMPlexGlobalToLocalBasis(), DMPlexLocalToGlobalBasis()
694: @*/
695: PetscErrorCode DMPlexCreateBasisRotation(DM dm, PetscReal alpha, PetscReal beta, PetscReal gamma)
696: {
697:   RotCtx        *rc;
698:   PetscInt       cdim;

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

716: /*@C
717:   DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector

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

731:   Output Parameter:
732: . locX   - A local vector to receives the boundary values

734:   Level: developer

736: .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
737: @*/
738: 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)
739: {
740:   PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
741:   void            **ctxs;
742:   PetscInt          numFields;
743:   PetscErrorCode    ierr;

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

755: /*@C
756:   DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector

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

771:   Output Parameter:
772: . locX   - A local vector to receives the boundary values

774:   Level: developer

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

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

804: /*@C
805:   DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector

807:   Input Parameters:
808: + dm     - The DM, with a PetscDS that matches the problem being constrained
809: . time   - The time
810: . faceGeometry - A vector with the FVM face geometry information
811: . cellGeometry - A vector with the FVM cell geometry information
812: . Grad         - A vector with the FVM cell gradient information
813: . field  - The field to constrain
814: . Nc     - The number of constrained field components, or 0 for all components
815: . comps  - An array of constrained component numbers, or NULL for all components
816: . label  - The DMLabel defining constrained points
817: . numids - The number of DMLabel ids for constrained points
818: . ids    - An array of ids for constrained points
819: . func   - A pointwise function giving boundary values
820: - ctx    - An optional user context for bcFunc

822:   Output Parameter:
823: . locX   - A local vector to receives the boundary values

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

827:   Level: developer

829: .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary()
830: @*/
831: 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[],
832:                                                  PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
833: {
834:   PetscDS            prob;
835:   PetscSF            sf;
836:   DM                 dmFace, dmCell, dmGrad;
837:   const PetscScalar *facegeom, *cellgeom = NULL, *grad;
838:   const PetscInt    *leaves;
839:   PetscScalar       *x, *fx;
840:   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
841:   PetscErrorCode     ierr, ierru = 0;

844:   DMGetPointSF(dm, &sf);
845:   PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);
846:   nleaves = PetscMax(0, nleaves);
847:   DMGetDimension(dm, &dim);
848:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
849:   DMGetDS(dm, &prob);
850:   VecGetDM(faceGeometry, &dmFace);
851:   VecGetArrayRead(faceGeometry, &facegeom);
852:   if (cellGeometry) {
853:     VecGetDM(cellGeometry, &dmCell);
854:     VecGetArrayRead(cellGeometry, &cellgeom);
855:   }
856:   if (Grad) {
857:     PetscFV fv;

859:     PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);
860:     VecGetDM(Grad, &dmGrad);
861:     VecGetArrayRead(Grad, &grad);
862:     PetscFVGetNumComponents(fv, &pdim);
863:     DMGetWorkArray(dm, pdim, MPIU_SCALAR, &fx);
864:   }
865:   VecGetArray(locX, &x);
866:   for (i = 0; i < numids; ++i) {
867:     IS              faceIS;
868:     const PetscInt *faces;
869:     PetscInt        numFaces, f;

871:     DMLabelGetStratumIS(label, ids[i], &faceIS);
872:     if (!faceIS) continue; /* No points with that id on this process */
873:     ISGetLocalSize(faceIS, &numFaces);
874:     ISGetIndices(faceIS, &faces);
875:     for (f = 0; f < numFaces; ++f) {
876:       const PetscInt         face = faces[f], *cells;
877:       PetscFVFaceGeom        *fg;

879:       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
880:       PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);
881:       if (loc >= 0) continue;
882:       DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
883:       DMPlexGetSupport(dm, face, &cells);
884:       if (Grad) {
885:         PetscFVCellGeom       *cg;
886:         PetscScalar           *cx, *cgrad;
887:         PetscScalar           *xG;
888:         PetscReal              dx[3];
889:         PetscInt               d;

891:         DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);
892:         DMPlexPointLocalRead(dm, cells[0], x, &cx);
893:         DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);
894:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
895:         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
896:         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
897:         ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);
898:         if (ierru) {
899:           ISRestoreIndices(faceIS, &faces);
900:           ISDestroy(&faceIS);
901:           goto cleanup;
902:         }
903:       } else {
904:         PetscScalar       *xI;
905:         PetscScalar       *xG;

907:         DMPlexPointLocalRead(dm, cells[0], x, &xI);
908:         DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);
909:         ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
910:         if (ierru) {
911:           ISRestoreIndices(faceIS, &faces);
912:           ISDestroy(&faceIS);
913:           goto cleanup;
914:         }
915:       }
916:     }
917:     ISRestoreIndices(faceIS, &faces);
918:     ISDestroy(&faceIS);
919:   }
920:   cleanup:
921:   VecRestoreArray(locX, &x);
922:   if (Grad) {
923:     DMRestoreWorkArray(dm, pdim, MPIU_SCALAR, &fx);
924:     VecRestoreArrayRead(Grad, &grad);
925:   }
926:   if (cellGeometry) {VecRestoreArrayRead(cellGeometry, &cellgeom);}
927:   VecRestoreArrayRead(faceGeometry, &facegeom);
928:   CHKERRQ(ierru);
929:   return(0);
930: }

932: PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
933: {
934:   PetscDS        prob;
935:   PetscInt       numBd, b;

939:   DMGetDS(dm, &prob);
940:   PetscDSGetNumBoundary(prob, &numBd);
941:   for (b = 0; b < numBd; ++b) {
942:     DMBoundaryConditionType type;
943:     const char             *name, *labelname;
944:     DMLabel                 label;
945:     PetscInt                field, Nc;
946:     const PetscInt         *comps;
947:     PetscObject             obj;
948:     PetscClassId            id;
949:     void                    (*func)(void);
950:     PetscInt                numids;
951:     const PetscInt         *ids;
952:     void                   *ctx;

954:     DMGetBoundary(dm, b, &type, &name, &labelname, &field, &Nc, &comps, &func, &numids, &ids, &ctx);
955:     if (insertEssential != (type & DM_BC_ESSENTIAL)) continue;
956:     DMGetLabel(dm, labelname, &label);
957:     if (!label) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "Label %s for boundary condition %s does not exist in the DM", labelname, name);
958:     DMGetField(dm, field, NULL, &obj);
959:     PetscObjectGetClassId(obj, &id);
960:     if (id == PETSCFE_CLASSID) {
961:       switch (type) {
962:         /* for FEM, there is no insertion to be done for non-essential boundary conditions */
963:       case DM_BC_ESSENTIAL:
964:         DMPlexLabelAddCells(dm,label);
965:         DMPlexInsertBoundaryValuesEssential(dm, time, field, Nc, comps, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);
966:         DMPlexLabelClearCells(dm,label);
967:         break;
968:       case DM_BC_ESSENTIAL_FIELD:
969:         DMPlexLabelAddCells(dm,label);
970:         DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, Nc, comps, label, numids, ids,
971:                                                         (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
972:                                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
973:                                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[])) func, ctx, locX);
974:         DMPlexLabelClearCells(dm,label);
975:         break;
976:       default: break;
977:       }
978:     } else if (id == PETSCFV_CLASSID) {
979:       if (!faceGeomFVM) continue;
980:       DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, Nc, comps, label, numids, ids,
981:                                                (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);
982:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
983:   }
984:   return(0);
985: }

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

990:   Input Parameters:
991: + dm - The DM
992: . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions
993: . time - The time
994: . faceGeomFVM - Face geometry data for FV discretizations
995: . cellGeomFVM - Cell geometry data for FV discretizations
996: - gradFVM - Gradient reconstruction data for FV discretizations

998:   Output Parameters:
999: . locX - Solution updated with boundary values

1001:   Level: developer

1003: .seealso: DMProjectFunctionLabelLocal()
1004: @*/
1005: PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
1006: {

1015:   PetscTryMethod(dm,"DMPlexInsertBoundaryValues_C",(DM,PetscBool,Vec,PetscReal,Vec,Vec,Vec),(dm,insertEssential,locX,time,faceGeomFVM,cellGeomFVM,gradFVM));
1016:   return(0);
1017: }

1019: PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1020: {
1021:   Vec              localX;
1022:   PetscErrorCode   ierr;

1025:   DMGetLocalVector(dm, &localX);
1026:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, time, NULL, NULL, NULL);
1027:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1028:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1029:   DMPlexComputeL2DiffLocal(dm, time, funcs, ctxs, localX, diff);
1030:   DMRestoreLocalVector(dm, &localX);
1031:   return(0);
1032: }

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

1037:   Collective on dm

1039:   Input Parameters:
1040: + dm     - The DM
1041: . time   - The time
1042: . funcs  - The functions to evaluate for each field component
1043: . ctxs   - Optional array of contexts to pass to each function, or NULL.
1044: - localX - The coefficient vector u_h, a local vector

1046:   Output Parameter:
1047: . diff - The diff ||u - u_h||_2

1049:   Level: developer

1051: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
1052: @*/
1053: PetscErrorCode DMPlexComputeL2DiffLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec localX, PetscReal *diff)
1054: {
1055:   const PetscInt   debug = ((DM_Plex*)dm->data)->printL2;
1056:   DM               tdm;
1057:   Vec              tv;
1058:   PetscSection     section;
1059:   PetscQuadrature  quad;
1060:   PetscFEGeom      fegeom;
1061:   PetscScalar     *funcVal, *interpolant;
1062:   PetscReal       *coords, *gcoords;
1063:   PetscReal        localDiff = 0.0;
1064:   const PetscReal *quadWeights;
1065:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cellHeight, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1066:   PetscBool        transform;
1067:   PetscErrorCode   ierr;

1070:   DMGetDimension(dm, &dim);
1071:   DMGetCoordinateDim(dm, &coordDim);
1072:   fegeom.dimEmbed = coordDim;
1073:   DMGetLocalSection(dm, &section);
1074:   PetscSectionGetNumFields(section, &numFields);
1075:   DMGetBasisTransformDM_Internal(dm, &tdm);
1076:   DMGetBasisTransformVec_Internal(dm, &tv);
1077:   DMHasBasisTransform(dm, &transform);
1078:   for (field = 0; field < numFields; ++field) {
1079:     PetscObject  obj;
1080:     PetscClassId id;
1081:     PetscInt     Nc;

1083:     DMGetField(dm, field, NULL, &obj);
1084:     PetscObjectGetClassId(obj, &id);
1085:     if (id == PETSCFE_CLASSID) {
1086:       PetscFE fe = (PetscFE) obj;

1088:       PetscFEGetQuadrature(fe, &quad);
1089:       PetscFEGetNumComponents(fe, &Nc);
1090:     } else if (id == PETSCFV_CLASSID) {
1091:       PetscFV fv = (PetscFV) obj;

1093:       PetscFVGetQuadrature(fv, &quad);
1094:       PetscFVGetNumComponents(fv, &Nc);
1095:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1096:     numComponents += Nc;
1097:   }
1098:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
1099:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1100:   PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1101:   DMPlexGetVTKCellHeight(dm, &cellHeight);
1102:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1103:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1104:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1105:   for (c = cStart; c < cEnd; ++c) {
1106:     PetscScalar *x = NULL;
1107:     PetscReal    elemDiff = 0.0;
1108:     PetscInt     qc = 0;

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

1113:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1114:       PetscObject  obj;
1115:       PetscClassId id;
1116:       void * const ctx = ctxs ? ctxs[field] : NULL;
1117:       PetscInt     Nb, Nc, q, fc;

1119:       DMGetField(dm, field, NULL, &obj);
1120:       PetscObjectGetClassId(obj, &id);
1121:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1122:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1123:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1124:       if (debug) {
1125:         char title[1024];
1126:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1127:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1128:       }
1129:       for (q = 0; q < Nq; ++q) {
1130:         PetscFEGeom qgeom;

1132:         qgeom.dimEmbed = fegeom.dimEmbed;
1133:         qgeom.J        = &fegeom.J[q*coordDim*coordDim];
1134:         qgeom.invJ     = &fegeom.invJ[q*coordDim*coordDim];
1135:         qgeom.detJ     = &fegeom.detJ[q];
1136:         if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, point %D", fegeom.detJ[q], c, q);
1137:         if (transform) {
1138:           gcoords = &coords[coordDim*Nq];
1139:           DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);
1140:         } else {
1141:           gcoords = &coords[coordDim*q];
1142:         }
1143:         (*funcs[field])(coordDim, time, gcoords, Nc, funcVal, ctx);
1144:         if (ierr) {
1145:           PetscErrorCode ierr2;
1146:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1147:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1148:           ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1149: 
1150:         }
1151:         if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1152:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1153:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1154:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1155:         for (fc = 0; fc < Nc; ++fc) {
1156:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1157:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %D field %D,%D diff %g\n", c, field, fc, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]);}
1158:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1159:         }
1160:       }
1161:       fieldOffset += Nb;
1162:       qc += Nc;
1163:     }
1164:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1165:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
1166:     localDiff += elemDiff;
1167:   }
1168:   PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1169:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1170:   *diff = PetscSqrtReal(*diff);
1171:   return(0);
1172: }

1174: 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)
1175: {
1176:   const PetscInt   debug = ((DM_Plex*)dm->data)->printL2;
1177:   DM               tdm;
1178:   PetscSection     section;
1179:   PetscQuadrature  quad;
1180:   Vec              localX, tv;
1181:   PetscScalar     *funcVal, *interpolant;
1182:   const PetscReal *quadWeights;
1183:   PetscFEGeom      fegeom;
1184:   PetscReal       *coords, *gcoords;
1185:   PetscReal        localDiff = 0.0;
1186:   PetscInt         dim, coordDim, qNc = 0, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1187:   PetscBool        transform;
1188:   PetscErrorCode   ierr;

1191:   DMGetDimension(dm, &dim);
1192:   DMGetCoordinateDim(dm, &coordDim);
1193:   fegeom.dimEmbed = coordDim;
1194:   DMGetLocalSection(dm, &section);
1195:   PetscSectionGetNumFields(section, &numFields);
1196:   DMGetLocalVector(dm, &localX);
1197:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1198:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1199:   DMGetBasisTransformDM_Internal(dm, &tdm);
1200:   DMGetBasisTransformVec_Internal(dm, &tv);
1201:   DMHasBasisTransform(dm, &transform);
1202:   for (field = 0; field < numFields; ++field) {
1203:     PetscFE  fe;
1204:     PetscInt Nc;

1206:     DMGetField(dm, field, NULL, (PetscObject *) &fe);
1207:     PetscFEGetQuadrature(fe, &quad);
1208:     PetscFEGetNumComponents(fe, &Nc);
1209:     numComponents += Nc;
1210:   }
1211:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights);
1212:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1213:   /* DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
1214:   PetscMalloc6(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ,numComponents*coordDim,&interpolant,Nq,&fegeom.detJ);
1215:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1216:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1217:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1218:   for (c = cStart; c < cEnd; ++c) {
1219:     PetscScalar *x = NULL;
1220:     PetscReal    elemDiff = 0.0;
1221:     PetscInt     qc = 0;

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

1226:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1227:       PetscFE          fe;
1228:       void * const     ctx = ctxs ? ctxs[field] : NULL;
1229:       PetscInt         Nb, Nc, q, fc;

1231:       DMGetField(dm, field, NULL, (PetscObject *) &fe);
1232:       PetscFEGetDimension(fe, &Nb);
1233:       PetscFEGetNumComponents(fe, &Nc);
1234:       if (debug) {
1235:         char title[1024];
1236:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1237:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1238:       }
1239:       for (q = 0; q < Nq; ++q) {
1240:         PetscFEGeom qgeom;

1242:         qgeom.dimEmbed = fegeom.dimEmbed;
1243:         qgeom.J        = &fegeom.J[q*coordDim*coordDim];
1244:         qgeom.invJ     = &fegeom.invJ[q*coordDim*coordDim];
1245:         qgeom.detJ     = &fegeom.detJ[q];
1246:         if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", fegeom.detJ[q], c, q);
1247:         if (transform) {
1248:           gcoords = &coords[coordDim*Nq];
1249:           DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);
1250:         } else {
1251:           gcoords = &coords[coordDim*q];
1252:         }
1253:         (*funcs[field])(coordDim, time, gcoords, n, Nc, funcVal, ctx);
1254:         if (ierr) {
1255:           PetscErrorCode ierr2;
1256:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1257:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1258:           ierr2 = PetscFree6(funcVal,coords,fegeom.J,fegeom.invJ,interpolant,fegeom.detJ);CHKERRQ(ierr2);
1259: 
1260:         }
1261:         if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1262:         PetscFEInterpolateGradient_Static(fe, &x[fieldOffset], &qgeom, q, interpolant);
1263:         /* Overwrite with the dot product if the normal is given */
1264:         if (n) {
1265:           for (fc = 0; fc < Nc; ++fc) {
1266:             PetscScalar sum = 0.0;
1267:             PetscInt    d;
1268:             for (d = 0; d < dim; ++d) sum += interpolant[fc*dim+d]*n[d];
1269:             interpolant[fc] = sum;
1270:           }
1271:         }
1272:         for (fc = 0; fc < Nc; ++fc) {
1273:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1274:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %D fieldDer %D,%D diff %g\n", c, field, fc, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]);}
1275:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1276:         }
1277:       }
1278:       fieldOffset += Nb;
1279:       qc          += Nc;
1280:     }
1281:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1282:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
1283:     localDiff += elemDiff;
1284:   }
1285:   PetscFree6(funcVal,coords,fegeom.J,fegeom.invJ,interpolant,fegeom.detJ);
1286:   DMRestoreLocalVector(dm, &localX);
1287:   MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1288:   *diff = PetscSqrtReal(*diff);
1289:   return(0);
1290: }

1292: PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
1293: {
1294:   const PetscInt   debug = ((DM_Plex*)dm->data)->printL2;
1295:   DM               tdm;
1296:   PetscSection     section;
1297:   PetscQuadrature  quad;
1298:   Vec              localX, tv;
1299:   PetscFEGeom      fegeom;
1300:   PetscScalar     *funcVal, *interpolant;
1301:   PetscReal       *coords, *gcoords;
1302:   PetscReal       *localDiff;
1303:   const PetscReal *quadPoints, *quadWeights;
1304:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1305:   PetscBool        transform;
1306:   PetscErrorCode   ierr;

1309:   DMGetDimension(dm, &dim);
1310:   DMGetCoordinateDim(dm, &coordDim);
1311:   DMGetLocalSection(dm, &section);
1312:   PetscSectionGetNumFields(section, &numFields);
1313:   DMGetLocalVector(dm, &localX);
1314:   VecSet(localX, 0.0);
1315:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1316:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1317:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1318:   DMGetBasisTransformDM_Internal(dm, &tdm);
1319:   DMGetBasisTransformVec_Internal(dm, &tv);
1320:   DMHasBasisTransform(dm, &transform);
1321:   for (field = 0; field < numFields; ++field) {
1322:     PetscObject  obj;
1323:     PetscClassId id;
1324:     PetscInt     Nc;

1326:     DMGetField(dm, field, NULL, &obj);
1327:     PetscObjectGetClassId(obj, &id);
1328:     if (id == PETSCFE_CLASSID) {
1329:       PetscFE fe = (PetscFE) obj;

1331:       PetscFEGetQuadrature(fe, &quad);
1332:       PetscFEGetNumComponents(fe, &Nc);
1333:     } else if (id == PETSCFV_CLASSID) {
1334:       PetscFV fv = (PetscFV) obj;

1336:       PetscFVGetQuadrature(fv, &quad);
1337:       PetscFVGetNumComponents(fv, &Nc);
1338:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1339:     numComponents += Nc;
1340:   }
1341:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1342:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1343:   PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*(Nq+1),&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1344:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1345:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1346:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1347:   for (c = cStart; c < cEnd; ++c) {
1348:     PetscScalar *x = NULL;
1349:     PetscInt     qc = 0;

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

1354:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1355:       PetscObject  obj;
1356:       PetscClassId id;
1357:       void * const ctx = ctxs ? ctxs[field] : NULL;
1358:       PetscInt     Nb, Nc, q, fc;

1360:       PetscReal       elemDiff = 0.0;

1362:       DMGetField(dm, field, NULL, &obj);
1363:       PetscObjectGetClassId(obj, &id);
1364:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1365:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1366:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1367:       if (debug) {
1368:         char title[1024];
1369:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
1370:         DMPrintCellVector(c, title, Nb, &x[fieldOffset]);
1371:       }
1372:       for (q = 0; q < Nq; ++q) {
1373:         PetscFEGeom qgeom;

1375:         qgeom.dimEmbed = fegeom.dimEmbed;
1376:         qgeom.J        = &fegeom.J[q*coordDim*coordDim];
1377:         qgeom.invJ     = &fegeom.invJ[q*coordDim*coordDim];
1378:         qgeom.detJ     = &fegeom.detJ[q];
1379:         if (fegeom.detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature point %D", fegeom.detJ, c, q);
1380:         if (transform) {
1381:           gcoords = &coords[coordDim*Nq];
1382:           DMPlexBasisTransformApplyReal_Internal(dm, &coords[coordDim*q], PETSC_TRUE, coordDim, &coords[coordDim*q], gcoords, dm->transformCtx);
1383:         } else {
1384:           gcoords = &coords[coordDim*q];
1385:         }
1386:         (*funcs[field])(coordDim, time, gcoords, Nc, funcVal, ctx);
1387:         if (ierr) {
1388:           PetscErrorCode ierr2;
1389:           ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1390:           ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1391:           ierr2 = PetscFree7(localDiff,funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1392: 
1393:         }
1394:         if (transform) {DMPlexBasisTransformApply_Internal(dm, &coords[coordDim*q], PETSC_FALSE, Nc, funcVal, funcVal, dm->transformCtx);}
1395:         if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1396:         else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1397:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1398:         for (fc = 0; fc < Nc; ++fc) {
1399:           const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1400:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %D field %D,%D point %g %g %g diff %g\n", c, field, fc, coordDim > 0 ? coords[coordDim*q] : 0., coordDim > 1 ? coords[coordDim*q+1] : 0., coordDim > 2 ? coords[coordDim*q+2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q]);}
1401:           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1402:         }
1403:       }
1404:       fieldOffset += Nb;
1405:       qc          += Nc;
1406:       localDiff[field] += elemDiff;
1407:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d field %d cum diff %g\n", c, field, localDiff[field]);}
1408:     }
1409:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1410:   }
1411:   DMRestoreLocalVector(dm, &localX);
1412:   MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1413:   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
1414:   PetscFree7(localDiff,funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1415:   return(0);
1416: }

1418: /*@C
1419:   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.

1421:   Collective on dm

1423:   Input Parameters:
1424: + dm    - The DM
1425: . time  - The time
1426: . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation
1427: . ctxs  - Optional array of contexts to pass to each function, or NULL.
1428: - X     - The coefficient vector u_h

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

1433:   Level: developer

1435: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1436: @*/
1437: PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D)
1438: {
1439:   PetscSection     section;
1440:   PetscQuadrature  quad;
1441:   Vec              localX;
1442:   PetscFEGeom      fegeom;
1443:   PetscScalar     *funcVal, *interpolant;
1444:   PetscReal       *coords;
1445:   const PetscReal *quadPoints, *quadWeights;
1446:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset;
1447:   PetscErrorCode   ierr;

1450:   VecSet(D, 0.0);
1451:   DMGetDimension(dm, &dim);
1452:   DMGetCoordinateDim(dm, &coordDim);
1453:   DMGetLocalSection(dm, &section);
1454:   PetscSectionGetNumFields(section, &numFields);
1455:   DMGetLocalVector(dm, &localX);
1456:   DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);
1457:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1458:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1459:   for (field = 0; field < numFields; ++field) {
1460:     PetscObject  obj;
1461:     PetscClassId id;
1462:     PetscInt     Nc;

1464:     DMGetField(dm, field, NULL, &obj);
1465:     PetscObjectGetClassId(obj, &id);
1466:     if (id == PETSCFE_CLASSID) {
1467:       PetscFE fe = (PetscFE) obj;

1469:       PetscFEGetQuadrature(fe, &quad);
1470:       PetscFEGetNumComponents(fe, &Nc);
1471:     } else if (id == PETSCFV_CLASSID) {
1472:       PetscFV fv = (PetscFV) obj;

1474:       PetscFVGetQuadrature(fv, &quad);
1475:       PetscFVGetNumComponents(fv, &Nc);
1476:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1477:     numComponents += Nc;
1478:   }
1479:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1480:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1481:   PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1482:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1483:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1484:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1485:   for (c = cStart; c < cEnd; ++c) {
1486:     PetscScalar *x = NULL;
1487:     PetscScalar  elemDiff = 0.0;
1488:     PetscInt     qc = 0;

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

1493:     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1494:       PetscObject  obj;
1495:       PetscClassId id;
1496:       void * const ctx = ctxs ? ctxs[field] : NULL;
1497:       PetscInt     Nb, Nc, q, fc;

1499:       DMGetField(dm, field, NULL, &obj);
1500:       PetscObjectGetClassId(obj, &id);
1501:       if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1502:       else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1503:       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1504:       if (funcs[field]) {
1505:         for (q = 0; q < Nq; ++q) {
1506:           PetscFEGeom qgeom;

1508:           qgeom.dimEmbed = fegeom.dimEmbed;
1509:           qgeom.J        = &fegeom.J[q*coordDim*coordDim];
1510:           qgeom.invJ     = &fegeom.invJ[q*coordDim*coordDim];
1511:           qgeom.detJ     = &fegeom.detJ[q];
1512:           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);
1513:           (*funcs[field])(coordDim, time, &coords[q*coordDim], Nc, funcVal, ctx);
1514:           if (ierr) {
1515:             PetscErrorCode ierr2;
1516:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2);
1517:             ierr2 = PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1518:             ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2);
1519: 
1520:           }
1521:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1522:           else if (id == PETSCFV_CLASSID) {PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);}
1523:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1524:           for (fc = 0; fc < Nc; ++fc) {
1525:             const PetscReal wt = quadWeights[q*qNc+(qNc == 1 ? 0 : qc+fc)];
1526:             elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*wt*fegeom.detJ[q];
1527:           }
1528:         }
1529:       }
1530:       fieldOffset += Nb;
1531:       qc          += Nc;
1532:     }
1533:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
1534:     VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);
1535:   }
1536:   PetscFree6(funcVal,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1537:   DMRestoreLocalVector(dm, &localX);
1538:   VecSqrtAbs(D);
1539:   return(0);
1540: }

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

1545:   Collective on dm

1547:   Input Parameters:
1548: + dm - The DM
1549: - LocX  - The coefficient vector u_h

1551:   Output Parameter:
1552: . locC - A Vec which holds the Clement interpolant of the gradient

1554:   Notes:
1555:     Add citation to (Clement, 1975) and definition of the interpolant
1556:   \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

1558:   Level: developer

1560: .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff()
1561: @*/
1562: PetscErrorCode DMPlexComputeGradientClementInterpolant(DM dm, Vec locX, Vec locC)
1563: {
1564:   DM_Plex         *mesh  = (DM_Plex *) dm->data;
1565:   PetscInt         debug = mesh->printFEM;
1566:   DM               dmC;
1567:   PetscSection     section;
1568:   PetscQuadrature  quad;
1569:   PetscScalar     *interpolant, *gradsum;
1570:   PetscFEGeom      fegeom;
1571:   PetscReal       *coords;
1572:   const PetscReal *quadPoints, *quadWeights;
1573:   PetscInt         dim, coordDim, numFields, numComponents = 0, qNc, Nq, cStart, cEnd, cEndInterior, vStart, vEnd, v, field, fieldOffset;
1574:   PetscErrorCode   ierr;

1577:   VecGetDM(locC, &dmC);
1578:   VecSet(locC, 0.0);
1579:   DMGetDimension(dm, &dim);
1580:   DMGetCoordinateDim(dm, &coordDim);
1581:   fegeom.dimEmbed = coordDim;
1582:   DMGetLocalSection(dm, &section);
1583:   PetscSectionGetNumFields(section, &numFields);
1584:   for (field = 0; field < numFields; ++field) {
1585:     PetscObject  obj;
1586:     PetscClassId id;
1587:     PetscInt     Nc;

1589:     DMGetField(dm, field, NULL, &obj);
1590:     PetscObjectGetClassId(obj, &id);
1591:     if (id == PETSCFE_CLASSID) {
1592:       PetscFE fe = (PetscFE) obj;

1594:       PetscFEGetQuadrature(fe, &quad);
1595:       PetscFEGetNumComponents(fe, &Nc);
1596:     } else if (id == PETSCFV_CLASSID) {
1597:       PetscFV fv = (PetscFV) obj;

1599:       PetscFVGetQuadrature(fv, &quad);
1600:       PetscFVGetNumComponents(fv, &Nc);
1601:     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1602:     numComponents += Nc;
1603:   }
1604:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1605:   if ((qNc != 1) && (qNc != numComponents)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_SIZ, "Quadrature components %D != %D field components", qNc, numComponents);
1606:   PetscMalloc6(coordDim*numComponents*2,&gradsum,coordDim*numComponents,&interpolant,coordDim*Nq,&coords,Nq,&fegeom.detJ,coordDim*coordDim*Nq,&fegeom.J,coordDim*coordDim*Nq,&fegeom.invJ);
1607:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1608:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1609:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1610:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1611:   for (v = vStart; v < vEnd; ++v) {
1612:     PetscScalar volsum = 0.0;
1613:     PetscInt   *star = NULL;
1614:     PetscInt    starSize, st, d, fc;

1616:     PetscArrayzero(gradsum, coordDim*numComponents);
1617:     DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1618:     for (st = 0; st < starSize*2; st += 2) {
1619:       const PetscInt cell = star[st];
1620:       PetscScalar   *grad = &gradsum[coordDim*numComponents];
1621:       PetscScalar   *x    = NULL;
1622:       PetscReal      vol  = 0.0;

1624:       if ((cell < cStart) || (cell >= cEnd)) continue;
1625:       DMPlexComputeCellGeometryFEM(dm, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ);
1626:       DMPlexVecGetClosure(dm, NULL, locX, cell, NULL, &x);
1627:       for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1628:         PetscObject  obj;
1629:         PetscClassId id;
1630:         PetscInt     Nb, Nc, q, qc = 0;

1632:         PetscArrayzero(grad, coordDim*numComponents);
1633:         DMGetField(dm, field, NULL, &obj);
1634:         PetscObjectGetClassId(obj, &id);
1635:         if (id == PETSCFE_CLASSID)      {PetscFEGetNumComponents((PetscFE) obj, &Nc);PetscFEGetDimension((PetscFE) obj, &Nb);}
1636:         else if (id == PETSCFV_CLASSID) {PetscFVGetNumComponents((PetscFV) obj, &Nc);Nb = 1;}
1637:         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1638:         for (q = 0; q < Nq; ++q) {
1639:           PetscFEGeom qgeom;

1641:           qgeom.dimEmbed = fegeom.dimEmbed;
1642:           qgeom.J        = &fegeom.J[q*coordDim*coordDim];
1643:           qgeom.invJ     = &fegeom.invJ[q*coordDim*coordDim];
1644:           qgeom.detJ     = &fegeom.detJ[q];
1645:           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);
1646:           if (ierr) {
1647:             PetscErrorCode ierr2;
1648:             ierr2 = DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);CHKERRQ(ierr2);
1649:             ierr2 = DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr2);
1650:             ierr2 = PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);CHKERRQ(ierr2);
1651: 
1652:           }
1653:           if (id == PETSCFE_CLASSID)      {PetscFEInterpolateGradient_Static((PetscFE) obj, &x[fieldOffset], &qgeom, q, interpolant);}
1654:           else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1655:           for (fc = 0; fc < Nc; ++fc) {
1656:             const PetscReal wt = quadWeights[q*qNc+qc+fc];

1658:             for (d = 0; d < coordDim; ++d) grad[fc*coordDim+d] += interpolant[fc*dim+d]*wt*fegeom.detJ[q];
1659:           }
1660:           vol += quadWeights[q*qNc]*fegeom.detJ[q];
1661:         }
1662:         fieldOffset += Nb;
1663:         qc          += Nc;
1664:       }
1665:       DMPlexVecRestoreClosure(dm, NULL, locX, cell, NULL, &x);
1666:       for (fc = 0; fc < numComponents; ++fc) {
1667:         for (d = 0; d < coordDim; ++d) {
1668:           gradsum[fc*coordDim+d] += grad[fc*coordDim+d];
1669:         }
1670:       }
1671:       volsum += vol;
1672:       if (debug) {
1673:         PetscPrintf(PETSC_COMM_SELF, "Cell %D gradient: [", cell);
1674:         for (fc = 0; fc < numComponents; ++fc) {
1675:           for (d = 0; d < coordDim; ++d) {
1676:             if (fc || d > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
1677:             PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(grad[fc*coordDim+d]));
1678:           }
1679:         }
1680:         PetscPrintf(PETSC_COMM_SELF, "]\n");
1681:       }
1682:     }
1683:     for (fc = 0; fc < numComponents; ++fc) {
1684:       for (d = 0; d < coordDim; ++d) gradsum[fc*coordDim+d] /= volsum;
1685:     }
1686:     DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &starSize, &star);
1687:     DMPlexVecSetClosure(dmC, NULL, locC, v, gradsum, INSERT_VALUES);
1688:   }
1689:   PetscFree6(gradsum,interpolant,coords,fegeom.detJ,fegeom.J,fegeom.invJ);
1690:   return(0);
1691: }

1693: static PetscErrorCode DMPlexComputeIntegral_Internal(DM dm, Vec X, PetscInt cStart, PetscInt cEnd, PetscScalar *cintegral, void *user)
1694: {
1695:   DM                 dmAux = NULL;
1696:   PetscDS            prob,    probAux = NULL;
1697:   PetscSection       section, sectionAux;
1698:   Vec                locX,    locA;
1699:   PetscInt           dim, numCells = cEnd - cStart, c, f;
1700:   PetscBool          useFVM = PETSC_FALSE;
1701:   /* DS */
1702:   PetscInt           Nf,    totDim,    *uOff, *uOff_x, numConstants;
1703:   PetscInt           NfAux, totDimAux, *aOff;
1704:   PetscScalar       *u, *a;
1705:   const PetscScalar *constants;
1706:   /* Geometry */
1707:   PetscFEGeom       *cgeomFEM;
1708:   DM                 dmGrad;
1709:   PetscQuadrature    affineQuad = NULL;
1710:   Vec                cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL;
1711:   PetscFVCellGeom   *cgeomFVM;
1712:   const PetscScalar *lgrad;
1713:   PetscInt           maxDegree;
1714:   DMField            coordField;
1715:   IS                 cellIS;
1716:   PetscErrorCode     ierr;

1719:   DMGetDS(dm, &prob);
1720:   DMGetDimension(dm, &dim);
1721:   DMGetLocalSection(dm, &section);
1722:   PetscSectionGetNumFields(section, &Nf);
1723:   /* Determine which discretizations we have */
1724:   for (f = 0; f < Nf; ++f) {
1725:     PetscObject  obj;
1726:     PetscClassId id;

1728:     PetscDSGetDiscretization(prob, f, &obj);
1729:     PetscObjectGetClassId(obj, &id);
1730:     if (id == PETSCFV_CLASSID) useFVM = PETSC_TRUE;
1731:   }
1732:   /* Get local solution with boundary values */
1733:   DMGetLocalVector(dm, &locX);
1734:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
1735:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
1736:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
1737:   /* Read DS information */
1738:   PetscDSGetTotalDimension(prob, &totDim);
1739:   PetscDSGetComponentOffsets(prob, &uOff);
1740:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
1741:   ISCreateStride(PETSC_COMM_SELF,numCells,cStart,1,&cellIS);
1742:   PetscDSGetConstants(prob, &numConstants, &constants);
1743:   /* Read Auxiliary DS information */
1744:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1745:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1746:   if (dmAux) {
1747:     DMGetDS(dmAux, &probAux);
1748:     PetscDSGetNumFields(probAux, &NfAux);
1749:     DMGetLocalSection(dmAux, &sectionAux);
1750:     PetscDSGetTotalDimension(probAux, &totDimAux);
1751:     PetscDSGetComponentOffsets(probAux, &aOff);
1752:   }
1753:   /* Allocate data  arrays */
1754:   PetscCalloc1(numCells*totDim, &u);
1755:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1756:   /* Read out geometry */
1757:   DMGetCoordinateField(dm,&coordField);
1758:   DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
1759:   if (maxDegree <= 1) {
1760:     DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
1761:     if (affineQuad) {
1762:       DMFieldCreateFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&cgeomFEM);
1763:     }
1764:   }
1765:   if (useFVM) {
1766:     PetscFV   fv = NULL;
1767:     Vec       grad;
1768:     PetscInt  fStart, fEnd;
1769:     PetscBool compGrad;

1771:     for (f = 0; f < Nf; ++f) {
1772:       PetscObject  obj;
1773:       PetscClassId id;

1775:       PetscDSGetDiscretization(prob, f, &obj);
1776:       PetscObjectGetClassId(obj, &id);
1777:       if (id == PETSCFV_CLASSID) {fv = (PetscFV) obj; break;}
1778:     }
1779:     PetscFVGetComputeGradients(fv, &compGrad);
1780:     PetscFVSetComputeGradients(fv, PETSC_TRUE);
1781:     DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);
1782:     DMPlexComputeGradientFVM(dm, fv, faceGeometryFVM, cellGeometryFVM, &dmGrad);
1783:     PetscFVSetComputeGradients(fv, compGrad);
1784:     VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1785:     /* Reconstruct and limit cell gradients */
1786:     DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1787:     DMGetGlobalVector(dmGrad, &grad);
1788:     DMPlexReconstructGradients_Internal(dm, fv, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1789:     /* Communicate gradient values */
1790:     DMGetLocalVector(dmGrad, &locGrad);
1791:     DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1792:     DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1793:     DMRestoreGlobalVector(dmGrad, &grad);
1794:     /* Handle non-essential (e.g. outflow) boundary values */
1795:     DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);
1796:     VecGetArrayRead(locGrad, &lgrad);
1797:   }
1798:   /* Read out data from inputs */
1799:   for (c = cStart; c < cEnd; ++c) {
1800:     PetscScalar *x = NULL;
1801:     PetscInt     i;

1803:     DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1804:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1805:     DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1806:     if (dmAux) {
1807:       DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1808:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1809:       DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1810:     }
1811:   }
1812:   /* Do integration for each field */
1813:   for (f = 0; f < Nf; ++f) {
1814:     PetscObject  obj;
1815:     PetscClassId id;
1816:     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

1818:     PetscDSGetDiscretization(prob, f, &obj);
1819:     PetscObjectGetClassId(obj, &id);
1820:     if (id == PETSCFE_CLASSID) {
1821:       PetscFE         fe = (PetscFE) obj;
1822:       PetscQuadrature q;
1823:       PetscFEGeom     *chunkGeom = NULL;
1824:       PetscInt        Nq, Nb;

1826:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1827:       PetscFEGetQuadrature(fe, &q);
1828:       PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1829:       PetscFEGetDimension(fe, &Nb);
1830:       blockSize = Nb*Nq;
1831:       batchSize = numBlocks * blockSize;
1832:       PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1833:       numChunks = numCells / (numBatches*batchSize);
1834:       Ne        = numChunks*numBatches*batchSize;
1835:       Nr        = numCells % (numBatches*batchSize);
1836:       offset    = numCells - Nr;
1837:       if (!affineQuad) {
1838:         DMFieldCreateFEGeom(coordField,cellIS,q,PETSC_FALSE,&cgeomFEM);
1839:       }
1840:       PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
1841:       PetscFEIntegrate(prob, f, Ne, chunkGeom, u, probAux, a, cintegral);
1842:       PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&chunkGeom);
1843:       PetscFEIntegrate(prob, f, Nr, chunkGeom, &u[offset*totDim], probAux, &a[offset*totDimAux], &cintegral[offset*Nf]);
1844:       PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&chunkGeom);
1845:       if (!affineQuad) {
1846:         PetscFEGeomDestroy(&cgeomFEM);
1847:       }
1848:     } else if (id == PETSCFV_CLASSID) {
1849:       PetscInt       foff;
1850:       PetscPointFunc obj_func;
1851:       PetscScalar    lint;

1853:       PetscDSGetObjective(prob, f, &obj_func);
1854:       PetscDSGetFieldOffset(prob, f, &foff);
1855:       if (obj_func) {
1856:         for (c = 0; c < numCells; ++c) {
1857:           PetscScalar *u_x;

1859:           DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);
1860:           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);
1861:           cintegral[c*Nf+f] += PetscRealPart(lint)*cgeomFVM[c].volume;
1862:         }
1863:       }
1864:     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1865:   }
1866:   /* Cleanup data arrays */
1867:   if (useFVM) {
1868:     VecRestoreArrayRead(locGrad, &lgrad);
1869:     VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1870:     DMRestoreLocalVector(dmGrad, &locGrad);
1871:     VecDestroy(&faceGeometryFVM);
1872:     VecDestroy(&cellGeometryFVM);
1873:     DMDestroy(&dmGrad);
1874:   }
1875:   if (dmAux) {PetscFree(a);}
1876:   PetscFree(u);
1877:   /* Cleanup */
1878:   if (affineQuad) {
1879:     PetscFEGeomDestroy(&cgeomFEM);
1880:   }
1881:   PetscQuadratureDestroy(&affineQuad);
1882:   ISDestroy(&cellIS);
1883:   DMRestoreLocalVector(dm, &locX);
1884:   return(0);
1885: }

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

1890:   Input Parameters:
1891: + dm - The mesh
1892: . X  - Global input vector
1893: - user - The user context

1895:   Output Parameter:
1896: . integral - Integral for each field

1898:   Level: developer

1900: .seealso: DMPlexComputeResidualFEM()
1901: @*/
1902: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
1903: {
1904:   DM_Plex       *mesh = (DM_Plex *) dm->data;
1905:   PetscScalar   *cintegral, *lintegral;
1906:   PetscInt       Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;

1913:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1914:   DMGetNumFields(dm, &Nf);
1915:   DMPlexGetVTKCellHeight(dm, &cellHeight);
1916:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1917:   DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1918:   cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1919:   /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1920:   PetscCalloc2(Nf, &lintegral, (cEnd-cStart)*Nf, &cintegral);
1921:   DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1922:   /* Sum up values */
1923:   for (cell = cStart; cell < cEnd; ++cell) {
1924:     const PetscInt c = cell - cStart;

1926:     if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1927:     for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c*Nf+f];
1928:   }
1929:   MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1930:   if (mesh->printFEM) {
1931:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Integral:");
1932:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", (double) PetscRealPart(integral[f]));}
1933:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
1934:   }
1935:   PetscFree2(lintegral, cintegral);
1936:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1937:   return(0);
1938: }

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

1943:   Input Parameters:
1944: + dm - The mesh
1945: . X  - Global input vector
1946: - user - The user context

1948:   Output Parameter:
1949: . integral - Cellwise integrals for each field

1951:   Level: developer

1953: .seealso: DMPlexComputeResidualFEM()
1954: @*/
1955: PetscErrorCode DMPlexComputeCellwiseIntegralFEM(DM dm, Vec X, Vec F, void *user)
1956: {
1957:   DM_Plex       *mesh = (DM_Plex *) dm->data;
1958:   DM             dmF;
1959:   PetscSection   sectionF;
1960:   PetscScalar   *cintegral, *af;
1961:   PetscInt       Nf, f, cellHeight, cStart, cEnd, cEndInterior[4], cell;

1968:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
1969:   DMGetNumFields(dm, &Nf);
1970:   DMPlexGetVTKCellHeight(dm, &cellHeight);
1971:   DMPlexGetHeightStratum(dm, cellHeight, &cStart, &cEnd);
1972:   DMPlexGetHybridBounds(dm, &cEndInterior[0], &cEndInterior[1], &cEndInterior[2], &cEndInterior[3]);
1973:   cEnd = cEndInterior[cellHeight] < 0 ? cEnd : cEndInterior[cellHeight];
1974:   /* TODO Introduce a loop over large chunks (right now this is a single chunk) */
1975:   PetscCalloc1((cEnd-cStart)*Nf, &cintegral);
1976:   DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user);
1977:   /* Put values in F*/
1978:   VecGetDM(F, &dmF);
1979:   DMGetLocalSection(dmF, &sectionF);
1980:   VecGetArray(F, &af);
1981:   for (cell = cStart; cell < cEnd; ++cell) {
1982:     const PetscInt c = cell - cStart;
1983:     PetscInt       dof, off;

1985:     if (mesh->printFEM > 1) {DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c*Nf]);}
1986:     PetscSectionGetDof(sectionF, cell, &dof);
1987:     PetscSectionGetOffset(sectionF, cell, &off);
1988:     if (dof != Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "The number of cell dofs %D != %D", dof, Nf);
1989:     for (f = 0; f < Nf; ++f) af[off+f] = cintegral[c*Nf+f];
1990:   }
1991:   VecRestoreArray(F, &af);
1992:   PetscFree(cintegral);
1993:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
1994:   return(0);
1995: }

1997: static PetscErrorCode DMPlexComputeBdIntegral_Internal(DM dm, Vec locX, IS pointIS,
1998:                                                        void (*func)(PetscInt, PetscInt, PetscInt,
1999:                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2000:                                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2001:                                                                     PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
2002:                                                        PetscScalar *fintegral, void *user)
2003: {
2004:   DM                 plex = NULL, plexA = NULL;
2005:   PetscDS            prob, probAux = NULL;
2006:   PetscSection       section, sectionAux = NULL;
2007:   Vec                locA = NULL;
2008:   DMField            coordField;
2009:   PetscInt           Nf,        totDim,        *uOff, *uOff_x;
2010:   PetscInt           NfAux = 0, totDimAux = 0, *aOff = NULL;
2011:   PetscScalar       *u, *a = NULL;
2012:   const PetscScalar *constants;
2013:   PetscInt           numConstants, f;
2014:   PetscErrorCode     ierr;

2017:   DMGetCoordinateField(dm, &coordField);
2018:   DMConvert(dm, DMPLEX, &plex);
2019:   DMGetDS(dm, &prob);
2020:   DMGetLocalSection(dm, &section);
2021:   PetscSectionGetNumFields(section, &Nf);
2022:   /* Determine which discretizations we have */
2023:   for (f = 0; f < Nf; ++f) {
2024:     PetscObject  obj;
2025:     PetscClassId id;

2027:     PetscDSGetDiscretization(prob, f, &obj);
2028:     PetscObjectGetClassId(obj, &id);
2029:     if (id == PETSCFV_CLASSID) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Not supported for FVM (field %D)", f);
2030:   }
2031:   /* Read DS information */
2032:   PetscDSGetTotalDimension(prob, &totDim);
2033:   PetscDSGetComponentOffsets(prob, &uOff);
2034:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
2035:   PetscDSGetConstants(prob, &numConstants, &constants);
2036:   /* Read Auxiliary DS information */
2037:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
2038:   if (locA) {
2039:     DM dmAux;

2041:     VecGetDM(locA, &dmAux);
2042:     DMConvert(dmAux, DMPLEX, &plexA);
2043:     DMGetDS(dmAux, &probAux);
2044:     PetscDSGetNumFields(probAux, &NfAux);
2045:     DMGetLocalSection(dmAux, &sectionAux);
2046:     PetscDSGetTotalDimension(probAux, &totDimAux);
2047:     PetscDSGetComponentOffsets(probAux, &aOff);
2048:   }
2049:   /* Integrate over points */
2050:   {
2051:     PetscFEGeom    *fgeom, *chunkGeom = NULL;
2052:     PetscInt        maxDegree;
2053:     PetscQuadrature qGeom = NULL;
2054:     const PetscInt *points;
2055:     PetscInt        numFaces, face, Nq, field;
2056:     PetscInt        numChunks, chunkSize, chunk, Nr, offset;

2058:     ISGetLocalSize(pointIS, &numFaces);
2059:     ISGetIndices(pointIS, &points);
2060:     PetscCalloc2(numFaces*totDim, &u, locA ? numFaces*totDimAux : 0, &a);
2061:     DMFieldGetDegree(coordField, pointIS, NULL, &maxDegree);
2062:     for (field = 0; field < Nf; ++field) {
2063:       PetscFE fe;

2065:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
2066:       if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, pointIS, &qGeom);}
2067:       if (!qGeom) {
2068:         PetscFEGetFaceQuadrature(fe, &qGeom);
2069:         PetscObjectReference((PetscObject) qGeom);
2070:       }
2071:       PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
2072:       DMPlexGetFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
2073:       for (face = 0; face < numFaces; ++face) {
2074:         const PetscInt point = points[face], *support, *cone;
2075:         PetscScalar    *x    = NULL;
2076:         PetscInt       i, coneSize, faceLoc;

2078:         DMPlexGetSupport(dm, point, &support);
2079:         DMPlexGetConeSize(dm, support[0], &coneSize);
2080:         DMPlexGetCone(dm, support[0], &cone);
2081:         for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break;
2082:         if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of support[0] %D", face, support[0]);
2083:         fgeom->face[face][0] = faceLoc;
2084:         DMPlexVecGetClosure(plex, section, locX, support[0], NULL, &x);
2085:         for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
2086:         DMPlexVecRestoreClosure(plex, section, locX, support[0], NULL, &x);
2087:         if (locA) {
2088:           PetscInt subp;
2089:           DMPlexGetSubpoint(plexA, support[0], &subp);
2090:           DMPlexVecGetClosure(plexA, sectionAux, locA, subp, NULL, &x);
2091:           for (i = 0; i < totDimAux; ++i) a[f*totDimAux+i] = x[i];
2092:           DMPlexVecRestoreClosure(plexA, sectionAux, locA, subp, NULL, &x);
2093:         }
2094:       }
2095:       /* Get blocking */
2096:       {
2097:         PetscQuadrature q;
2098:         PetscInt        numBatches, batchSize, numBlocks, blockSize;
2099:         PetscInt        Nq, Nb;

2101:         PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2102:         PetscFEGetQuadrature(fe, &q);
2103:         PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
2104:         PetscFEGetDimension(fe, &Nb);
2105:         blockSize = Nb*Nq;
2106:         batchSize = numBlocks * blockSize;
2107:         chunkSize = numBatches*batchSize;
2108:         PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2109:         numChunks = numFaces / chunkSize;
2110:         Nr        = numFaces % chunkSize;
2111:         offset    = numFaces - Nr;
2112:       }
2113:       /* Do integration for each field */
2114:       for (chunk = 0; chunk < numChunks; ++chunk) {
2115:         PetscFEGeomGetChunk(fgeom, chunk*chunkSize, (chunk+1)*chunkSize, &chunkGeom);
2116:         PetscFEIntegrateBd(prob, field, func, chunkSize, chunkGeom, u, probAux, a, fintegral);
2117:         PetscFEGeomRestoreChunk(fgeom, 0, offset, &chunkGeom);
2118:       }
2119:       PetscFEGeomGetChunk(fgeom, offset, numFaces, &chunkGeom);
2120:       PetscFEIntegrateBd(prob, field, func, Nr, chunkGeom, &u[offset*totDim], probAux, a ? &a[offset*totDimAux] : NULL, &fintegral[offset*Nf]);
2121:       PetscFEGeomRestoreChunk(fgeom, offset, numFaces, &chunkGeom);
2122:       /* Cleanup data arrays */
2123:       DMPlexRestoreFEGeom(coordField, pointIS, qGeom, PETSC_TRUE, &fgeom);
2124:       PetscQuadratureDestroy(&qGeom);
2125:       PetscFree2(u, a);
2126:       ISRestoreIndices(pointIS, &points);
2127:     }
2128:   }
2129:   if (plex)  {DMDestroy(&plex);}
2130:   if (plexA) {DMDestroy(&plexA);}
2131:   return(0);
2132: }

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

2137:   Input Parameters:
2138: + dm      - The mesh
2139: . X       - Global input vector
2140: . label   - The boundary DMLabel
2141: . numVals - The number of label values to use, or PETSC_DETERMINE for all values
2142: . vals    - The label values to use, or PETSC_NULL for all values
2143: . func    = The function to integrate along the boundary
2144: - user    - The user context

2146:   Output Parameter:
2147: . integral - Integral for each field

2149:   Level: developer

2151: .seealso: DMPlexComputeIntegralFEM(), DMPlexComputeBdResidualFEM()
2152: @*/
2153: PetscErrorCode DMPlexComputeBdIntegral(DM dm, Vec X, DMLabel label, PetscInt numVals, const PetscInt vals[],
2154:                                        void (*func)(PetscInt, PetscInt, PetscInt,
2155:                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2156:                                                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
2157:                                                     PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
2158:                                        PetscScalar *integral, void *user)
2159: {
2160:   Vec            locX;
2161:   PetscSection   section;
2162:   DMLabel        depthLabel;
2163:   IS             facetIS;
2164:   PetscInt       dim, Nf, f, v;

2173:   PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);
2174:   DMPlexGetDepthLabel(dm, &depthLabel);
2175:   DMGetDimension(dm, &dim);
2176:   DMLabelGetStratumIS(depthLabel, dim-1, &facetIS);
2177:   DMGetLocalSection(dm, &section);
2178:   PetscSectionGetNumFields(section, &Nf);
2179:   /* Get local solution with boundary values */
2180:   DMGetLocalVector(dm, &locX);
2181:   DMPlexInsertBoundaryValues(dm, PETSC_TRUE, locX, 0.0, NULL, NULL, NULL);
2182:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
2183:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
2184:   /* Loop over label values */
2185:   PetscArrayzero(integral, Nf);
2186:   for (v = 0; v < numVals; ++v) {
2187:     IS           pointIS;
2188:     PetscInt     numFaces, face;
2189:     PetscScalar *fintegral;

2191:     DMLabelGetStratumIS(label, vals[v], &pointIS);
2192:     if (!pointIS) continue; /* No points with that id on this process */
2193:     {
2194:       IS isectIS;

2196:       /* TODO: Special cases of ISIntersect where it is quick to check a priori if one is a superset of the other */
2197:       ISIntersect_Caching_Internal(facetIS, pointIS, &isectIS);
2198:       ISDestroy(&pointIS);
2199:       pointIS = isectIS;
2200:     }
2201:     ISGetLocalSize(pointIS, &numFaces);
2202:     PetscCalloc1(numFaces*Nf, &fintegral);
2203:     DMPlexComputeBdIntegral_Internal(dm, locX, pointIS, func, fintegral, user);
2204:     /* Sum point contributions into integral */
2205:     for (f = 0; f < Nf; ++f) for (face = 0; face < numFaces; ++face) integral[f] += fintegral[face*Nf+f];
2206:     PetscFree(fintegral);
2207:     ISDestroy(&pointIS);
2208:   }
2209:   DMRestoreLocalVector(dm, &locX);
2210:   ISDestroy(&facetIS);
2211:   PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);
2212:   return(0);
2213: }

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

2218:   Input Parameters:
2219: + dmf  - The fine mesh
2220: . dmc  - The coarse mesh
2221: - user - The user context

2223:   Output Parameter:
2224: . In  - The interpolation matrix

2226:   Level: developer

2228: .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2229: @*/
2230: PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user)
2231: {
2232:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
2233:   const char       *name  = "Interpolator";
2234:   PetscDS           prob;
2235:   PetscFE          *feRef;
2236:   PetscFV          *fvRef;
2237:   PetscSection      fsection, fglobalSection;
2238:   PetscSection      csection, cglobalSection;
2239:   PetscScalar      *elemMat;
2240:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
2241:   PetscInt          cTotDim, rTotDim = 0;
2242:   PetscErrorCode    ierr;

2245:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2246:   DMGetDimension(dmf, &dim);
2247:   DMGetLocalSection(dmf, &fsection);
2248:   DMGetGlobalSection(dmf, &fglobalSection);
2249:   DMGetLocalSection(dmc, &csection);
2250:   DMGetGlobalSection(dmc, &cglobalSection);
2251:   PetscSectionGetNumFields(fsection, &Nf);
2252:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
2253:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
2254:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2255:   DMGetDS(dmf, &prob);
2256:   PetscCalloc2(Nf,&feRef,Nf,&fvRef);
2257:   for (f = 0; f < Nf; ++f) {
2258:     PetscObject  obj;
2259:     PetscClassId id;
2260:     PetscInt     rNb = 0, Nc = 0;

2262:     PetscDSGetDiscretization(prob, f, &obj);
2263:     PetscObjectGetClassId(obj, &id);
2264:     if (id == PETSCFE_CLASSID) {
2265:       PetscFE fe = (PetscFE) obj;

2267:       PetscFERefine(fe, &feRef[f]);
2268:       PetscFEGetDimension(feRef[f], &rNb);
2269:       PetscFEGetNumComponents(fe, &Nc);
2270:     } else if (id == PETSCFV_CLASSID) {
2271:       PetscFV        fv = (PetscFV) obj;
2272:       PetscDualSpace Q;

2274:       PetscFVRefine(fv, &fvRef[f]);
2275:       PetscFVGetDualSpace(fvRef[f], &Q);
2276:       PetscDualSpaceGetDimension(Q, &rNb);
2277:       PetscFVGetNumComponents(fv, &Nc);
2278:     }
2279:     rTotDim += rNb;
2280:   }
2281:   PetscDSGetTotalDimension(prob, &cTotDim);
2282:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
2283:   PetscArrayzero(elemMat, rTotDim*cTotDim);
2284:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
2285:     PetscDualSpace   Qref;
2286:     PetscQuadrature  f;
2287:     const PetscReal *qpoints, *qweights;
2288:     PetscReal       *points;
2289:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

2291:     /* Compose points from all dual basis functionals */
2292:     if (feRef[fieldI]) {
2293:       PetscFEGetDualSpace(feRef[fieldI], &Qref);
2294:       PetscFEGetNumComponents(feRef[fieldI], &Nc);
2295:     } else {
2296:       PetscFVGetDualSpace(fvRef[fieldI], &Qref);
2297:       PetscFVGetNumComponents(fvRef[fieldI], &Nc);
2298:     }
2299:     PetscDualSpaceGetDimension(Qref, &fpdim);
2300:     for (i = 0; i < fpdim; ++i) {
2301:       PetscDualSpaceGetFunctional(Qref, i, &f);
2302:       PetscQuadratureGetData(f, NULL, NULL, &Np, NULL, NULL);
2303:       npoints += Np;
2304:     }
2305:     PetscMalloc1(npoints*dim,&points);
2306:     for (i = 0, k = 0; i < fpdim; ++i) {
2307:       PetscDualSpaceGetFunctional(Qref, i, &f);
2308:       PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2309:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
2310:     }

2312:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
2313:       PetscObject  obj;
2314:       PetscClassId id;
2315:       PetscReal   *B;
2316:       PetscInt     NcJ = 0, cpdim = 0, j, qNc;

2318:       PetscDSGetDiscretization(prob, fieldJ, &obj);
2319:       PetscObjectGetClassId(obj, &id);
2320:       if (id == PETSCFE_CLASSID) {
2321:         PetscFE fe = (PetscFE) obj;

2323:         /* Evaluate basis at points */
2324:         PetscFEGetNumComponents(fe, &NcJ);
2325:         PetscFEGetDimension(fe, &cpdim);
2326:         /* For now, fields only interpolate themselves */
2327:         if (fieldI == fieldJ) {
2328:           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);
2329:           PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
2330:           for (i = 0, k = 0; i < fpdim; ++i) {
2331:             PetscDualSpaceGetFunctional(Qref, i, &f);
2332:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
2333:             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);
2334:             for (p = 0; p < Np; ++p, ++k) {
2335:               for (j = 0; j < cpdim; ++j) {
2336:                 /*
2337:                    cTotDim:            Total columns in element interpolation matrix, sum of number of dual basis functionals in each field
2338:                    offsetI, offsetJ:   Offsets into the larger element interpolation matrix for different fields
2339:                    fpdim, i, cpdim, j: Dofs for fine and coarse grids, correspond to dual space basis functionals
2340:                    qNC, Nc, Ncj, c:    Number of components in this field
2341:                    Np, p:              Number of quad points in the fine grid functional i
2342:                    k:                  i*Np + p, overall point number for the interpolation
2343:                 */
2344:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p*qNc+c];
2345:               }
2346:             }
2347:           }
2348:           PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
2349:         }
2350:       } else if (id == PETSCFV_CLASSID) {
2351:         PetscFV        fv = (PetscFV) obj;

2353:         /* Evaluate constant function at points */
2354:         PetscFVGetNumComponents(fv, &NcJ);
2355:         cpdim = 1;
2356:         /* For now, fields only interpolate themselves */
2357:         if (fieldI == fieldJ) {
2358:           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);
2359:           for (i = 0, k = 0; i < fpdim; ++i) {
2360:             PetscDualSpaceGetFunctional(Qref, i, &f);
2361:             PetscQuadratureGetData(f, NULL, &qNc, &Np, NULL, &qweights);
2362:             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);
2363:             for (p = 0; p < Np; ++p, ++k) {
2364:               for (j = 0; j < cpdim; ++j) {
2365:                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i)*cTotDim + offsetJ + j] += 1.0*qweights[p*qNc+c];
2366:               }
2367:             }
2368:           }
2369:         }
2370:       }
2371:       offsetJ += cpdim;
2372:     }
2373:     offsetI += fpdim;
2374:     PetscFree(points);
2375:   }
2376:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
2377:   /* Preallocate matrix */
2378:   {
2379:     Mat          preallocator;
2380:     PetscScalar *vals;
2381:     PetscInt    *cellCIndices, *cellFIndices;
2382:     PetscInt     locRows, locCols, cell;

2384:     MatGetLocalSize(In, &locRows, &locCols);
2385:     MatCreate(PetscObjectComm((PetscObject) In), &preallocator);
2386:     MatSetType(preallocator, MATPREALLOCATOR);
2387:     MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);
2388:     MatSetUp(preallocator);
2389:     PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
2390:     for (cell = cStart; cell < cEnd; ++cell) {
2391:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
2392:       MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);
2393:     }
2394:     PetscFree3(vals,cellCIndices,cellFIndices);
2395:     MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);
2396:     MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);
2397:     MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);
2398:     MatDestroy(&preallocator);
2399:   }
2400:   /* Fill matrix */
2401:   MatZeroEntries(In);
2402:   for (c = cStart; c < cEnd; ++c) {
2403:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
2404:   }
2405:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
2406:   PetscFree2(feRef,fvRef);
2407:   PetscFree(elemMat);
2408:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2409:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2410:   if (mesh->printFEM) {
2411:     PetscPrintf(PetscObjectComm((PetscObject)In), "%s:\n", name);
2412:     MatChop(In, 1.0e-10);
2413:     MatView(In, NULL);
2414:   }
2415:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2416:   return(0);
2417: }

2419: PetscErrorCode DMPlexComputeMassMatrixNested(DM dmc, DM dmf, Mat mass, void *user)
2420: {
2421:   SETERRQ(PetscObjectComm((PetscObject) dmc), PETSC_ERR_SUP, "Laziness");
2422: }

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

2427:   Input Parameters:
2428: + dmf  - The fine mesh
2429: . dmc  - The coarse mesh
2430: - user - The user context

2432:   Output Parameter:
2433: . In  - The interpolation matrix

2435:   Level: developer

2437: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2438: @*/
2439: PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user)
2440: {
2441:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
2442:   const char    *name = "Interpolator";
2443:   PetscDS        prob;
2444:   PetscSection   fsection, csection, globalFSection, globalCSection;
2445:   PetscHSetIJ    ht;
2446:   PetscLayout    rLayout;
2447:   PetscInt      *dnz, *onz;
2448:   PetscInt       locRows, rStart, rEnd;
2449:   PetscReal     *x, *v0, *J, *invJ, detJ;
2450:   PetscReal     *v0c, *Jc, *invJc, detJc;
2451:   PetscScalar   *elemMat;
2452:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

2456:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2457:   DMGetCoordinateDim(dmc, &dim);
2458:   DMGetDS(dmc, &prob);
2459:   PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);
2460:   PetscDSGetNumFields(prob, &Nf);
2461:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2462:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2463:   DMGetLocalSection(dmf, &fsection);
2464:   DMGetGlobalSection(dmf, &globalFSection);
2465:   DMGetLocalSection(dmc, &csection);
2466:   DMGetGlobalSection(dmc, &globalCSection);
2467:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2468:   PetscDSGetTotalDimension(prob, &totDim);
2469:   PetscMalloc1(totDim, &elemMat);

2471:   MatGetLocalSize(In, &locRows, NULL);
2472:   PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
2473:   PetscLayoutSetLocalSize(rLayout, locRows);
2474:   PetscLayoutSetBlockSize(rLayout, 1);
2475:   PetscLayoutSetUp(rLayout);
2476:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2477:   PetscLayoutDestroy(&rLayout);
2478:   PetscCalloc2(locRows,&dnz,locRows,&onz);
2479:   PetscHSetIJCreate(&ht);
2480:   for (field = 0; field < Nf; ++field) {
2481:     PetscObject      obj;
2482:     PetscClassId     id;
2483:     PetscDualSpace   Q = NULL;
2484:     PetscQuadrature  f;
2485:     const PetscReal *qpoints;
2486:     PetscInt         Nc, Np, fpdim, i, d;

2488:     PetscDSGetDiscretization(prob, field, &obj);
2489:     PetscObjectGetClassId(obj, &id);
2490:     if (id == PETSCFE_CLASSID) {
2491:       PetscFE fe = (PetscFE) obj;

2493:       PetscFEGetDualSpace(fe, &Q);
2494:       PetscFEGetNumComponents(fe, &Nc);
2495:     } else if (id == PETSCFV_CLASSID) {
2496:       PetscFV fv = (PetscFV) obj;

2498:       PetscFVGetDualSpace(fv, &Q);
2499:       Nc   = 1;
2500:     }
2501:     PetscDualSpaceGetDimension(Q, &fpdim);
2502:     /* For each fine grid cell */
2503:     for (cell = cStart; cell < cEnd; ++cell) {
2504:       PetscInt *findices,   *cindices;
2505:       PetscInt  numFIndices, numCIndices;

2507:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2508:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2509:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
2510:       for (i = 0; i < fpdim; ++i) {
2511:         Vec             pointVec;
2512:         PetscScalar    *pV;
2513:         PetscSF         coarseCellSF = NULL;
2514:         const PetscSFNode *coarseCells;
2515:         PetscInt        numCoarseCells, q, c;

2517:         /* Get points from the dual basis functional quadrature */
2518:         PetscDualSpaceGetFunctional(Q, i, &f);
2519:         PetscQuadratureGetData(f, NULL, NULL, &Np, &qpoints, NULL);
2520:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2521:         VecSetBlockSize(pointVec, dim);
2522:         VecGetArray(pointVec, &pV);
2523:         for (q = 0; q < Np; ++q) {
2524:           const PetscReal xi0[3] = {-1., -1., -1.};

2526:           /* Transform point to real space */
2527:           CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2528:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2529:         }
2530:         VecRestoreArray(pointVec, &pV);
2531:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2532:         /* OPT: Pack all quad points from fine cell */
2533:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2534:         PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2535:         /* Update preallocation info */
2536:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2537:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2538:         {
2539:           PetscHashIJKey key;
2540:           PetscBool      missing;

2542:           key.i = findices[i];
2543:           if (key.i >= 0) {
2544:             /* Get indices for coarse elements */
2545:             for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2546:               DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2547:               for (c = 0; c < numCIndices; ++c) {
2548:                 key.j = cindices[c];
2549:                 if (key.j < 0) continue;
2550:                 PetscHSetIJQueryAdd(ht, key, &missing);
2551:                 if (missing) {
2552:                   if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2553:                   else                                     ++onz[key.i-rStart];
2554:                 }
2555:               }
2556:               DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2557:             }
2558:           }
2559:         }
2560:         PetscSFDestroy(&coarseCellSF);
2561:         VecDestroy(&pointVec);
2562:       }
2563:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2564:     }
2565:   }
2566:   PetscHSetIJDestroy(&ht);
2567:   MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
2568:   MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2569:   PetscFree2(dnz,onz);
2570:   for (field = 0; field < Nf; ++field) {
2571:     PetscObject      obj;
2572:     PetscClassId     id;
2573:     PetscDualSpace   Q = NULL;
2574:     PetscQuadrature  f;
2575:     const PetscReal *qpoints, *qweights;
2576:     PetscInt         Nc, qNc, Np, fpdim, i, d;

2578:     PetscDSGetDiscretization(prob, field, &obj);
2579:     PetscObjectGetClassId(obj, &id);
2580:     if (id == PETSCFE_CLASSID) {
2581:       PetscFE fe = (PetscFE) obj;

2583:       PetscFEGetDualSpace(fe, &Q);
2584:       PetscFEGetNumComponents(fe, &Nc);
2585:     } else if (id == PETSCFV_CLASSID) {
2586:       PetscFV fv = (PetscFV) obj;

2588:       PetscFVGetDualSpace(fv, &Q);
2589:       Nc   = 1;
2590:     } else SETERRQ1(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONG,"Unknown discretization type for field %d",field);
2591:     PetscDualSpaceGetDimension(Q, &fpdim);
2592:     /* For each fine grid cell */
2593:     for (cell = cStart; cell < cEnd; ++cell) {
2594:       PetscInt *findices,   *cindices;
2595:       PetscInt  numFIndices, numCIndices;

2597:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2598:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2599:       if (numFIndices != fpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim);
2600:       for (i = 0; i < fpdim; ++i) {
2601:         Vec             pointVec;
2602:         PetscScalar    *pV;
2603:         PetscSF         coarseCellSF = NULL;
2604:         const PetscSFNode *coarseCells;
2605:         PetscInt        numCoarseCells, cpdim, q, c, j;

2607:         /* Get points from the dual basis functional quadrature */
2608:         PetscDualSpaceGetFunctional(Q, i, &f);
2609:         PetscQuadratureGetData(f, NULL, &qNc, &Np, &qpoints, &qweights);
2610:         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);
2611:         VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);
2612:         VecSetBlockSize(pointVec, dim);
2613:         VecGetArray(pointVec, &pV);
2614:         for (q = 0; q < Np; ++q) {
2615:           const PetscReal xi0[3] = {-1., -1., -1.};

2617:           /* Transform point to real space */
2618:           CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2619:           for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2620:         }
2621:         VecRestoreArray(pointVec, &pV);
2622:         /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2623:         /* OPT: Read this out from preallocation information */
2624:         DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2625:         /* Update preallocation info */
2626:         PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2627:         if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2628:         VecGetArray(pointVec, &pV);
2629:         for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2630:           PetscReal pVReal[3];
2631:           const PetscReal xi0[3] = {-1., -1., -1.};

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

2639:           if (id == PETSCFE_CLASSID) {
2640:             PetscFE    fe = (PetscFE) obj;
2641:             PetscReal *B;

2643:             /* Evaluate coarse basis on contained point */
2644:             PetscFEGetDimension(fe, &cpdim);
2645:             PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2646:             PetscArrayzero(elemMat, cpdim);
2647:             /* Get elemMat entries by multiplying by weight */
2648:             for (j = 0; j < cpdim; ++j) {
2649:               for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*qweights[ccell*qNc + c];
2650:             }
2651:             PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2652:           } else {
2653:             cpdim = 1;
2654:             for (j = 0; j < cpdim; ++j) {
2655:               for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*qweights[ccell*qNc + c];
2656:             }
2657:           }
2658:           /* Update interpolator */
2659:           if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2660:           if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2661:           MatSetValues(In, 1, &findices[i], numCIndices, cindices, elemMat, INSERT_VALUES);
2662:           DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2663:         }
2664:         VecRestoreArray(pointVec, &pV);
2665:         PetscSFDestroy(&coarseCellSF);
2666:         VecDestroy(&pointVec);
2667:       }
2668:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2669:     }
2670:   }
2671:   PetscFree3(v0,J,invJ);
2672:   PetscFree3(v0c,Jc,invJc);
2673:   PetscFree(elemMat);
2674:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
2675:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
2676:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
2677:   return(0);
2678: }

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

2683:   Input Parameters:
2684: + dmf  - The fine mesh
2685: . dmc  - The coarse mesh
2686: - user - The user context

2688:   Output Parameter:
2689: . mass  - The mass matrix

2691:   Level: developer

2693: .seealso: DMPlexComputeMassMatrixNested(), DMPlexComputeInterpolatorNested(), DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM()
2694: @*/
2695: PetscErrorCode DMPlexComputeMassMatrixGeneral(DM dmc, DM dmf, Mat mass, void *user)
2696: {
2697:   DM_Plex       *mesh = (DM_Plex *) dmf->data;
2698:   const char    *name = "Mass Matrix";
2699:   PetscDS        prob;
2700:   PetscSection   fsection, csection, globalFSection, globalCSection;
2701:   PetscHSetIJ    ht;
2702:   PetscLayout    rLayout;
2703:   PetscInt      *dnz, *onz;
2704:   PetscInt       locRows, rStart, rEnd;
2705:   PetscReal     *x, *v0, *J, *invJ, detJ;
2706:   PetscReal     *v0c, *Jc, *invJc, detJc;
2707:   PetscScalar   *elemMat;
2708:   PetscInt       dim, Nf, field, totDim, cStart, cEnd, cell, ccell;

2712:   DMGetCoordinateDim(dmc, &dim);
2713:   DMGetDS(dmc, &prob);
2714:   PetscDSGetWorkspace(prob, &x, NULL, NULL, NULL, NULL);
2715:   PetscDSGetNumFields(prob, &Nf);
2716:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
2717:   PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);
2718:   DMGetLocalSection(dmf, &fsection);
2719:   DMGetGlobalSection(dmf, &globalFSection);
2720:   DMGetLocalSection(dmc, &csection);
2721:   DMGetGlobalSection(dmc, &globalCSection);
2722:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
2723:   PetscDSGetTotalDimension(prob, &totDim);
2724:   PetscMalloc1(totDim, &elemMat);

2726:   MatGetLocalSize(mass, &locRows, NULL);
2727:   PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);
2728:   PetscLayoutSetLocalSize(rLayout, locRows);
2729:   PetscLayoutSetBlockSize(rLayout, 1);
2730:   PetscLayoutSetUp(rLayout);
2731:   PetscLayoutGetRange(rLayout, &rStart, &rEnd);
2732:   PetscLayoutDestroy(&rLayout);
2733:   PetscCalloc2(locRows,&dnz,locRows,&onz);
2734:   PetscHSetIJCreate(&ht);
2735:   for (field = 0; field < Nf; ++field) {
2736:     PetscObject      obj;
2737:     PetscClassId     id;
2738:     PetscQuadrature  quad;
2739:     const PetscReal *qpoints;
2740:     PetscInt         Nq, Nc, i, d;

2742:     PetscDSGetDiscretization(prob, field, &obj);
2743:     PetscObjectGetClassId(obj, &id);
2744:     if (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, &quad);}
2745:     else                       {PetscFVGetQuadrature((PetscFV) obj, &quad);}
2746:     PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);
2747:     /* For each fine grid cell */
2748:     for (cell = cStart; cell < cEnd; ++cell) {
2749:       Vec                pointVec;
2750:       PetscScalar       *pV;
2751:       PetscSF            coarseCellSF = NULL;
2752:       const PetscSFNode *coarseCells;
2753:       PetscInt           numCoarseCells, q, c;
2754:       PetscInt          *findices,   *cindices;
2755:       PetscInt           numFIndices, numCIndices;

2757:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2758:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
2759:       /* Get points from the quadrature */
2760:       VecCreateSeq(PETSC_COMM_SELF, Nq*dim, &pointVec);
2761:       VecSetBlockSize(pointVec, dim);
2762:       VecGetArray(pointVec, &pV);
2763:       for (q = 0; q < Nq; ++q) {
2764:         const PetscReal xi0[3] = {-1., -1., -1.};

2766:         /* Transform point to real space */
2767:         CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2768:         for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2769:       }
2770:       VecRestoreArray(pointVec, &pV);
2771:       /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2772:       DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2773:       PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");
2774:       /* Update preallocation info */
2775:       PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2776:       if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2777:       {
2778:         PetscHashIJKey key;
2779:         PetscBool      missing;

2781:         for (i = 0; i < numFIndices; ++i) {
2782:           key.i = findices[i];
2783:           if (key.i >= 0) {
2784:             /* Get indices for coarse elements */
2785:             for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2786:               DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2787:               for (c = 0; c < numCIndices; ++c) {
2788:                 key.j = cindices[c];
2789:                 if (key.j < 0) continue;
2790:                 PetscHSetIJQueryAdd(ht, key, &missing);
2791:                 if (missing) {
2792:                   if ((key.j >= rStart) && (key.j < rEnd)) ++dnz[key.i-rStart];
2793:                   else                                     ++onz[key.i-rStart];
2794:                 }
2795:               }
2796:               DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2797:             }
2798:           }
2799:         }
2800:       }
2801:       PetscSFDestroy(&coarseCellSF);
2802:       VecDestroy(&pointVec);
2803:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2804:     }
2805:   }
2806:   PetscHSetIJDestroy(&ht);
2807:   MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
2808:   MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
2809:   PetscFree2(dnz,onz);
2810:   for (field = 0; field < Nf; ++field) {
2811:     PetscObject      obj;
2812:     PetscClassId     id;
2813:     PetscQuadrature  quad;
2814:     PetscReal       *Bfine;
2815:     const PetscReal *qpoints, *qweights;
2816:     PetscInt         Nq, Nc, i, d;

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

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

2842:         /* Transform point to real space */
2843:         CoordinatesRefToReal(dim, dim, xi0, v0, J, &qpoints[q*dim], x);
2844:         for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d];
2845:       }
2846:       VecRestoreArray(pointVec, &pV);
2847:       /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */
2848:       DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);
2849:       /* Update matrix */
2850:       PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);
2851:       if (numCoarseCells != Nq) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located");
2852:       VecGetArray(pointVec, &pV);
2853:       for (ccell = 0; ccell < numCoarseCells; ++ccell) {
2854:         PetscReal pVReal[3];
2855:         const PetscReal xi0[3] = {-1., -1., -1.};


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

2864:         if (id == PETSCFE_CLASSID) {
2865:           PetscFE    fe = (PetscFE) obj;
2866:           PetscReal *B;

2868:           /* Evaluate coarse basis on contained point */
2869:           PetscFEGetDimension(fe, &cpdim);
2870:           PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);
2871:           /* Get elemMat entries by multiplying by weight */
2872:           for (i = 0; i < numFIndices; ++i) {
2873:             PetscArrayzero(elemMat, cpdim);
2874:             for (j = 0; j < cpdim; ++j) {
2875:               for (c = 0; c < Nc; ++c) elemMat[j] += B[j*Nc + c]*Bfine[(ccell*numFIndices + i)*Nc + c]*qweights[ccell*Nc + c]*detJ;
2876:             }
2877:             /* Update interpolator */
2878:             if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2879:             if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2880:             MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2881:           }
2882:           PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);
2883:         } else {
2884:           cpdim = 1;
2885:           for (i = 0; i < numFIndices; ++i) {
2886:             PetscArrayzero(elemMat, cpdim);
2887:             for (j = 0; j < cpdim; ++j) {
2888:               for (c = 0; c < Nc; ++c) elemMat[j] += 1.0*1.0*qweights[ccell*Nc + c]*detJ;
2889:             }
2890:             /* Update interpolator */
2891:             if (mesh->printFEM > 1) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
2892:             PetscPrintf(PETSC_COMM_SELF, "Nq: %d %d Nf: %d %d Nc: %d %d\n", ccell, Nq, i, numFIndices, j, numCIndices);
2893:             if (numCIndices != cpdim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D", numCIndices, cpdim);
2894:             MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
2895:           }
2896:         }
2897:         DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);
2898:       }
2899:       VecRestoreArray(pointVec, &pV);
2900:       PetscSFDestroy(&coarseCellSF);
2901:       VecDestroy(&pointVec);
2902:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
2903:     }
2904:   }
2905:   PetscFree3(v0,J,invJ);
2906:   PetscFree3(v0c,Jc,invJc);
2907:   PetscFree(elemMat);
2908:   MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
2909:   MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
2910:   return(0);
2911: }

2913: /*@
2914:   DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns

2916:   Input Parameters:
2917: + dmc  - The coarse mesh
2918: - dmf  - The fine mesh
2919: - user - The user context

2921:   Output Parameter:
2922: . sc   - The mapping

2924:   Level: developer

2926: .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM()
2927: @*/
2928: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
2929: {
2930:   PetscDS        prob;
2931:   PetscFE       *feRef;
2932:   PetscFV       *fvRef;
2933:   Vec            fv, cv;
2934:   IS             fis, cis;
2935:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
2936:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
2937:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m;
2938:   PetscBool     *needAvg;

2942:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
2943:   DMGetDimension(dmf, &dim);
2944:   DMGetLocalSection(dmf, &fsection);
2945:   DMGetGlobalSection(dmf, &fglobalSection);
2946:   DMGetLocalSection(dmc, &csection);
2947:   DMGetGlobalSection(dmc, &cglobalSection);
2948:   PetscSectionGetNumFields(fsection, &Nf);
2949:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
2950:   DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);
2951:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2952:   DMGetDS(dmc, &prob);
2953:   PetscCalloc3(Nf,&feRef,Nf,&fvRef,Nf,&needAvg);
2954:   for (f = 0; f < Nf; ++f) {
2955:     PetscObject  obj;
2956:     PetscClassId id;
2957:     PetscInt     fNb = 0, Nc = 0;

2959:     PetscDSGetDiscretization(prob, f, &obj);
2960:     PetscObjectGetClassId(obj, &id);
2961:     if (id == PETSCFE_CLASSID) {
2962:       PetscFE    fe = (PetscFE) obj;
2963:       PetscSpace sp;
2964:       PetscInt   maxDegree;

2966:       PetscFERefine(fe, &feRef[f]);
2967:       PetscFEGetDimension(feRef[f], &fNb);
2968:       PetscFEGetNumComponents(fe, &Nc);
2969:       PetscFEGetBasisSpace(fe, &sp);
2970:       PetscSpaceGetDegree(sp, NULL, &maxDegree);
2971:       if (!maxDegree) needAvg[f] = PETSC_TRUE;
2972:     } else if (id == PETSCFV_CLASSID) {
2973:       PetscFV        fv = (PetscFV) obj;
2974:       PetscDualSpace Q;

2976:       PetscFVRefine(fv, &fvRef[f]);
2977:       PetscFVGetDualSpace(fvRef[f], &Q);
2978:       PetscDualSpaceGetDimension(Q, &fNb);
2979:       PetscFVGetNumComponents(fv, &Nc);
2980:       needAvg[f] = PETSC_TRUE;
2981:     }
2982:     fTotDim += fNb;
2983:   }
2984:   PetscDSGetTotalDimension(prob, &cTotDim);
2985:   PetscMalloc1(cTotDim,&cmap);
2986:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
2987:     PetscFE        feC;
2988:     PetscFV        fvC;
2989:     PetscDualSpace QF, QC;
2990:     PetscInt       order = -1, NcF, NcC, fpdim, cpdim;

2992:     if (feRef[field]) {
2993:       PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
2994:       PetscFEGetNumComponents(feC, &NcC);
2995:       PetscFEGetNumComponents(feRef[field], &NcF);
2996:       PetscFEGetDualSpace(feRef[field], &QF);
2997:       PetscDualSpaceGetOrder(QF, &order);
2998:       PetscDualSpaceGetDimension(QF, &fpdim);
2999:       PetscFEGetDualSpace(feC, &QC);
3000:       PetscDualSpaceGetDimension(QC, &cpdim);
3001:     } else {
3002:       PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);
3003:       PetscFVGetNumComponents(fvC, &NcC);
3004:       PetscFVGetNumComponents(fvRef[field], &NcF);
3005:       PetscFVGetDualSpace(fvRef[field], &QF);
3006:       PetscDualSpaceGetDimension(QF, &fpdim);
3007:       PetscFVGetDualSpace(fvC, &QC);
3008:       PetscDualSpaceGetDimension(QC, &cpdim);
3009:     }
3010:     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);
3011:     for (c = 0; c < cpdim; ++c) {
3012:       PetscQuadrature  cfunc;
3013:       const PetscReal *cqpoints, *cqweights;
3014:       PetscInt         NqcC, NpC;
3015:       PetscBool        found = PETSC_FALSE;

3017:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
3018:       PetscQuadratureGetData(cfunc, NULL, &NqcC, &NpC, &cqpoints, &cqweights);
3019:       if (NqcC != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcC, NcC);
3020:       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
3021:       for (f = 0; f < fpdim; ++f) {
3022:         PetscQuadrature  ffunc;
3023:         const PetscReal *fqpoints, *fqweights;
3024:         PetscReal        sum = 0.0;
3025:         PetscInt         NqcF, NpF;

3027:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
3028:         PetscQuadratureGetData(ffunc, NULL, &NqcF, &NpF, &fqpoints, &fqweights);
3029:         if (NqcF != NcF) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of quadrature components %D must match number of field components", NqcF, NcF);
3030:         if (NpC != NpF) continue;
3031:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
3032:         if (sum > 1.0e-9) continue;
3033:         for (d = 0; d < NcC; ++d) sum += PetscAbsReal(cqweights[d]*fqweights[d]);
3034:         if (sum < 1.0e-9) continue;
3035:         cmap[offsetC+c] = offsetF+f;
3036:         found = PETSC_TRUE;
3037:         break;
3038:       }
3039:       if (!found) {
3040:         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
3041:         if (fvRef[field] || (feRef[field] && order == 0)) {
3042:           cmap[offsetC+c] = offsetF+0;
3043:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
3044:       }
3045:     }
3046:     offsetC += cpdim;
3047:     offsetF += fpdim;
3048:   }
3049:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);PetscFVDestroy(&fvRef[f]);}
3050:   PetscFree3(feRef,fvRef,needAvg);

3052:   DMGetGlobalVector(dmf, &fv);
3053:   DMGetGlobalVector(dmc, &cv);
3054:   VecGetOwnershipRange(cv, &startC, &endC);
3055:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
3056:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
3057:   PetscMalloc1(m,&cindices);
3058:   PetscMalloc1(m,&findices);
3059:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
3060:   for (c = cStart; c < cEnd; ++c) {
3061:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
3062:     for (d = 0; d < cTotDim; ++d) {
3063:       if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue;
3064:       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]]);
3065:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
3066:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
3067:     }
3068:   }
3069:   PetscFree(cmap);
3070:   PetscFree2(cellCIndices,cellFIndices);

3072:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
3073:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
3074:   VecScatterCreate(cv, cis, fv, fis, sc);
3075:   ISDestroy(&cis);
3076:   ISDestroy(&fis);
3077:   DMRestoreGlobalVector(dmf, &fv);
3078:   DMRestoreGlobalVector(dmc, &cv);
3079:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
3080:   return(0);
3081: }

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

3086:   Input Parameters:
3087: + dm     - The DM
3088: . cellIS - The cells to include
3089: . locX   - A local vector with the solution fields
3090: . locX_t - A local vector with solution field time derivatives, or NULL
3091: - locA   - A local vector with auxiliary fields, or NULL

3093:   Output Parameters:
3094: + u   - The field coefficients
3095: . u_t - The fields derivative coefficients
3096: - a   - The auxiliary field coefficients

3098:   Level: developer

3100: .seealso: DMPlexGetFaceFields()
3101: @*/
3102: PetscErrorCode DMPlexGetCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
3103: {
3104:   DM              plex, plexA = NULL;
3105:   PetscSection    section, sectionAux;
3106:   PetscDS         prob;
3107:   const PetscInt *cells;
3108:   PetscInt        cStart, cEnd, numCells, totDim, totDimAux, c;
3109:   PetscErrorCode  ierr;

3119:   DMPlexConvertPlex(dm, &plex, PETSC_FALSE);
3120:   ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3121:   DMGetLocalSection(dm, &section);
3122:   DMGetCellDS(dm, cStart, &prob);
3123:   PetscDSGetTotalDimension(prob, &totDim);
3124:   if (locA) {
3125:     DM      dmAux;
3126:     PetscDS probAux;

3128:     VecGetDM(locA, &dmAux);
3129:     DMPlexConvertPlex(dmAux, &plexA, PETSC_FALSE);
3130:     DMGetLocalSection(dmAux, &sectionAux);
3131:     DMGetDS(dmAux, &probAux);
3132:     PetscDSGetTotalDimension(probAux, &totDimAux);
3133:   }
3134:   numCells = cEnd - cStart;
3135:   DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u);
3136:   if (locX_t) {DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u_t);} else {*u_t = NULL;}
3137:   if (locA)   {DMGetWorkArray(dm, numCells*totDimAux, MPIU_SCALAR, a);} else {*a = NULL;}
3138:   for (c = cStart; c < cEnd; ++c) {
3139:     const PetscInt cell = cells ? cells[c] : c;
3140:     const PetscInt cind = c - cStart;
3141:     PetscScalar   *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
3142:     PetscInt       i;

3144:     DMPlexVecGetClosure(plex, section, locX, cell, NULL, &x);
3145:     for (i = 0; i < totDim; ++i) ul[cind*totDim+i] = x[i];
3146:     DMPlexVecRestoreClosure(plex, section, locX, cell, NULL, &x);
3147:     if (locX_t) {
3148:       DMPlexVecGetClosure(plex, section, locX_t, cell, NULL, &x_t);
3149:       for (i = 0; i < totDim; ++i) ul_t[cind*totDim+i] = x_t[i];
3150:       DMPlexVecRestoreClosure(plex, section, locX_t, cell, NULL, &x_t);
3151:     }
3152:     if (locA) {
3153:       PetscInt subcell;
3154:       DMPlexGetAuxiliaryPoint(plex, plexA, cell, &subcell);
3155:       DMPlexVecGetClosure(plexA, sectionAux, locA, subcell, NULL, &x);
3156:       for (i = 0; i < totDimAux; ++i) al[cind*totDimAux+i] = x[i];
3157:       DMPlexVecRestoreClosure(plexA, sectionAux, locA, subcell, NULL, &x);
3158:     }
3159:   }
3160:   DMDestroy(&plex);
3161:   if (locA) {DMDestroy(&plexA);}
3162:   ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
3163:   return(0);
3164: }

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

3169:   Input Parameters:
3170: + dm     - The DM
3171: . cellIS - The cells to include
3172: . locX   - A local vector with the solution fields
3173: . locX_t - A local vector with solution field time derivatives, or NULL
3174: - locA   - A local vector with auxiliary fields, or NULL

3176:   Output Parameters:
3177: + u   - The field coefficients
3178: . u_t - The fields derivative coefficients
3179: - a   - The auxiliary field coefficients

3181:   Level: developer

3183: .seealso: DMPlexGetFaceFields()
3184: @*/
3185: PetscErrorCode DMPlexRestoreCellFields(DM dm, IS cellIS, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
3186: {

3190:   DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u);
3191:   if (locX_t) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u_t);}
3192:   if (locA)   {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, a);}
3193:   return(0);
3194: }

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

3199:   Input Parameters:
3200: + dm     - The DM
3201: . fStart - The first face to include
3202: . fEnd   - The first face to exclude
3203: . locX   - A local vector with the solution fields
3204: . locX_t - A local vector with solution field time derivatives, or NULL
3205: . faceGeometry - A local vector with face geometry
3206: . cellGeometry - A local vector with cell geometry
3207: - locaGrad - A local vector with field gradients, or NULL

3209:   Output Parameters:
3210: + Nface - The number of faces with field values
3211: . uL - The field values at the left side of the face
3212: - uR - The field values at the right side of the face

3214:   Level: developer

3216: .seealso: DMPlexGetCellFields()
3217: @*/
3218: 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)
3219: {
3220:   DM                 dmFace, dmCell, dmGrad = NULL;
3221:   PetscSection       section;
3222:   PetscDS            prob;
3223:   DMLabel            ghostLabel;
3224:   const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
3225:   PetscBool         *isFE;
3226:   PetscInt           dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
3227:   PetscErrorCode     ierr;

3238:   DMGetDimension(dm, &dim);
3239:   DMGetDS(dm, &prob);
3240:   DMGetLocalSection(dm, &section);
3241:   PetscDSGetNumFields(prob, &Nf);
3242:   PetscDSGetTotalComponents(prob, &Nc);
3243:   PetscMalloc1(Nf, &isFE);
3244:   for (f = 0; f < Nf; ++f) {
3245:     PetscObject  obj;
3246:     PetscClassId id;

3248:     PetscDSGetDiscretization(prob, f, &obj);
3249:     PetscObjectGetClassId(obj, &id);
3250:     if (id == PETSCFE_CLASSID)      {isFE[f] = PETSC_TRUE;}
3251:     else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
3252:     else                            {isFE[f] = PETSC_FALSE;}
3253:   }
3254:   DMGetLabel(dm, "ghost", &ghostLabel);
3255:   VecGetArrayRead(locX, &x);
3256:   VecGetDM(faceGeometry, &dmFace);
3257:   VecGetArrayRead(faceGeometry, &facegeom);
3258:   VecGetDM(cellGeometry, &dmCell);
3259:   VecGetArrayRead(cellGeometry, &cellgeom);
3260:   if (locGrad) {
3261:     VecGetDM(locGrad, &dmGrad);
3262:     VecGetArrayRead(locGrad, &lgrad);
3263:   }
3264:   DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uL);
3265:   DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uR);
3266:   /* Right now just eat the extra work for FE (could make a cell loop) */
3267:   for (face = fStart, iface = 0; face < fEnd; ++face) {
3268:     const PetscInt        *cells;
3269:     PetscFVFaceGeom       *fg;
3270:     PetscFVCellGeom       *cgL, *cgR;
3271:     PetscScalar           *xL, *xR, *gL, *gR;
3272:     PetscScalar           *uLl = *uL, *uRl = *uR;
3273:     PetscInt               ghost, nsupp, nchild;

3275:     DMLabelGetValue(ghostLabel, face, &ghost);
3276:     DMPlexGetSupportSize(dm, face, &nsupp);
3277:     DMPlexGetTreeChildren(dm, face, &nchild, NULL);
3278:     if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
3279:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
3280:     DMPlexGetSupport(dm, face, &cells);
3281:     DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
3282:     DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
3283:     for (f = 0; f < Nf; ++f) {
3284:       PetscInt off;

3286:       PetscDSGetComponentOffset(prob, f, &off);
3287:       if (isFE[f]) {
3288:         const PetscInt *cone;
3289:         PetscInt        comp, coneSizeL, coneSizeR, faceLocL, faceLocR, ldof, rdof, d;

3291:         xL = xR = NULL;
3292:         PetscSectionGetFieldComponents(section, f, &comp);
3293:         DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
3294:         DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
3295:         DMPlexGetCone(dm, cells[0], &cone);
3296:         DMPlexGetConeSize(dm, cells[0], &coneSizeL);
3297:         for (faceLocL = 0; faceLocL < coneSizeL; ++faceLocL) if (cone[faceLocL] == face) break;
3298:         DMPlexGetCone(dm, cells[1], &cone);
3299:         DMPlexGetConeSize(dm, cells[1], &coneSizeR);
3300:         for (faceLocR = 0; faceLocR < coneSizeR; ++faceLocR) if (cone[faceLocR] == face) break;
3301:         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]);
3302:         /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
3303:         /* TODO: this is a hack that might not be right for nonconforming */
3304:         if (faceLocL < coneSizeL) {
3305:           PetscFEEvaluateFaceFields_Internal(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
3306:           if (rdof == ldof && faceLocR < coneSizeR) {PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
3307:           else              {for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
3308:         }
3309:         else {
3310:           PetscFEEvaluateFaceFields_Internal(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);
3311:           PetscSectionGetFieldComponents(section, f, &comp);
3312:           for(d = 0; d < comp; ++d) uLl[iface*Nc+off+d] = uRl[iface*Nc+off+d];
3313:         }
3314:         DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
3315:         DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
3316:       } else {
3317:         PetscFV  fv;
3318:         PetscInt numComp, c;

3320:         PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
3321:         PetscFVGetNumComponents(fv, &numComp);
3322:         DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
3323:         DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
3324:         if (dmGrad) {
3325:           PetscReal dxL[3], dxR[3];

3327:           DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
3328:           DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
3329:           DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
3330:           DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
3331:           for (c = 0; c < numComp; ++c) {
3332:             uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
3333:             uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
3334:           }
3335:         } else {
3336:           for (c = 0; c < numComp; ++c) {
3337:             uLl[iface*Nc+off+c] = xL[c];
3338:             uRl[iface*Nc+off+c] = xR[c];
3339:           }
3340:         }
3341:       }
3342:     }
3343:     ++iface;
3344:   }
3345:   *Nface = iface;
3346:   VecRestoreArrayRead(locX, &x);
3347:   VecRestoreArrayRead(faceGeometry, &facegeom);
3348:   VecRestoreArrayRead(cellGeometry, &cellgeom);
3349:   if (locGrad) {
3350:     VecRestoreArrayRead(locGrad, &lgrad);
3351:   }
3352:   PetscFree(isFE);
3353:   return(0);
3354: }

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

3359:   Input Parameters:
3360: + dm     - The DM
3361: . fStart - The first face to include
3362: . fEnd   - The first face to exclude
3363: . locX   - A local vector with the solution fields
3364: . locX_t - A local vector with solution field time derivatives, or NULL
3365: . faceGeometry - A local vector with face geometry
3366: . cellGeometry - A local vector with cell geometry
3367: - locaGrad - A local vector with field gradients, or NULL

3369:   Output Parameters:
3370: + Nface - The number of faces with field values
3371: . uL - The field values at the left side of the face
3372: - uR - The field values at the right side of the face

3374:   Level: developer

3376: .seealso: DMPlexGetFaceFields()
3377: @*/
3378: 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)
3379: {

3383:   DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uL);
3384:   DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uR);
3385:   return(0);
3386: }

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

3391:   Input Parameters:
3392: + dm     - The DM
3393: . fStart - The first face to include
3394: . fEnd   - The first face to exclude
3395: . faceGeometry - A local vector with face geometry
3396: - cellGeometry - A local vector with cell geometry

3398:   Output Parameters:
3399: + Nface - The number of faces with field values
3400: . fgeom - The extract the face centroid and normal
3401: - vol   - The cell volume

3403:   Level: developer

3405: .seealso: DMPlexGetCellFields()
3406: @*/
3407: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
3408: {
3409:   DM                 dmFace, dmCell;
3410:   DMLabel            ghostLabel;
3411:   const PetscScalar *facegeom, *cellgeom;
3412:   PetscInt           dim, numFaces = fEnd - fStart, iface, face;
3413:   PetscErrorCode     ierr;

3421:   DMGetDimension(dm, &dim);
3422:   DMGetLabel(dm, "ghost", &ghostLabel);
3423:   VecGetDM(faceGeometry, &dmFace);
3424:   VecGetArrayRead(faceGeometry, &facegeom);
3425:   VecGetDM(cellGeometry, &dmCell);
3426:   VecGetArrayRead(cellGeometry, &cellgeom);
3427:   PetscMalloc1(numFaces, fgeom);
3428:   DMGetWorkArray(dm, numFaces*2, MPIU_SCALAR, vol);
3429:   for (face = fStart, iface = 0; face < fEnd; ++face) {
3430:     const PetscInt        *cells;
3431:     PetscFVFaceGeom       *fg;
3432:     PetscFVCellGeom       *cgL, *cgR;
3433:     PetscFVFaceGeom       *fgeoml = *fgeom;
3434:     PetscReal             *voll   = *vol;
3435:     PetscInt               ghost, d, nchild, nsupp;

3437:     DMLabelGetValue(ghostLabel, face, &ghost);
3438:     DMPlexGetSupportSize(dm, face, &nsupp);
3439:     DMPlexGetTreeChildren(dm, face, &nchild, NULL);
3440:     if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
3441:     DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
3442:     DMPlexGetSupport(dm, face, &cells);
3443:     DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
3444:     DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
3445:     for (d = 0; d < dim; ++d) {
3446:       fgeoml[iface].centroid[d] = fg->centroid[d];
3447:       fgeoml[iface].normal[d]   = fg->normal[d];
3448:     }
3449:     voll[iface*2+0] = cgL->volume;
3450:     voll[iface*2+1] = cgR->volume;
3451:     ++iface;
3452:   }
3453:   *Nface = iface;
3454:   VecRestoreArrayRead(faceGeometry, &facegeom);
3455:   VecRestoreArrayRead(cellGeometry, &cellgeom);
3456:   return(0);
3457: }

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

3462:   Input Parameters:
3463: + dm     - The DM
3464: . fStart - The first face to include
3465: . fEnd   - The first face to exclude
3466: . faceGeometry - A local vector with face geometry
3467: - cellGeometry - A local vector with cell geometry

3469:   Output Parameters:
3470: + Nface - The number of faces with field values
3471: . fgeom - The extract the face centroid and normal
3472: - vol   - The cell volume

3474:   Level: developer

3476: .seealso: DMPlexGetFaceFields()
3477: @*/
3478: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
3479: {

3483:   PetscFree(*fgeom);
3484:   DMRestoreWorkArray(dm, 0, MPIU_REAL, vol);
3485:   return(0);
3486: }

3488: PetscErrorCode DMSNESGetFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3489: {
3490:   char            composeStr[33] = {0};
3491:   PetscObjectId   id;
3492:   PetscContainer  container;
3493:   PetscErrorCode  ierr;

3496:   PetscObjectGetId((PetscObject)quad,&id);
3497:   PetscSNPrintf(composeStr, 32, "DMSNESGetFEGeom_%x\n", id);
3498:   PetscObjectQuery((PetscObject) pointIS, composeStr, (PetscObject *) &container);
3499:   if (container) {
3500:     PetscContainerGetPointer(container, (void **) geom);
3501:   } else {
3502:     DMFieldCreateFEGeom(coordField, pointIS, quad, faceData, geom);
3503:     PetscContainerCreate(PETSC_COMM_SELF,&container);
3504:     PetscContainerSetPointer(container, (void *) *geom);
3505:     PetscContainerSetUserDestroy(container, PetscContainerUserDestroy_PetscFEGeom);
3506:     PetscObjectCompose((PetscObject) pointIS, composeStr, (PetscObject) container);
3507:     PetscContainerDestroy(&container);
3508:   }
3509:   return(0);
3510: }

3512: PetscErrorCode DMSNESRestoreFEGeom(DMField coordField, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
3513: {
3515:   *geom = NULL;
3516:   return(0);
3517: }

3519: PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM dm, PetscSection section, IS cellIS, PetscReal t, Vec locX, Vec locX_t, Vec locF, void *user)
3520: {
3521:   DM_Plex         *mesh       = (DM_Plex *) dm->data;
3522:   const char      *name       = "Residual";
3523:   DM               dmAux      = NULL;
3524:   DMLabel          ghostLabel = NULL;
3525:   PetscDS          prob       = NULL;
3526:   PetscDS          probAux    = NULL;
3527:   PetscBool        useFEM     = PETSC_FALSE;
3528:   PetscBool        isImplicit = (locX_t || t == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
3529:   DMField          coordField = NULL;
3530:   Vec              locA;
3531:   PetscScalar     *u = NULL, *u_t, *a, *uL = NULL, *uR = NULL;
3532:   IS               chunkIS;
3533:   const PetscInt  *cells;
3534:   PetscInt         cStart, cEnd, numCells;
3535:   PetscInt         Nf, f, totDim, totDimAux, numChunks, cellChunkSize, chunk, fStart, fEnd;
3536:   PetscInt         maxDegree = PETSC_MAX_INT;
3537:   PetscQuadrature  affineQuad = NULL, *quads = NULL;
3538:   PetscFEGeom     *affineGeom = NULL, **geoms = NULL;
3539:   PetscErrorCode   ierr;

3542:   PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
3543:   /* FEM+FVM */
3544:   /* 1: Get sizes from dm and dmAux */
3545:   DMGetLabel(dm, "ghost", &ghostLabel);
3546:   DMGetDS(dm, &prob);
3547:   PetscDSGetNumFields(prob, &Nf);
3548:   PetscDSGetTotalDimension(prob, &totDim);
3549:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
3550:   if (locA) {
3551:     VecGetDM(locA, &dmAux);
3552:     DMGetDS(dmAux, &probAux);
3553:     PetscDSGetTotalDimension(probAux, &totDimAux);
3554:   }
3555:   /* 2: Get geometric data */
3556:   for (f = 0; f < Nf; ++f) {
3557:     PetscObject  obj;
3558:     PetscClassId id;
3559:     PetscBool    fimp;

3561:     PetscDSGetImplicit(prob, f, &fimp);
3562:     if (isImplicit != fimp) continue;
3563:     PetscDSGetDiscretization(prob, f, &obj);
3564:     PetscObjectGetClassId(obj, &id);
3565:     if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
3566:     if (id == PETSCFV_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Use of FVM with PCPATCH not yet implemented");
3567:   }
3568:   if (useFEM) {
3569:     DMGetCoordinateField(dm, &coordField);
3570:     DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
3571:     if (maxDegree <= 1) {
3572:       DMFieldCreateDefaultQuadrature(coordField,cellIS,&affineQuad);
3573:       if (affineQuad) {
3574:         DMSNESGetFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3575:       }
3576:     } else {
3577:       PetscCalloc2(Nf,&quads,Nf,&geoms);
3578:       for (f = 0; f < Nf; ++f) {
3579:         PetscObject  obj;
3580:         PetscClassId id;
3581:         PetscBool    fimp;

3583:         PetscDSGetImplicit(prob, f, &fimp);
3584:         if (isImplicit != fimp) continue;
3585:         PetscDSGetDiscretization(prob, f, &obj);
3586:         PetscObjectGetClassId(obj, &id);
3587:         if (id == PETSCFE_CLASSID) {
3588:           PetscFE fe = (PetscFE) obj;

3590:           PetscFEGetQuadrature(fe, &quads[f]);
3591:           PetscObjectReference((PetscObject)quads[f]);
3592:           DMSNESGetFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3593:         }
3594:       }
3595:     }
3596:   }
3597:   /* Loop over chunks */
3598:   ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3599:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
3600:   if (useFEM) {ISCreate(PETSC_COMM_SELF, &chunkIS);}
3601:   numCells      = cEnd - cStart;
3602:   numChunks     = 1;
3603:   cellChunkSize = numCells/numChunks;
3604:   numChunks     = PetscMin(1,numCells);
3605:   for (chunk = 0; chunk < numChunks; ++chunk) {
3606:     PetscScalar     *elemVec, *fluxL = NULL, *fluxR = NULL;
3607:     PetscReal       *vol = NULL;
3608:     PetscFVFaceGeom *fgeom = NULL;
3609:     PetscInt         cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, c;
3610:     PetscInt         numFaces = 0;

3612:     /* Extract field coefficients */
3613:     if (useFEM) {
3614:       ISGetPointSubrange(chunkIS, cS, cE, cells);
3615:       DMPlexGetCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3616:       DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3617:       PetscArrayzero(elemVec, numCells*totDim);
3618:     }
3619:     /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
3620:     /* Loop over fields */
3621:     for (f = 0; f < Nf; ++f) {
3622:       PetscObject  obj;
3623:       PetscClassId id;
3624:       PetscBool    fimp;
3625:       PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;

3627:       PetscDSGetImplicit(prob, f, &fimp);
3628:       if (isImplicit != fimp) continue;
3629:       PetscDSGetDiscretization(prob, f, &obj);
3630:       PetscObjectGetClassId(obj, &id);
3631:       if (id == PETSCFE_CLASSID) {
3632:         PetscFE         fe = (PetscFE) obj;
3633:         PetscFEGeom    *geom = affineGeom ? affineGeom : geoms[f];
3634:         PetscFEGeom    *chunkGeom = NULL;
3635:         PetscQuadrature quad = affineQuad ? affineQuad : quads[f];
3636:         PetscInt        Nq, Nb;

3638:         PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3639:         PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
3640:         PetscFEGetDimension(fe, &Nb);
3641:         blockSize = Nb;
3642:         batchSize = numBlocks * blockSize;
3643:         PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3644:         numChunks = numCells / (numBatches*batchSize);
3645:         Ne        = numChunks*numBatches*batchSize;
3646:         Nr        = numCells % (numBatches*batchSize);
3647:         offset    = numCells - Nr;
3648:         /* Integrate FE residual to get elemVec (need fields at quadrature points) */
3649:         /*   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) */
3650:         PetscFEGeomGetChunk(geom,0,offset,&chunkGeom);
3651:         PetscFEIntegrateResidual(prob, f, Ne, chunkGeom, u, u_t, probAux, a, t, elemVec);
3652:         PetscFEGeomGetChunk(geom,offset,numCells,&chunkGeom);
3653:         PetscFEIntegrateResidual(prob, f, Nr, chunkGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);
3654:         PetscFEGeomRestoreChunk(geom,offset,numCells,&chunkGeom);
3655:       } else if (id == PETSCFV_CLASSID) {
3656:         PetscFV fv = (PetscFV) obj;

3658:         Ne = numFaces;
3659:         /* Riemann solve over faces (need fields at face centroids) */
3660:         /*   We need to evaluate FE fields at those coordinates */
3661:         PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
3662:       } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
3663:     }
3664:     /* Loop over domain */
3665:     if (useFEM) {
3666:       /* Add elemVec to locX */
3667:       for (c = cS; c < cE; ++c) {
3668:         const PetscInt cell = cells ? cells[c] : c;
3669:         const PetscInt cind = c - cStart;

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

3675:           DMLabelGetValue(ghostLabel,cell,&ghostVal);
3676:           if (ghostVal > 0) continue;
3677:         }
3678:         DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cind*totDim], ADD_ALL_VALUES);
3679:       }
3680:     }
3681:     /* Handle time derivative */
3682:     if (locX_t) {
3683:       PetscScalar *x_t, *fa;

3685:       VecGetArray(locF, &fa);
3686:       VecGetArray(locX_t, &x_t);
3687:       for (f = 0; f < Nf; ++f) {
3688:         PetscFV      fv;
3689:         PetscObject  obj;
3690:         PetscClassId id;
3691:         PetscInt     pdim, d;

3693:         PetscDSGetDiscretization(prob, f, &obj);
3694:         PetscObjectGetClassId(obj, &id);
3695:         if (id != PETSCFV_CLASSID) continue;
3696:         fv   = (PetscFV) obj;
3697:         PetscFVGetNumComponents(fv, &pdim);
3698:         for (c = cS; c < cE; ++c) {
3699:           const PetscInt cell = cells ? cells[c] : c;
3700:           PetscScalar   *u_t, *r;

3702:           if (ghostLabel) {
3703:             PetscInt ghostVal;

3705:             DMLabelGetValue(ghostLabel, cell, &ghostVal);
3706:             if (ghostVal > 0) continue;
3707:           }
3708:           DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
3709:           DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
3710:           for (d = 0; d < pdim; ++d) r[d] += u_t[d];
3711:         }
3712:       }
3713:       VecRestoreArray(locX_t, &x_t);
3714:       VecRestoreArray(locF, &fa);
3715:     }
3716:     if (useFEM) {
3717:       DMPlexRestoreCellFields(dm, chunkIS, locX, locX_t, locA, &u, &u_t, &a);
3718:       DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
3719:     }
3720:   }
3721:   if (useFEM) {ISDestroy(&chunkIS);}
3722:   ISRestorePointRange(cellIS, &cStart, &cEnd, &cells);
3723:   /* TODO Could include boundary residual here (see DMPlexComputeResidual_Internal) */
3724:   if (useFEM) {
3725:     if (maxDegree <= 1) {
3726:       DMSNESRestoreFEGeom(coordField,cellIS,affineQuad,PETSC_FALSE,&affineGeom);
3727:       PetscQuadratureDestroy(&affineQuad);
3728:     } else {
3729:       for (f = 0; f < Nf; ++f) {
3730:         DMSNESRestoreFEGeom(coordField,cellIS,quads[f],PETSC_FALSE,&geoms[f]);
3731:         PetscQuadratureDestroy(&quads[f]);
3732:       }
3733:       PetscFree2(quads,geoms);
3734:     }
3735:   }
3736:   PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
3737:   return(0);
3738: }

3740: /*
3741:   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

3743:   X   - The local solution vector
3744:   X_t - The local solution time derviative vector, or NULL
3745: */
3746: PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM dm, PetscSection section, PetscSection globalSection, IS cellIS,
3747:                                                     PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP, void *ctx)
3748: {
3749:   DM_Plex         *mesh  = (DM_Plex *) dm->data;
3750:   const char      *name = "Jacobian", *nameP = "JacobianPre";
3751:   DM               dmAux = NULL;
3752:   PetscDS          prob,   probAux = NULL;
3753:   PetscSection     sectionAux = NULL;
3754:   Vec              A;
3755:   DMField          coordField;
3756:   PetscFEGeom     *cgeomFEM;
3757:   PetscQuadrature  qGeom = NULL;
3758:   Mat              J = Jac, JP = JacP;
3759:   PetscScalar     *work, *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL, *elemMatP = NULL, *elemMatD = NULL;
3760:   PetscBool        hasJac, hasPrec, hasDyn, assembleJac, isMatIS, isMatISP, *isFE, hasFV = PETSC_FALSE;
3761:   const PetscInt  *cells;
3762:   PetscInt         Nf, fieldI, fieldJ, maxDegree, numCells, cStart, cEnd, numChunks, chunkSize, chunk, totDim, totDimAux = 0, sz, wsz, off = 0, offCell = 0;
3763:   PetscErrorCode   ierr;

3766:   CHKMEMQ;
3767:   ISGetLocalSize(cellIS, &numCells);
3768:   ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
3769:   PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
3770:   DMGetDS(dm, &prob);
3771:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3772:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3773:   if (dmAux) {
3774:     DMGetLocalSection(dmAux, &sectionAux);
3775:     DMGetDS(dmAux, &probAux);
3776:   }
3777:   /* Get flags */
3778:   PetscDSGetNumFields(prob, &Nf);
3779:   DMGetWorkArray(dm, Nf, MPIU_BOOL, &isFE);
3780:   for (fieldI = 0; fieldI < Nf; ++fieldI) {
3781:     PetscObject  disc;
3782:     PetscClassId id;
3783:     PetscDSGetDiscretization(prob, fieldI, &disc);
3784:     PetscObjectGetClassId(disc, &id);
3785:     if (id == PETSCFE_CLASSID)      {isFE[fieldI] = PETSC_TRUE;}
3786:     else if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; isFE[fieldI] = PETSC_FALSE;}
3787:   }
3788:   PetscDSHasJacobian(prob, &hasJac);
3789:   PetscDSHasJacobianPreconditioner(prob, &hasPrec);
3790:   PetscDSHasDynamicJacobian(prob, &hasDyn);
3791:   assembleJac = hasJac && hasPrec && (Jac != JacP) ? PETSC_TRUE : PETSC_FALSE;
3792:   hasDyn      = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
3793:   PetscObjectTypeCompare((PetscObject) Jac,  MATIS, &isMatIS);
3794:   PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
3795:   /* Setup input data and temp arrays (should be DMGetWorkArray) */
3796:   if (isMatISP || isMatISP) {DMPlexGetSubdomainSection(dm, &globalSection);}
3797:   if (isMatIS)  {MatISGetLocalMat(Jac,  &J);}
3798:   if (isMatISP) {MatISGetLocalMat(JacP, &JP);}
3799:   if (hasFV)    {MatSetOption(JP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);} /* No allocated space for FV stuff, so ignore the zero entries */
3800:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
3801:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
3802:   PetscDSGetTotalDimension(prob, &totDim);
3803:   if (probAux) {PetscDSGetTotalDimension(probAux, &totDimAux);}
3804:   CHKMEMQ;
3805:   /* Compute batch sizes */
3806:   if (isFE[0]) {
3807:     PetscFE         fe;
3808:     PetscQuadrature q;
3809:     PetscInt        numQuadPoints, numBatches, batchSize, numBlocks, blockSize, Nb;

3811:     PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3812:     PetscFEGetQuadrature(fe, &q);
3813:     PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
3814:     PetscFEGetDimension(fe, &Nb);
3815:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
3816:     blockSize = Nb*numQuadPoints;
3817:     batchSize = numBlocks  * blockSize;
3818:     chunkSize = numBatches * batchSize;
3819:     numChunks = numCells / chunkSize + numCells % chunkSize;
3820:     PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
3821:   } else {
3822:     chunkSize = numCells;
3823:     numChunks = 1;
3824:   }
3825:   /* Get work space */
3826:   wsz  = (((X?1:0) + (X_t?1:0) + (dmAux?1:0))*totDim + ((hasJac?1:0) + (hasPrec?1:0) + (hasDyn?1:0))*totDim*totDim)*chunkSize;
3827:   DMGetWorkArray(dm, wsz, MPIU_SCALAR, &work);
3828:   PetscArrayzero(work, wsz);
3829:   off      = 0;
3830:   u        = X       ? (sz = chunkSize*totDim,        off += sz, work+off-sz) : NULL;
3831:   u_t      = X_t     ? (sz = chunkSize*totDim,        off += sz, work+off-sz) : NULL;
3832:   a        = dmAux   ? (sz = chunkSize*totDimAux,     off += sz, work+off-sz) : NULL;
3833:   elemMat  = hasJac  ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3834:   elemMatP = hasPrec ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3835:   elemMatD = hasDyn  ? (sz = chunkSize*totDim*totDim, off += sz, work+off-sz) : NULL;
3836:   if (off != wsz) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error is workspace size %D should be %D", off, wsz);
3837:   /* Setup geometry */
3838:   DMGetCoordinateField(dm, &coordField);
3839:   DMFieldGetDegree(coordField, cellIS, NULL, &maxDegree);
3840:   if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField, cellIS, &qGeom);}
3841:   if (!qGeom) {
3842:     PetscFE fe;

3844:     PetscDSGetDiscretization(prob, 0, (PetscObject *) &fe);
3845:     PetscFEGetQuadrature(fe, &qGeom);
3846:     PetscObjectReference((PetscObject) qGeom);
3847:   }
3848:   DMSNESGetFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
3849:   /* Compute volume integrals */
3850:   if (assembleJac) {MatZeroEntries(J);}
3851:   MatZeroEntries(JP);
3852:   for (chunk = 0; chunk < numChunks; ++chunk, offCell += chunkSize) {
3853:     const PetscInt   Ncell = PetscMin(chunkSize, numCells - offCell);
3854:     PetscInt         c;

3856:     /* Extract values */
3857:     for (c = 0; c < Ncell; ++c) {
3858:       const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
3859:       PetscScalar   *x = NULL,  *x_t = NULL;
3860:       PetscInt       i;

3862:       if (X) {
3863:         DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);
3864:         for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
3865:         DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);
3866:       }
3867:       if (X_t) {
3868:         DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);
3869:         for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
3870:         DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);
3871:       }
3872:       if (dmAux) {
3873:         DMPlexVecGetClosure(dmAux, sectionAux, A, cell, NULL, &x);
3874:         for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
3875:         DMPlexVecRestoreClosure(dmAux, sectionAux, A, cell, NULL, &x);
3876:       }
3877:     }
3878:     CHKMEMQ;
3879:     for (fieldI = 0; fieldI < Nf; ++fieldI) {
3880:       PetscFE fe;
3881:       PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3882:       for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
3883:         if (hasJac)  {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN,     fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMat);}
3884:         if (hasPrec) {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatP);}
3885:         if (hasDyn)  {PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ncell, cgeomFEM, u, u_t, probAux, a, t, X_tShift, elemMatD);}
3886:       }
3887:       /* For finite volume, add the identity */
3888:       if (!isFE[fieldI]) {
3889:         PetscFV  fv;
3890:         PetscInt eOffset = 0, Nc, fc, foff;

3892:         PetscDSGetFieldOffset(prob, fieldI, &foff);
3893:         PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);
3894:         PetscFVGetNumComponents(fv, &Nc);
3895:         for (c = 0; c < chunkSize; ++c, eOffset += totDim*totDim) {
3896:           for (fc = 0; fc < Nc; ++fc) {
3897:             const PetscInt i = foff + fc;
3898:             if (hasJac)  {elemMat [eOffset+i*totDim+i] = 1.0;}
3899:             if (hasPrec) {elemMatP[eOffset+i*totDim+i] = 1.0;}
3900:           }
3901:         }
3902:       }
3903:     }
3904:     CHKMEMQ;
3905:     /*   Add contribution from X_t */
3906:     if (hasDyn) {for (c = 0; c < chunkSize*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];}
3907:     /* Insert values into matrix */
3908:     for (c = 0; c < Ncell; ++c) {
3909:       const PetscInt cell = cells ? cells[c+offCell] : c+offCell;
3910:       if (mesh->printFEM > 1) {
3911:         if (hasJac)  {DMPrintCellMatrix(cell, name,  totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
3912:         if (hasPrec) {DMPrintCellMatrix(cell, nameP, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
3913:       }
3914:       if (assembleJac) {DMPlexMatSetClosure(dm, section, globalSection, Jac, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);}
3915:       DMPlexMatSetClosure(dm, section, globalSection, JP, cell, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
3916:     }
3917:     CHKMEMQ;
3918:   }
3919:   /* Cleanup */
3920:   DMSNESRestoreFEGeom(coordField, cellIS, qGeom, PETSC_FALSE, &cgeomFEM);
3921:   PetscQuadratureDestroy(&qGeom);
3922:   if (hasFV) {MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);}
3923:   DMRestoreWorkArray(dm, Nf, MPIU_BOOL, &isFE);
3924:   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);
3925:   /* Compute boundary integrals */
3926:   /* DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, ctx); */
3927:   /* Assemble matrix */
3928:   if (assembleJac) {MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);}
3929:   MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
3930:   PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
3931:   CHKMEMQ;
3932:   return(0);
3933: }