Actual source code: plexfem.c

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

  3: #include <petsc-private/petscfeimpl.h>
  4: #include <petscfv.h>

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

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

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

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

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

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

 91:   Collective on DM

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

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

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

103:   Level: advanced

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

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

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

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

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

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

186: 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)
187: {
188:   PetscDualSpace *sp;
189:   PetscSection    section;
190:   PetscScalar    *values;
191:   PetscReal      *v0, *J, detJ;
192:   PetscBool      *fieldActive;
193:   PetscInt        numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, f, d, v, i, comp;
194:   PetscErrorCode  ierr;

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

218:     DMLabelGetStratumIS(label, ids[i], &pointIS);
219:     ISGetLocalSize(pointIS, &n);
220:     ISGetIndices(pointIS, &points);
221:     for (p = 0; p < n; ++p) {
222:       const PetscInt    point = points[p];
223:       PetscCellGeometry geom;

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

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

266:   DMGetDefaultSection(dm, &section);
267:   PetscSectionGetNumFields(section, &numFields);
268:   PetscMalloc1(numFields, &sp);
269:   for (f = 0; f < numFields; ++f) {
270:     PetscFE fe;

272:     DMGetField(dm, f, (PetscObject *) &fe);
273:     PetscFEGetDualSpace(fe, &sp[f]);
274:     PetscFEGetNumComponents(fe, &numComp);
275:     PetscDualSpaceGetDimension(sp[f], &spDim);
276:     totDim += spDim*numComp;
277:   }
278:   DMPlexGetDimension(dm, &dim);
279:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
280:   DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
281:   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
282:   DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
283:   PetscMalloc2(dim,&v0,dim*dim,&J);
284:   for (c = cStart; c < cEnd; ++c) {
285:     PetscCellGeometry geom;

287:     DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);
288:     geom.v0   = v0;
289:     geom.J    = J;
290:     geom.detJ = &detJ;
291:     for (f = 0, v = 0; f < numFields; ++f) {
292:       PetscFE      fe;
293:       void * const ctx = ctxs ? ctxs[f] : NULL;

295:       DMGetField(dm, f, (PetscObject *) &fe);
296:       PetscFEGetNumComponents(fe, &numComp);
297:       PetscDualSpaceGetDimension(sp[f], &spDim);
298:       for (d = 0; d < spDim; ++d) {
299:         if (funcs[f]) {
300:           PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);
301:         } else {
302:           for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0;
303:         }
304:         v += numComp;
305:       }
306:     }
307:     DMPlexVecSetClosure(dm, section, localX, c, values, mode);
308:   }
309:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
310:   PetscFree2(v0,J);
311:   PetscFree(sp);
312:   return(0);
313: }

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

320:   Input Parameters:
321: + dm      - The DM
322: . funcs   - The coordinate functions to evaluate, one per field
323: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
324: - mode    - The insertion mode for values

326:   Output Parameter:
327: . X - vector

329:   Level: developer

331: .seealso: DMPlexComputeL2Diff()
332: @*/
333: PetscErrorCode DMPlexProjectFunction(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
334: {
335:   Vec            localX;

340:   DMGetLocalVector(dm, &localX);
341:   DMPlexProjectFunctionLocal(dm, funcs, ctxs, mode, localX);
342:   DMLocalToGlobalBegin(dm, localX, mode, X);
343:   DMLocalToGlobalEnd(dm, localX, mode, X);
344:   DMRestoreLocalVector(dm, &localX);
345:   return(0);
346: }

350: PetscErrorCode DMPlexProjectFieldLocal(DM dm, Vec localU, void (**funcs)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal [], PetscScalar []), InsertMode mode, Vec localX)
351: {
352:   DM                dmAux;
353:   PetscDS           prob, probAux;
354:   Vec               A;
355:   PetscSection      section, sectionAux;
356:   PetscScalar      *values, *u, *u_x, *a, *a_x;
357:   PetscReal        *x, *v0, *J, *invJ, detJ, **basisField, **basisFieldDer, **basisFieldAux, **basisFieldDerAux;
358:   PetscInt          Nf, dim, spDim, totDim, numValues, cStart, cEnd, c, f, d, v, comp;
359:   PetscErrorCode    ierr;

362:   DMGetDS(dm, &prob);
363:   DMPlexGetDimension(dm, &dim);
364:   DMGetDefaultSection(dm, &section);
365:   PetscSectionGetNumFields(section, &Nf);
366:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
367:   PetscDSGetTotalDimension(prob, &totDim);
368:   PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
369:   PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
370:   PetscDSGetRefCoordArrays(prob, &x, NULL);
371:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
372:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
373:   if (dmAux) {
374:     DMGetDS(dmAux, &probAux);
375:     DMGetDefaultSection(dmAux, &sectionAux);
376:     PetscDSGetTabulation(prob, &basisFieldAux, &basisFieldDerAux);
377:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
378:   }
379:   DMPlexInsertBoundaryValuesFEM(dm, localU);
380:   DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);
381:   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
382:   DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);
383:   PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
384:   for (c = cStart; c < cEnd; ++c) {
385:     PetscScalar *coefficients = NULL, *coefficientsAux = NULL;

387:     DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
388:     DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);
389:     if (dmAux) {DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
390:     for (f = 0, v = 0; f < Nf; ++f) {
391:       PetscFE        fe;
392:       PetscDualSpace sp;
393:       PetscInt       Ncf;

395:       PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
396:       PetscFEGetDualSpace(fe, &sp);
397:       PetscFEGetNumComponents(fe, &Ncf);
398:       PetscDualSpaceGetDimension(sp, &spDim);
399:       for (d = 0; d < spDim; ++d) {
400:         PetscQuadrature  quad;
401:         const PetscReal *points, *weights;
402:         PetscInt         numPoints, q;

404:         if (funcs[f]) {
405:           PetscDualSpaceGetFunctional(sp, d, &quad);
406:           PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);
407:           PetscFEGetTabulation(fe, numPoints, points, &basisField[f], &basisFieldDer[f], NULL);
408:           for (q = 0; q < numPoints; ++q) {
409:             CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x);
410:             EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, coefficients,    NULL, u, u_x, NULL);
411:             EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);
412:             (*funcs[f])(u, NULL, u_x, a, NULL, a_x, x, &values[v]);
413:           }
414:           PetscFERestoreTabulation(fe, numPoints, points, &basisField[f], &basisFieldDer[f], NULL);
415:         } else {
416:           for (comp = 0; comp < Ncf; ++comp) values[v+comp] = 0.0;
417:         }
418:         v += Ncf;
419:       }
420:     }
421:     DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);
422:     if (dmAux) {DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);}
423:     DMPlexVecSetClosure(dm, section, localX, c, values, mode);
424:   }
425:   DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);
426:   PetscFree3(v0,J,invJ);
427:   return(0);
428: }

432: /*@C
433:   DMPlexProjectField - This projects the given function of the fields into the function space provided.

435:   Input Parameters:
436: + dm      - The DM
437: . U       - The input field vector
438: . funcs   - The functions to evaluate, one per field
439: - mode    - The insertion mode for values

441:   Output Parameter:
442: . X       - The output vector

444:   Level: developer

446: .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
447: @*/
448: PetscErrorCode DMPlexProjectField(DM dm, Vec U, void (**funcs)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal [], PetscScalar []), InsertMode mode, Vec X)
449: {
450:   Vec            localX, localU;

455:   DMGetLocalVector(dm, &localX);
456:   DMGetLocalVector(dm, &localU);
457:   DMGlobalToLocalBegin(dm, U, INSERT_VALUES, localU);
458:   DMGlobalToLocalEnd(dm, U, INSERT_VALUES, localU);
459:   DMPlexProjectFieldLocal(dm, localU, funcs, mode, localX);
460:   DMLocalToGlobalBegin(dm, localX, mode, X);
461:   DMLocalToGlobalEnd(dm, localX, mode, X);
462:   DMRestoreLocalVector(dm, &localX);
463:   DMRestoreLocalVector(dm, &localU);
464:   return(0);
465: }

