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