xref: /petsc/src/dm/impls/plex/plexfem.c (revision 0a6ece4efbe0fc0d64b36d89d274ad0f89f90441) 
1 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
3 #undef __FUNCT__
4 #define __FUNCT__ "DMPlexGetScale"
5 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
6 {
7   DM_Plex *mesh = (DM_Plex*) dm->data;
8 
9   PetscFunctionBegin;
10   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
11   PetscValidPointer(scale, 3);
12   *scale = mesh->scale[unit];
13   PetscFunctionReturn(0);
14 }
15 
16 #undef __FUNCT__
17 #define __FUNCT__ "DMPlexSetScale"
18 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
19 {
20   DM_Plex *mesh = (DM_Plex*) dm->data;
21 
22   PetscFunctionBegin;
23   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
24   mesh->scale[unit] = scale;
25   PetscFunctionReturn(0);
26 }
27 
28 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
29 {
30   switch (i) {
31   case 0:
32     switch (j) {
33     case 0: return 0;
34     case 1:
35       switch (k) {
36       case 0: return 0;
37       case 1: return 0;
38       case 2: return 1;
39       }
40     case 2:
41       switch (k) {
42       case 0: return 0;
43       case 1: return -1;
44       case 2: return 0;
45       }
46     }
47   case 1:
48     switch (j) {
49     case 0:
50       switch (k) {
51       case 0: return 0;
52       case 1: return 0;
53       case 2: return -1;
54       }
55     case 1: return 0;
56     case 2:
57       switch (k) {
58       case 0: return 1;
59       case 1: return 0;
60       case 2: return 0;
61       }
62     }
63   case 2:
64     switch (j) {
65     case 0:
66       switch (k) {
67       case 0: return 0;
68       case 1: return 1;
69       case 2: return 0;
70       }
71     case 1:
72       switch (k) {
73       case 0: return -1;
74       case 1: return 0;
75       case 2: return 0;
76       }
77     case 2: return 0;
78     }
79   }
80   return 0;
81 }
82 
83 #undef __FUNCT__
84 #define __FUNCT__ "DMPlexCreateRigidBody"
85 /*@C
86   DMPlexCreateRigidBody - create rigid body modes from coordinates
87 
88   Collective on DM
89 
90   Input Arguments:
91 + dm - the DM
92 . section - the local section associated with the rigid field, or NULL for the default section
93 - globalSection - the global section associated with the rigid field, or NULL for the default section
94 
95   Output Argument:
96 . sp - the null space
97 
98   Note: This is necessary to take account of Dirichlet conditions on the displacements
99 
100   Level: advanced
101 
102 .seealso: MatNullSpaceCreate()
103 @*/
104 PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp)
105 {
106   MPI_Comm       comm;
107   Vec            coordinates, localMode, mode[6];
108   PetscSection   coordSection;
109   PetscScalar   *coords;
110   PetscInt       dim, vStart, vEnd, v, n, m, d, i, j;
111   PetscErrorCode ierr;
112 
113   PetscFunctionBegin;
114   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
115   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
116   if (dim == 1) {
117     ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr);
118     PetscFunctionReturn(0);
119   }
120   if (!section)       {ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);}
121   if (!globalSection) {ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);}
122   ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr);
123   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
124   ierr = DMPlexGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
125   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
126   m    = (dim*(dim+1))/2;
127   ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr);
128   ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr);
129   ierr = VecSetUp(mode[0]);CHKERRQ(ierr);
130   for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);}
131   /* Assume P1 */
132   ierr = DMGetLocalVector(dm, &localMode);CHKERRQ(ierr);
133   for (d = 0; d < dim; ++d) {
134     PetscScalar values[3] = {0.0, 0.0, 0.0};
135 
136     values[d] = 1.0;
137     ierr      = VecSet(localMode, 0.0);CHKERRQ(ierr);
138     for (v = vStart; v < vEnd; ++v) {
139       ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
140     }
141     ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
142     ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
143   }
144   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
145   for (d = dim; d < dim*(dim+1)/2; ++d) {
146     PetscInt i, j, k = dim > 2 ? d - dim : d;
147 
148     ierr = VecSet(localMode, 0.0);CHKERRQ(ierr);
149     for (v = vStart; v < vEnd; ++v) {
150       PetscScalar values[3] = {0.0, 0.0, 0.0};
151       PetscInt    off;
152 
153       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
154       for (i = 0; i < dim; ++i) {
155         for (j = 0; j < dim; ++j) {
156           values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]);
157         }
158       }
159       ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
160     }
161     ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
162     ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
163   }
164   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
165   ierr = DMRestoreLocalVector(dm, &localMode);CHKERRQ(ierr);
166   for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);}
167   /* Orthonormalize system */
168   for (i = dim; i < m; ++i) {
169     PetscScalar dots[6];
170 
171     ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr);
172     for (j = 0; j < i; ++j) dots[j] *= -1.0;
173     ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr);
174     ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);
175   }
176   ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr);
177   for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);}
178   PetscFunctionReturn(0);
179 }
180 /*******************************************************************************
181 This should be in a separate Discretization object, but I am not sure how to lay
182 it out yet, so I am stuffing things here while I experiment.
183 *******************************************************************************/
184 #undef __FUNCT__
185 #define __FUNCT__ "DMPlexSetFEMIntegration"
186 PetscErrorCode DMPlexSetFEMIntegration(DM dm,
187                                           PetscErrorCode (*integrateResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
188                                                                                  const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
189                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
190                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
191                                           PetscErrorCode (*integrateBdResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
192                                                                                    const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
193                                                                                    void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]),
194                                                                                    void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]), PetscScalar[]),
195                                           PetscErrorCode (*integrateJacobianActionFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], const PetscScalar[],
196                                                                                        const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
197                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
198                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
199                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
200                                                                                        void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
201                                           PetscErrorCode (*integrateJacobianFEM)(PetscInt, PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
202                                                                                  const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
203                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
204                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
205                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
206                                                                                  void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]))
207 {
208   DM_Plex *mesh = (DM_Plex*) dm->data;
209 
210   PetscFunctionBegin;
211   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
212   mesh->integrateResidualFEM       = integrateResidualFEM;
213   mesh->integrateBdResidualFEM     = integrateBdResidualFEM;
214   mesh->integrateJacobianActionFEM = integrateJacobianActionFEM;
215   mesh->integrateJacobianFEM       = integrateJacobianFEM;
216   PetscFunctionReturn(0);
217 }
218 
219 #undef __FUNCT__
220 #define __FUNCT__ "DMPlexProjectFunctionLocal"
221 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscInt numComp, void (**funcs)(const PetscReal [], PetscScalar *), InsertMode mode, Vec localX)
222 {
223   Vec            coordinates;
224   PetscSection   section, cSection;
225   PetscInt       dim, vStart, vEnd, v, c, d;
226   PetscScalar   *values, *cArray;
227   PetscReal     *coords;
228   PetscErrorCode ierr;
229 
230   PetscFunctionBegin;
231   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
232   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
233   ierr = DMPlexGetCoordinateSection(dm, &cSection);CHKERRQ(ierr);
234   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
235   ierr = PetscMalloc(numComp * sizeof(PetscScalar), &values);CHKERRQ(ierr);
236   ierr = VecGetArray(coordinates, &cArray);CHKERRQ(ierr);
237   ierr = PetscSectionGetDof(cSection, vStart, &dim);CHKERRQ(ierr);
238   ierr = PetscMalloc(dim * sizeof(PetscReal),&coords);CHKERRQ(ierr);
239   for (v = vStart; v < vEnd; ++v) {
240     PetscInt dof, off;
241 
242     ierr = PetscSectionGetDof(cSection, v, &dof);CHKERRQ(ierr);
243     ierr = PetscSectionGetOffset(cSection, v, &off);CHKERRQ(ierr);
244     if (dof > dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot have more coordinates %d then dimensions %d", dof, dim);
245     for (d = 0; d < dof; ++d) coords[d] = PetscRealPart(cArray[off+d]);
246     for (c = 0; c < numComp; ++c) (*funcs[c])(coords, &values[c]);
247     ierr = VecSetValuesSection(localX, section, v, values, mode);CHKERRQ(ierr);
248   }
249   ierr = VecRestoreArray(coordinates, &cArray);CHKERRQ(ierr);
250   /* Temporary, must be replaced by a projection on the finite element basis */
251   {
252     PetscInt eStart = 0, eEnd = 0, e, depth;
253 
254     ierr = DMPlexGetLabelSize(dm, "depth", &depth);CHKERRQ(ierr);
255     --depth;
256     if (depth > 1) {ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);}
257     for (e = eStart; e < eEnd; ++e) {
258       const PetscInt *cone = NULL;
259       PetscInt        coneSize, d;
260       PetscScalar    *coordsA, *coordsB;
261 
262       ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr);
263       ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr);
264       if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Cone size %d for point %d should be 2", coneSize, e);
265       ierr = VecGetValuesSection(coordinates, cSection, cone[0], &coordsA);CHKERRQ(ierr);
266       ierr = VecGetValuesSection(coordinates, cSection, cone[1], &coordsB);CHKERRQ(ierr);
267       for (d = 0; d < dim; ++d) {
268         coords[d] = 0.5*(PetscRealPart(coordsA[d]) + PetscRealPart(coordsB[d]));
269       }
270       for (c = 0; c < numComp; ++c) (*funcs[c])(coords, &values[c]);
271       ierr = VecSetValuesSection(localX, section, e, values, mode);CHKERRQ(ierr);
272     }
273   }
274 
275   ierr = PetscFree(coords);CHKERRQ(ierr);
276   ierr = PetscFree(values);CHKERRQ(ierr);
277 #if 0
278   const PetscInt localDof = this->_mesh->sizeWithBC(s, *cells->begin());
279   PetscReal      detJ;
280 
281   ierr = PetscMalloc(localDof * sizeof(PetscScalar), &values);CHKERRQ(ierr);
282   ierr = PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&J);CHKERRQ(ierr);
283   ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV(PetscPowInt(this->_mesh->getSieve()->getMaxConeSize(),dim+1), true);
284 
285   for (PetscInt c = cStart; c < cEnd; ++c) {
286     ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*this->_mesh->getSieve(), c, pV);
287     const PETSC_MESH_TYPE::point_type *oPoints = pV.getPoints();
288     const int                          oSize   = pV.getSize();
289     int                                v       = 0;
290 
291     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr);
292     for (PetscInt cl = 0; cl < oSize; ++cl) {
293       const PetscInt fDim;
294 
295       ierr = PetscSectionGetDof(oPoints[cl], &fDim);CHKERRQ(ierr);
296       if (pointDim) {
297         for (PetscInt d = 0; d < fDim; ++d, ++v) {
298           values[v] = (*this->_options.integrate)(v0, J, v, initFunc);
299         }
300       }
301     }
302     ierr = DMPlexVecSetClosure(dm, NULL, localX, c, values);CHKERRQ(ierr);
303     pV.clear();
304   }
305   ierr = PetscFree2(v0,J);CHKERRQ(ierr);
306   ierr = PetscFree(values);CHKERRQ(ierr);
307 #endif
308   PetscFunctionReturn(0);
309 }
310 
311 #undef __FUNCT__
312 #define __FUNCT__ "DMPlexProjectFunction"
313 /*@C
314   DMPlexProjectFunction - This projects the given function into the function space provided.
315 
316   Input Parameters:
317 + dm      - The DM
318 . numComp - The number of components (functions)
319 . funcs   - The coordinate functions to evaluate
320 - mode    - The insertion mode for values
321 
322   Output Parameter:
323 . X - vector
324 
325   Level: developer
326 
327   Note:
328   This currently just calls the function with the coordinates of each vertex and edge midpoint, and stores the result in a vector.
329   We will eventually fix it.
330 
331 .seealso: DMPlexComputeL2Diff()
332 @*/
333 PetscErrorCode DMPlexProjectFunction(DM dm, PetscInt numComp, void (**funcs)(const PetscReal [], PetscScalar *), InsertMode mode, Vec X)
334 {
335   Vec            localX;
336   PetscErrorCode ierr;
337 
338   PetscFunctionBegin;
339   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
340   ierr = DMPlexProjectFunctionLocal(dm, numComp, funcs, mode, localX);CHKERRQ(ierr);
341   ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr);
342   ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr);
343   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
344   PetscFunctionReturn(0);
345 }
346 
347 #undef __FUNCT__
348 #define __FUNCT__ "DMPlexComputeL2Diff"
349 /*@C
350   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
351 
352   Input Parameters:
353 + dm    - The DM
354 . quad  - The PetscQuadrature object for each field
355 . funcs - The functions to evaluate for each field component
356 - X     - The coefficient vector u_h
357 
358   Output Parameter:
359 . diff - The diff ||u - u_h||_2
360 
361   Level: developer
362 
363 .seealso: DMPlexProjectFunction()
364 @*/
365 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscQuadrature quad[], void (**funcs)(const PetscReal [], PetscScalar *), Vec X, PetscReal *diff)
366 {
367   const PetscInt debug = 0;
368   PetscSection   section;
369   Vec            localX;
370   PetscReal     *coords, *v0, *J, *invJ, detJ;
371   PetscReal      localDiff = 0.0;
372   PetscInt       dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
373   PetscErrorCode ierr;
374 
375   PetscFunctionBegin;
376   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
377   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
378   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
379   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
380   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
381   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
382   for (field = 0; field < numFields; ++field) {
383     numComponents += quad[field].numComponents;
384   }
385   ierr = DMPlexProjectFunctionLocal(dm, numComponents, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
386   ierr = PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);CHKERRQ(ierr);
387   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
388   for (c = cStart; c < cEnd; ++c) {
389     PetscScalar *x;
390     PetscReal    elemDiff = 0.0;
391 
392     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
393     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
394     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
395 
396     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
397       const PetscInt   numQuadPoints = quad[field].numQuadPoints;
398       const PetscReal *quadPoints    = quad[field].quadPoints;
399       const PetscReal *quadWeights   = quad[field].quadWeights;
400       const PetscInt   numBasisFuncs = quad[field].numBasisFuncs;
401       const PetscInt   numBasisComps = quad[field].numComponents;
402       const PetscReal *basis         = quad[field].basis;
403       PetscInt         q, d, e, fc, f;
404 
405       if (debug) {
406         char title[1024];
407         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
408         ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr);
409       }
410       for (q = 0; q < numQuadPoints; ++q) {
411         for (d = 0; d < dim; d++) {
412           coords[d] = v0[d];
413           for (e = 0; e < dim; e++) {
414             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
415           }
416         }
417         for (fc = 0; fc < numBasisComps; ++fc) {
418           PetscScalar funcVal;
419           PetscScalar interpolant = 0.0;
420 
421           (*funcs[comp+fc])(coords, &funcVal);
422           for (f = 0; f < numBasisFuncs; ++f) {
423             const PetscInt fidx = f*numBasisComps+fc;
424             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
425           }
426           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);}
427           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal))*quadWeights[q]*detJ;
428         }
429       }
430       comp        += numBasisComps;
431       fieldOffset += numBasisFuncs*numBasisComps;
432     }
433     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
434     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
435     localDiff += elemDiff;
436   }
437   ierr  = PetscFree4(coords,v0,J,invJ);CHKERRQ(ierr);
438   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
439   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
440   *diff = PetscSqrtReal(*diff);
441   PetscFunctionReturn(0);
442 }
443 
444 #undef __FUNCT__
445 #define __FUNCT__ "DMPlexComputeResidualFEM"
446 /*@
447   DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
448 
449   Input Parameters:
450 + dm - The mesh
451 . X  - Local input vector
452 - user - The user context
453 
454   Output Parameter:
455 . F  - Local output vector
456 
457   Note:
458   The second member of the user context must be an FEMContext.
459 
460   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
461   like a GPU, or vectorize on a multicore machine.
462 
463   Level: developer
464 
465 .seealso: DMPlexComputeJacobianActionFEM()
466 @*/
467 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
468 {
469   DM_Plex         *mesh   = (DM_Plex*) dm->data;
470   PetscFEM        *fem    = (PetscFEM*) &((DM*) user)[1];
471   PetscQuadrature *quad   = fem->quad;
472   PetscQuadrature *quadBd = fem->quadBd;
473   PetscSection     section;
474   PetscReal       *v0, *n, *J, *invJ, *detJ;
475   PetscScalar     *elemVec, *u;
476   PetscInt         dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
477   PetscInt         cellDof, numComponents;
478   PetscBool        has;
479   PetscErrorCode   ierr;
480 
481   PetscFunctionBegin;
482   /* ierr = PetscLogEventBegin(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */
483   ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
484   ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
485   ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
486   ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
487   numCells = cEnd - cStart;
488   for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) {
489     cellDof       += quad[field].numBasisFuncs*quad[field].numComponents;
490     numComponents += quad[field].numComponents;
491   }
492   ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
493   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
494   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);
495   for (c = cStart; c < cEnd; ++c) {
496     PetscScalar *x;
497     PetscInt     i;
498 
499     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
500     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
501     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
502 
503     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
504     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
505   }
506   for (field = 0; field < numFields; ++field) {
507     const PetscInt numQuadPoints = quad[field].numQuadPoints;
508     const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
509     void           (*f0)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[field];
510     void           (*f1)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[field];
511     /* Conforming batches */
512     PetscInt blockSize  = numBasisFuncs*numQuadPoints;
513     PetscInt numBlocks  = 1;
514     PetscInt batchSize  = numBlocks * blockSize;
515     PetscInt numBatches = numBatchesTmp;
516     PetscInt numChunks  = numCells / (numBatches*batchSize);
517     /* Remainder */
518     PetscInt numRemainder = numCells % (numBatches * batchSize);
519     PetscInt offset       = numCells - numRemainder;
520 
521     ierr = (*mesh->integrateResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, v0, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr);
522     ierr = (*mesh->integrateResidualFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
523                                          f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
524   }
525   for (c = cStart; c < cEnd; ++c) {
526     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Residual", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
527     ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
528   }
529   ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
530   /* Integration over the boundary:
531      - This can probably be generalized to integration over a set of labels, however
532        the idea here is to do integration where we need the cell normal
533      - We can replace hardcoding with a registration process, and this is how we hook
534        up the system to something like FEniCS
535   */
536   ierr = DMPlexHasLabel(dm, "boundary", &has);CHKERRQ(ierr);
537   if (has && quadBd) {
538     DMLabel         label;
539     IS              pointIS;
540     const PetscInt *points;
541     PetscInt        numPoints, p;
542 
543     ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
544     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
545     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
546     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
547     for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) {
548       cellDof       += quadBd[field].numBasisFuncs*quadBd[field].numComponents;
549       numComponents += quadBd[field].numComponents;
550     }
551     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);
552     for (p = 0; p < numPoints; ++p) {
553       const PetscInt point = points[p];
554       PetscScalar   *x;
555       PetscInt       i;
556 
557       /* TODO: Add normal determination here */
558       ierr = DMPlexComputeCellGeometry(dm, point, &v0[p*dim], &J[p*dim*dim], &invJ[p*dim*dim], &detJ[p]);CHKERRQ(ierr);
559       if (detJ[p] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[p], point);
560       ierr = DMPlexVecGetClosure(dm, NULL, X, point, NULL, &x);CHKERRQ(ierr);
561 
562       for (i = 0; i < cellDof; ++i) u[p*cellDof+i] = x[i];
563       ierr = DMPlexVecRestoreClosure(dm, NULL, X, point, NULL, &x);CHKERRQ(ierr);
564     }
565     for (field = 0; field < numFields; ++field) {
566       const PetscInt numQuadPoints = quadBd[field].numQuadPoints;
567       const PetscInt numBasisFuncs = quadBd[field].numBasisFuncs;
568       void           (*f0)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[field];
569       void           (*f1)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[field];
570       /* Conforming batches */
571       PetscInt blockSize  = numBasisFuncs*numQuadPoints;
572       PetscInt numBlocks  = 1;
573       PetscInt batchSize  = numBlocks * blockSize;
574       PetscInt numBatches = numBatchesTmp;
575       PetscInt numChunks  = numPoints / (numBatches*batchSize);
576       /* Remainder */
577       PetscInt numRemainder = numPoints % (numBatches * batchSize);
578       PetscInt offset       = numPoints - numRemainder;
579 
580       ierr = (*mesh->integrateBdResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quadBd, u, v0, n, J, invJ, detJ, f0, f1, elemVec);CHKERRQ(ierr);
581       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],
582                                              f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
583     }
584     for (p = 0; p < numPoints; ++p) {
585       const PetscInt point = points[p];
586 
587       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "Residual", cellDof, &elemVec[p*cellDof]);CHKERRQ(ierr);}
588       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[p*cellDof], ADD_VALUES);CHKERRQ(ierr);
589     }
590     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
591     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
592     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
593   }
594   if (mesh->printFEM) {
595     PetscMPIInt rank, numProcs;
596     PetscInt    p;
597 
598     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
599     ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
600     ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Residual:\n");CHKERRQ(ierr);
601     for (p = 0; p < numProcs; ++p) {
602       if (p == rank) {
603         Vec f;
604 
605         ierr = VecDuplicate(F, &f);CHKERRQ(ierr);
606         ierr = VecCopy(F, f);CHKERRQ(ierr);
607         ierr = VecChop(f, 1.0e-10);CHKERRQ(ierr);
608         ierr = VecView(f, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
609         ierr = VecDestroy(&f);CHKERRQ(ierr);
610         ierr = PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
611       }
612       ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
613     }
614   }
615   /* ierr = PetscLogEventEnd(ResidualFEMEvent,0,0,0,0);CHKERRQ(ierr); */
616   PetscFunctionReturn(0);
617 }
618 
619 #undef __FUNCT__
620 #define __FUNCT__ "DMPlexComputeJacobianActionFEM"
621 /*@C
622   DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user
623 
624   Input Parameters:
625 + dm - The mesh
626 . J  - The Jacobian shell matrix
627 . X  - Local input vector
628 - user - The user context
629 
630   Output Parameter:
631 . F  - Local output vector
632 
633   Note:
634   The second member of the user context must be an FEMContext.
635 
636   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
637   like a GPU, or vectorize on a multicore machine.
638 
639   Level: developer
640 
641 .seealso: DMPlexComputeResidualFEM()
642 @*/
643 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
644 {
645   DM_Plex         *mesh = (DM_Plex*) dm->data;
646   PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
647   PetscQuadrature *quad = fem->quad;
648   PetscSection     section;
649   JacActionCtx    *jctx;
650   PetscReal       *v0, *J, *invJ, *detJ;
651   PetscScalar     *elemVec, *u, *a;
652   PetscInt         dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
653   PetscInt         cellDof = 0;
654   PetscErrorCode   ierr;
655 
656   PetscFunctionBegin;
657   /* ierr = PetscLogEventBegin(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */
658   ierr     = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
659   ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
660   ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
661   ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
662   ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
663   numCells = cEnd - cStart;
664   for (field = 0; field < numFields; ++field) {
665     cellDof += quad[field].numBasisFuncs*quad[field].numComponents;
666   }
667   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
668   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);
669   for (c = cStart; c < cEnd; ++c) {
670     PetscScalar *x;
671     PetscInt     i;
672 
673     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
674     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
675     ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
676     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
677     ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
678     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
679     for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
680     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
681   }
682   for (field = 0; field < numFields; ++field) {
683     const PetscInt numQuadPoints = quad[field].numQuadPoints;
684     const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
685     /* Conforming batches */
686     PetscInt blockSize  = numBasisFuncs*numQuadPoints;
687     PetscInt numBlocks  = 1;
688     PetscInt batchSize  = numBlocks * blockSize;
689     PetscInt numBatches = numBatchesTmp;
690     PetscInt numChunks  = numCells / (numBatches*batchSize);
691     /* Remainder */
692     PetscInt numRemainder = numCells % (numBatches * batchSize);
693     PetscInt offset       = numCells - numRemainder;
694 
695     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);
696     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],
697                                                fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr);
698   }
699   for (c = cStart; c < cEnd; ++c) {
700     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
701     ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
702   }
703   ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
704   if (mesh->printFEM) {
705     PetscMPIInt rank, numProcs;
706     PetscInt    p;
707 
708     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
709     ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
710     ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");CHKERRQ(ierr);
711     for (p = 0; p < numProcs; ++p) {
712       if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}
713       ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
714     }
715   }
716   /* ierr = PetscLogEventEnd(JacobianActionFEMEvent,0,0,0,0);CHKERRQ(ierr); */
717   PetscFunctionReturn(0);
718 }
719 
720 #undef __FUNCT__
721 #define __FUNCT__ "DMPlexComputeJacobianFEM"
722 /*@
723   DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
724 
725   Input Parameters:
726 + dm - The mesh
727 . X  - Local input vector
728 - user - The user context
729 
730   Output Parameter:
731 . Jac  - Jacobian matrix
732 
733   Note:
734   The second member of the user context must be an FEMContext.
735 
736   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
737   like a GPU, or vectorize on a multicore machine.
738 
739   Level: developer
740 
741 .seealso: FormFunctionLocal()
742 @*/
743 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, MatStructure *str,void *user)
744 {
745   DM_Plex         *mesh = (DM_Plex*) dm->data;
746   PetscFEM        *fem  = (PetscFEM*) &((DM*) user)[1];
747   PetscQuadrature *quad = fem->quad;
748   PetscSection     section;
749   PetscReal       *v0, *J, *invJ, *detJ;
750   PetscScalar     *elemMat, *u;
751   PetscInt         dim, numFields, field, fieldI, numBatchesTmp = 1, numCells, cStart, cEnd, c;
752   PetscInt         cellDof = 0, numComponents = 0;
753   PetscBool        isShell;
754   PetscErrorCode   ierr;
755 
756   PetscFunctionBegin;
757   /* ierr = PetscLogEventBegin(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */
758   ierr     = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
759   ierr     = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
760   ierr     = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
761   ierr     = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
762   numCells = cEnd - cStart;
763   for (field = 0; field < numFields; ++field) {
764     cellDof       += quad[field].numBasisFuncs*quad[field].numComponents;
765     numComponents += quad[field].numComponents;
766   }
767   ierr = DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
768   ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
769   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);
770   for (c = cStart; c < cEnd; ++c) {
771     PetscScalar *x;
772     PetscInt     i;
773 
774     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
775     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
776     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
777 
778     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
779     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
780   }
781   ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr);
782   for (fieldI = 0; fieldI < numFields; ++fieldI) {
783     const PetscInt numQuadPoints = quad[fieldI].numQuadPoints;
784     const PetscInt numBasisFuncs = quad[fieldI].numBasisFuncs;
785     PetscInt       fieldJ;
786 
787     for (fieldJ = 0; fieldJ < numFields; ++fieldJ) {
788       void (*g0)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*numFields+fieldJ];
789       void (*g1)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*numFields+fieldJ];
790       void (*g2)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*numFields+fieldJ];
791       void (*g3)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*numFields+fieldJ];
792       /* Conforming batches */
793       PetscInt blockSize  = numBasisFuncs*numQuadPoints;
794       PetscInt numBlocks  = 1;
795       PetscInt batchSize  = numBlocks * blockSize;
796       PetscInt numBatches = numBatchesTmp;
797       PetscInt numChunks  = numCells / (numBatches*batchSize);
798       /* Remainder */
799       PetscInt numRemainder = numCells % (numBatches * batchSize);
800       PetscInt offset       = numCells - numRemainder;
801 
802       ierr = (*mesh->integrateJacobianFEM)(numChunks*numBatches*batchSize, numFields, fieldI, fieldJ, quad, u, v0, J, invJ, detJ, g0, g1, g2, g3, elemMat);CHKERRQ(ierr);
803       ierr = (*mesh->integrateJacobianFEM)(numRemainder, numFields, fieldI, fieldJ, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
804                                            g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);CHKERRQ(ierr);
805     }
806   }
807   for (c = cStart; c < cEnd; ++c) {
808     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, "Jacobian", cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);}
809     ierr = DMPlexMatSetClosure(dm, NULL, NULL, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr);
810   }
811   ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
812 
813   /* Assemble matrix, using the 2-step process:
814        MatAssemblyBegin(), MatAssemblyEnd(). */
815   ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
816   ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
817 
818   if (mesh->printFEM) {
819     ierr = PetscPrintf(PETSC_COMM_WORLD, "Jacobian:\n");CHKERRQ(ierr);
820     ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
821     ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
822   }
823   /* ierr = PetscLogEventEnd(JacobianFEMEvent,0,0,0,0);CHKERRQ(ierr); */
824   ierr = PetscObjectTypeCompare((PetscObject)Jac, MATSHELL, &isShell);CHKERRQ(ierr);
825   if (isShell) {
826     JacActionCtx *jctx;
827 
828     ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
829     ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
830   }
831   *str = SAME_NONZERO_PATTERN;
832   PetscFunctionReturn(0);
833 }
834