469: PetscErrorCode DMPlexInsertBoundaryValuesFEM(DM dm, Vec localX)
470: {
471:   void        (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx);
472:   void         **ctxs;
473:   PetscFE       *fe;
474:   PetscInt       numFields, f, numBd, b;

480:   DMGetNumFields(dm, &numFields);
481:   PetscMalloc3(numFields,&fe,numFields,&funcs,numFields,&ctxs);
482:   for (f = 0; f < numFields; ++f) {DMGetField(dm, f, (PetscObject *) &fe[f]);}
483:   /* OPT: Could attempt to do multiple BCs at once */
484:   DMPlexGetNumBoundary(dm, &numBd);
485:   for (b = 0; b < numBd; ++b) {
486:     DMLabel         label;
487:     const PetscInt *ids;
488:     const char     *labelname;
489:     PetscInt        numids, field;
490:     PetscBool       isEssential;
491:     void          (*func)();
492:     void           *ctx;

494:     DMPlexGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, &func, &numids, &ids, &ctx);
495:     DMPlexGetLabel(dm, labelname, &label);
496:     for (f = 0; f < numFields; ++f) {
497:       funcs[f] = field == f ? (void (*)(const PetscReal[], PetscScalar *, void *)) func : NULL;
498:       ctxs[f]  = field == f ? ctx : NULL;
499:     }
500:     DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, fe, funcs, ctxs, INSERT_BC_VALUES, localX);
501:   }
502:   PetscFree3(fe,funcs,ctxs);
503:   return(0);
504: }

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

511:   Input Parameters:
512: + dm    - The DM
513: . funcs - The functions to evaluate for each field component
514: . ctxs  - Optional array of contexts to pass to each function, or NULL.
515: - X     - The coefficient vector u_h

517:   Output Parameter:
518: . diff - The diff ||u - u_h||_2

520:   Level: developer

522: .seealso: DMPlexProjectFunction(), DMPlexComputeL2GradientDiff()
523: @*/
524: PetscErrorCode DMPlexComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
525: {
526:   const PetscInt  debug = 0;
527:   PetscSection    section;
528:   PetscQuadrature quad;
529:   Vec             localX;
530:   PetscScalar    *funcVal;
531:   PetscReal      *coords, *v0, *J, *invJ, detJ;
532:   PetscReal       localDiff = 0.0;
533:   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
534:   PetscErrorCode  ierr;

537:   DMPlexGetDimension(dm, &dim);
538:   DMGetDefaultSection(dm, &section);
539:   PetscSectionGetNumFields(section, &numFields);
540:   DMGetLocalVector(dm, &localX);
541:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
542:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
543:   for (field = 0; field < numFields; ++field) {
544:     PetscFE  fe;
545:     PetscInt Nc;

547:     DMGetField(dm, field, (PetscObject *) &fe);
548:     PetscFEGetQuadrature(fe, &quad);
549:     PetscFEGetNumComponents(fe, &Nc);
550:     numComponents += Nc;
551:   }
552:   DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);
553:   PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
554:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
555:   for (c = cStart; c < cEnd; ++c) {
556:     PetscScalar *x = NULL;
557:     PetscReal    elemDiff = 0.0;

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

563:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
564:       PetscFE          fe;
565:       void * const     ctx = ctxs ? ctxs[field] : NULL;
566:       const PetscReal *quadPoints, *quadWeights;
567:       PetscReal       *basis;
568:       PetscInt         numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;

570:       DMGetField(dm, field, (PetscObject *) &fe);
571:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
572:       PetscFEGetDimension(fe, &numBasisFuncs);
573:       PetscFEGetNumComponents(fe, &numBasisComps);
574:       PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
575:       if (debug) {
576:         char title[1024];
577:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
578:         DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);
579:       }
580:       for (q = 0; q < numQuadPoints; ++q) {
581:         for (d = 0; d < dim; d++) {
582:           coords[d] = v0[d];
583:           for (e = 0; e < dim; e++) {
584:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
585:           }
586:         }
587:         (*funcs[field])(coords, funcVal, ctx);
588:         for (fc = 0; fc < numBasisComps; ++fc) {
589:           PetscScalar interpolant = 0.0;

591:           for (f = 0; f < numBasisFuncs; ++f) {
592:             const PetscInt fidx = f*numBasisComps+fc;
593:             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
594:           }
595:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
596:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
597:         }
598:       }
599:       comp        += numBasisComps;
600:       fieldOffset += numBasisFuncs*numBasisComps;
601:     }
602:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
603:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
604:     localDiff += elemDiff;
605:   }
606:   PetscFree5(funcVal,coords,v0,J,invJ);
607:   DMRestoreLocalVector(dm, &localX);
608:   MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
609:   *diff = PetscSqrtReal(*diff);
610:   return(0);
611: }

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

618:   Input Parameters:
619: + dm    - The DM
620: . funcs - The gradient functions to evaluate for each field component
621: . ctxs  - Optional array of contexts to pass to each function, or NULL.
622: . X     - The coefficient vector u_h
623: - n     - The vector to project along

625:   Output Parameter:
626: . diff - The diff ||(grad u - grad u_h) . n||_2

628:   Level: developer

630: .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
631: @*/
632: PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
633: {
634:   const PetscInt  debug = 0;
635:   PetscSection    section;
636:   PetscQuadrature quad;
637:   Vec             localX;
638:   PetscScalar    *funcVal, *interpolantVec;
639:   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
640:   PetscReal       localDiff = 0.0;
641:   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
642:   PetscErrorCode  ierr;

645:   DMPlexGetDimension(dm, &dim);
646:   DMGetDefaultSection(dm, &section);
647:   PetscSectionGetNumFields(section, &numFields);
648:   DMGetLocalVector(dm, &localX);
649:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
650:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
651:   for (field = 0; field < numFields; ++field) {
652:     PetscFE  fe;
653:     PetscInt Nc;

655:     DMGetField(dm, field, (PetscObject *) &fe);
656:     PetscFEGetQuadrature(fe, &quad);
657:     PetscFEGetNumComponents(fe, &Nc);
658:     numComponents += Nc;
659:   }
660:   /* DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX); */
661:   PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);
662:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
663:   for (c = cStart; c < cEnd; ++c) {
664:     PetscScalar *x = NULL;
665:     PetscReal    elemDiff = 0.0;

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

671:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
672:       PetscFE          fe;
673:       void * const     ctx = ctxs ? ctxs[field] : NULL;
674:       const PetscReal *quadPoints, *quadWeights;
675:       PetscReal       *basisDer;
676:       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;

678:       DMGetField(dm, field, (PetscObject *) &fe);
679:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
680:       PetscFEGetDimension(fe, &Nb);
681:       PetscFEGetNumComponents(fe, &Ncomp);
682:       PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
683:       if (debug) {
684:         char title[1024];
685:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
686:         DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);
687:       }
688:       for (q = 0; q < numQuadPoints; ++q) {
689:         for (d = 0; d < dim; d++) {
690:           coords[d] = v0[d];
691:           for (e = 0; e < dim; e++) {
692:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
693:           }
694:         }
695:         (*funcs[field])(coords, n, funcVal, ctx);
696:         for (fc = 0; fc < Ncomp; ++fc) {
697:           PetscScalar interpolant = 0.0;

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

703:             for (d = 0; d < dim; ++d) {
704:               realSpaceDer[d] = 0.0;
705:               for (g = 0; g < dim; ++g) {
706:                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
707:               }
708:               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
709:             }
710:           }
711:           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
712:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
713:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
714:         }
715:       }
716:       comp        += Ncomp;
717:       fieldOffset += Nb*Ncomp;
718:     }
719:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
720:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);}
721:     localDiff += elemDiff;
722:   }
723:   PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);
724:   DMRestoreLocalVector(dm, &localX);
725:   MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
726:   *diff = PetscSqrtReal(*diff);
727:   return(0);
728: }

