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