Actual source code: plexfem.c

petsc-dev 2014-04-22
Report Typos and Errors
  1: #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/

  3: #include <petscfe.h>

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

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

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

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

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

 87: /*@C
 88:   DMPlexCreateRigidBody - create rigid body modes from coordinates

 90:   Collective on DM

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

 97:   Output Argument:
 98: . sp - the null space

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

102:   Level: advanced

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

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

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

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

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

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

185: PetscErrorCode DMPlexProjectFunctionLabelLocal(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
186: {
187:   PetscDualSpace *sp;
188:   PetscSection    section;
189:   PetscScalar    *values;
190:   PetscReal      *v0, *J, detJ;
191:   PetscInt        numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, f, d, v, i, comp;
192:   PetscErrorCode  ierr;

195:   DMPlexGetDimension(dm, &dim);
196:   DMGetDefaultSection(dm, &section);
197:   PetscSectionGetNumFields(section, &numFields);
198:   PetscMalloc3(numFields,&sp,dim,&v0,dim*dim,&J);
199:   for (f = 0; f < numFields; ++f) {
200:     PetscFEGetDualSpace(fe[f], &sp[f]);
201:     PetscFEGetNumComponents(fe[f], &numComp);
202:     PetscDualSpaceGetDimension(sp[f], &spDim);
203:     totDim += spDim*numComp;
204:   }
205:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
206:   DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
207:   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
208:   DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
209:   for (i = 0; i < numIds; ++i) {
210:     IS              pointIS;
211:     const PetscInt *points;
212:     PetscInt        n, p;

214:     DMLabelGetStratumIS(label, ids[i], &pointIS);
215:     ISGetLocalSize(pointIS, &n);
216:     ISGetIndices(pointIS, &points);
217:     for (p = 0; p < n; ++p) {
218:       const PetscInt    point = points[p];
219:       PetscCellGeometry geom;

221:       if ((point < cStart) || (point >= cEnd)) continue;
222:       DMPlexComputeCellGeometry(dm, point, v0, J, NULL, &detJ);
223:       geom.v0   = v0;
224:       geom.J    = J;
225:       geom.detJ = &detJ;
226:       for (f = 0, v = 0; f < numFields; ++f) {
227:         void * const ctx = ctxs ? ctxs[f] : NULL;
228:         PetscFEGetNumComponents(fe[f], &numComp);
229:         PetscDualSpaceGetDimension(sp[f], &spDim);
230:         for (d = 0; d < spDim; ++d) {
231:           if (funcs[f]) {
232:             PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);
233:           } else {
234:             for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0;
235:           }
236:           v += numComp;
237:         }
238:       }
239:       DMPlexVecSetClosure(dm, section, localX, point, values, mode);
240:     }
241:     ISRestoreIndices(pointIS, &points);
242:     ISDestroy(&pointIS);
243:   }
244:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
245:   PetscFree3(sp,v0,J);
246:   return(0);
247: }

251: PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
252: {
253:   PetscDualSpace *sp;
254:   PetscSection    section;
255:   PetscScalar    *values;
256:   PetscReal      *v0, *J, detJ;
257:   PetscInt        numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, c, f, d, v, comp;
258:   PetscErrorCode  ierr;

261:   DMGetDefaultSection(dm, &section);
262:   PetscSectionGetNumFields(section, &numFields);
263:   PetscMalloc1(numFields, &sp);
264:   for (f = 0; f < numFields; ++f) {
265:     PetscFEGetDualSpace(fe[f], &sp[f]);
266:     PetscFEGetNumComponents(fe[f], &numComp);
267:     PetscDualSpaceGetDimension(sp[f], &spDim);
268:     totDim += spDim*numComp;
269:   }
270:   DMPlexGetDimension(dm, &dim);
271:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
272:   DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
273:   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
274:   DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
275:   PetscMalloc2(dim,&v0,dim*dim,&J);
276:   for (c = cStart; c < cEnd; ++c) {
277:     PetscCellGeometry geom;

279:     DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);
280:     geom.v0   = v0;
281:     geom.J    = J;
282:     geom.detJ = &detJ;
283:     for (f = 0, v = 0; f < numFields; ++f) {
284:       void * const ctx = ctxs ? ctxs[f] : NULL;
285:       PetscFEGetNumComponents(fe[f], &numComp);
286:       PetscDualSpaceGetDimension(sp[f], &spDim);
287:       for (d = 0; d < spDim; ++d) {
288:         if (funcs[f]) {
289:           PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);
290:         } else {
291:           for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0;
292:         }
293:         v += numComp;
294:       }
295:     }
296:     DMPlexVecSetClosure(dm, section, localX, c, values, mode);
297:   }
298:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
299:   PetscFree2(v0,J);
300:   PetscFree(sp);
301:   return(0);
302: }

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

309:   Input Parameters:
310: + dm      - The DM
311: . fe      - The PetscFE associated with the field
312: . funcs   - The coordinate functions to evaluate, one per field
313: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
314: - mode    - The insertion mode for values

316:   Output Parameter:
317: . X - vector

319:   Level: developer

321: .seealso: DMPlexComputeL2Diff()
322: @*/
323: PetscErrorCode DMPlexProjectFunction(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
324: {
325:   Vec            localX;

330:   DMGetLocalVector(dm, &localX);
331:   DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, mode, localX);
332:   DMLocalToGlobalBegin(dm, localX, mode, X);
333:   DMLocalToGlobalEnd(dm, localX, mode, X);
334:   DMRestoreLocalVector(dm, &localX);
335:   return(0);
336: }

340: PetscErrorCode DMPlexInsertBoundaryValuesFEM(DM dm, Vec localX)
341: {
342:   void        (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx);
343:   void         **ctxs;
344:   PetscFE       *fe;
345:   PetscInt       numFields, f, numBd, b;

351:   DMGetNumFields(dm, &numFields);
352:   PetscMalloc3(numFields,&fe,numFields,&funcs,numFields,&ctxs);
353:   for (f = 0; f < numFields; ++f) {DMGetField(dm, f, (PetscObject *) &fe[f]);}
354:   /* OPT: Could attempt to do multiple BCs at once */
355:   DMPlexGetNumBoundary(dm, &numBd);
356:   for (b = 0; b < numBd; ++b) {
357:     DMLabel         label;
358:     const PetscInt *ids;
359:     const char     *name;
360:     PetscInt        numids, field;
361:     PetscBool       isEssential;
362:     void          (*func)();
363:     void           *ctx;

365:     /* TODO: We need to set only the part indicated by the ids */
366:     DMPlexGetBoundary(dm, b, &isEssential, &name, &field, &func, &numids, &ids, &ctx);
367:     DMPlexGetLabel(dm, name, &label);
368:     for (f = 0; f < numFields; ++f) {
369:       funcs[f] = field == f ? (void (*)(const PetscReal[], PetscScalar *, void *)) func : NULL;
370:       ctxs[f]  = field == f ? ctx : NULL;
371:     }
372:     DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, fe, funcs, ctxs, INSERT_BC_VALUES, localX);
373:   }
374:   PetscFree3(fe,funcs,ctxs);
375:   return(0);
376: }