732: PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
733: {
734:   const PetscInt  debug = 0;
735:   PetscSection    section;
736:   PetscQuadrature quad;
737:   Vec             localX;
738:   PetscScalar    *funcVal;
739:   PetscReal      *coords, *v0, *J, *invJ, detJ;
740:   PetscReal      *localDiff;
741:   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
742:   PetscErrorCode  ierr;

745:   DMPlexGetDimension(dm, &dim);
746:   DMGetDefaultSection(dm, &section);
747:   PetscSectionGetNumFields(section, &numFields);
748:   DMGetLocalVector(dm, &localX);
749:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
750:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
751:   for (field = 0; field < numFields; ++field) {
752:     PetscFE  fe;
753:     PetscInt Nc;

755:     DMGetField(dm, field, (PetscObject *) &fe);
756:     PetscFEGetQuadrature(fe, &quad);
757:     PetscFEGetNumComponents(fe, &Nc);
758:     numComponents += Nc;
759:   }
760:   DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);
761:   PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);
762:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
763:   for (c = cStart; c < cEnd; ++c) {
764:     PetscScalar *x = NULL;

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

770:     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
771:       PetscFE          fe;
772:       void * const     ctx = ctxs ? ctxs[field] : NULL;
773:       const PetscReal *quadPoints, *quadWeights;
774:       PetscReal       *basis, elemDiff = 0.0;
775:       PetscInt         numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;

777:       DMGetField(dm, field, (PetscObject *) &fe);
778:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);
779:       PetscFEGetDimension(fe, &numBasisFuncs);
780:       PetscFEGetNumComponents(fe, &numBasisComps);
781:       PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
782:       if (debug) {
783:         char title[1024];
784:         PetscSNPrintf(title, 1023, "Solution for Field %d", field);
785:         DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);
786:       }
787:       for (q = 0; q < numQuadPoints; ++q) {
788:         for (d = 0; d < dim; d++) {
789:           coords[d] = v0[d];
790:           for (e = 0; e < dim; e++) {
791:             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
792:           }
793:         }
794:         (*funcs[field])(coords, funcVal, ctx);
795:         for (fc = 0; fc < numBasisComps; ++fc) {
796:           PetscScalar interpolant = 0.0;

798:           for (f = 0; f < numBasisFuncs; ++f) {
799:             const PetscInt fidx = f*numBasisComps+fc;
800:             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
801:           }
802:           if (debug) {PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);}
803:           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
804:         }
805:       }
806:       comp        += numBasisComps;
807:       fieldOffset += numBasisFuncs*numBasisComps;
808:       localDiff[field] += elemDiff;
809:     }
810:     DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
811:   }
812:   DMRestoreLocalVector(dm, &localX);
813:   MPI_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
814:   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
815:   PetscFree6(localDiff,funcVal,coords,v0,J,invJ);
816:   return(0);
817: }

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

824:   Input Parameters:
825: + dm - The mesh
826: . X  - Local input vector
827: - user - The user context

829:   Output Parameter:
830: . integral - Local integral for each field

832:   Level: developer

834: .seealso: DMPlexComputeResidualFEM()
835: @*/
836: PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
837: {
838:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
839:   DM                dmAux;
840:   Vec               localX, A;
841:   PetscDS           prob, probAux = NULL;
842:   PetscQuadrature   q;
843:   PetscCellGeometry geom;
844:   PetscSection      section, sectionAux;
845:   PetscReal        *v0, *J, *invJ, *detJ;
846:   PetscScalar      *u, *a = NULL;
847:   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c;
848:   PetscInt          totDim, totDimAux;
849:   PetscErrorCode    ierr;

852:   /*PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);*/
853:   DMGetLocalVector(dm, &localX);
854:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
855:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
856:   DMPlexGetDimension(dm, &dim);
857:   DMGetDefaultSection(dm, &section);
858:   DMGetDS(dm, &prob);
859:   PetscDSGetTotalDimension(prob, &totDim);
860:   PetscSectionGetNumFields(section, &Nf);
861:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
862:   numCells = cEnd - cStart;
863:   for (f = 0; f < Nf; ++f) {integral[f]    = 0.0;}
864:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
865:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
866:   if (dmAux) {
867:     DMGetDefaultSection(dmAux, &sectionAux);
868:     DMGetDS(dmAux, &probAux);
869:     PetscDSGetTotalDimension(probAux, &totDimAux);
870:   }
871:   DMPlexInsertBoundaryValuesFEM(dm, localX);
872:   PetscMalloc5(numCells*totDim,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ);
873:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
874:   for (c = cStart; c < cEnd; ++c) {
875:     PetscScalar *x = NULL;
876:     PetscInt     i;

878:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
879:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
880:     DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);
881:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
882:     DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);
883:     if (dmAux) {
884:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
885:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
886:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
887:     }
888:   }
889:   for (f = 0; f < Nf; ++f) {
890:     PetscFE  fe;
891:     PetscInt numQuadPoints, Nb;
892:     /* Conforming batches */
893:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
894:     /* Remainder */
895:     PetscInt Nr, offset;

897:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
898:     PetscFEGetQuadrature(fe, &q);
899:     PetscFEGetDimension(fe, &Nb);
900:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
901:     PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
902:     blockSize = Nb*numQuadPoints;
903:     batchSize = numBlocks * blockSize;
904:      PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
905:     numChunks = numCells / (numBatches*batchSize);
906:     Ne        = numChunks*numBatches*batchSize;
907:     Nr        = numCells % (numBatches*batchSize);
908:     offset    = numCells - Nr;
909:     geom.v0   = v0;
910:     geom.J    = J;
911:     geom.invJ = invJ;
912:     geom.detJ = detJ;
913:     PetscFEIntegrate(fe, prob, f, Ne, geom, u, probAux, a, integral);
914:     geom.v0   = &v0[offset*dim];
915:     geom.J    = &J[offset*dim*dim];
916:     geom.invJ = &invJ[offset*dim*dim];
917:     geom.detJ = &detJ[offset];
918:     PetscFEIntegrate(fe, prob, f, Nr, geom, &u[offset*totDim], probAux, &a[offset*totDimAux], integral);
919:   }
920:   PetscFree5(u,v0,J,invJ,detJ);
921:   if (dmAux) {PetscFree(a);}
922:   if (mesh->printFEM) {
923:     PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");
924:     for (f = 0; f < Nf; ++f) {PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", integral[f]);}
925:     PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");
926:   }
927:   DMRestoreLocalVector(dm, &localX);
928:   /* TODO: Allreduce for integral */
929:   /*PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);*/
930:   return(0);
931: }

