xref: /petsc/src/dm/impls/plex/plexfem.c (revision c0cd03017aa709fd99d1ab926cbd6798e2380b6f) 
1 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
3 #include <petscfe.h>
4 
5 #undef __FUNCT__
6 #define __FUNCT__ "DMPlexGetScale"
7 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
8 {
9   DM_Plex *mesh = (DM_Plex*) dm->data;
10 
11   PetscFunctionBegin;
12   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
13   PetscValidPointer(scale, 3);
14   *scale = mesh->scale[unit];
15   PetscFunctionReturn(0);
16 }
17 
18 #undef __FUNCT__
19 #define __FUNCT__ "DMPlexSetScale"
20 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
21 {
22   DM_Plex *mesh = (DM_Plex*) dm->data;
23 
24   PetscFunctionBegin;
25   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
26   mesh->scale[unit] = scale;
27   PetscFunctionReturn(0);
28 }
29 
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 }
84 
85 #undef __FUNCT__
86 #define __FUNCT__ "DMPlexCreateRigidBody"
87 /*@C
88   DMPlexCreateRigidBody - create rigid body modes from coordinates
89 
90   Collective on DM
91 
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
96 
97   Output Argument:
98 . sp - the null space
99 
100   Note: This is necessary to take account of Dirichlet conditions on the displacements
101 
102   Level: advanced
103 
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;
113   PetscErrorCode ierr;
114 
115   PetscFunctionBegin;
116   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
117   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
118   if (dim == 1) {
119     ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr);
120     PetscFunctionReturn(0);
121   }
122   if (!section)       {ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);}
123   if (!globalSection) {ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);}
124   ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr);
125   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
126   ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
127   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
128   m    = (dim*(dim+1))/2;
129   ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr);
130   ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr);
131   ierr = VecSetUp(mode[0]);CHKERRQ(ierr);
132   for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);}
133   /* Assume P1 */
134   ierr = DMGetLocalVector(dm, &localMode);CHKERRQ(ierr);
135   for (d = 0; d < dim; ++d) {
136     PetscScalar values[3] = {0.0, 0.0, 0.0};
137 
138     values[d] = 1.0;
139     ierr      = VecSet(localMode, 0.0);CHKERRQ(ierr);
140     for (v = vStart; v < vEnd; ++v) {
141       ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
142     }
143     ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
144     ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
145   }
146   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
147   for (d = dim; d < dim*(dim+1)/2; ++d) {
148     PetscInt i, j, k = dim > 2 ? d - dim : d;
149 
150     ierr = VecSet(localMode, 0.0);CHKERRQ(ierr);
151     for (v = vStart; v < vEnd; ++v) {
152       PetscScalar values[3] = {0.0, 0.0, 0.0};
153       PetscInt    off;
154 
155       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
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       ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
162     }
163     ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
164     ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
165   }
166   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
167   ierr = DMRestoreLocalVector(dm, &localMode);CHKERRQ(ierr);
168   for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);}
169   /* Orthonormalize system */
170   for (i = dim; i < m; ++i) {
171     PetscScalar dots[6];
172 
173     ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr);
174     for (j = 0; j < i; ++j) dots[j] *= -1.0;
175     ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr);
176     ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);
177   }
178   ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr);
179   for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);}
180   PetscFunctionReturn(0);
181 }
182 
183 #undef __FUNCT__
184 #define __FUNCT__ "DMPlexProjectFunctionLocal"
185 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, 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, c, f, d, v, comp;
192   PetscErrorCode  ierr;
193 
194   PetscFunctionBegin;
195   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
196   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
197   ierr = PetscMalloc1(numFields, &sp);CHKERRQ(ierr);
198   for (f = 0; f < numFields; ++f) {
199     ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr);
200     ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr);
201     ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
202     totDim += spDim*numComp;
203   }
204   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
205   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
206   ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr);
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   ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
209   ierr = PetscMalloc2(dim,&v0,dim*dim,&J);CHKERRQ(ierr);
210   for (c = cStart; c < cEnd; ++c) {
211     PetscCellGeometry geom;
212 
213     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr);
214     geom.v0   = v0;
215     geom.J    = J;
216     geom.detJ = &detJ;
217     for (f = 0, v = 0; f < numFields; ++f) {
218       void * const ctx = ctxs ? ctxs[f] : NULL;
219       ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr);
220       ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
221       for (d = 0; d < spDim; ++d) {
222         if (funcs[f]) {
223           ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr);
224         } else {
225           for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0;
226         }
227         v += numComp;
228       }
229     }
230     ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr);
231   }
232   ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
233   ierr = PetscFree2(v0,J);CHKERRQ(ierr);
234   ierr = PetscFree(sp);CHKERRQ(ierr);
235   PetscFunctionReturn(0);
236 }
237 
238 #undef __FUNCT__
239 #define __FUNCT__ "DMPlexProjectFunction"
240 /*@C
241   DMPlexProjectFunction - This projects the given function into the function space provided.
242 
243   Input Parameters:
244 + dm      - The DM
245 . fe      - The PetscFE associated with the field
246 . funcs   - The coordinate functions to evaluate, one per field
247 . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
248 - mode    - The insertion mode for values
249 
250   Output Parameter:
251 . X - vector
252 
253   Level: developer
254 
255 .seealso: DMPlexComputeL2Diff()
256 @*/
257 PetscErrorCode DMPlexProjectFunction(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
258 {
259   Vec            localX;
260   PetscErrorCode ierr;
261 
262   PetscFunctionBegin;
263   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
264   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
265   ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, mode, localX);CHKERRQ(ierr);
266   ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr);
267   ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr);
268   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
269   PetscFunctionReturn(0);
270 }
271 
272 
273 #undef __FUNCT__
274 #define __FUNCT__ "DMPlexInsertBoundaryValues"
275 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec localX)
276 {
277   void        (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx);
278   void         **ctxs;
279   PetscFE       *fe;
280   PetscInt       numFields, f, numBd, b;
281   PetscErrorCode ierr;
282 
283   PetscFunctionBegin;
284   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
285   PetscValidHeaderSpecific(localX, VEC_CLASSID, 2);
286   ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr);
287   ierr = PetscMalloc3(numFields,&fe,numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr);
288   for (f = 0; f < numFields; ++f) {ierr = DMGetField(dm, f, (PetscObject *) &fe[f]);CHKERRQ(ierr);}
289   /* OPT: Could attempt to do multiple BCs at once */
290   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
291   for (b = 0; b < numBd; ++b) {
292     const PetscInt *ids;
293     PetscInt        numids, field;
294     PetscBool       isEssential;
295     void          (*func)();
296     void           *ctx;
297 
298     /* TODO: We need to set only the part indicated by the ids */
299     ierr = DMPlexGetBoundary(dm, b, &isEssential, NULL, &field, &func, &numids, &ids, &ctx);CHKERRQ(ierr);
300     for (f = 0; f < numFields; ++f) {
301       funcs[f] = field == f ? (void (*)(const PetscReal[], PetscScalar *, void *)) func : NULL;
302       ctxs[f]  = field == f ? ctx : NULL;
303     }
304     ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
305   }
306   ierr = PetscFree3(fe,funcs,ctxs);CHKERRQ(ierr);
307   PetscFunctionReturn(0);
308 }
309 
310 #undef __FUNCT__
311 #define __FUNCT__ "DMPlexComputeL2Diff"
312 /*@C
313   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
314 
315   Input Parameters:
316 + dm    - The DM
317 . fe    - The PetscFE object for each field
318 . funcs - The functions to evaluate for each field component
319 . ctxs  - Optional array of contexts to pass to each function, or NULL.
320 - X     - The coefficient vector u_h
321 
322   Output Parameter:
323 . diff - The diff ||u - u_h||_2
324 
325   Level: developer
326 
327 .seealso: DMPlexProjectFunction(), DMPlexComputeL2GradientDiff()
328 @*/
329 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
330 {
331   const PetscInt  debug = 0;
332   PetscSection    section;
333   PetscQuadrature quad;
334   Vec             localX;
335   PetscScalar    *funcVal;
336   PetscReal      *coords, *v0, *J, *invJ, detJ;
337   PetscReal       localDiff = 0.0;
338   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
339   PetscErrorCode  ierr;
340 
341   PetscFunctionBegin;
342   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
343   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
344   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
345   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
346   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
347   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
348   for (field = 0; field < numFields; ++field) {
349     PetscInt Nc;
350 
351     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
352     numComponents += Nc;
353   }
354   ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
355   ierr = PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
356   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
357   ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr);
358   for (c = cStart; c < cEnd; ++c) {
359     PetscScalar *x = NULL;
360     PetscReal    elemDiff = 0.0;
361 
362     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
363     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
364     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
365 
366     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
367       void * const     ctx = ctxs ? ctxs[field] : NULL;
368       const PetscReal *quadPoints, *quadWeights;
369       PetscReal       *basis;
370       PetscInt         numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;
371 
372       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
373       ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr);
374       ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr);
375       ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr);
376       if (debug) {
377         char title[1024];
378         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
379         ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr);
380       }
381       for (q = 0; q < numQuadPoints; ++q) {
382         for (d = 0; d < dim; d++) {
383           coords[d] = v0[d];
384           for (e = 0; e < dim; e++) {
385             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
386           }
387         }
388         (*funcs[field])(coords, funcVal, ctx);
389         for (fc = 0; fc < numBasisComps; ++fc) {
390           PetscScalar interpolant = 0.0;
391 
392           for (f = 0; f < numBasisFuncs; ++f) {
393             const PetscInt fidx = f*numBasisComps+fc;
394             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
395           }
396           if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);}
397           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
398         }
399       }
400       comp        += numBasisComps;
401       fieldOffset += numBasisFuncs*numBasisComps;
402     }
403     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
404     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
405     localDiff += elemDiff;
406   }
407   ierr  = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr);
408   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
409   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
410   *diff = PetscSqrtReal(*diff);
411   PetscFunctionReturn(0);
412 }
413 
414 #undef __FUNCT__
415 #define __FUNCT__ "DMPlexComputeL2GradientDiff"
416 /*@C
417   DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
418 
419   Input Parameters:
420 + dm    - The DM
421 . fe    - The PetscFE object for each field
422 . funcs - The gradient functions to evaluate for each field component
423 . ctxs  - Optional array of contexts to pass to each function, or NULL.
424 . X     - The coefficient vector u_h
425 - n     - The vector to project along
426 
427   Output Parameter:
428 . diff - The diff ||(grad u - grad u_h) . n||_2
429 
430   Level: developer
431 
432 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
433 @*/
434 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
435 {
436   const PetscInt  debug = 0;
437   PetscSection    section;
438   PetscQuadrature quad;
439   Vec             localX;
440   PetscScalar    *funcVal, *interpolantVec;
441   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
442   PetscReal       localDiff = 0.0;
443   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
444   PetscErrorCode  ierr;
445 
446   PetscFunctionBegin;
447   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
448   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
449   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
450   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
451   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
452   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
453   for (field = 0; field < numFields; ++field) {
454     PetscInt Nc;
455 
456     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
457     numComponents += Nc;
458   }
459   /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */
460   ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr);
461   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
462   ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr);
463   for (c = cStart; c < cEnd; ++c) {
464     PetscScalar *x = NULL;
465     PetscReal    elemDiff = 0.0;
466 
467     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
468     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
469     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
470 
471     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
472       void * const     ctx = ctxs ? ctxs[field] : NULL;
473       const PetscReal *quadPoints, *quadWeights;
474       PetscReal       *basisDer;
475       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;
476 
477       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
478       ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
479       ierr = PetscFEGetNumComponents(fe[field], &Ncomp);CHKERRQ(ierr);
480       ierr = PetscFEGetDefaultTabulation(fe[field], NULL, &basisDer, NULL);CHKERRQ(ierr);
481       if (debug) {
482         char title[1024];
483         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
484         ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr);
485       }
486       for (q = 0; q < numQuadPoints; ++q) {
487         for (d = 0; d < dim; d++) {
488           coords[d] = v0[d];
489           for (e = 0; e < dim; e++) {
490             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
491           }
492         }
493         (*funcs[field])(coords, n, funcVal, ctx);
494         for (fc = 0; fc < Ncomp; ++fc) {
495           PetscScalar interpolant = 0.0;
496 
497           for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
498           for (f = 0; f < Nb; ++f) {
499             const PetscInt fidx = f*Ncomp+fc;
500 
501             for (d = 0; d < dim; ++d) {
502               realSpaceDer[d] = 0.0;
503               for (g = 0; g < dim; ++g) {
504                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
505               }
506               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
507             }
508           }
509           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
510           if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "    elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);}
511           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
512         }
513       }
514       comp        += Ncomp;
515       fieldOffset += Nb*Ncomp;
516     }
517     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
518     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
519     localDiff += elemDiff;
520   }
521   ierr  = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr);
522   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
523   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
524   *diff = PetscSqrtReal(*diff);
525   PetscFunctionReturn(0);
526 }
527 
528 #undef __FUNCT__
529 #define __FUNCT__ "DMPlexComputeResidualFEM"
530 /*@
531   DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
532 
533   Input Parameters:
534 + dm - The mesh
535 . X  - Local input vector
536 - user - The user context
537 
538   Output Parameter:
539 . F  - Local output vector
540 
541   Note:
542   The first member of the user context must be an FEMContext.
543 
544   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
545   like a GPU, or vectorize on a multicore machine.
546 
547   Level: developer
548 
549 .seealso: DMPlexComputeJacobianActionFEM()
550 @*/
551 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
552 {
553   DM_Plex          *mesh  = (DM_Plex *) dm->data;
554   PetscFEM         *fem   = (PetscFEM *) user;
555   PetscFE          *fe    = fem->fe;
556   PetscFE          *feAux = fem->feAux;
557   PetscFE          *feBd  = fem->feBd;
558   const char       *name  = "Residual";
559   DM                dmAux;
560   Vec               A;
561   PetscQuadrature   q;
562   PetscCellGeometry geom;
563   PetscSection      section, sectionAux;
564   PetscReal        *v0, *J, *invJ, *detJ;
565   PetscScalar      *elemVec, *u, *a = NULL;
566   PetscInt          dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c;
567   PetscInt          cellDof = 0, numComponents = 0;
568   PetscInt          cellDofAux = 0, numComponentsAux = 0;
569   PetscErrorCode    ierr;
570 
571   PetscFunctionBegin;
572   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
573   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
574   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
575   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
576   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
577   numCells = cEnd - cStart;
578   for (f = 0; f < Nf; ++f) {
579     PetscInt Nb, Nc;
580 
581     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
582     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
583     cellDof       += Nb*Nc;
584     numComponents += Nc;
585   }
586   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
587   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
588   if (dmAux) {
589     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
590     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
591   }
592   for (f = 0; f < NfAux; ++f) {
593     PetscInt Nb, Nc;
594 
595     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
596     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
597     cellDofAux       += Nb*Nc;
598     numComponentsAux += Nc;
599   }
600   ierr = DMPlexProjectFunctionLocal(dm, fe, fem->bcFuncs, fem->bcCtxs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
601   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
602   ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr);
603   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
604   for (c = cStart; c < cEnd; ++c) {
605     PetscScalar *x = NULL;
606     PetscInt     i;
607 
608     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
609     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
610     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
611     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
612     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
613     if (dmAux) {
614       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
615       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
616       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
617     }
618   }
619   for (f = 0; f < Nf; ++f) {
620     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[f];
621     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[f];
622     PetscInt numQuadPoints, Nb;
623     /* Conforming batches */
624     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
625     /* Remainder */
626     PetscInt Nr, offset;
627 
628     ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr);
629     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
630     ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
631     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
632     blockSize = Nb*numQuadPoints;
633     batchSize = numBlocks * blockSize;
634     ierr =  PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
635     numChunks = numCells / (numBatches*batchSize);
636     Ne        = numChunks*numBatches*batchSize;
637     Nr        = numCells % (numBatches*batchSize);
638     offset    = numCells - Nr;
639     geom.v0   = v0;
640     geom.J    = J;
641     geom.invJ = invJ;
642     geom.detJ = detJ;
643     ierr = PetscFEIntegrateResidual(fe[f], Ne, Nf, fe, f, geom, u, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr);
644     geom.v0   = &v0[offset*dim];
645     geom.J    = &J[offset*dim*dim];
646     geom.invJ = &invJ[offset*dim*dim];
647     geom.detJ = &detJ[offset];
648     ierr = PetscFEIntegrateResidual(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
649   }
650   for (c = cStart; c < cEnd; ++c) {
651     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
652     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
653   }
654   ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
655   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
656   if (feBd) {
657     DMLabel         label, depth;
658     IS              pointIS;
659     const PetscInt *points;
660     PetscInt        dep, numPoints, p, numFaces;
661     PetscReal      *n;
662 
663     ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
664     ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
665     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
666     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
667     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
668     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
669       PetscInt Nb, Nc;
670 
671       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
672       ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr);
673       cellDof       += Nb*Nc;
674       numComponents += Nc;
675     }
676     for (p = 0, numFaces = 0; p < numPoints; ++p) {
677       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
678       if (dep == dim-1) ++numFaces;
679     }
680     ierr = PetscMalloc7(numFaces*cellDof,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*cellDof,&elemVec);CHKERRQ(ierr);
681     for (p = 0, f = 0; p < numPoints; ++p) {
682       const PetscInt point = points[p];
683       PetscScalar   *x     = NULL;
684       PetscInt       i;
685 
686       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
687       if (dep != dim-1) continue;
688       ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
689       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
690       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
691       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
692       for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
693       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
694       ++f;
695     }
696     for (f = 0; f < Nf; ++f) {
697       void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[f];
698       void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[f];
699       PetscInt numQuadPoints, Nb;
700       /* Conforming batches */
701       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
702       /* Remainder */
703       PetscInt Nr, offset;
704 
705       ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr);
706       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
707       ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
708       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
709       blockSize = Nb*numQuadPoints;
710       batchSize = numBlocks * blockSize;
711       ierr =  PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
712       numChunks = numFaces / (numBatches*batchSize);
713       Ne        = numChunks*numBatches*batchSize;
714       Nr        = numFaces % (numBatches*batchSize);
715       offset    = numFaces - Nr;
716       geom.v0   = v0;
717       geom.n    = n;
718       geom.J    = J;
719       geom.invJ = invJ;
720       geom.detJ = detJ;
721       ierr = PetscFEIntegrateBdResidual(feBd[f], Ne, Nf, feBd, f, geom, u, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr);
722       geom.v0   = &v0[offset*dim];
723       geom.n    = &n[offset*dim];
724       geom.J    = &J[offset*dim*dim];
725       geom.invJ = &invJ[offset*dim*dim];
726       geom.detJ = &detJ[offset];
727       ierr = PetscFEIntegrateBdResidual(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
728     }
729     for (p = 0, f = 0; p < numPoints; ++p) {
730       const PetscInt point = points[p];
731 
732       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
733       if (dep != dim-1) continue;
734       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);}
735       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr);
736       ++f;
737     }
738     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
739     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
740     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
741   }
742   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
743   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
744   PetscFunctionReturn(0);
745 }
746 
747 #undef __FUNCT__
748 #define __FUNCT__ "DMPlexComputeIFunctionFEM"
749 /*@
750   DMPlexComputeIFunctionFEM - Form the local implicit function F from the local input X, X_t using pointwise functions specified by the user
751 
752   Input Parameters:
753 + dm - The mesh
754 . X  - Local input vector
755 . X_t  - Time derivative of the local input vector
756 - user - The user context
757 
758   Output Parameter:
759 . F  - Local output vector
760 
761   Note:
762   The first member of the user context must be an FEMContext.
763 
764   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
765   like a GPU, or vectorize on a multicore machine.
766 
767   Level: developer
768 
769 .seealso: DMPlexComputeResidualFEM()
770 @*/
771 PetscErrorCode DMPlexComputeIFunctionFEM(DM dm, Vec X, Vec X_t, Vec F, void *user)
772 {
773   DM_Plex          *mesh  = (DM_Plex *) dm->data;
774   PetscFEM         *fem   = (PetscFEM *) user;
775   PetscFE          *fe    = fem->fe;
776   PetscFE          *feAux = fem->feAux;
777   PetscFE          *feBd  = fem->feBd;
778   const char       *name  = "Residual";
779   DM                dmAux;
780   Vec               A;
781   PetscQuadrature   q;
782   PetscCellGeometry geom;
783   PetscSection      section, sectionAux;
784   PetscReal        *v0, *J, *invJ, *detJ;
785   PetscScalar      *elemVec, *u, *u_t, *a = NULL;
786   PetscInt          dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c;
787   PetscInt          cellDof = 0, numComponents = 0;
788   PetscInt          cellDofAux = 0, numComponentsAux = 0;
789   PetscErrorCode    ierr;
790 
791   PetscFunctionBegin;
792   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
793   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
794   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
795   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
796   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
797   numCells = cEnd - cStart;
798   for (f = 0; f < Nf; ++f) {
799     PetscInt Nb, Nc;
800 
801     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
802     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
803     cellDof       += Nb*Nc;
804     numComponents += Nc;
805   }
806   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
807   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
808   if (dmAux) {
809     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
810     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
811   }
812   for (f = 0; f < NfAux; ++f) {
813     PetscInt Nb, Nc;
814 
815     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
816     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
817     cellDofAux       += Nb*Nc;
818     numComponentsAux += Nc;
819   }
820   ierr = DMPlexProjectFunctionLocal(dm, fe, fem->bcFuncs, fem->bcCtxs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
821   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
822   ierr = PetscMalloc7(numCells*cellDof,&u,numCells*cellDof,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr);
823   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
824   for (c = cStart; c < cEnd; ++c) {
825     PetscScalar *x = NULL, *x_t = NULL;
826     PetscInt     i;
827 
828     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
829     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
830     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
831     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
832     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
833     ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
834     for (i = 0; i < cellDof; ++i) u_t[c*cellDof+i] = x_t[i];
835     ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
836     if (dmAux) {
837       PetscScalar *x_a = NULL;
838       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr);
839       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x_a[i];
840       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr);
841     }
842   }
843   for (f = 0; f < Nf; ++f) {
844     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0IFuncs[f];
845     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1IFuncs[f];
846     PetscInt numQuadPoints, Nb;
847     /* Conforming batches */
848     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
849     /* Remainder */
850     PetscInt Nr, offset;
851 
852     ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr);
853     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
854     ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
855     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
856     blockSize = Nb*numQuadPoints;
857     batchSize = numBlocks * blockSize;
858     ierr =  PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
859     numChunks = numCells / (numBatches*batchSize);
860     Ne        = numChunks*numBatches*batchSize;
861     Nr        = numCells % (numBatches*batchSize);
862     offset    = numCells - Nr;
863     geom.v0   = v0;
864     geom.J    = J;
865     geom.invJ = invJ;
866     geom.detJ = detJ;
867     ierr = PetscFEIntegrateIFunction(fe[f], Ne, Nf, fe, f, geom, u, u_t, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr);
868     geom.v0   = &v0[offset*dim];
869     geom.J    = &J[offset*dim*dim];
870     geom.invJ = &invJ[offset*dim*dim];
871     geom.detJ = &detJ[offset];
872     ierr = 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]);CHKERRQ(ierr);
873   }
874   for (c = cStart; c < cEnd; ++c) {
875     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
876     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
877   }
878   ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
879   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
880   if (feBd) {
881     DMLabel         label, depth;
882     IS              pointIS;
883     const PetscInt *points;
884     PetscInt        dep, numPoints, p, numFaces;
885     PetscReal      *n;
886 
887     ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
888     ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
889     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
890     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
891     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
892     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
893       PetscInt Nb, Nc;
894 
895       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
896       ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr);
897       cellDof       += Nb*Nc;
898       numComponents += Nc;
899     }
900     for (p = 0, numFaces = 0; p < numPoints; ++p) {
901       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
902       if (dep == dim-1) ++numFaces;
903     }
904     ierr = PetscMalloc7(numFaces*cellDof,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*cellDof,&elemVec);CHKERRQ(ierr);
905     ierr = PetscMalloc1(numFaces*cellDof,&u_t);CHKERRQ(ierr);
906     for (p = 0, f = 0; p < numPoints; ++p) {
907       const PetscInt point = points[p];
908       PetscScalar   *x     = NULL;
909       PetscInt       i;
910 
911       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
912       if (dep != dim-1) continue;
913       ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
914       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
915       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
916       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
917       for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
918       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
919       ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
920       for (i = 0; i < cellDof; ++i) u_t[f*cellDof+i] = x[i];
921       ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
922       ++f;
923     }
924     for (f = 0; f < Nf; ++f) {
925       void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdIFuncs[f];
926       void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdIFuncs[f];
927       PetscInt numQuadPoints, Nb;
928       /* Conforming batches */
929       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
930       /* Remainder */
931       PetscInt Nr, offset;
932 
933       ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr);
934       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
935       ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
936       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
937       blockSize = Nb*numQuadPoints;
938       batchSize = numBlocks * blockSize;
939       ierr =  PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
940       numChunks = numFaces / (numBatches*batchSize);
941       Ne        = numChunks*numBatches*batchSize;
942       Nr        = numFaces % (numBatches*batchSize);
943       offset    = numFaces - Nr;
944       geom.v0   = v0;
945       geom.n    = n;
946       geom.J    = J;
947       geom.invJ = invJ;
948       geom.detJ = detJ;
949       ierr = PetscFEIntegrateBdIFunction(feBd[f], Ne, Nf, feBd, f, geom, u, u_t, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr);
950       geom.v0   = &v0[offset*dim];
951       geom.n    = &n[offset*dim];
952       geom.J    = &J[offset*dim*dim];
953       geom.invJ = &invJ[offset*dim*dim];
954       geom.detJ = &detJ[offset];
955       ierr = PetscFEIntegrateBdIFunction(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], &u_t[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
956     }
957     for (p = 0, f = 0; p < numPoints; ++p) {
958       const PetscInt point = points[p];
959 
960       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
961       if (dep != dim-1) continue;
962       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);}
963       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr);
964       ++f;
965     }
966     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
967     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
968     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
969     ierr = PetscFree(u_t);CHKERRQ(ierr);
970   }
971   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
972   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
973   PetscFunctionReturn(0);
974 }
975 
976 #undef __FUNCT__
977 #define __FUNCT__ "DMPlexComputeJacobianActionFEM"
978 /*@C
979   DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user
980 
981   Input Parameters:
982 + dm - The mesh
983 . J  - The Jacobian shell matrix
984 . X  - Local input vector
985 - user - The user context
986 
987   Output Parameter:
988 . F  - Local output vector
989 
990   Note:
991   The first member of the user context must be an FEMContext.
992 
993   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
994   like a GPU, or vectorize on a multicore machine.
995 
996   Level: developer
997 
998 .seealso: DMPlexComputeResidualFEM()
999 @*/
1000 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
1001 {
1002   DM_Plex          *mesh = (DM_Plex *) dm->data;
1003   PetscFEM         *fem  = (PetscFEM *) user;
1004   PetscFE          *fe   = fem->fe;
1005   PetscQuadrature   quad;
1006   PetscCellGeometry geom;
1007   PetscSection      section;
1008   JacActionCtx     *jctx;
1009   PetscReal        *v0, *J, *invJ, *detJ;
1010   PetscScalar      *elemVec, *u, *a;
1011   PetscInt          dim, numFields, field, numCells, cStart, cEnd, c;
1012   PetscInt          cellDof = 0;
1013   PetscErrorCode    ierr;
1014 
1015   PetscFunctionBegin;
1016   /* ierr = PetscLogEventBegin(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */
1017   ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
1018   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
1019   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1020   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
1021   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1022   numCells = cEnd - cStart;
1023   for (field = 0; field < numFields; ++field) {
1024     PetscInt Nb, Nc;
1025 
1026     ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
1027     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
1028     cellDof += Nb*Nc;
1029   }
1030   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
1031   ierr = PetscMalloc7(numCells*cellDof,&u,numCells*cellDof,&a,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr);
1032   for (c = cStart; c < cEnd; ++c) {
1033     PetscScalar *x = NULL;
1034     PetscInt     i;
1035 
1036     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1037     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1038     ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
1039     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
1040     ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
1041     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
1042     for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
1043     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
1044   }
1045   for (field = 0; field < numFields; ++field) {
1046     PetscInt numQuadPoints, Nb;
1047     /* Conforming batches */
1048     PetscInt numBlocks  = 1;
1049     PetscInt numBatches = 1;
1050     PetscInt numChunks, Ne, blockSize, batchSize;
1051     /* Remainder */
1052     PetscInt Nr, offset;
1053 
1054     ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr);
1055     ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
1056     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1057     blockSize = Nb*numQuadPoints;
1058     batchSize = numBlocks * blockSize;
1059     numChunks = numCells / (numBatches*batchSize);
1060     Ne        = numChunks*numBatches*batchSize;
1061     Nr        = numCells % (numBatches*batchSize);
1062     offset    = numCells - Nr;
1063     geom.v0   = v0;
1064     geom.J    = J;
1065     geom.invJ = invJ;
1066     geom.detJ = detJ;
1067     ierr = PetscFEIntegrateJacobianAction(fe[field], Ne, numFields, fe, field, geom, u, a, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr);
1068     geom.v0   = &v0[offset*dim];
1069     geom.J    = &J[offset*dim*dim];
1070     geom.invJ = &invJ[offset*dim*dim];
1071     geom.detJ = &detJ[offset];
1072     ierr = PetscFEIntegrateJacobianAction(fe[field], Nr, numFields, fe, field, geom, &u[offset*cellDof], &a[offset*cellDof],
1073                                           fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr);
1074   }
1075   for (c = cStart; c < cEnd; ++c) {
1076     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
1077     ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
1078   }
1079   ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
1080   if (mesh->printFEM) {
1081     PetscMPIInt rank, numProcs;
1082     PetscInt    p;
1083 
1084     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
1085     ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
1086     ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");CHKERRQ(ierr);
1087     for (p = 0; p < numProcs; ++p) {
1088       if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}
1089       ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
1090     }
1091   }
1092   /* ierr = PetscLogEventEnd(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */
1093   PetscFunctionReturn(0);
1094 }
1095 
1096 #undef __FUNCT__
1097 #define __FUNCT__ "DMPlexComputeJacobianFEM"
1098 /*@
1099   DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
1100 
1101   Input Parameters:
1102 + dm - The mesh
1103 . X  - Local input vector
1104 - user - The user context
1105 
1106   Output Parameter:
1107 . Jac  - Jacobian matrix
1108 
1109   Note:
1110   The first member of the user context must be an FEMContext.
1111 
1112   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1113   like a GPU, or vectorize on a multicore machine.
1114 
1115   Level: developer
1116 
1117 .seealso: FormFunctionLocal()
1118 @*/
1119 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1120 {
1121   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1122   PetscFEM         *fem   = (PetscFEM *) user;
1123   PetscFE          *fe    = fem->fe;
1124   PetscFE          *feAux = fem->feAux;
1125   const char       *name  = "Jacobian";
1126   DM                dmAux;
1127   Vec               A;
1128   PetscQuadrature   quad;
1129   PetscCellGeometry geom;
1130   PetscSection      section, globalSection, sectionAux;
1131   PetscReal        *v0, *J, *invJ, *detJ;
1132   PetscScalar      *elemMat, *u, *a;
1133   PetscInt          dim, Nf, NfAux = 0, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1134   PetscInt          cellDof = 0, numComponents = 0;
1135   PetscInt          cellDofAux = 0, numComponentsAux = 0;
1136   PetscBool         isShell;
1137   PetscErrorCode    ierr;
1138 
1139   PetscFunctionBegin;
1140   ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1141   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
1142   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1143   ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);
1144   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1145   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1146   numCells = cEnd - cStart;
1147   for (f = 0; f < Nf; ++f) {
1148     PetscInt Nb, Nc;
1149 
1150     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
1151     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1152     cellDof       += Nb*Nc;
1153     numComponents += Nc;
1154   }
1155   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1156   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1157   if (dmAux) {
1158     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1159     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
1160   }
1161   for (f = 0; f < NfAux; ++f) {
1162     PetscInt Nb, Nc;
1163 
1164     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
1165     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
1166     cellDofAux       += Nb*Nc;
1167     numComponentsAux += Nc;
1168   }
1169   ierr = DMPlexProjectFunctionLocal(dm, fe, fem->bcFuncs, fem->bcCtxs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
1170   ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
1171   ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof*cellDof,&elemMat);CHKERRQ(ierr);
1172   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
1173   for (c = cStart; c < cEnd; ++c) {
1174     PetscScalar *x = NULL;
1175     PetscInt     i;
1176 
1177     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1178     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1179     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1180     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
1181     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1182     if (dmAux) {
1183       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1184       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
1185       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1186     }
1187   }
1188   ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1189   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1190     PetscInt numQuadPoints, Nb;
1191     /* Conforming batches */
1192     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1193     /* Remainder */
1194     PetscInt Nr, offset;
1195 
1196     ierr = PetscFEGetQuadrature(fe[fieldI], &quad);CHKERRQ(ierr);
1197     ierr = PetscFEGetDimension(fe[fieldI], &Nb);CHKERRQ(ierr);
1198     ierr = PetscFEGetTileSizes(fe[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1199     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1200     blockSize = Nb*numQuadPoints;
1201     batchSize = numBlocks * blockSize;
1202     ierr = PetscFESetTileSizes(fe[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1203     numChunks = numCells / (numBatches*batchSize);
1204     Ne        = numChunks*numBatches*batchSize;
1205     Nr        = numCells % (numBatches*batchSize);
1206     offset    = numCells - Nr;
1207     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1208       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ];
1209       void   (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*Nf+fieldJ];
1210       void   (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*Nf+fieldJ];
1211       void   (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*Nf+fieldJ];
1212 
1213       geom.v0   = v0;
1214       geom.J    = J;
1215       geom.invJ = invJ;
1216       geom.detJ = detJ;
1217       ierr = PetscFEIntegrateJacobian(fe[fieldI], Ne, Nf, fe, fieldI, fieldJ, geom, u, NfAux, feAux, a, g0, g1, g2, g3, elemMat);CHKERRQ(ierr);
1218       geom.v0   = &v0[offset*dim];
1219       geom.J    = &J[offset*dim*dim];
1220       geom.invJ = &invJ[offset*dim*dim];
1221       geom.detJ = &detJ[offset];
1222       ierr = 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]);CHKERRQ(ierr);
1223     }
1224   }
1225   for (c = cStart; c < cEnd; ++c) {
1226     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);}
1227     ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr);
1228   }
1229   ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1230   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1231   ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1232   ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1233   if (mesh->printFEM) {
1234     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1235     ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
1236     ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1237   }
1238   ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1239   ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr);
1240   if (isShell) {
1241     JacActionCtx *jctx;
1242 
1243     ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
1244     ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
1245   }
1246   PetscFunctionReturn(0);
1247 }
1248 
1249 #if 0
1250 
1251 static void g0_identity_1d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1252 {
1253   g3[0] = 1.0;
1254 }
1255 
1256 static void g0_identity_2d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1257 {
1258   g3[0] = g3[3] = 1.0;
1259 }
1260 
1261 static void g0_identity_3d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1262 {
1263   g3[0] = g3[4] = g3[8] = 1.0;
1264 }
1265 
1266 #undef __FUNCT__
1267 #define __FUNCT__ "DMPlexComputeInterpolatorFEMBroken"
1268 PetscErrorCode DMPlexComputeInterpolatorFEMBroken(DM dmc, DM dmf, Mat I, void *user)
1269 {
1270   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1271   PetscFEM         *fem   = (PetscFEM *) user;
1272   PetscFE          *fe    = fem->fe;
1273   const char       *name  = "Interpolator";
1274   PetscFE          *feRef;
1275   PetscQuadrature   quad, quadOld;
1276   PetscCellGeometry geom;
1277   PetscSection      fsection, fglobalSection;
1278   PetscSection      csection, cglobalSection;
1279   PetscReal        *v0, *J, *invJ, *detJ;
1280   PetscScalar      *elemMat;
1281   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1282   PetscInt          rCellDof = 0, cCellDof = 0, numComponents = 0;
1283   PetscErrorCode    ierr;
1284 
1285   PetscFunctionBegin;
1286 #if 0
1287   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1288 #endif
1289   ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr);
1290   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1291   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1292   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1293   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1294   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1295   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1296   numCells = cEnd - cStart;
1297   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1298   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1299     ierr = PetscFERefine(fe[fieldI], &feRef[fieldI]);CHKERRQ(ierr);
1300   }
1301   for (f = 0; f < Nf; ++f) {
1302     PetscInt rNb, cNb, Nc;
1303 
1304     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1305     ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr);
1306     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1307     numComponents += Nc;
1308     rCellDof += rNb*Nc;
1309     cCellDof += cNb*Nc;
1310   }
1311   ierr = MatZeroEntries(I);CHKERRQ(ierr);
1312   ierr = PetscMalloc5(numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*rCellDof*cCellDof,&elemMat);CHKERRQ(ierr);
1313   for (c = cStart; c < cEnd; ++c) {
1314     ierr = DMPlexComputeCellGeometry(dmc, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1315     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1316   }
1317   ierr = PetscMemzero(elemMat, numCells*rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1318   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1319     PetscInt numQuadPoints, Nb;
1320     /* Conforming batches */
1321     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1322     /* Remainder */
1323     PetscInt Nr, offset;
1324 
1325     /* Make new fine FE which refines the ref cell and the quadrature rule */
1326     ierr = PetscFEGetQuadrature(feRef[fieldI], &quad);CHKERRQ(ierr);
1327     ierr = PetscFEGetDimension(feRef[fieldI], &Nb);CHKERRQ(ierr);
1328     ierr = PetscFEGetTileSizes(feRef[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1329     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1330     blockSize = Nb*numQuadPoints;
1331     batchSize = numBlocks * blockSize;
1332     ierr = PetscFESetTileSizes(feRef[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1333     numChunks = numCells / (numBatches*batchSize);
1334     Ne        = numChunks*numBatches*batchSize;
1335     Nr        = numCells % (numBatches*batchSize);
1336     offset    = numCells - Nr;
1337 
1338     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1339       /* void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ]; */
1340       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = g0_identity_2d_static;
1341 
1342       /* Replace quadrature in coarse FE with refined quadrature */
1343       ierr = PetscFEGetQuadrature(fe[fieldJ], &quadOld);CHKERRQ(ierr);
1344       ierr = PetscObjectReference((PetscObject) quadOld);CHKERRQ(ierr);
1345       ierr = PetscFESetQuadrature(fe[fieldJ], quad);CHKERRQ(ierr);
1346       geom.v0   = v0;
1347       geom.J    = J;
1348       geom.invJ = invJ;
1349       geom.detJ = detJ;
1350       ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Ne, Nf, feRef, fieldI, fe, fieldJ, geom, g0, elemMat);CHKERRQ(ierr);
1351       geom.v0   = &v0[offset*dim];
1352       geom.J    = &J[offset*dim*dim];
1353       geom.invJ = &invJ[offset*dim*dim];
1354       geom.detJ = &detJ[offset];
1355       ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Nr, Nf, feRef, fieldI, fe, fieldJ, geom, g0, &elemMat[offset*rCellDof*cCellDof]);CHKERRQ(ierr);
1356       ierr = PetscFESetQuadrature(fe[fieldJ], quadOld);CHKERRQ(ierr);
1357     }
1358   }
1359   for (c = cStart; c < cEnd; ++c) {
1360     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, rCellDof, cCellDof, &elemMat[c*rCellDof*cCellDof]);CHKERRQ(ierr);}
1361     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, I, c, &elemMat[c*rCellDof*cCellDof], ADD_VALUES);CHKERRQ(ierr);
1362   }
1363   ierr = PetscFree5(v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1364   ierr = PetscFree(feRef);CHKERRQ(ierr);
1365   ierr = MatAssemblyBegin(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1366   ierr = MatAssemblyEnd(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1367   if (mesh->printFEM) {
1368     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1369     ierr = MatChop(I, 1.0e-10);CHKERRQ(ierr);
1370     ierr = MatView(I, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1371   }
1372 #if 0
1373   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1374 #endif
1375   PetscFunctionReturn(0);
1376 }
1377 #endif
1378 
1379 #undef __FUNCT__
1380 #define __FUNCT__ "DMPlexComputeInterpolatorFEM"
1381 /*@
1382   DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1383 
1384   Input Parameters:
1385 + dmf  - The fine mesh
1386 . dmc  - The coarse mesh
1387 - user - The user context
1388 
1389   Output Parameter:
1390 . In  - The interpolation matrix
1391 
1392   Note:
1393   The first member of the user context must be an FEMContext.
1394 
1395   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1396   like a GPU, or vectorize on a multicore machine.
1397 
1398   Level: developer
1399 
1400 .seealso: DMPlexComputeJacobianFEM()
1401 @*/
1402 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1403 {
1404   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1405   PetscFEM         *fem   = (PetscFEM *) user;
1406   PetscFE          *fe    = fem->fe;
1407   const char       *name  = "Interpolator";
1408   PetscFE          *feRef;
1409   PetscSection      fsection, fglobalSection;
1410   PetscSection      csection, cglobalSection;
1411   PetscScalar      *elemMat;
1412   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c;
1413   PetscInt          rCellDof = 0, cCellDof = 0;
1414   PetscErrorCode    ierr;
1415 
1416   PetscFunctionBegin;
1417 #if 0
1418   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1419 #endif
1420   ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr);
1421   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1422   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1423   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1424   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1425   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1426   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1427   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1428   for (f = 0; f < Nf; ++f) {
1429     PetscInt rNb, cNb, Nc;
1430 
1431     ierr = PetscFERefine(fe[f], &feRef[f]);CHKERRQ(ierr);
1432     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1433     ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr);
1434     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1435     rCellDof += rNb*Nc;
1436     cCellDof += cNb*Nc;
1437   }
1438   ierr = MatZeroEntries(In);CHKERRQ(ierr);
1439   ierr = PetscMalloc1(rCellDof*cCellDof,&elemMat);CHKERRQ(ierr);
1440   ierr = PetscMemzero(elemMat, rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1441   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1442     PetscDualSpace   Qref;
1443     PetscQuadrature  f;
1444     const PetscReal *qpoints, *qweights;
1445     PetscReal       *points;
1446     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;
1447 
1448     /* Compose points from all dual basis functionals */
1449     ierr = PetscFEGetNumComponents(fe[fieldI], &Nc);CHKERRQ(ierr);
1450     ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr);
1451     ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr);
1452     for (i = 0; i < fpdim; ++i) {
1453       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1454       ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr);
1455       npoints += Np;
1456     }
1457     ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr);
1458     for (i = 0, k = 0; i < fpdim; ++i) {
1459       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1460       ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr);
1461       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1462     }
1463 
1464     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1465       PetscReal *B;
1466       PetscInt   NcJ, cpdim, j;
1467 
1468       /* Evaluate basis at points */
1469       ierr = PetscFEGetNumComponents(fe[fieldJ], &NcJ);CHKERRQ(ierr);
1470       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);
1471       ierr = PetscFEGetDimension(fe[fieldJ], &cpdim);CHKERRQ(ierr);
1472       /* For now, fields only interpolate themselves */
1473       if (fieldI == fieldJ) {
1474         ierr = PetscFEGetTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);
1475         for (i = 0, k = 0; i < fpdim; ++i) {
1476           ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1477           ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr);
1478           for (p = 0; p < Np; ++p, ++k) {
1479             for (j = 0; j < cpdim; ++j) {
1480               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];
1481             }
1482           }
1483         }
1484         ierr = PetscFERestoreTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr);
1485       }
1486       offsetJ += cpdim*NcJ;
1487     }
1488     offsetI += fpdim*Nc;
1489     ierr = PetscFree(points);CHKERRQ(ierr);
1490   }
1491   if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rCellDof, cCellDof, elemMat);CHKERRQ(ierr);}
1492   for (c = cStart; c < cEnd; ++c) {
1493     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr);
1494   }
1495   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);}
1496   ierr = PetscFree(feRef);CHKERRQ(ierr);
1497   ierr = PetscFree(elemMat);CHKERRQ(ierr);
1498   ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1499   ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1500   if (mesh->printFEM) {
1501     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1502     ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr);
1503     ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1504   }
1505 #if 0
1506   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1507 #endif
1508   PetscFunctionReturn(0);
1509 }
1510 
1511 #undef __FUNCT__
1512 #define __FUNCT__ "DMPlexAddBoundary"
1513 /* The ids can be overridden by the command line option -bc_<boundary name> */
1514 PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], PetscInt field, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx)
1515 {
1516   DM_Plex       *mesh = (DM_Plex *) dm->data;
1517   DMBoundary     b;
1518   PetscErrorCode ierr;
1519 
1520   PetscFunctionBegin;
1521   ierr = PetscNew(&b);CHKERRQ(ierr);
1522   ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr);
1523   ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr);
1524   ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr);
1525   b->essential   = isEssential;
1526   b->field       = field;
1527   b->func        = bcFunc;
1528   b->numids      = numids;
1529   b->ctx         = ctx;
1530   b->next        = mesh->boundary;
1531   mesh->boundary = b;
1532   PetscFunctionReturn(0);
1533 }
1534 
1535 #undef __FUNCT__
1536 #define __FUNCT__ "DMPlexGetNumBoundary"
1537 PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd)
1538 {
1539   DM_Plex   *mesh = (DM_Plex *) dm->data;
1540   DMBoundary b    = mesh->boundary;
1541 
1542   PetscFunctionBegin;
1543   *numBd = 0;
1544   while (b) {++(*numBd); b = b->next;}
1545   PetscFunctionReturn(0);
1546 }
1547 
1548 #undef __FUNCT__
1549 #define __FUNCT__ "DMPlexGetBoundary"
1550 PetscErrorCode DMPlexGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, PetscInt *field, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx)
1551 {
1552   DM_Plex   *mesh = (DM_Plex *) dm->data;
1553   DMBoundary b    = mesh->boundary;
1554   PetscInt   n    = 0;
1555 
1556   PetscFunctionBegin;
1557   while (b) {
1558     if (n == bd) break;
1559     b = b->next;
1560     ++n;
1561   }
1562   if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
1563   if (isEssential) {
1564     PetscValidPointer(isEssential, 3);
1565     *isEssential = b->essential;
1566   }
1567   if (name) {
1568     PetscValidPointer(name, 4);
1569     *name = b->name;
1570   }
1571   if (field) {
1572     PetscValidPointer(field, 5);
1573     *field = b->field;
1574   }
1575   if (func) {
1576     PetscValidPointer(func, 6);
1577     *func = b->func;
1578   }
1579   if (numids) {
1580     PetscValidPointer(numids, 7);
1581     *numids = b->numids;
1582   }
1583   if (ids) {
1584     PetscValidPointer(ids, 8);
1585     *ids = b->ids;
1586   }
1587   if (ctx) {
1588     PetscValidPointer(ctx, 9);
1589     *ctx = b->ctx;
1590   }
1591   PetscFunctionReturn(0);
1592 }
1593