378: /* Assuming dim == 3 */
379: typedef struct {
380:   PetscScalar normal[3];   /* Area-scaled normals */
381:   PetscScalar centroid[3]; /* Location of centroid (quadrature point) */
382:   PetscScalar grad[2][3];  /* Face contribution to gradient in left and right cell */
383: } FaceGeom;

387: PetscErrorCode DMPlexInsertBoundaryValuesFVM(DM dm, PetscReal time, Vec locX)
388: {
389:   DM                 dmFace;
390:   Vec                faceGeometry;
391:   DMLabel            label;
392:   const PetscScalar *facegeom;
393:   PetscScalar       *x;
394:   PetscInt           numBd, b, fStart, fEnd;
395:   PetscErrorCode     ierr;

398:   /* TODO Pull this geometry calculation into the library */
399:   PetscObjectQuery((PetscObject) dm, "FaceGeometry", (PetscObject *) &faceGeometry);
400:   DMPlexGetLabel(dm, "Face Sets", &label);
401:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
402:   DMPlexGetNumBoundary(dm, &numBd);
403:   VecGetDM(faceGeometry, &dmFace);
404:   VecGetArrayRead(faceGeometry, &facegeom);
405:   VecGetArray(locX, &x);
406:   for (b = 0; b < numBd; ++b) {
407:     PetscErrorCode (*func)(PetscReal,const PetscScalar*,const PetscScalar*,const PetscScalar*,PetscScalar*,void*);
408:     const PetscInt  *ids;
409:     PetscInt         numids, i;
410:     void            *ctx;

412:     DMPlexGetBoundary(dm, b, NULL, NULL, NULL, (void (**)()) &func, &numids, &ids, &ctx);
413:     for (i = 0; i < numids; ++i) {
414:       IS              faceIS;
415:       const PetscInt *faces;
416:       PetscInt        numFaces, f;

418:       DMLabelGetStratumIS(label, ids[i], &faceIS);
419:       if (!faceIS) continue; /* No points with that id on this process */
420:       ISGetLocalSize(faceIS, &numFaces);
421:       ISGetIndices(faceIS, &faces);
422:       for (f = 0; f < numFaces; ++f) {
423:         const PetscInt     face = faces[f], *cells;
424:         const PetscScalar *xI;
425:         PetscScalar       *xG;
426:         const FaceGeom    *fg;

428:         if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
429:         DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
430:         DMPlexGetSupport(dm, face, &cells);
431:         DMPlexPointLocalRead(dm, cells[0], x, &xI);
432:         DMPlexPointLocalRef(dm, cells[1], x, &xG);
433:         (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);
434:       }
435:       ISRestoreIndices(faceIS, &faces);
436:       ISDestroy(&faceIS);
437:     }
438:   }
439:   VecRestoreArrayRead(faceGeometry, &facegeom);
440:   VecRestoreArray(locX, &x);
441:   return(0);
442: }

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

449:   Input Parameters:
450: + dm    - The DM
451: . fe    - The PetscFE object for each field
452: . funcs - The functions to evaluate for each field component
453: . ctxs  - Optional array of contexts to pass to each function, or NULL.
454: - X     - The coefficient vector u_h

456:   Output Parameter:
457: . diff - The diff ||u - u_h||_2

459:   Level: developer

461: .seealso: DMPlexProjectFunction(), DMPlexComputeL2GradientDiff()
462: @*/
463: PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
464: {
465:   const PetscInt  debug = 0;
466:   PetscSection    section;
467:   PetscQuadrature quad;
468:   Vec             localX;
469:   PetscScalar    *funcVal;
470:   PetscReal      *coords, *v0, *J, *invJ, detJ;
471:   PetscReal       localDiff = 0.0;
472:   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
473:   PetscErrorCode  ierr;

476:   DMPlexGetDimension(dm, &dim);
477:   DMGetDefaultSection(dm, &section);
478:   PetscSectionGetNumFields(section, &numFields);
479:   DMGetLocalVector(dm, &localX);
480:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
481:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
482:   for (field = 0; field < numFields; ++field) {
483:     PetscInt Nc;

485:     PetscFEGetNumComponents(fe[field], &Nc);
486:     numComponents += Nc;
487:   }
488:   DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, INSERT_BC_VALUES, localX);
489:   PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
490:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
491:   PetscFEGetQuadrature(fe[0], &quad);
492:   for (c = cStart; c < cEnd; ++c) {
493:     PetscScalar *x = NULL;
494:     PetscReal    elemDiff = 0.0;

496:     DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
497:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
498:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

500:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
501:       void * const     ctx = ctxs ? ctxs[field] : NULL;
502:       const PetscReal *quadPoints, *quadWeights;
503:       PetscReal       *basis;
504:       PetscInt         numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;

506:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
507:       PetscFEGetDimension(fe[field], &numBasisFuncs);
508:       PetscFEGetNumComponents(fe[field], &numBasisComps);
509:       PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);
510:       if (debug) {
511:         char title[1024];
512:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
513:         DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);
514:       }
515:       for (q = 0; q < numQuadPoints; ++q) {
516:         for (d = 0; d < dim; d++) {
517:           coords[d] = v0[d];
518:           for (e = 0; e < dim; e++) {
519:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
520:           }
521:         }
522:         (*funcs[field])(coords, funcVal, ctx);
523:         for (fc = 0; fc < numBasisComps; ++fc) {
524:           PetscScalar interpolant = 0.0;

526:           for (f = 0; f < numBasisFuncs; ++f) {
527:             const PetscInt fidx = f*numBasisComps+fc;
528:             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
529:           }
530:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
531:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
532:         }
533:       }
534:       comp        += numBasisComps;
535:       fieldOffset += numBasisFuncs*numBasisComps;
536:     }
537:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
538:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
539:     localDiff += elemDiff;
540:   }
541:   PetscFree5(funcVal,coords,v0,J,invJ);
542:   DMRestoreLocalVector(dm, &localX);
543:   MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
544:   *diff = PetscSqrtReal(*diff);
545:   return(0);
546: }

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

553:   Input Parameters:
554: + dm    - The DM
555: . fe    - The PetscFE object for each field
556: . funcs - The gradient functions to evaluate for each field component
557: . ctxs  - Optional array of contexts to pass to each function, or NULL.
558: . X     - The coefficient vector u_h
559: - n     - The vector to project along

561:   Output Parameter:
562: . diff - The diff ||(grad u - grad u_h) . n||_2

564:   Level: developer