935: PetscErrorCode DMPlexComputeResidualFEM_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
936: {
937:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
938:   const char       *name  = "Residual";
939:   DM                dmAux;
940:   DMLabel           depth;
941:   Vec               A;
942:   PetscDS           prob, probAux = NULL;
943:   PetscQuadrature   q;
944:   PetscCellGeometry geom;
945:   PetscSection      section, sectionAux;
946:   PetscReal        *v0, *J, *invJ, *detJ;
947:   PetscScalar      *elemVec, *u, *u_t, *a = NULL;
948:   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c, numBd, bd;
949:   PetscInt          totDim, totDimBd, totDimAux;
950:   PetscErrorCode    ierr;

953:   PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
954:   DMPlexGetDimension(dm, &dim);
955:   DMGetDefaultSection(dm, &section);
956:   DMGetDS(dm, &prob);
957:   PetscDSGetTotalDimension(prob, &totDim);
958:   PetscDSGetTotalBdDimension(prob, &totDimBd);
959:   PetscSectionGetNumFields(section, &Nf);
960:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
961:   numCells = cEnd - cStart;
962:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
963:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
964:   if (dmAux) {
965:     DMGetDefaultSection(dmAux, &sectionAux);
966:     DMGetDS(dmAux, &probAux);
967:     PetscDSGetTotalDimension(probAux, &totDimAux);
968:   }
969:   DMPlexInsertBoundaryValuesFEM(dm, X);
970:   VecSet(F, 0.0);
971:   PetscMalloc7(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*totDim,&elemVec);
972:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
973:   for (c = cStart; c < cEnd; ++c) {
974:     PetscScalar *x = NULL, *x_t = NULL;
975:     PetscInt     i;

977:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
978:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
979:     DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
980:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
981:     DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
982:     if (X_t) {
983:       DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
984:       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
985:       DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
986:     }
987:     if (dmAux) {
988:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
989:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
990:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
991:     }
992:   }
993:   for (f = 0; f < Nf; ++f) {
994:     PetscFE  fe;
995:     PetscInt numQuadPoints, Nb;
996:     /* Conforming batches */
997:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
998:     /* Remainder */
999:     PetscInt Nr, offset;

1001:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1002:     PetscFEGetQuadrature(fe, &q);
1003:     PetscFEGetDimension(fe, &Nb);
1004:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1005:     PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1006:     blockSize = Nb*numQuadPoints;
1007:     batchSize = numBlocks * blockSize;
1008:      PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1009:     numChunks = numCells / (numBatches*batchSize);
1010:     Ne        = numChunks*numBatches*batchSize;
1011:     Nr        = numCells % (numBatches*batchSize);
1012:     offset    = numCells - Nr;
1013:     geom.v0   = v0;
1014:     geom.J    = J;
1015:     geom.invJ = invJ;
1016:     geom.detJ = detJ;
1017:     PetscFEIntegrateResidual(fe, prob, f, Ne, geom, u, u_t, probAux, a, elemVec);
1018:     geom.v0   = &v0[offset*dim];
1019:     geom.J    = &J[offset*dim*dim];
1020:     geom.invJ = &invJ[offset*dim*dim];
1021:     geom.detJ = &detJ[offset];
1022:     PetscFEIntegrateResidual(fe, prob, f, Nr, geom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);
1023:   }
1024:   for (c = cStart; c < cEnd; ++c) {
1025:     if (mesh->printFEM > 1) {DMPrintCellVector(c, name, totDim, &elemVec[c*totDim]);}
1026:     DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);
1027:   }
1028:   PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);
1029:   if (dmAux) {PetscFree(a);}
1030:   DMPlexGetDepthLabel(dm, &depth);
1031:   DMPlexGetNumBoundary(dm, &numBd);
1032:   for (bd = 0; bd < numBd; ++bd) {
1033:     const char     *bdLabel;
1034:     DMLabel         label;
1035:     IS              pointIS;
1036:     const PetscInt *points;
1037:     const PetscInt *values;
1038:     PetscReal      *n;
1039:     PetscInt        field, numValues, numPoints, p, dep, numFaces;
1040:     PetscBool       isEssential;

1042:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);
1043:     if (isEssential) continue;
1044:     if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
1045:     DMPlexGetLabel(dm, bdLabel, &label);
1046:     DMLabelGetStratumSize(label, 1, &numPoints);
1047:     DMLabelGetStratumIS(label, 1, &pointIS);
1048:     ISGetIndices(pointIS, &points);
1049:     for (p = 0, numFaces = 0; p < numPoints; ++p) {
1050:       DMLabelGetValue(depth, points[p], &dep);
1051:       if (dep == dim-1) ++numFaces;
1052:     }
1053:     PetscMalloc7(numFaces*totDimBd,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*totDimBd,&elemVec);
1054:     if (X_t) {PetscMalloc1(numFaces*totDimBd,&u_t);}
1055:     for (p = 0, f = 0; p < numPoints; ++p) {
1056:       const PetscInt point = points[p];
1057:       PetscScalar   *x     = NULL;
1058:       PetscInt       i;

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

1083:       PetscDSGetBdDiscretization(prob, f, (PetscObject *) &fe);
1084:       PetscFEGetQuadrature(fe, &q);
1085:       PetscFEGetDimension(fe, &Nb);
1086:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1087:       PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1088:       blockSize = Nb*numQuadPoints;
1089:       batchSize = numBlocks * blockSize;
1090:        PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1091:       numChunks = numFaces / (numBatches*batchSize);
1092:       Ne        = numChunks*numBatches*batchSize;
1093:       Nr        = numFaces % (numBatches*batchSize);
1094:       offset    = numFaces - Nr;
1095:       geom.v0   = v0;
1096:       geom.n    = n;
1097:       geom.J    = J;
1098:       geom.invJ = invJ;
1099:       geom.detJ = detJ;
1100:       PetscFEIntegrateBdResidual(fe, prob, f, Ne, geom, u, u_t, NULL, NULL, elemVec);
1101:       geom.v0   = &v0[offset*dim];
1102:       geom.n    = &n[offset*dim];
1103:       geom.J    = &J[offset*dim*dim];
1104:       geom.invJ = &invJ[offset*dim*dim];
1105:       geom.detJ = &detJ[offset];
1106:       PetscFEIntegrateBdResidual(fe, prob, f, Nr, geom, &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemVec[offset*totDimBd]);
1107:     }
1108:     for (p = 0, f = 0; p < numPoints; ++p) {
1109:       const PetscInt point = points[p];

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

1129: static PetscErrorCode DMPlexComputeResidualFEM_Check(DM dm, Vec X, Vec X_t, Vec F, void *user)
1130: {
1131:   DM                dmCh, dmAux;
1132:   Vec               A;
1133:   PetscDS           prob, probCh, probAux = NULL;
1134:   PetscQuadrature   q;
1135:   PetscCellGeometry geom;
1136:   PetscSection      section, sectionAux;
1137:   PetscReal        *v0, *J, *invJ, *detJ;
1138:   PetscScalar      *elemVec, *elemVecCh, *u, *u_t, *a = NULL;
1139:   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c;
1140:   PetscInt          totDim, totDimAux, diffCell = 0;
1141:   PetscErrorCode    ierr;

1144:   DMPlexGetDimension(dm, &dim);
1145:   DMGetDefaultSection(dm, &section);
1146:   DMGetDS(dm, &prob);
1147:   PetscDSGetTotalDimension(prob, &totDim);
1148:   PetscSectionGetNumFields(section, &Nf);
1149:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1150:   numCells = cEnd - cStart;
1151:   PetscObjectQuery((PetscObject) dm, "dmCh", (PetscObject *) &dmCh);
1152:   DMGetDS(dmCh, &probCh);
1153:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1154:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1155:   if (dmAux) {
1156:     DMGetDefaultSection(dmAux, &sectionAux);
1157:     DMGetDS(dmAux, &probAux);
1158:     PetscDSGetTotalDimension(probAux, &totDimAux);
1159:   }
1160:   DMPlexInsertBoundaryValuesFEM(dm, X);
1161:   VecSet(F, 0.0);
1162:   PetscMalloc7(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*totDim,&elemVec);
1163:   PetscMalloc1(numCells*totDim,&elemVecCh);
1164:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1165:   for (c = cStart; c < cEnd; ++c) {
1166:     PetscScalar *x = NULL, *x_t = NULL;
1167:     PetscInt     i;

1169:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
1170:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1171:     DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
1172:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1173:     DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
1174:     if (X_t) {
1175:       DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
1176:       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1177:       DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
1178:     }
1179:     if (dmAux) {
1180:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1181:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1182:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1183:     }
1184:   }
1185:   for (f = 0; f < Nf; ++f) {
1186:     PetscFE  fe, feCh;
1187:     PetscInt numQuadPoints, Nb;
1188:     /* Conforming batches */
1189:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1190:     /* Remainder */
1191:     PetscInt Nr, offset;

1193:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1194:     PetscDSGetDiscretization(probCh, f, (PetscObject *) &feCh);
1195:     PetscFEGetQuadrature(fe, &q);
1196:     PetscFEGetDimension(fe, &Nb);
1197:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1198:     PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1199:     blockSize = Nb*numQuadPoints;
1200:     batchSize = numBlocks * blockSize;
1201:      PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1202:     numChunks = numCells / (numBatches*batchSize);
1203:     Ne        = numChunks*numBatches*batchSize;
1204:     Nr        = numCells % (numBatches*batchSize);
1205:     offset    = numCells - Nr;
1206:     geom.v0   = v0;
1207:     geom.J    = J;
1208:     geom.invJ = invJ;
1209:     geom.detJ = detJ;
1210:     PetscFEIntegrateResidual(fe, prob, f, Ne, geom, u, u_t, probAux, a, elemVec);
1211:     PetscFEIntegrateResidual(feCh, prob, f, Ne, geom, u, u_t, probAux, a, elemVecCh);
1212:     geom.v0   = &v0[offset*dim];
1213:     geom.J    = &J[offset*dim*dim];
1214:     geom.invJ = &invJ[offset*dim*dim];
1215:     geom.detJ = &detJ[offset];
1216:     PetscFEIntegrateResidual(fe, prob, f, Nr, geom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);
1217:     PetscFEIntegrateResidual(feCh, prob, f, Nr, geom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVecCh[offset*totDim]);
1218:   }
1219:   for (c = cStart; c < cEnd; ++c) {
1220:     PetscBool diff = PETSC_FALSE;
1221:     PetscInt  d;

1223:     for (d = 0; d < totDim; ++d) if (PetscAbsScalar(elemVec[c*totDim+d] - elemVecCh[c*totDim+d]) > 1.0e-7) {diff = PETSC_TRUE;break;}
1224:     if (diff) {
1225:       PetscPrintf(PetscObjectComm((PetscObject) dm), "Different cell %d\n", c);
1226:       DMPrintCellVector(c, "Residual", totDim, &elemVec[c*totDim]);
1227:       DMPrintCellVector(c, "Check Residual", totDim, &elemVecCh[c*totDim]);
1228:       ++diffCell;
1229:     }
1230:     if (diffCell > 9) break;
1231:     DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);
1232:   }
1233:   PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);
1234:   PetscFree(elemVecCh);
1235:   if (dmAux) {PetscFree(a);}
1236:   return(0);
1237: }

