xref: /petsc/src/dm/impls/plex/plexfem.c (revision a0845e3a928e8c3de76a3c8bfba8b69f6cc922fe) 
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 = DMPlexGetCoordinateSection(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 This should be in a separate Discretization object, but I am not sure how to lay
184 it out yet, so I am stuffing things here while I experiment.
185 *******************************************************************************/
186 #undef __FUNCT__
187 #define __FUNCT__ "DMPlexSetFEMIntegration"
188 PetscErrorCode DMPlexSetFEMIntegration(DM dm,
189                                           PetscErrorCode (*integrateResidualFEM)(PetscInt, PetscInt, PetscFE[], PetscInt, PetscCellGeometry, const PetscScalar[],
190                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
191                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
192                                                                                  PetscScalar[]),
193                                           PetscErrorCode (*integrateBdResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
194                                                                                    const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
195                                                                                    void (*)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]),
196                                                                                    void (*)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]), PetscScalar[]),
197                                           PetscErrorCode (*integrateJacobianActionFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], const PetscScalar[],
198                                                                                        const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
199                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
200                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
201                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
202                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
203                                           PetscErrorCode (*integrateJacobianFEM)(PetscInt, PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
204                                                                                  const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
205                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
206                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
207                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
208                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]))
209 {
210   DM_Plex *mesh = (DM_Plex*) dm->data;
211 
212   PetscFunctionBegin;
213   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
214   mesh->integrateResidualFEM       = integrateResidualFEM;
215   mesh->integrateBdResidualFEM     = integrateBdResidualFEM;
216   mesh->integrateJacobianActionFEM = integrateJacobianActionFEM;
217   mesh->integrateJacobianFEM       = integrateJacobianFEM;
218   PetscFunctionReturn(0);
219 }
220 
221 #undef __FUNCT__
222 #define __FUNCT__ "DMPlexProjectFunctionLocal"
223 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscInt numComp, void (**funcs)(const PetscReal [], PetscScalar *), InsertMode mode, Vec localX)
224 {
225   Vec            coordinates;
226   PetscSection   section, cSection;
227   PetscInt       dim, vStart, vEnd, v, c, d;
228   PetscScalar   *values, *cArray;
229   PetscReal     *coords;
230   PetscErrorCode ierr;
231 
232   PetscFunctionBegin;
233   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
234   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
235   ierr = DMPlexGetCoordinateSection(dm, &cSection);CHKERRQ(ierr);
236   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
237   ierr = PetscMalloc(numComp * sizeof(PetscScalar), &values);CHKERRQ(ierr);
238   ierr = VecGetArray(coordinates, &cArray);CHKERRQ(ierr);
239   ierr = PetscSectionGetDof(cSection, vStart, &dim);CHKERRQ(ierr);
240   ierr = PetscMalloc(dim * sizeof(PetscReal),&coords);CHKERRQ(ierr);
241   for (v = vStart; v < vEnd; ++v) {
242     PetscInt dof, off;
243 
244     ierr = PetscSectionGetDof(cSection, v, &dof);CHKERRQ(ierr);
245     ierr = PetscSectionGetOffset(cSection, v, &off);CHKERRQ(ierr);
246     if (dof > dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot have more coordinates %d then dimensions %d", dof, dim);
247     for (d = 0; d < dof; ++d) coords[d] = PetscRealPart(cArray[off+d]);
248     for (c = 0; c < numComp; ++c) (*funcs[c])(coords, &values[c]);
249     ierr = VecSetValuesSection(localX, section, v, values, mode);CHKERRQ(ierr);
250   }
251   ierr = VecRestoreArray(coordinates, &cArray);CHKERRQ(ierr);
252   /* Temporary, must be replaced by a projection on the finite element basis */
253   {
254     PetscInt eStart = 0, eEnd = 0, e, depth;
255 
256     ierr = DMPlexGetLabelSize(dm, "depth", &depth);CHKERRQ(ierr);
257     --depth;
258     if (depth > 1) {ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);}
259     for (e = eStart; e < eEnd; ++e) {
260       const PetscInt *cone = NULL;
261       PetscInt        coneSize, d;
262       PetscScalar    *coordsA, *coordsB;
263 
264       ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr);
265       ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr);
266       if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Cone size %d for point %d should be 2", coneSize, e);
267       ierr = VecGetValuesSection(coordinates, cSection, cone[0], &coordsA);CHKERRQ(ierr);
268       ierr = VecGetValuesSection(coordinates, cSection, cone[1], &coordsB);CHKERRQ(ierr);
269       for (d = 0; d < dim; ++d) {
270         coords[d] = 0.5*(PetscRealPart(coordsA[d]) + PetscRealPart(coordsB[d]));
271       }
272       for (c = 0; c < numComp; ++c) (*funcs[c])(coords, &values[c]);
273       ierr = VecSetValuesSection(localX, section, e, values, mode);CHKERRQ(ierr);
274     }
275   }
276 
277   ierr = PetscFree(coords);CHKERRQ(ierr);
278   ierr = PetscFree(values);CHKERRQ(ierr);
279 #if 0
280   const PetscInt localDof = this->_mesh->sizeWithBC(s, *cells->begin());
281   PetscReal      detJ;
282 
283   ierr = PetscMalloc(localDof * sizeof(PetscScalar), &values);CHKERRQ(ierr);
284   ierr = PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&J);CHKERRQ(ierr);
285   ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV(PetscPowInt(this->_mesh->getSieve()->getMaxConeSize(),dim+1), true);
286 
287   for (PetscInt c = cStart; c < cEnd; ++c) {
288     ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*this->_mesh->getSieve(), c, pV);
289     const PETSC_MESH_TYPE::point_type *oPoints = pV.getPoints();
290     const int                          oSize   = pV.getSize();
291     int                                v       = 0;
292 
293     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr);
294     for (PetscInt cl = 0; cl < oSize; ++cl) {
295       const PetscInt fDim;
296 
297       ierr = PetscSectionGetDof(oPoints[cl], &fDim);CHKERRQ(ierr);
298       if (pointDim) {
299         for (PetscInt d = 0; d < fDim; ++d, ++v) {
300           values[v] = (*this->_options.integrate)(v0, J, v, initFunc);
301         }
302       }
303     }
304     ierr = DMPlexVecSetClosure(dm, NULL, localX, c, values);CHKERRQ(ierr);
305     pV.clear();
306   }
307   ierr = PetscFree2(v0,J);CHKERRQ(ierr);
308   ierr = PetscFree(values);CHKERRQ(ierr);
309 #endif
310   PetscFunctionReturn(0);
311 }
312 
313 #undef __FUNCT__
314 #define __FUNCT__ "DMPlexProjectFunction"
315 /*@C
316   DMPlexProjectFunction - This projects the given function into the function space provided.
317 
318   Input Parameters:
319 + dm      - The DM
320 . numComp - The number of components (functions)
321 . funcs   - The coordinate functions to evaluate
322 - mode    - The insertion mode for values
323 
324   Output Parameter:
325 . X - vector
326 
327   Level: developer
328 
329   Note:
330   This currently just calls the function with the coordinates of each vertex and edge midpoint, and stores the result in a vector.
331   We will eventually fix it.
332 
333 .seealso: DMPlexComputeL2Diff()
334 @*/
335 PetscErrorCode DMPlexProjectFunction(DM dm, PetscInt numComp, void (**funcs)(const PetscReal [], PetscScalar *), InsertMode mode, Vec X)
336 {
337   Vec            localX;
338   PetscErrorCode ierr;
339 
340   PetscFunctionBegin;
341   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
342   ierr = DMPlexProjectFunctionLocal(dm, numComp, funcs, mode, localX);CHKERRQ(ierr);
343   ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr);
344   ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr);
345   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
346   PetscFunctionReturn(0);
347 }
348 
349 #undef __FUNCT__
350 #define __FUNCT__ "DMPlexComputeL2Diff"
351 /*@C
352   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
353 
354   Input Parameters:
355 + dm    - The DM
356 . quad  - The PetscQuadrature object for each field
357 . funcs - The functions to evaluate for each field component
358 - X     - The coefficient vector u_h
359 
360   Output Parameter:
361 . diff - The diff ||u - u_h||_2
362 
363   Level: developer
364 
365 .seealso: DMPlexProjectFunction()
366 @*/
367 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscQuadrature quad[], void (**funcs)(const PetscReal [], PetscScalar *), Vec X, PetscReal *diff)
368 {
369   const PetscInt debug = 0;
370   PetscSection   section;
371   Vec            localX;
372   PetscReal     *coords, *v0, *J, *invJ, detJ;
373   PetscReal      localDiff = 0.0;
374   PetscInt       dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
375   PetscErrorCode ierr;
376 
377   PetscFunctionBegin;
378   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
379   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
380   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
381   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
382   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
383   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
384   for (field = 0; field < numFields; ++field) {
385     numComponents += quad[field].numComponents;
386   }
387   ierr = DMPlexProjectFunctionLocal(dm, numComponents, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
388   ierr = PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);CHKERRQ(ierr);
389   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
390   for (c = cStart; c < cEnd; ++c) {
391     PetscScalar *x;
392     PetscReal    elemDiff = 0.0;
393 
394     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
395     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
396     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
397 
398     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
399       const PetscInt   numQuadPoints = quad[field].numQuadPoints;
400       const PetscReal *quadPoints    = quad[field].quadPoints;
401       const PetscReal *quadWeights   = quad[field].quadWeights;
402       const PetscInt   numBasisFuncs = quad[field].numBasisFuncs;
403       const PetscInt   numBasisComps = quad[field].numComponents;
404       const PetscReal *basis         = quad[field].basis;
405       PetscInt         q, d, e, fc, f;
406 
407       if (debug) {
408         char title[1024];
409         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
410         ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr);
411       }
412       for (q = 0; q < numQuadPoints; ++q) {
413         for (d = 0; d < dim; d++) {
414           coords[d] = v0[d];
415           for (e = 0; e < dim; e++) {
416             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
417           }
418         }
419         for (fc = 0; fc < numBasisComps; ++fc) {
420           PetscScalar funcVal;
421           PetscScalar interpolant = 0.0;
422 
423           (*funcs[comp+fc])(coords, &funcVal);
424           for (f = 0; f < numBasisFuncs; ++f) {
425             const PetscInt fidx = f*numBasisComps+fc;
426             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
427           }
428           if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "    elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal))*quadWeights[q]*detJ);CHKERRQ(ierr);}
429           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal))*quadWeights[q]*detJ;
430         }
431       }
432       comp        += numBasisComps;
433       fieldOffset += numBasisFuncs*numBasisComps;
434     }
435     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
436     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
437     localDiff += elemDiff;
438   }
439   ierr  = PetscFree4(coords,v0,J,invJ);CHKERRQ(ierr);
440   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
441   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
442   *diff = PetscSqrtReal(*diff);
443   PetscFunctionReturn(0);
444 }
445 
446 #if 0
447 
448 #undef __FUNCT__
449 #define __FUNCT__ "DMPlexComputeResidualFEM"
450 /*@
451   DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
452 
453   Input Parameters:
454 + dm - The mesh
455 . X  - Local input vector
456 - user - The user context
457 
458   Output Parameter:
459 . F  - Local output vector
460 
461   Note:
462   The second member of the user context must be an FEMContext.
463 
464   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
465   like a GPU, or vectorize on a multicore machine.
466 
467   Level: developer
468 
469 .seealso: DMPlexComputeJacobianActionFEM()
470 @*/
471 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
472 {
473   DM_Plex         *mesh   = (DM_Plex*) dm->data;
474   PetscFEM        *fem    = (PetscFEM*) &((DM*) user)[1];
475   PetscQuadrature *quad   = fem->quad;
476   PetscQuadrature *quadBd = fem->quadBd;
477   PetscSection     section;
478   PetscReal       *v0, *n, *J, *invJ, *detJ;
479   PetscScalar     *elemVec, *u;
480   PetscInt         dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
481   PetscInt         cellDof, numComponents;
482   PetscBool        has;
483   PetscErrorCode   ierr;
484 
485   PetscFunctionBegin;
486   /* ierr = PetscLogEventBegin(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */
487   ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
488   ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
489   ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
490   ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
491   numCells = cEnd - cStart;
492   for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) {
493     cellDof       += quad[field].numBasisFuncs*quad[field].numComponents;
494     numComponents += quad[field].numComponents;
495   }
496   ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
497   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
498   ierr = PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr);
499   for (c = cStart; c < cEnd; ++c) {
500     PetscScalar *x;
501     PetscInt     i;
502 
503     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
504     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
505     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
506 
507     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
508     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
509   }
510   for (field = 0; field < numFields; ++field) {
511     const PetscInt numQuadPoints = quad[field].numQuadPoints;
512     const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
513     void           (*f0)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[field];
514     void           (*f1)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[field];
515     /* Conforming batches */
516     PetscInt blockSize  = numBasisFuncs*numQuadPoints;
517     PetscInt numBlocks  = 1;
518     PetscInt batchSize  = numBlocks * blockSize;
519     PetscInt numBatches = numBatchesTmp;
520     PetscInt numChunks  = numCells / (numBatches*batchSize);
521     /* Remainder */
522     PetscInt numRemainder = numCells % (numBatches * batchSize);
523     PetscInt offset       = numCells - numRemainder;
524 
525     ierr = (*mesh->integrateResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, v0, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr);
526     ierr = (*mesh->integrateResidualFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
527                                          f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
528   }
529   for (c = cStart; c < cEnd; ++c) {
530     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Residual", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
531     ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
532   }
533   ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
534   /* Integration over the boundary:
535      - This can probably be generalized to integration over a set of labels, however
536        the idea here is to do integration where we need the cell normal
537      - We can replace hardcoding with a registration process, and this is how we hook
538        up the system to something like FEniCS
539   */
540   ierr = DMPlexHasLabel(dm, "boundary", &has);CHKERRQ(ierr);
541   if (has && quadBd) {
542     DMLabel         label;
543     IS              pointIS;
544     const PetscInt *points;
545     PetscInt        numPoints, p;
546 
547     ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
548     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
549     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
550     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
551     for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) {
552       cellDof       += quadBd[field].numBasisFuncs*quadBd[field].numComponents;
553       numComponents += quadBd[field].numComponents;
554     }
555     ierr = PetscMalloc7(numPoints*cellDof,PetscScalar,&u,numPoints*dim,PetscReal,&v0,numPoints*dim,PetscReal,&n,numPoints*dim*dim,PetscReal,&J,numPoints*dim*dim,PetscReal,&invJ,numPoints,PetscReal,&detJ,numPoints*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr);
556     for (p = 0; p < numPoints; ++p) {
557       const PetscInt point = points[p];
558       PetscScalar   *x;
559       PetscInt       i;
560 
561       /* TODO: Add normal determination here */
562       ierr = DMPlexComputeCellGeometry(dm, point, &v0[p*dim], &J[p*dim*dim], &invJ[p*dim*dim], &detJ[p]);CHKERRQ(ierr);
563       if (detJ[p] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[p], point);
564       ierr = DMPlexVecGetClosure(dm, NULL, X, point, NULL, &x);CHKERRQ(ierr);
565 
566       for (i = 0; i < cellDof; ++i) u[p*cellDof+i] = x[i];
567       ierr = DMPlexVecRestoreClosure(dm, NULL, X, point, NULL, &x);CHKERRQ(ierr);
568     }
569     for (field = 0; field < numFields; ++field) {
570       const PetscInt numQuadPoints = quadBd[field].numQuadPoints;
571       const PetscInt numBasisFuncs = quadBd[field].numBasisFuncs;
572       void           (*f0)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[field];
573       void           (*f1)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[field];
574       /* Conforming batches */
575       PetscInt blockSize  = numBasisFuncs*numQuadPoints;
576       PetscInt numBlocks  = 1;
577       PetscInt batchSize  = numBlocks * blockSize;
578       PetscInt numBatches = numBatchesTmp;
579       PetscInt numChunks  = numPoints / (numBatches*batchSize);
580       /* Remainder */
581       PetscInt numRemainder = numPoints % (numBatches * batchSize);
582       PetscInt offset       = numPoints - numRemainder;
583 
584       ierr = (*mesh->integrateBdResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quadBd, u, v0, n, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr);
585       ierr = (*mesh->integrateBdResidualFEM)(numRemainder, numFields, field, quadBd, &u[offset*cellDof], &v0[offset*dim], &n[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
586                                              f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
587     }
588     for (p = 0; p < numPoints; ++p) {
589       const PetscInt point = points[p];
590 
591       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "Residual", cellDof, &elemVec[p*cellDof]);CHKERRQ(ierr);}
592       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[p*cellDof], ADD_VALUES);CHKERRQ(ierr);
593     }
594     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
595     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
596     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
597   }
598   if (mesh->printFEM) {
599     PetscMPIInt rank, numProcs;
600     PetscInt    p;
601 
602     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
603     ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
604     ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Residual:\n");CHKERRQ(ierr);
605     for (p = 0; p < numProcs; ++p) {
606       if (p == rank) {
607         Vec f;
608 
609         ierr = VecDuplicate(F, &f);CHKERRQ(ierr);
610         ierr = VecCopy(F, f);CHKERRQ(ierr);
611         ierr = VecChop(f, 1.0e-10);CHKERRQ(ierr);
612         ierr = VecView(f, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
613         ierr = VecDestroy(&f);CHKERRQ(ierr);
614         ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
615       }
616       ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
617     }
618   }
619   /* ierr = PetscLogEventEnd(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */
620   PetscFunctionReturn(0);
621 }
622 
623 #else
624 
625 #undef __FUNCT__
626 #define __FUNCT__ "DMPlexComputeResidualFEM"
627 /*@
628   DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
629 
630   Input Parameters:
631 + dm - The mesh
632 . X  - Local input vector
633 - user - The user context
634 
635   Output Parameter:
636 . F  - Local output vector
637 
638   Note:
639   The second member of the user context must be an FEMContext.
640 
641   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
642   like a GPU, or vectorize on a multicore machine.
643 
644   Level: developer
645 
646 .seealso: DMPlexComputeJacobianActionFEM()
647 @*/
648 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
649 {
650   DM_Plex         *mesh = (DM_Plex*) dm->data;
651   PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
652   PetscFE         *fe   = fem->fe;
653   const char      *name = "Residual";
654   PetscQuadrature  q;
655   PetscCellGeometry geom;
656   PetscSection     section;
657   PetscReal       *v0, *J, *invJ, *detJ;
658   PetscScalar     *elemVec, *u;
659   PetscInt         dim, numFields, f, numCells, cStart, cEnd, c;
660   PetscInt         cellDof = 0, numComponents = 0;
661   PetscErrorCode   ierr;
662 
663   PetscFunctionBegin;
664   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
665   ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
666   ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
667   ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
668   ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
669   numCells = cEnd - cStart;
670   for (f = 0; f < numFields; ++f) {
671     PetscInt Nb, Nc;
672 
673     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
674     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
675     cellDof       += Nb*Nc;
676     numComponents += Nc;
677   }
678   ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
679   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
680   ierr = PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr);
681   for (c = cStart; c < cEnd; ++c) {
682     PetscScalar *x = NULL;
683     PetscInt     i;
684 
685     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
686     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
687     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
688     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
689     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
690   }
691   for (f = 0; f < numFields; ++f) {
692     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscScalar[]) = fem->f0Funcs[f];
693     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscScalar[]) = fem->f1Funcs[f];
694     PetscInt Nb;
695     /* Conforming batches */
696     PetscInt numBlocks  = 1;
697     PetscInt numBatches = 1;
698     PetscInt numChunks, Ne, blockSize, batchSize;
699     /* Remainder */
700     PetscInt Nr, offset;
701 
702     ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr);
703     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
704     blockSize = Nb*q.numQuadPoints;
705     batchSize = numBlocks * blockSize;
706     numChunks = numCells / (numBatches*batchSize);
707     Ne        = numChunks*numBatches*batchSize;
708     Nr        = numCells % (numBatches*batchSize);
709     offset    = numCells - Nr;
710     geom.v0   = v0;
711     geom.J    = J;
712     geom.invJ = invJ;
713     geom.detJ = detJ;
714     ierr = (*mesh->integrateResidualFEM)(Ne, numFields, fe, f, geom, u, f0, f1, elemVec);CHKERRQ(ierr);
715     geom.v0   = &v0[offset*dim];
716     geom.J    = &J[offset*dim*dim];
717     geom.invJ = &invJ[offset*dim*dim];
718     geom.detJ = &detJ[offset];
719     ierr = (*mesh->integrateResidualFEM)(Nr, numFields, fe, f, geom, &u[offset*cellDof], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
720   }
721   for (c = cStart; c < cEnd; ++c) {
722     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
723     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
724   }
725   ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
726   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, F);CHKERRQ(ierr);}
727   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
728   PetscFunctionReturn(0);
729 }
730 
731 #endif
732 
733 #undef __FUNCT__
734 #define __FUNCT__ "DMPlexComputeJacobianActionFEM"
735 /*@C
736   DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user
737 
738   Input Parameters:
739 + dm - The mesh
740 . J  - The Jacobian shell matrix
741 . X  - Local input vector
742 - user - The user context
743 
744   Output Parameter:
745 . F  - Local output vector
746 
747   Note:
748   The second member of the user context must be an FEMContext.
749 
750   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
751   like a GPU, or vectorize on a multicore machine.
752 
753   Level: developer
754 
755 .seealso: DMPlexComputeResidualFEM()
756 @*/
757 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
758 {
759   DM_Plex         *mesh = (DM_Plex*) dm->data;
760   PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
761   PetscQuadrature *quad = fem->quad;
762   PetscSection     section;
763   JacActionCtx    *jctx;
764   PetscReal       *v0, *J, *invJ, *detJ;
765   PetscScalar     *elemVec, *u, *a;
766   PetscInt         dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
767   PetscInt         cellDof = 0;
768   PetscErrorCode   ierr;
769 
770   PetscFunctionBegin;
771   /* ierr = PetscLogEventBegin(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */
772   ierr     = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
773   ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
774   ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
775   ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
776   ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
777   numCells = cEnd - cStart;
778   for (field = 0; field < numFields; ++field) {
779     cellDof += quad[field].numBasisFuncs*quad[field].numComponents;
780   }
781   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
782   ierr = PetscMalloc7(numCells*cellDof,PetscScalar,&u,numCells*cellDof,PetscScalar,&a,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr);
783   for (c = cStart; c < cEnd; ++c) {
784     PetscScalar *x;
785     PetscInt     i;
786 
787     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
788     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
789     ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
790     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
791     ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
792     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
793     for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
794     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
795   }
796   for (field = 0; field < numFields; ++field) {
797     const PetscInt numQuadPoints = quad[field].numQuadPoints;
798     const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
799     /* Conforming batches */
800     PetscInt blockSize  = numBasisFuncs*numQuadPoints;
801     PetscInt numBlocks  = 1;
802     PetscInt batchSize  = numBlocks * blockSize;
803     PetscInt numBatches = numBatchesTmp;
804     PetscInt numChunks  = numCells / (numBatches*batchSize);
805     /* Remainder */
806     PetscInt numRemainder = numCells % (numBatches * batchSize);
807     PetscInt offset       = numCells - numRemainder;
808 
809     ierr = (*mesh->integrateJacobianActionFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, a, v0, J, invJ, detJ, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr);
810     ierr = (*mesh->integrateJacobianActionFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &a[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
811                                                fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr);
812   }
813   for (c = cStart; c < cEnd; ++c) {
814     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
815     ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
816   }
817   ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
818   if (mesh->printFEM) {
819     PetscMPIInt rank, numProcs;
820     PetscInt    p;
821 
822     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
823     ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
824     ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");CHKERRQ(ierr);
825     for (p = 0; p < numProcs; ++p) {
826       if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}
827       ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
828     }
829   }
830   /* ierr = PetscLogEventEnd(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */
831   PetscFunctionReturn(0);
832 }
833 
834 #undef __FUNCT__
835 #define __FUNCT__ "DMPlexComputeJacobianFEM"
836 /*@
837   DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
838 
839   Input Parameters:
840 + dm - The mesh
841 . X  - Local input vector
842 - user - The user context
843 
844   Output Parameter:
845 . Jac  - Jacobian matrix
846 
847   Note:
848   The second member of the user context must be an FEMContext.
849 
850   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
851   like a GPU, or vectorize on a multicore machine.
852 
853   Level: developer
854 
855 .seealso: FormFunctionLocal()
856 @*/
857 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, MatStructure *str,void *user)
858 {
859   DM_Plex         *mesh = (DM_Plex*) dm->data;
860   PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
861   PetscQuadrature *quad = fem->quad;
862   PetscSection     section;
863   PetscReal       *v0, *J, *invJ, *detJ;
864   PetscScalar     *elemMat, *u;
865   PetscInt         dim, numFields, field, fieldI, numBatchesTmp = 1, numCells, cStart, cEnd, c;
866   PetscInt         cellDof = 0, numComponents = 0;
867   PetscBool        isShell;
868   PetscErrorCode   ierr;
869 
870   PetscFunctionBegin;
871   /* ierr = PetscLogEventBegin(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */
872   ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
873   ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
874   ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
875   ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
876   numCells = cEnd - cStart;
877   for (field = 0; field < numFields; ++field) {
878     cellDof       += quad[field].numBasisFuncs*quad[field].numComponents;
879     numComponents += quad[field].numComponents;
880   }
881   ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
882   ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
883   ierr = PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof*cellDof,PetscScalar,&elemMat);CHKERRQ(ierr);
884   for (c = cStart; c < cEnd; ++c) {
885     PetscScalar *x;
886     PetscInt     i;
887 
888     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
889     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
890     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
891 
892     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
893     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
894   }
895   ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr);
896   for (fieldI = 0; fieldI < numFields; ++fieldI) {
897     const PetscInt numQuadPoints = quad[fieldI].numQuadPoints;
898     const PetscInt numBasisFuncs = quad[fieldI].numBasisFuncs;
899     PetscInt       fieldJ;
900 
901     for (fieldJ = 0; fieldJ < numFields; ++fieldJ) {
902       void (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*numFields+fieldJ];
903       void (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*numFields+fieldJ];
904       void (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*numFields+fieldJ];
905       void (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*numFields+fieldJ];
906       /* Conforming batches */
907       PetscInt blockSize  = numBasisFuncs*numQuadPoints;
908       PetscInt numBlocks  = 1;
909       PetscInt batchSize  = numBlocks * blockSize;
910       PetscInt numBatches = numBatchesTmp;
911       PetscInt numChunks  = numCells / (numBatches*batchSize);
912       /* Remainder */
913       PetscInt numRemainder = numCells % (numBatches * batchSize);
914       PetscInt offset       = numCells - numRemainder;
915 
916       ierr = (*mesh->integrateJacobianFEM)(numChunks*numBatches*batchSize, numFields, fieldI, fieldJ, quad, u, v0, J, invJ, detJ, g0, g1, g2, g3, elemMat);CHKERRQ(ierr);
917       ierr = (*mesh->integrateJacobianFEM)(numRemainder, numFields, fieldI, fieldJ, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
918                                            g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);CHKERRQ(ierr);
919     }
920   }
921   for (c = cStart; c < cEnd; ++c) {
922     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, "Jacobian", cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);}
923     ierr = DMPlexMatSetClosure(dm, NULL, NULL, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr);
924   }
925   ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
926 
927   /* Assemble matrix, using the 2-step process:
928        MatAssemblyBegin(), MatAssemblyEnd(). */
929   ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
930   ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
931 
932   if (mesh->printFEM) {
933     ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian:\n");CHKERRQ(ierr);
934     ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
935     ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
936   }
937   /* ierr = PetscLogEventEnd(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */
938   ierr = PetscObjectTypeCompare((PetscObject)Jac, MATSHELL, &isShell);CHKERRQ(ierr);
939   if (isShell) {
940     JacActionCtx *jctx;
941 
942     ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
943     ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
944   }
945   *str = SAME_NONZERO_PATTERN;
946   PetscFunctionReturn(0);
947 }
948