566: .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
567: @*/
568: PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
569: {
570:   const PetscInt  debug = 0;
571:   PetscSection    section;
572:   PetscQuadrature quad;
573:   Vec             localX;
574:   PetscScalar    *funcVal, *interpolantVec;
575:   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
576:   PetscReal       localDiff = 0.0;
577:   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
578:   PetscErrorCode  ierr;

581:   DMPlexGetDimension(dm, &dim);
582:   DMGetDefaultSection(dm, &section);
583:   PetscSectionGetNumFields(section, &numFields);
584:   DMGetLocalVector(dm, &localX);
585:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
586:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
587:   for (field = 0; field < numFields; ++field) {
588:     PetscInt Nc;

590:     PetscFEGetNumComponents(fe[field], &Nc);
591:     numComponents += Nc;
592:   }
593:   /* DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
594:   PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
595:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
596:   PetscFEGetQuadrature(fe[0], &quad);
597:   for (c = cStart; c < cEnd; ++c) {
598:     PetscScalar *x = NULL;
599:     PetscReal    elemDiff = 0.0;

601:     DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
602:     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
603:     DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);

605:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
606:       void * const     ctx = ctxs ? ctxs[field] : NULL;
607:       const PetscReal *quadPoints, *quadWeights;
608:       PetscReal       *basisDer;
609:       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;

611:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
612:       PetscFEGetDimension(fe[field], &Nb);
613:       PetscFEGetNumComponents(fe[field], &Ncomp);
614:       PetscFEGetDefaultTabulation(fe[field], NULL, &basisDer, NULL);
615:       if (debug) {
616:         char title[1024];
617:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
618:         DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);
619:       }
620:       for (q = 0; q < numQuadPoints; ++q) {
621:         for (d = 0; d < dim; d++) {
622:           coords[d] = v0[d];
623:           for (e = 0; e < dim; e++) {
624:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
625:           }
626:         }
627:         (*funcs[field])(coords, n, funcVal, ctx);
628:         for (fc = 0; fc < Ncomp; ++fc) {
629:           PetscScalar interpolant = 0.0;

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

635:             for (d = 0; d < dim; ++d) {
636:               realSpaceDer[d] = 0.0;
637:               for (g = 0; g < dim; ++g) {
638:                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
639:               }
640:               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
641:             }
642:           }
643:           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
644:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
645:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
646:         }
647:       }
648:       comp        += Ncomp;
649:       fieldOffset += Nb*Ncomp;
650:     }
651:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
652:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
653:     localDiff += elemDiff;
654:   }
655:   PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
656:   DMRestoreLocalVector(dm, &localX);
657:   MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
658:   *diff = PetscSqrtReal(*diff);
659:   return(0);
660: }

664: /*@
665:   DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user

667:   Input Parameters:
668: + dm - The mesh
669: . X  - Local input vector
670: - user - The user context

672:   Output Parameter:
673: . F  - Local output vector

675:   Note:
676:   The first member of the user context must be an FEMContext.

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

681:   Level: developer

683: .seealso: DMPlexComputeJacobianActionFEM()
684: @*/
685: PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
686: {
687:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
688:   PetscFEM         *fem   = (PetscFEM *) user;
689:   PetscFE          *fe    = fem->fe;
690:   PetscFE          *feAux = fem->feAux;
691:   PetscFE          *feBd  = fem->feBd;
692:   const char       *name  = "Residual";
693:   DM                dmAux;
694:   Vec               A;
695:   PetscQuadrature   q;
696:   PetscCellGeometry geom;
697:   PetscSection      section, sectionAux;
698:   PetscReal        *v0, *J, *invJ, *detJ;
699:   PetscScalar      *elemVec, *u, *a = NULL;
700:   PetscInt          dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c;
701:   PetscInt          cellDof = 0, numComponents = 0;
702:   PetscInt          cellDofAux = 0, numComponentsAux = 0;
703:   PetscErrorCode    ierr;

706:   PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
707:   DMPlexGetDimension(dm, &dim);
708:   DMGetDefaultSection(dm, &section);
709:   PetscSectionGetNumFields(section, &Nf);
710:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
711:   numCells = cEnd - cStart;
712:   for (f = 0; f < Nf; ++f) {
713:     PetscInt Nb, Nc;

715:     PetscFEGetDimension(fe[f], &Nb);
716:     PetscFEGetNumComponents(fe[f], &Nc);
717:     cellDof       += Nb*Nc;
718:     numComponents += Nc;
719:   }
720:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
721:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
722:   if (dmAux) {
723:     DMGetDefaultSection(dmAux, &sectionAux);
724:     PetscSectionGetNumFields(sectionAux, &NfAux);
725:   }
726:   for (f = 0; f < NfAux; ++f) {
727:     PetscInt Nb, Nc;

729:     PetscFEGetDimension(feAux[f], &Nb);
730:     PetscFEGetNumComponents(feAux[f], &Nc);
731:     cellDofAux       += Nb*Nc;
732:     numComponentsAux += Nc;
733:   }
734:   DMPlexInsertBoundaryValuesFEM(dm, X);
735:   VecSet(F, 0.0);
736:   PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);
737:   if (dmAux) {PetscMalloc1(numCells*cellDofAux, &a);}
738:   for (c = cStart; c < cEnd; ++c) {
739:     PetscScalar *x = NULL;
740:     PetscInt     i;

742:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
743:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
744:     DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
745:     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
746:     DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
747:     if (dmAux) {
748:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
749:       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
750:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
751:     }
752:   }
753:   for (f = 0; f < Nf; ++f) {
754:     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[f];
755:     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[f];
756:     PetscInt numQuadPoints, Nb;
757:     /* Conforming batches */
758:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
759:     /* Remainder */
760:     PetscInt Nr, offset;

762:     PetscFEGetQuadrature(fe[f], &q);
763:     PetscFEGetDimension(fe[f], &Nb);
764:     PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);
765:     PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
766:     blockSize = Nb*numQuadPoints;
767:     batchSize = numBlocks * blockSize;
768:      PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);
769:     numChunks = numCells / (numBatches*batchSize);
770:     Ne        = numChunks*numBatches*batchSize;
771:     Nr        = numCells % (numBatches*batchSize);
772:     offset    = numCells - Nr;
773:     geom.v0   = v0;
774:     geom.J    = J;
775:     geom.invJ = invJ;
776:     geom.detJ = detJ;
777:     PetscFEIntegrateResidual(fe[f], Ne, Nf, fe, f, geom, u, NfAux, feAux, a, f0, f1, elemVec);
778:     geom.v0   = &v0[offset*dim];
779:     geom.J    = &J[offset*dim*dim];
780:     geom.invJ = &invJ[offset*dim*dim];
781:     geom.detJ = &detJ[offset];
782:     PetscFEIntegrateResidual(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);
783:   }
784:   for (c = cStart; c < cEnd; ++c) {
785:     if (mesh->printFEM > 1) {DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);}
786:     DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);
787:   }
788:   PetscFree6(u,v0,J,invJ,detJ,elemVec);
789:   if (dmAux) {PetscFree(a);}
790:   if (feBd) {
791:     DMLabel  depth;
792:     PetscInt numBd, bd;

794:     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
795:       PetscInt Nb, Nc;

797:       PetscFEGetDimension(feBd[f], &Nb);
798:       PetscFEGetNumComponents(feBd[f], &Nc);
799:       cellDof       += Nb*Nc;
800:       numComponents += Nc;
801:     }
802:     DMPlexGetDepthLabel(dm, &depth);
803:     DMPlexGetNumBoundary(dm, &numBd);
804:     for (bd = 0; bd < numBd; ++bd) {
805:       const char     *bdLabel;
806:       DMLabel         label;
807:       IS              pointIS;
808:       const PetscInt *points;
809:       const PetscInt *values;
810:       PetscReal      *n;
811:       PetscInt        field, numValues, numPoints, p, dep, numFaces;
812:       PetscBool       isEssential;

814:       DMPlexGetBoundary(dm, bd, &isEssential, &bdLabel, &field, NULL, &numValues, &values, NULL);
815:       if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
816:       DMPlexGetLabel(dm, bdLabel, &label);
817:       DMLabelGetStratumSize(label, 1, &numPoints);
818:       DMLabelGetStratumIS(label, 1, &pointIS);
819:       ISGetIndices(pointIS, &points);
820:       for (p = 0, numFaces = 0; p < numPoints; ++p) {
821:         DMLabelGetValue(depth, points[p], &dep);
822:         if (dep == dim-1) ++numFaces;
823:       }
824:       PetscMalloc7(numFaces*cellDof,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*cellDof,&elemVec);
825:       for (p = 0, f = 0; p < numPoints; ++p) {
826:         const PetscInt point = points[p];
827:         PetscScalar   *x     = NULL;
828:         PetscInt       i;

830:         DMLabelGetValue(depth, points[p], &dep);
831:         if (dep != dim-1) continue;
832:         DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);
833:         DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
834:         if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
835:         DMPlexVecGetClosure(dm, section, X, point, NULL, &x);
836:         for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
837:         DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);
838:         ++f;
839:       }
840:       for (f = 0; f < Nf; ++f) {
841:         void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[f];
842:         void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[f];
843:         PetscInt numQuadPoints, Nb;
844:         /* Conforming batches */
845:         PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
846:         /* Remainder */
847:         PetscInt Nr, offset;

849:         PetscFEGetQuadrature(feBd[f], &q);
850:         PetscFEGetDimension(feBd[f], &Nb);
851:         PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);
852:         PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
853:         blockSize = Nb*numQuadPoints;
854:         batchSize = numBlocks * blockSize;
855:          PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);
856:         numChunks = numFaces / (numBatches*batchSize);
857:         Ne        = numChunks*numBatches*batchSize;
858:         Nr        = numFaces % (numBatches*batchSize);
859:         offset    = numFaces - Nr;
860:         geom.v0   = v0;
861:         geom.n    = n;
862:         geom.J    = J;
863:         geom.invJ = invJ;
864:         geom.detJ = detJ;
865:         PetscFEIntegrateBdResidual(feBd[f], Ne, Nf, feBd, f, geom, u, 0, NULL, NULL, f0, f1, elemVec);
866:         geom.v0   = &v0[offset*dim];
867:         geom.n    = &n[offset*dim];
868:         geom.J    = &J[offset*dim*dim];
869:         geom.invJ = &invJ[offset*dim*dim];
870:         geom.detJ = &detJ[offset];
871:         PetscFEIntegrateBdResidual(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);
872:       }
873:       for (p = 0, f = 0; p < numPoints; ++p) {
874:         const PetscInt point = points[p];

876:         DMLabelGetValue(depth, point, &dep);
877:         if (dep != dim-1) continue;
878:         if (mesh->printFEM > 1) {DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);}
879:         DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);
880:         ++f;
881:       }
882:       ISRestoreIndices(pointIS, &points);
883:       ISDestroy(&pointIS);
884:       PetscFree7(u,v0,n,J,invJ,detJ,elemVec);
885:     }
886:   }
887:   if (mesh->printFEM) {DMPrintLocalVec(dm, name, mesh->printTol, F);}
888:   PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
889:   return(0);
890: }