1241: /*@
1242:   DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user

1244:   Input Parameters:
1245: + dm - The mesh
1246: . X  - Local solution
1247: - user - The user context

1249:   Output Parameter:
1250: . F  - Local output vector

1252:   Level: developer

1254: .seealso: DMPlexComputeJacobianActionFEM()
1255: @*/
1256: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1257: {
1258:   PetscObject    check;

1262:   PetscObjectQuery((PetscObject) dm, "dmCh", &check);
1263:   if (check) {DMPlexComputeResidualFEM_Check(dm, X, NULL, F, user);}
1264:   else       {DMPlexComputeResidualFEM_Internal(dm, X, NULL, F, user);}
1265:   return(0);
1266: }

1270: /*@
1271:   DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user

1273:   Input Parameters:
1274: + dm - The mesh
1275: . t - The time
1276: . X  - Local solution
1277: . X_t - Local solution time derivative, or NULL
1278: - user - The user context

1280:   Output Parameter:
1281: . F  - Local output vector

1283:   Level: developer

1285: .seealso: DMPlexComputeJacobianActionFEM()
1286: @*/
1287: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec X, Vec X_t, Vec F, void *user)
1288: {

1292:   DMPlexComputeResidualFEM_Internal(dm, X, X_t, F, user);
1293:   return(0);
1294: }

1298: PetscErrorCode DMPlexComputeJacobianFEM_Internal(DM dm, Vec X, Vec X_t, Mat Jac, Mat JacP,void *user)
1299: {
1300:   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1301:   const char       *name  = "Jacobian";
1302:   DM                dmAux;
1303:   DMLabel           depth;
1304:   Vec               A;
1305:   PetscDS           prob, probAux = NULL;
1306:   PetscQuadrature   quad;
1307:   PetscCellGeometry geom;
1308:   PetscSection      section, globalSection, sectionAux;
1309:   PetscReal        *v0, *J, *invJ, *detJ;
1310:   PetscScalar      *elemMat, *u, *u_t, *a = NULL;
1311:   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1312:   PetscInt          totDim, totDimBd, totDimAux, numBd, bd;
1313:   PetscBool         isShell;
1314:   PetscErrorCode    ierr;

1317:   PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
1318:   DMPlexGetDimension(dm, &dim);
1319:   DMGetDefaultSection(dm, &section);
1320:   DMGetDefaultGlobalSection(dm, &globalSection);
1321:   DMGetDS(dm, &prob);
1322:   PetscDSGetTotalDimension(prob, &totDim);
1323:   PetscDSGetTotalBdDimension(prob, &totDimBd);
1324:   PetscSectionGetNumFields(section, &Nf);
1325:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1326:   numCells = cEnd - cStart;
1327:   PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1328:   PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1329:   if (dmAux) {
1330:     DMGetDefaultSection(dmAux, &sectionAux);
1331:     DMGetDS(dmAux, &probAux);
1332:     PetscDSGetTotalDimension(probAux, &totDimAux);
1333:   }
1334:   DMPlexInsertBoundaryValuesFEM(dm, X);
1335:   MatZeroEntries(JacP);
1336:   PetscMalloc7(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*totDim*totDim,&elemMat);
1337:   if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1338:   for (c = cStart; c < cEnd; ++c) {
1339:     PetscScalar *x = NULL,  *x_t = NULL;
1340:     PetscInt     i;

1342:     DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
1343:     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1344:     DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
1345:     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1346:     DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
1347:     if (X_t) {
1348:       DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
1349:       for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1350:       DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
1351:     }
1352:     if (dmAux) {
1353:       DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1354:       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1355:       DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1356:     }
1357:   }
1358:   PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));
1359:   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1360:     PetscFE  fe;
1361:     PetscInt numQuadPoints, Nb;
1362:     /* Conforming batches */
1363:     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1364:     /* Remainder */
1365:     PetscInt Nr, offset;

