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