894: /*@
895:   DMPlexComputeIFunctionFEM - Form the local implicit function F from the local input X, X_t using pointwise functions specified by the user

897:   Input Parameters:
898: + dm - The mesh
899: . time - The current time
900: . X  - Local input vector
901: . X_t  - Time derivative of the local input vector
902: - user - The user context

904:   Output Parameter:
905: . F  - Local output vector

907:   Note:
908:   The first member of the user context must be an FEMContext.

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

913:   Level: developer

915: .seealso: DMPlexComputeResidualFEM()
916: @*/
917: PetscErrorCode DMPlexComputeIFunctionFEM(DM dm, PetscReal time, Vec X, Vec X_t, Vec F, void *user)
918: {
919:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
920:   PetscFEM         *fem   = (PetscFEM *) user;
921:   PetscFE          *fe    = fem->fe;
922:   PetscFE          *feAux = fem->feAux;
923:   PetscFE          *feBd  = fem->feBd;
924:   const char       *name  = "Residual";
925:   DM                dmAux;
926:   Vec               A;
927:   PetscQuadrature   q;
928:   PetscCellGeometry geom;
929:   PetscSection      section, sectionAux;
930:   PetscReal        *v0, *J, *invJ, *detJ;
931:   PetscScalar      *elemVec, *u, *u_t, *a = NULL;
932:   PetscInt          dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c;
933:   PetscInt          cellDof = 0, numComponents = 0;
934:   PetscInt          cellDofAux = 0, numComponentsAux = 0;
935:   PetscErrorCode    ierr;

938:   PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
939:   DMPlexGetDimension(dm, &dim);
940:   DMGetDefaultSection(dm, &section);
941:   PetscSectionGetNumFields(section, &Nf);
942:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
943:   numCells = cEnd - cStart;
944:   for (f = 0; f < Nf; ++f) {
945:     PetscInt Nb, Nc;

947:     PetscFEGetDimension(fe[f], &Nb);
948:     PetscFEGetNumComponents(fe[f], &Nc);
949:     cellDof       += Nb*Nc;
950:     numComponents += Nc;
951:   }
952:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
953:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
954:   if (dmAux) {
955:     DMGetDefaultSection(dmAux, &sectionAux);
956:     PetscSectionGetNumFields(sectionAux, &NfAux);
957:   }
958:   for (f = 0; f < NfAux; ++f) {
959:     PetscInt Nb, Nc;

961:     PetscFEGetDimension(feAux[f], &Nb);
962:     PetscFEGetNumComponents(feAux[f], &Nc);
963:     cellDofAux       += Nb*Nc;
964:     numComponentsAux += Nc;
965:   }
966:   DMPlexInsertBoundaryValuesFEM(dm, X);
967:   VecSet(F, 0.0);
968:   PetscMalloc7(numCells*cellDof,&u,numCells*cellDof,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);
969:   if (dmAux) {PetscMalloc1(numCells*cellDofAux, &a);}
970:   for (c = cStart; c < cEnd; ++c) {
971:     PetscScalar *x = NULL, *x_t = NULL;
972:     PetscInt     i;

974:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
975:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
976:     DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
977:     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
978:     DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
979:     DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
980:     for (i = 0; i < cellDof; ++i) u_t[c*cellDof+i] = x_t[i];
981:     DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
982:     if (dmAux) {
983:       PetscScalar *x_a = NULL;
984:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x_a);
985:       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x_a[i];
986:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x_a);
987:     }
988:   }
989:   for (f = 0; f < Nf; ++f) {
990:     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0IFuncs[f];
991:     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1IFuncs[f];
992:     PetscInt numQuadPoints, Nb;
993:     /* Conforming batches */
994:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
995:     /* Remainder */
996:     PetscInt Nr, offset;

998:     PetscFEGetQuadrature(fe[f], &q);
999:     PetscFEGetDimension(fe[f], &Nb);
1000:     PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);
1001:     PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1002:     blockSize = Nb*numQuadPoints;
1003:     batchSize = numBlocks * blockSize;
1004:      PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);
1005:     numChunks = numCells / (numBatches*batchSize);
1006:     Ne        = numChunks*numBatches*batchSize;
1007:     Nr        = numCells % (numBatches*batchSize);
1008:     offset    = numCells - Nr;
1009:     geom.v0   = v0;
1010:     geom.J    = J;
1011:     geom.invJ = invJ;
1012:     geom.detJ = detJ;
1013:     PetscFEIntegrateIFunction(fe[f], Ne, Nf, fe, f, geom, u, u_t, NfAux, feAux, a, f0, f1, elemVec);
1014:     geom.v0   = &v0[offset*dim];
1015:     geom.J    = &J[offset*dim*dim];
1016:     geom.invJ = &invJ[offset*dim*dim];
1017:     geom.detJ = &detJ[offset];
1018:     PetscFEIntegrateIFunction(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], &u_t[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);
1019:   }
1020:   for (c = cStart; c < cEnd; ++c) {
1021:     if (mesh->printFEM > 1) {DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);}
1022:     DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);
1023:   }
1024:   PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);
1025:   if (dmAux) {PetscFree(a);}
1026:   if (feBd) {
1027:     DMLabel         label, depth;
1028:     IS              pointIS;
1029:     const PetscInt *points;
1030:     PetscInt        dep, numPoints, p, numFaces;
1031:     PetscReal      *n;

1033:     DMPlexGetLabel(dm, "boundary", &label);
1034:     DMPlexGetDepthLabel(dm, &depth);
1035:     DMLabelGetStratumSize(label, 1, &numPoints);
1036:     DMLabelGetStratumIS(label, 1, &pointIS);
1037:     ISGetIndices(pointIS, &points);
1038:     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
1039:       PetscInt Nb, Nc;

1041:       PetscFEGetDimension(feBd[f], &Nb);
1042:       PetscFEGetNumComponents(feBd[f], &Nc);
1043:       cellDof       += Nb*Nc;
1044:       numComponents += Nc;
1045:     }
1046:     for (p = 0, numFaces = 0; p < numPoints; ++p) {
1047:       DMLabelGetValue(depth, points[p], &dep);
1048:       if (dep == dim-1) ++numFaces;
1049:     }
1050:     PetscMalloc7(numFaces*cellDof,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*cellDof,&elemVec);
1051:     PetscMalloc1(numFaces*cellDof,&u_t);
1052:     for (p = 0, f = 0; p < numPoints; ++p) {
1053:       const PetscInt point = points[p];
1054:       PetscScalar   *x     = NULL;
1055:       PetscInt       i;

1057:       DMLabelGetValue(depth, points[p], &dep);
1058:       if (dep != dim-1) continue;
1059:       DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);
1060:       DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
1061:       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
1062:       DMPlexVecGetClosure(dm, section, X, point, NULL, &x);
1063:       for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
1064:       DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);
1065:       DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);
1066:       for (i = 0; i < cellDof; ++i) u_t[f*cellDof+i] = x[i];
1067:       DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);
1068:       ++f;
1069:     }
1070:     for (f = 0; f < Nf; ++f) {
1071:       void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdIFuncs[f];
1072:       void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdIFuncs[f];
1073:       PetscInt numQuadPoints, Nb;
1074:       /* Conforming batches */
1075:       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1076:       /* Remainder */
1077:       PetscInt Nr, offset;

1079:       PetscFEGetQuadrature(feBd[f], &q);
1080:       PetscFEGetDimension(feBd[f], &Nb);
1081:       PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);
1082:       PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1083:       blockSize = Nb*numQuadPoints;
1084:       batchSize = numBlocks * blockSize;
1085:        PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);
1086:       numChunks = numFaces / (numBatches*batchSize);
1087:       Ne        = numChunks*numBatches*batchSize;
1088:       Nr        = numFaces % (numBatches*batchSize);
1089:       offset    = numFaces - Nr;
1090:       geom.v0   = v0;
1091:       geom.n    = n;
1092:       geom.J    = J;
1093:       geom.invJ = invJ;
1094:       geom.detJ = detJ;
1095:       PetscFEIntegrateBdIFunction(feBd[f], Ne, Nf, feBd, f, geom, u, u_t, 0, NULL, NULL, f0, f1, elemVec);
1096:       geom.v0   = &v0[offset*dim];
1097:       geom.n    = &n[offset*dim];
1098:       geom.J    = &J[offset*dim*dim];
1099:       geom.invJ = &invJ[offset*dim*dim];
1100:       geom.detJ = &detJ[offset];
1101:       PetscFEIntegrateBdIFunction(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], &u_t[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);
1102:     }
1103:     for (p = 0, f = 0; p < numPoints; ++p) {
1104:       const PetscInt point = points[p];

1106:       DMLabelGetValue(depth, point, &dep);
1107:       if (dep != dim-1) continue;
1108:       if (mesh->printFEM > 1) {DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);}
1109:       DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);
1110:       ++f;
1111:     }
1112:     ISRestoreIndices(pointIS, &points);
1113:     ISDestroy(&pointIS);
1114:     PetscFree7(u,v0,n,J,invJ,detJ,elemVec);
1115:     PetscFree(u_t);
1116:   }
1117:   if (mesh->printFEM) {DMPrintLocalVec(dm, name, mesh->printTol, F);}
1118:   PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
1119:   return(0);
1120: }