1367:     PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1368:     PetscFEGetQuadrature(fe, &quad);
1369:     PetscFEGetDimension(fe, &Nb);
1370:     PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1371:     PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
1372:     blockSize = Nb*numQuadPoints;
1373:     batchSize = numBlocks * blockSize;
1374:     PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1375:     numChunks = numCells / (numBatches*batchSize);
1376:     Ne        = numChunks*numBatches*batchSize;
1377:     Nr        = numCells % (numBatches*batchSize);
1378:     offset    = numCells - Nr;
1379:     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1380:       geom.v0   = v0;
1381:       geom.J    = J;
1382:       geom.invJ = invJ;
1383:       geom.detJ = detJ;
1384:       PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Ne, geom, u, u_t, probAux, a, elemMat);
1385:       geom.v0   = &v0[offset*dim];
1386:       geom.J    = &J[offset*dim*dim];
1387:       geom.invJ = &invJ[offset*dim*dim];
1388:       geom.detJ = &detJ[offset];
1389:       PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Nr, geom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemMat[offset*totDim*totDim]);
1390:     }
1391:   }
1392:   for (c = cStart; c < cEnd; ++c) {
1393:     if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[c*totDim*totDim]);}
1394:     DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*totDim*totDim], ADD_VALUES);
1395:   }
1396:   PetscFree7(u,u_t,v0,J,invJ,detJ,elemMat);
1397:   if (dmAux) {PetscFree(a);}
1398:   DMPlexGetDepthLabel(dm, &depth);
1399:   DMPlexGetNumBoundary(dm, &numBd);
1400:   DMPlexGetDepthLabel(dm, &depth);
1401:   DMPlexGetNumBoundary(dm, &numBd);
1402:   for (bd = 0; bd < numBd; ++bd) {
1403:     const char     *bdLabel;
1404:     DMLabel         label;
1405:     IS              pointIS;
1406:     const PetscInt *points;
1407:     const PetscInt *values;
1408:     PetscReal      *n;
1409:     PetscInt        field, numValues, numPoints, p, dep, numFaces;
1410:     PetscBool       isEssential;

1412:     DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);
1413:     if (isEssential) continue;
1414:     if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
1415:     DMPlexGetLabel(dm, bdLabel, &label);
1416:     DMLabelGetStratumSize(label, 1, &numPoints);
1417:     DMLabelGetStratumIS(label, 1, &pointIS);
1418:     ISGetIndices(pointIS, &points);
1419:     for (p = 0, numFaces = 0; p < numPoints; ++p) {
1420:       DMLabelGetValue(depth, points[p], &dep);
1421:       if (dep == dim-1) ++numFaces;
1422:     }
1423:     PetscMalloc7(numFaces*totDimBd,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*totDimBd*totDimBd,&elemMat);
1424:     if (X_t) {PetscMalloc1(numFaces*totDimBd,&u_t);}
1425:     for (p = 0, f = 0; p < numPoints; ++p) {
1426:       const PetscInt point = points[p];
1427:       PetscScalar   *x     = NULL;
1428:       PetscInt       i;

1430:       DMLabelGetValue(depth, points[p], &dep);
1431:       if (dep != dim-1) continue;
1432:       DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);
1433:       DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
1434:       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
1435:       DMPlexVecGetClosure(dm, section, X, point, NULL, &x);
1436:       for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
1437:       DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);
1438:       if (X_t) {
1439:         DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);
1440:         for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
1441:         DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);
1442:       }
1443:       ++f;
1444:     }
1445:     PetscMemzero(elemMat, numFaces*totDimBd*totDimBd * sizeof(PetscScalar));
1446:     for (fieldI = 0; fieldI < Nf; ++fieldI) {
1447:       PetscFE  fe;
1448:       PetscInt numQuadPoints, Nb;
1449:       /* Conforming batches */
1450:       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1451:       /* Remainder */
1452:       PetscInt Nr, offset;

1454:       PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
1455:       PetscFEGetQuadrature(fe, &quad);
1456:       PetscFEGetDimension(fe, &Nb);
1457:       PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1458:       PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
1459:       blockSize = Nb*numQuadPoints;
1460:       batchSize = numBlocks * blockSize;
1461:        PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1462:       numChunks = numFaces / (numBatches*batchSize);
1463:       Ne        = numChunks*numBatches*batchSize;
1464:       Nr        = numFaces % (numBatches*batchSize);
1465:       offset    = numFaces - Nr;
1466:       for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1467:         geom.v0   = v0;
1468:         geom.n    = n;
1469:         geom.J    = J;
1470:         geom.invJ = invJ;
1471:         geom.detJ = detJ;
1472:         PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, geom, u, u_t, NULL, NULL, elemMat);
1473:         geom.v0   = &v0[offset*dim];
1474:         geom.n    = &n[offset*dim];
1475:         geom.J    = &J[offset*dim*dim];
1476:         geom.invJ = &invJ[offset*dim*dim];
1477:         geom.detJ = &detJ[offset];
1478:         PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Nr, geom, &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemMat[offset*totDimBd*totDimBd]);
1479:       }
1480:     }
1481:     for (p = 0, f = 0; p < numPoints; ++p) {
1482:       const PetscInt point = points[p];

1484:       DMLabelGetValue(depth, point, &dep);
1485:       if (dep != dim-1) continue;
1486:       if (mesh->printFEM > 1) {DMPrintCellMatrix(point, "BdJacobian", totDimBd, totDimBd, &elemMat[f*totDimBd*totDimBd]);}
1487:       DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*totDimBd*totDimBd], ADD_VALUES);
1488:       ++f;
1489:     }
1490:     ISRestoreIndices(pointIS, &points);
1491:     ISDestroy(&pointIS);
1492:     PetscFree7(u,v0,n,J,invJ,detJ,elemMat);
1493:     if (X_t) {PetscFree(u_t);}
1494:   }
1495:   MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
1496:   MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
1497:   if (mesh->printFEM) {
1498:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1499:     MatChop(JacP, 1.0e-10);
1500:     MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);
1501:   }
1502:   PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
1503:   PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);
1504:   if (isShell) {
1505:     JacActionCtx *jctx;

1507:     MatShellGetContext(Jac, &jctx);
1508:     VecCopy(X, jctx->u);
1509:   }
1510:   return(0);
1511: }

1515: /*@
1516:   DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.

1518:   Input Parameters:
1519: + dm - The mesh
1520: . X  - Local input vector
1521: - user - The user context

1523:   Output Parameter:
1524: . Jac  - Jacobian matrix

1526:   Note:
1527:   The first member of the user context must be an FEMContext.

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

1532:   Level: developer

1534: .seealso: FormFunctionLocal()
1535: @*/
1536: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1537: {

1541:   DMPlexComputeJacobianFEM_Internal(dm, X, NULL, Jac, JacP, user);
1542:   return(0);
1543: }

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

1550:   Input Parameters:
1551: + dmf  - The fine mesh
1552: . dmc  - The coarse mesh
1553: - user - The user context

1555:   Output Parameter:
1556: . In  - The interpolation matrix

1558:   Note:
1559:   The first member of the user context must be an FEMContext.

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

1564:   Level: developer

1566: .seealso: DMPlexComputeJacobianFEM()
1567: @*/
1568: PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1569: {
1570:   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1571:   const char       *name  = "Interpolator";
1572:   PetscDS           prob;
1573:   PetscFE          *feRef;
1574:   PetscSection      fsection, fglobalSection;
1575:   PetscSection      csection, cglobalSection;
1576:   PetscScalar      *elemMat;
1577:   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c;
1578:   PetscInt          cTotDim, rTotDim = 0;
1579:   PetscErrorCode    ierr;

1582:   PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1583:   DMPlexGetDimension(dmf, &dim);
1584:   DMGetDefaultSection(dmf, &fsection);
1585:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1586:   DMGetDefaultSection(dmc, &csection);
1587:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1588:   PetscSectionGetNumFields(fsection, &Nf);
1589:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1590:   DMGetDS(dmf, &prob);
1591:   PetscMalloc1(Nf,&feRef);
1592:   for (f = 0; f < Nf; ++f) {
1593:     PetscFE  fe;
1594:     PetscInt rNb, Nc;

1596:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1597:     PetscFERefine(fe, &feRef[f]);
1598:     PetscFEGetDimension(feRef[f], &rNb);
1599:     PetscFEGetNumComponents(fe, &Nc);
1600:     rTotDim += rNb*Nc;
1601:   }
1602:   PetscDSGetTotalDimension(prob, &cTotDim);
1603:   PetscMalloc1(rTotDim*cTotDim,&elemMat);
1604:   PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));
1605:   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1606:     PetscDualSpace   Qref;
1607:     PetscQuadrature  f;
1608:     const PetscReal *qpoints, *qweights;
1609:     PetscReal       *points;
1610:     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;

1612:     /* Compose points from all dual basis functionals */
1613:     PetscFEGetDualSpace(feRef[fieldI], &Qref);
1614:     PetscFEGetNumComponents(feRef[fieldI], &Nc);
1615:     PetscDualSpaceGetDimension(Qref, &fpdim);
1616:     for (i = 0; i < fpdim; ++i) {
1617:       PetscDualSpaceGetFunctional(Qref, i, &f);
1618:       PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);
1619:       npoints += Np;
1620:     }
1621:     PetscMalloc1(npoints*dim,&points);
1622:     for (i = 0, k = 0; i < fpdim; ++i) {
1623:       PetscDualSpaceGetFunctional(Qref, i, &f);
1624:       PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);
1625:       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1626:     }

1628:     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1629:       PetscFE    fe;
1630:       PetscReal *B;
1631:       PetscInt   NcJ, cpdim, j;

1633:       /* Evaluate basis at points */
1634:       PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
1635:       PetscFEGetNumComponents(fe, &NcJ);
1636:       PetscFEGetDimension(fe, &cpdim);
1637:       /* For now, fields only interpolate themselves */
1638:       if (fieldI == fieldJ) {
1639:         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);
1640:         PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);
1641:         for (i = 0, k = 0; i < fpdim; ++i) {
1642:           PetscDualSpaceGetFunctional(Qref, i, &f);
1643:           PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);
1644:           for (p = 0; p < Np; ++p, ++k) {
1645:             for (j = 0; j < cpdim; ++j) {
1646:               for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p];
1647:             }
1648:           }
1649:         }
1650:         PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);
1651:       }
1652:       offsetJ += cpdim*NcJ;
1653:     }
1654:     offsetI += fpdim*Nc;
1655:     PetscFree(points);
1656:   }
1657:   if (mesh->printFEM > 1) {DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);}
1658:   /* Preallocate matrix */
1659:   {
1660:     PetscHashJK ht;
1661:     PetscLayout rLayout;
1662:     PetscInt   *dnz, *onz, *cellCIndices, *cellFIndices;
1663:     PetscInt    locRows, rStart, rEnd, cell, r;

1665:     MatGetLocalSize(In, &locRows, NULL);
1666:     PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);
1667:     PetscLayoutSetLocalSize(rLayout, locRows);
1668:     PetscLayoutSetBlockSize(rLayout, 1);
1669:     PetscLayoutSetUp(rLayout);
1670:     PetscLayoutGetRange(rLayout, &rStart, &rEnd);
1671:     PetscLayoutDestroy(&rLayout);
1672:     PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);
1673:     PetscHashJKCreate(&ht);
1674:     for (cell = cStart; cell < cEnd; ++cell) {
1675:       DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);
1676:       for (r = 0; r < rTotDim; ++r) {
1677:         PetscHashJKKey  key;
1678:         PetscHashJKIter missing, iter;

1680:         key.j = cellFIndices[r];
1681:         if (key.j < 0) continue;
1682:         for (c = 0; c < cTotDim; ++c) {
1683:           key.k = cellCIndices[c];
1684:           if (key.k < 0) continue;
1685:           PetscHashJKPut(ht, key, &missing, &iter);
1686:           if (missing) {
1687:             PetscHashJKSet(ht, iter, 1);
1688:             if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1689:             else                                     ++onz[key.j-rStart];
1690:           }
1691:         }
1692:       }
1693:     }
1694:     PetscHashJKDestroy(&ht);
1695:     MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);
1696:     MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
1697:     PetscFree4(dnz,onz,cellCIndices,cellFIndices);
1698:   }
1699:   /* Fill matrix */
1700:   MatZeroEntries(In);
1701:   for (c = cStart; c < cEnd; ++c) {
1702:     DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);
1703:   }
1704:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1705:   PetscFree(feRef);
1706:   PetscFree(elemMat);
1707:   MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);
1708:   MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);
1709:   if (mesh->printFEM) {
1710:     PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
1711:     MatChop(In, 1.0e-10);
1712:     MatView(In, PETSC_VIEWER_STDOUT_WORLD);
1713:   }
1714:   PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);
1715:   return(0);
1716: }

1720: PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1721: {
1722:   PetscDS        prob;
1723:   PetscFE       *feRef;
1724:   Vec            fv, cv;
1725:   IS             fis, cis;
1726:   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1727:   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1728:   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, c, dim, d, startC, offsetC, offsetF, m;

1732:   PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1733:   DMPlexGetDimension(dmf, &dim);
1734:   DMGetDefaultSection(dmf, &fsection);
1735:   DMGetDefaultGlobalSection(dmf, &fglobalSection);
1736:   DMGetDefaultSection(dmc, &csection);
1737:   DMGetDefaultGlobalSection(dmc, &cglobalSection);
1738:   PetscSectionGetNumFields(fsection, &Nf);
1739:   DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);
1740:   DMGetDS(dmc, &prob);
1741:   PetscMalloc1(Nf,&feRef);
1742:   for (f = 0; f < Nf; ++f) {
1743:     PetscFE  fe;
1744:     PetscInt fNb, Nc;

1746:     PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1747:     PetscFERefine(fe, &feRef[f]);
1748:     PetscFEGetDimension(feRef[f], &fNb);
1749:     PetscFEGetNumComponents(fe, &Nc);
1750:     fTotDim += fNb*Nc;
1751:   }
1752:   PetscDSGetTotalDimension(prob, &cTotDim);
1753:   PetscMalloc1(cTotDim,&cmap);
1754:   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1755:     PetscFE        feC;
1756:     PetscDualSpace QF, QC;
1757:     PetscInt       NcF, NcC, fpdim, cpdim;