1124: /*@C
1125:   DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user

1127:   Input Parameters:
1128: + dm - The mesh
1129: . J  - The Jacobian shell matrix
1130: . X  - Local input vector
1131: - user - The user context

1133:   Output Parameter:
1134: . F  - Local output vector

1136:   Note:
1137:   The first member of the user context must be an FEMContext.

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

1142:   Level: developer

1144: .seealso: DMPlexComputeResidualFEM()
1145: @*/
1146: PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
1147: {
1148:   DM_Plex          *mesh = (DM_Plex *) dm->data;
1149:   PetscFEM         *fem  = (PetscFEM *) user;
1150:   PetscFE          *fe   = fem->fe;
1151:   PetscQuadrature   quad;
1152:   PetscCellGeometry geom;
1153:   PetscSection      section;
1154:   JacActionCtx     *jctx;
1155:   PetscReal        *v0, *J, *invJ, *detJ;
1156:   PetscScalar      *elemVec, *u, *a;
1157:   PetscInt          dim, numFields, field, numCells, cStart, cEnd, c;
1158:   PetscInt          cellDof = 0;
1159:   PetscErrorCode    ierr;

1162:   /* PetscLogEventBegin(DMPLEX_JacobianActionFEM,dm,0,0,0); */
1163:   MatShellGetContext(Jac, &jctx);
1164:   DMPlexGetDimension(dm, &dim);
1165:   DMGetDefaultSection(dm, &section);
1166:   PetscSectionGetNumFields(section, &numFields);
1167:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1168:   numCells = cEnd - cStart;
1169:   for (field = 0; field < numFields; ++field) {
1170:     PetscInt Nb, Nc;

1172:     PetscFEGetDimension(fe[field], &Nb);
1173:     PetscFEGetNumComponents(fe[field], &Nc);
1174:     cellDof += Nb*Nc;
1175:   }
1176:   VecSet(F, 0.0);
1177:   PetscMalloc7(numCells*cellDof,&u,numCells*cellDof,&a,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);
1178:   for (c = cStart; c < cEnd; ++c) {
1179:     PetscScalar *x = NULL;
1180:     PetscInt     i;

1182:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
1183:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1184:     DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);
1185:     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
1186:     DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);
1187:     DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);
1188:     for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
1189:     DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);
1190:   }
1191:   for (field = 0; field < numFields; ++field) {
1192:     PetscInt numQuadPoints, Nb;
1193:     /* Conforming batches */
1194:     PetscInt numBlocks  = 1;
1195:     PetscInt numBatches = 1;
1196:     PetscInt numChunks, Ne, blockSize, batchSize;
1197:     /* Remainder */
1198:     PetscInt Nr, offset;

1200:     PetscFEGetQuadrature(fe[field], &quad);
1201:     PetscFEGetDimension(fe[field], &Nb);
1202:     PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
1203:     blockSize = Nb*numQuadPoints;
1204:     batchSize = numBlocks * blockSize;
1205:     numChunks = numCells / (numBatches*batchSize);
1206:     Ne        = numChunks*numBatches*batchSize;
1207:     Nr        = numCells % (numBatches*batchSize);
1208:     offset    = numCells - Nr;
1209:     geom.v0   = v0;
1210:     geom.J    = J;
1211:     geom.invJ = invJ;
1212:     geom.detJ = detJ;
1213:     PetscFEIntegrateJacobianAction(fe[field], Ne, numFields, fe, field, geom, u, a, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);
1214:     geom.v0   = &v0[offset*dim];
1215:     geom.J    = &J[offset*dim*dim];
1216:     geom.invJ = &invJ[offset*dim*dim];
1217:     geom.detJ = &detJ[offset];
1218:     PetscFEIntegrateJacobianAction(fe[field], Nr, numFields, fe, field, geom, &u[offset*cellDof], &a[offset*cellDof],
1219:                                           fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);
1220:   }
1221:   for (c = cStart; c < cEnd; ++c) {
1222:     if (mesh->printFEM > 1) {DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);}
1223:     DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);
1224:   }
1225:   PetscFree7(u,a,v0,J,invJ,detJ,elemVec);
1226:   if (mesh->printFEM) {
1227:     PetscMPIInt rank, numProcs;
1228:     PetscInt    p;

1230:     MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1231:     MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
1232:     PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");
1233:     for (p = 0; p < numProcs; ++p) {
1234:       if (p == rank) {VecView(F, PETSC_VIEWER_STDOUT_SELF);}
1235:       PetscBarrier((PetscObject) dm);
1236:     }
1237:   }
1238:   /* PetscLogEventEnd(DMPLEX_JacobianActionFEM,dm,0,0,0); */
1239:   return(0);
1240: }