1759:     PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);
1760:     PetscFEGetNumComponents(feC, &NcC);
1761:     PetscFEGetNumComponents(feRef[field], &NcF);
1762:     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);
1763:     PetscFEGetDualSpace(feRef[field], &QF);
1764:     PetscDualSpaceGetDimension(QF, &fpdim);
1765:     PetscFEGetDualSpace(feC, &QC);
1766:     PetscDualSpaceGetDimension(QC, &cpdim);
1767:     for (c = 0; c < cpdim; ++c) {
1768:       PetscQuadrature  cfunc;
1769:       const PetscReal *cqpoints;
1770:       PetscInt         NpC;

1772:       PetscDualSpaceGetFunctional(QC, c, &cfunc);
1773:       PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);
1774:       if (NpC != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1775:       for (f = 0; f < fpdim; ++f) {
1776:         PetscQuadrature  ffunc;
1777:         const PetscReal *fqpoints;
1778:         PetscReal        sum = 0.0;
1779:         PetscInt         NpF, comp;

1781:         PetscDualSpaceGetFunctional(QF, f, &ffunc);
1782:         PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);
1783:         if (NpC != NpF) continue;
1784:         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1785:         if (sum > 1.0e-9) continue;
1786:         for (comp = 0; comp < NcC; ++comp) {
1787:           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1788:         }
1789:         break;
1790:       }
1791:     }
1792:     offsetC += cpdim*NcC;
1793:     offsetF += fpdim*NcF;
1794:   }
1795:   for (f = 0; f < Nf; ++f) {PetscFEDestroy(&feRef[f]);}
1796:   PetscFree(feRef);

1798:   DMGetGlobalVector(dmf, &fv);
1799:   DMGetGlobalVector(dmc, &cv);
1800:   VecGetOwnershipRange(cv, &startC, NULL);
1801:   PetscSectionGetConstrainedStorageSize(cglobalSection, &m);
1802:   PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);
1803:   PetscMalloc1(m,&cindices);
1804:   PetscMalloc1(m,&findices);
1805:   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1806:   for (c = cStart; c < cEnd; ++c) {
1807:     DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);
1808:     for (d = 0; d < cTotDim; ++d) {
1809:       if (cellCIndices[d] < 0) continue;
1810:       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]]);
1811:       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1812:       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1813:     }
1814:   }
1815:   PetscFree(cmap);
1816:   PetscFree2(cellCIndices,cellFIndices);

1818:   ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);
1819:   ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);
1820:   VecScatterCreate(cv, cis, fv, fis, sc);
1821:   ISDestroy(&cis);
1822:   ISDestroy(&fis);
1823:   DMRestoreGlobalVector(dmf, &fv);
1824:   DMRestoreGlobalVector(dmc, &cv);
1825:   PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);
1826:   return(0);
1827: }

1831: static PetscErrorCode BoundaryDuplicate(DMBoundary bd, DMBoundary *boundary)
1832: {
1833:   DMBoundary     b = bd, b2, bold = NULL;

1837:   *boundary = NULL;
1838:   for (; b; b = b->next, bold = b2) {
1839:     PetscNew(&b2);
1840:     PetscStrallocpy(b->name, (char **) &b2->name);
1841:     PetscStrallocpy(b->labelname, (char **) &b2->labelname);
1842:     PetscMalloc1(b->numids, &b2->ids);
1843:     PetscMemcpy(b2->ids, b->ids, b->numids*sizeof(PetscInt));
1844:     b2->label     = NULL;
1845:     b2->essential = b->essential;
1846:     b2->field     = b->field;
1847:     b2->func      = b->func;
1848:     b2->numids    = b->numids;
1849:     b2->ctx       = b->ctx;
1850:     b2->next      = NULL;
1851:     if (!*boundary) *boundary   = b2;
1852:     if (bold)        bold->next = b2;
1853:   }
1854:   return(0);
1855: }

1859: PetscErrorCode DMPlexCopyBoundary(DM dm, DM dmNew)
1860: {
1861:   DM_Plex       *mesh    = (DM_Plex *) dm->data;
1862:   DM_Plex       *meshNew = (DM_Plex *) dmNew->data;
1863:   DMBoundary     b;

1867:   BoundaryDuplicate(mesh->boundary, &meshNew->boundary);
1868:   for (b = meshNew->boundary; b; b = b->next) {
1869:     if (b->labelname) {
1870:       DMPlexGetLabel(dmNew, b->labelname, &b->label);
1871:       if (!b->label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s does not exist in this DM", b->labelname);
1872:     }
1873:   }
1874:   return(0);
1875: }

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

1888:   PetscNew(&b);
1889:   PetscStrallocpy(name, (char **) &b->name);
1890:   PetscStrallocpy(labelname, (char **) &b->labelname);
1891:   PetscMalloc1(numids, &b->ids);
1892:   PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));
1893:   if (b->labelname) {
1894:     DMPlexGetLabel(dm, b->labelname, &b->label);
1895:     if (!b->label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s does not exist in this DM", b->labelname);
1896:   }
1897:   b->essential   = isEssential;
1898:   b->field       = field;
1899:   b->func        = bcFunc;
1900:   b->numids      = numids;
1901:   b->ctx         = ctx;
1902:   b->next        = mesh->boundary;
1903:   mesh->boundary = b;
1904:   return(0);
1905: }

1909: PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd)
1910: {
1911:   DM_Plex   *mesh = (DM_Plex *) dm->data;
1912:   DMBoundary b    = mesh->boundary;

1917:   *numBd = 0;
1918:   while (b) {++(*numBd); b = b->next;}
1919:   return(0);
1920: }

1924: PetscErrorCode DMPlexGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, const char **labelname, PetscInt *field, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx)
1925: {
1926:   DM_Plex   *mesh = (DM_Plex *) dm->data;
1927:   DMBoundary b    = mesh->boundary;
1928:   PetscInt   n    = 0;

1932:   while (b) {
1933:     if (n == bd) break;
1934:     b = b->next;
1935:     ++n;
1936:   }
1937:   if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
1938:   if (isEssential) {
1940:     *isEssential = b->essential;
1941:   }
1942:   if (name) {
1944:     *name = b->name;
1945:   }
1946:   if (labelname) {
1948:     *labelname = b->labelname;
1949:   }
1950:   if (field) {
1952:     *field = b->field;
1953:   }
1954:   if (func) {
1956:     *func = b->func;
1957:   }
1958:   if (numids) {
1960:     *numids = b->numids;
1961:   }
1962:   if (ids) {
1964:     *ids = b->ids;
1965:   }
1966:   if (ctx) {
1968:     *ctx = b->ctx;
1969:   }
1970:   return(0);
1971: }

1975: PetscErrorCode DMPlexIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
1976: {
1977:   DM_Plex       *mesh = (DM_Plex *) dm->data;
1978:   DMBoundary     b    = mesh->boundary;

1984:   *isBd = PETSC_FALSE;
1985:   while (b && !(*isBd)) {
1986:     if (b->label) {
1987:       PetscInt i;

1989:       for (i = 0; i < b->numids && !(*isBd); ++i) {
1990:         DMLabelStratumHasPoint(b->label, b->ids[i], point, isBd);
1991:       }
1992:     }
1993:     b = b->next;
1994:   }
1995:   return(0);
1996: }