1244: /*@
1245:   DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.

1247:   Input Parameters:
1248: + dm - The mesh
1249: . X  - Local input vector
1250: - user - The user context

1252:   Output Parameter:
1253: . Jac  - Jacobian matrix

1255:   Note:
1256:   The first member of the user context must be an FEMContext.

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

1261:   Level: developer

1263: .seealso: FormFunctionLocal()
1264: @*/
1265: PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1266: {
1267:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1268:   PetscFEM         *fem   = (PetscFEM *) user;
1269:   PetscFE          *fe    = fem->fe;
1270:   PetscFE          *feAux = fem->feAux;
1271:   PetscFE          *feBd  = fem->feBd;
1272:   const char       *name  = "Jacobian";
1273:   DM                dmAux;
1274:   Vec               A;
1275:   PetscQuadrature   quad;
1276:   PetscCellGeometry geom;
1277:   PetscSection      section, globalSection, sectionAux;
1278:   PetscReal        *v0, *J, *invJ, *detJ;
1279:   PetscScalar      *elemMat, *u, *a;
1280:   PetscInt          dim, Nf, NfAux = 0, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1281:   PetscInt          cellDof = 0, numComponents = 0;
1282:   PetscInt          cellDofAux = 0, numComponentsAux = 0;
1283:   PetscBool         isShell;
1284:   PetscErrorCode    ierr;

1287:   PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
1288:   DMPlexGetDimension(dm, &dim);
1289:   DMGetDefaultSection(dm, &section);
1290:   DMGetDefaultGlobalSection(dm, &globalSection);
1291:   PetscSectionGetNumFields(section, &Nf);
1292:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1293:   numCells = cEnd - cStart;
1294:   for (f = 0; f < Nf; ++f) {
1295:     PetscInt Nb, Nc;

1297:     PetscFEGetDimension(fe[f], &Nb);
1298:     PetscFEGetNumComponents(fe[f], &Nc);
1299:     cellDof       += Nb*Nc;
1300:     numComponents += Nc;
1301:   }
1302:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1303:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1304:   if (dmAux) {
1305:     DMGetDefaultSection(dmAux, &sectionAux);
1306:     PetscSectionGetNumFields(sectionAux, &NfAux);
1307:   }
1308:   for (f = 0; f < NfAux; ++f) {
1309:     PetscInt Nb, Nc;

1311:     PetscFEGetDimension(feAux[f], &Nb);
1312:     PetscFEGetNumComponents(feAux[f], &Nc);
1313:     cellDofAux       += Nb*Nc;
1314:     numComponentsAux += Nc;
1315:   }
1316:   DMPlexInsertBoundaryValuesFEM(dm, X);
1317:   MatZeroEntries(JacP);
1318:   PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof*cellDof,&elemMat);
1319:   if (dmAux) {PetscMalloc1(numCells*cellDofAux, &a);}
1320:   for (c = cStart; c < cEnd; ++c) {
1321:     PetscScalar *x = NULL;
1322:     PetscInt     i;

1324:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
1325:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1326:     DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
1327:     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
1328:     DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
1329:     if (dmAux) {
1330:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1331:       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
1332:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1333:     }
1334:   }
1335:   PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));
1336:   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1337:     PetscInt numQuadPoints, Nb;
1338:     /* Conforming batches */
1339:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1340:     /* Remainder */
1341:     PetscInt Nr, offset;

1343:     PetscFEGetQuadrature(fe[fieldI], &quad);
1344:     PetscFEGetDimension(fe[fieldI], &Nb);
1345:     PetscFEGetTileSizes(fe[fieldI], NULL, &numBlocks, NULL, &numBatches);
1346:     PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
1347:     blockSize = Nb*numQuadPoints;
1348:     batchSize = numBlocks * blockSize;
1349:     PetscFESetTileSizes(fe[fieldI], blockSize, numBlocks, batchSize, numBatches);
1350:     numChunks = numCells / (numBatches*batchSize);
1351:     Ne        = numChunks*numBatches*batchSize;
1352:     Nr        = numCells % (numBatches*batchSize);
1353:     offset    = numCells - Nr;
1354:     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1355:       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ];
1356:       void   (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*Nf+fieldJ];
1357:       void   (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*Nf+fieldJ];
1358:       void   (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*Nf+fieldJ];

1360:       geom.v0   = v0;
1361:       geom.J    = J;
1362:       geom.invJ = invJ;
1363:       geom.detJ = detJ;
1364:       PetscFEIntegrateJacobian(fe[fieldI], Ne, Nf, fe, fieldI, fieldJ, geom, u, NfAux, feAux, a, g0, g1, g2, g3, elemMat);
1365:       geom.v0   = &v0[offset*dim];
1366:       geom.J    = &J[offset*dim*dim];
1367:       geom.invJ = &invJ[offset*dim*dim];
1368:       geom.detJ = &detJ[offset];
1369:       PetscFEIntegrateJacobian(fe[fieldI], Nr, Nf, fe, fieldI, fieldJ, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);
1370:     }
1371:   }
1372:   for (c = cStart; c < cEnd; ++c) {
1373:     if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, cellDof, cellDof, &elemMat[c*cellDof*cellDof]);}
1374:     DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);
1375:   }
1376:   PetscFree6(u,v0,J,invJ,detJ,elemMat);
1377:   if (dmAux) {PetscFree(a);}
1378:   if (feBd) {
1379:     DMLabel  depth;
1380:     PetscInt numBd, bd;

1382:     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
1383:       PetscInt Nb, Nc;

1385:       PetscFEGetDimension(feBd[f], &Nb);
1386:       PetscFEGetNumComponents(feBd[f], &Nc);
1387:       cellDof       += Nb*Nc;
1388:       numComponents += Nc;
1389:     }
1390:     DMPlexGetDepthLabel(dm, &depth);
1391:     DMPlexGetNumBoundary(dm, &numBd);
1392:     for (bd = 0; bd < numBd; ++bd) {
1393:       const char     *bdLabel;
1394:       DMLabel         label;
1395:       IS              pointIS;
1396:       const PetscInt *points;
1397:       const PetscInt *values;
1398:       PetscReal      *n;
1399:       PetscInt        field, numValues, numPoints, p, dep, numFaces;
1400:       PetscBool       isEssential;

1402:       DMPlexGetBoundary(dm, bd, &isEssential, &bdLabel, &field, NULL, &numValues, &values, NULL);
1403:       if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
1404:       DMPlexGetLabel(dm, bdLabel, &label);
1405:       DMLabelGetStratumSize(label, 1, &numPoints);
1406:       DMLabelGetStratumIS(label, 1, &pointIS);
1407:       ISGetIndices(pointIS, &points);
1408:       for (p = 0, numFaces = 0; p < numPoints; ++p) {
1409:         DMLabelGetValue(depth, points[p], &dep);
1410:         if (dep == dim-1) ++numFaces;
1411:       }
1412:       PetscMalloc7(numFaces*cellDof,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*cellDof*cellDof,&elemMat);
1413:       for (p = 0, f = 0; p < numPoints; ++p) {
1414:         const PetscInt point = points[p];
1415:         PetscScalar   *x     = NULL;
1416:         PetscInt       i;

1418:         DMLabelGetValue(depth, points[p], &dep);
1419:         if (dep != dim-1) continue;
1420:         DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);
1421:         DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
1422:         if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
1423:         DMPlexVecGetClosure(dm, section, X, point, NULL, &x);
1424:         for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
1425:         DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);
1426:         ++f;
1427:       }
1428:       PetscMemzero(elemMat, numFaces*cellDof*cellDof * sizeof(PetscScalar));
1429:       for (fieldI = 0; fieldI < Nf; ++fieldI) {
1430:         PetscInt numQuadPoints, Nb;
1431:         /* Conforming batches */
1432:         PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1433:         /* Remainder */
1434:         PetscInt Nr, offset;

1436:         PetscFEGetQuadrature(feBd[fieldI], &quad);
1437:         PetscFEGetDimension(feBd[fieldI], &Nb);
1438:         PetscFEGetTileSizes(feBd[fieldI], NULL, &numBlocks, NULL, &numBatches);
1439:         PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
1440:         blockSize = Nb*numQuadPoints;
1441:         batchSize = numBlocks * blockSize;
1442:          PetscFESetTileSizes(feBd[fieldI], blockSize, numBlocks, batchSize, numBatches);
1443:         numChunks = numFaces / (numBatches*batchSize);
1444:         Ne        = numChunks*numBatches*batchSize;
1445:         Nr        = numFaces % (numBatches*batchSize);
1446:         offset    = numFaces - Nr;
1447:         for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1448:           void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g0BdFuncs[fieldI*Nf+fieldJ];
1449:           void   (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g1BdFuncs[fieldI*Nf+fieldJ];
1450:           void   (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g2BdFuncs[fieldI*Nf+fieldJ];
1451:           void   (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g3BdFuncs[fieldI*Nf+fieldJ];

1453:           geom.v0   = v0;
1454:           geom.n    = n;
1455:           geom.J    = J;
1456:           geom.invJ = invJ;
1457:           geom.detJ = detJ;
1458:           PetscFEIntegrateBdJacobian(feBd[fieldI], Ne, Nf, feBd, fieldI, fieldJ, geom, u, 0, NULL, NULL, g0, g1, g2, g3, elemMat);
1459:           geom.v0   = &v0[offset*dim];
1460:           geom.n    = &n[offset*dim];
1461:           geom.J    = &J[offset*dim*dim];
1462:           geom.invJ = &invJ[offset*dim*dim];
1463:           geom.detJ = &detJ[offset];
1464:           PetscFEIntegrateBdJacobian(feBd[fieldI], Nr, Nf, feBd, fieldI, fieldJ, geom, &u[offset*cellDof], 0, NULL, NULL, g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);
1465:         }
1466:       }
1467:       for (p = 0, f = 0; p < numPoints; ++p) {
1468:         const PetscInt point = points[p];

1470:         DMLabelGetValue(depth, point, &dep);
1471:         if (dep != dim-1) continue;
1472:         if (mesh->printFEM > 1) {DMPrintCellMatrix(point, "BdJacobian", cellDof, cellDof, &elemMat[f*cellDof*cellDof]);}
1473:         DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*cellDof*cellDof], ADD_VALUES);
1474:         ++f;
1475:       }
1476:       ISRestoreIndices(pointIS, &points);
1477:       ISDestroy(&pointIS);
1478:       PetscFree7(u,v0,n,J,invJ,detJ,elemMat);
1479:     }
1480:   }
1481:   MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
1482:   MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
1483:   if (mesh->printFEM) {
1484:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1485:     MatChop(JacP, 1.0e-10);
1486:     MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);
1487:   }
1488:   PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
1489:   PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);
1490:   if (isShell) {
1491:     JacActionCtx *jctx;

1493:     MatShellGetContext(Jac, &jctx);
1494:     VecCopy(X, jctx->u);
1495:   }
1496:   return(0);
1497: }

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

1504:   Input Parameters:
1505: + dmf  - The fine mesh
1506: . dmc  - The coarse mesh
1507: - user - The user context

1509:   Output Parameter:
1510: . In  - The interpolation matrix

1512:   Note:
1513:   The first member of the user context must be an FEMContext.

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

1518:   Level: developer

1520: .seealso: DMPlexComputeJacobianFEM()
1521: @*/
1522: PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1523: {
1524:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1525:   PetscFEM         *fem   = (PetscFEM *) user;
1526:   PetscFE          *fe    = fem->fe;
1527:   const char       *name  = "Interpolator";
1528:   PetscFE          *feRef;
1529:   PetscSection      fsection, fglobalSection;
1530:   PetscSection      csection, cglobalSection;
1531:   PetscScalar      *elemMat;
1532:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c;
1533:   PetscInt          rCellDof = 0, cCellDof = 0;
1534:   PetscErrorCode    ierr;

1537: #if 0
1538:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1539: #endif
1540:   DMPlexGetDimension(dmf, &dim);
1541:   DMGetDefaultSection(dmf, &fsection);
1542:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1543:   DMGetDefaultSection(dmc, &csection);
1544:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1545:   PetscSectionGetNumFields(fsection, &Nf);
1546:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1547:   PetscMalloc1(Nf,&feRef);
1548:   for (f = 0; f < Nf; ++f) {
1549:     PetscInt rNb, cNb, Nc;

1551:     PetscFERefine(fe[f], &feRef[f]);
1552:     PetscFEGetDimension(feRef[f], &rNb);
1553:     PetscFEGetDimension(fe[f], &cNb);
1554:     PetscFEGetNumComponents(fe[f], &Nc);
1555:     rCellDof += rNb*Nc;
1556:     cCellDof += cNb*Nc;
1557:   }
1558:   MatZeroEntries(In);
1559:   PetscMalloc1(rCellDof*cCellDof,&elemMat);
1560:   PetscMemzero(elemMat, rCellDof*cCellDof * sizeof(PetscScalar));
1561:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1562:     PetscDualSpace   Qref;
1563:     PetscQuadrature  f;
1564:     const PetscReal *qpoints, *qweights;
1565:     PetscReal       *points;
1566:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1568:     /* Compose points from all dual basis functionals */
1569:     PetscFEGetNumComponents(fe[fieldI], &Nc);
1570:     PetscFEGetDualSpace(feRef[fieldI], &Qref);
1571:     PetscDualSpaceGetDimension(Qref, &fpdim);
1572:     for (i = 0; i < fpdim; ++i) {
1573:       PetscDualSpaceGetFunctional(Qref, i, &f);
1574:       PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1575:       npoints += Np;
1576:     }
1577:     PetscMalloc1(npoints*dim,&points);
1578:     for (i = 0, k = 0; i < fpdim; ++i) {
1579:       PetscDualSpaceGetFunctional(Qref, i, &f);
1580:       PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1581:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1582:     }

1584:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1585:       PetscReal *B;
1586:       PetscInt   NcJ, cpdim, j;

1588:       /* Evaluate basis at points */
1589:       PetscFEGetNumComponents(fe[fieldJ], &NcJ);
1590:       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);
1591:       PetscFEGetDimension(fe[fieldJ], &cpdim);
1592:       /* For now, fields only interpolate themselves */
1593:       if (fieldI == fieldJ) {
1594:         PetscFEGetTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);
1595:         for (i = 0, k = 0; i < fpdim; ++i) {
1596:           PetscDualSpaceGetFunctional(Qref, i, &f);
1597:           PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1598:           for (p = 0; p < Np; ++p, ++k) {
1599:             for (j = 0; j < cpdim; ++j) {
1600:               for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cCellDof + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p];
1601:             }
1602:           }
1603:         }
1604:         PetscFERestoreTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);
1605:       }
1606:       offsetJ += cpdim*NcJ;
1607:     }
1608:     offsetI += fpdim*Nc;
1609:     PetscFree(points);
1610:   }
1611:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rCellDof, cCellDof, elemMat);}
1612:   for (c = cStart; c < cEnd; ++c) {
1613:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1614:   }
1615:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1616:   PetscFree(feRef);
1617:   PetscFree(elemMat);
1618:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1619:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1620:   if (mesh->printFEM) {
1621:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1622:     MatChop(In, 1.0e-10);
1623:     MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1624:   }
1625: #if 0
1626:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1627: #endif
1628:   return(0);
1629: }

1633: /* The ids can be overridden by the command line option -bc_<boundary name> */
1634: PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], PetscInt field, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx)
1635: {
1636:   DM_Plex       *mesh = (DM_Plex *) dm->data;
1637:   DMBoundary     b;

1641:   PetscNew(&b);
1642:   PetscStrallocpy(name, (char **) &b->name);
1643:   PetscMalloc1(numids, &b->ids);
1644:   PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));
1645:   b->essential   = isEssential;
1646:   b->field       = field;
1647:   b->func        = bcFunc;
1648:   b->numids      = numids;
1649:   b->ctx         = ctx;
1650:   b->next        = mesh->boundary;
1651:   mesh->boundary = b;
1652:   return(0);
1653: }

1657: PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd)
1658: {
1659:   DM_Plex   *mesh = (DM_Plex *) dm->data;
1660:   DMBoundary b    = mesh->boundary;

1663:   *numBd = 0;
1664:   while (b) {++(*numBd); b = b->next;}
1665:   return(0);
1666: }

1670: PetscErrorCode DMPlexGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, PetscInt *field, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx)
1671: {
1672:   DM_Plex   *mesh = (DM_Plex *) dm->data;
1673:   DMBoundary b    = mesh->boundary;
1674:   PetscInt   n    = 0;

1677:   while (b) {
1678:     if (n == bd) break;
1679:     b = b->next;
1680:     ++n;
1681:   }
1682:   if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
1683:   if (isEssential) {
1685:     *isEssential = b->essential;
1686:   }
1687:   if (name) {
1689:     *name = b->name;
1690:   }
1691:   if (field) {
1693:     *field = b->field;
1694:   }
1695:   if (func) {
1697:     *func = b->func;
1698:   }
1699:   if (numids) {
1701:     *numids = b->numids;
1702:   }
1703:   if (ids) {
1705:     *ids = b->ids;
1706:   }
1707:   if (ctx) {
1709:     *ctx = b->ctx;
1710:   }
1711:   return(0);
1712: }