xref: /petsc/src/dm/impls/plex/plexfem.c (revision 1c41a8ca76d8605dc7d579627879ffaa2f37c0cf) 
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__ "DMPlexProjectFunctionLabelLocal"
185 PetscErrorCode DMPlexProjectFunctionLabelLocal(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, 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, f, d, v, i, comp;
192   PetscErrorCode  ierr;
193 
194   PetscFunctionBegin;
195   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
196   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
197   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
198   ierr = PetscMalloc3(numFields,&sp,dim,&v0,dim*dim,&J);CHKERRQ(ierr);
199   for (f = 0; f < numFields; ++f) {
200     ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr);
201     ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr);
202     ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
203     totDim += spDim*numComp;
204   }
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   for (i = 0; i < numIds; ++i) {
210     IS              pointIS;
211     const PetscInt *points;
212     PetscInt        n, p;
213 
214     ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr);
215     ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr);
216     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
217     for (p = 0; p < n; ++p) {
218       const PetscInt    point = points[p];
219       PetscCellGeometry geom;
220 
221       if ((point < cStart) || (point >= cEnd)) continue;
222       ierr = DMPlexComputeCellGeometry(dm, point, v0, J, NULL, &detJ);CHKERRQ(ierr);
223       geom.v0   = v0;
224       geom.J    = J;
225       geom.detJ = &detJ;
226       for (f = 0, v = 0; f < numFields; ++f) {
227         void * const ctx = ctxs ? ctxs[f] : NULL;
228         ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr);
229         ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
230         for (d = 0; d < spDim; ++d) {
231           if (funcs[f]) {
232             ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr);
233           } else {
234             for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0;
235           }
236           v += numComp;
237         }
238       }
239       ierr = DMPlexVecSetClosure(dm, section, localX, point, values, mode);CHKERRQ(ierr);
240     }
241     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
242     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
243   }
244   ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
245   ierr = PetscFree3(sp,v0,J);CHKERRQ(ierr);
246   PetscFunctionReturn(0);
247 }
248 
249 #undef __FUNCT__
250 #define __FUNCT__ "DMPlexProjectFunctionLocal"
251 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
252 {
253   PetscDualSpace *sp;
254   PetscSection    section;
255   PetscScalar    *values;
256   PetscReal      *v0, *J, detJ;
257   PetscInt        numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, c, f, d, v, comp;
258   PetscErrorCode  ierr;
259 
260   PetscFunctionBegin;
261   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
262   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
263   ierr = PetscMalloc1(numFields, &sp);CHKERRQ(ierr);
264   for (f = 0; f < numFields; ++f) {
265     ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr);
266     ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr);
267     ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
268     totDim += spDim*numComp;
269   }
270   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
271   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
272   ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr);
273   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
274   ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
275   ierr = PetscMalloc2(dim,&v0,dim*dim,&J);CHKERRQ(ierr);
276   for (c = cStart; c < cEnd; ++c) {
277     PetscCellGeometry geom;
278 
279     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr);
280     geom.v0   = v0;
281     geom.J    = J;
282     geom.detJ = &detJ;
283     for (f = 0, v = 0; f < numFields; ++f) {
284       void * const ctx = ctxs ? ctxs[f] : NULL;
285       ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr);
286       ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
287       for (d = 0; d < spDim; ++d) {
288         if (funcs[f]) {
289           ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr);
290         } else {
291           for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0;
292         }
293         v += numComp;
294       }
295     }
296     ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr);
297   }
298   ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
299   ierr = PetscFree2(v0,J);CHKERRQ(ierr);
300   ierr = PetscFree(sp);CHKERRQ(ierr);
301   PetscFunctionReturn(0);
302 }
303 
304 #undef __FUNCT__
305 #define __FUNCT__ "DMPlexProjectFunction"
306 /*@C
307   DMPlexProjectFunction - This projects the given function into the function space provided.
308 
309   Input Parameters:
310 + dm      - The DM
311 . fe      - The PetscFE associated with the field
312 . funcs   - The coordinate functions to evaluate, one per field
313 . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
314 - mode    - The insertion mode for values
315 
316   Output Parameter:
317 . X - vector
318 
319   Level: developer
320 
321 .seealso: DMPlexComputeL2Diff()
322 @*/
323 PetscErrorCode DMPlexProjectFunction(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
324 {
325   Vec            localX;
326   PetscErrorCode ierr;
327 
328   PetscFunctionBegin;
329   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
330   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
331   ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, mode, localX);CHKERRQ(ierr);
332   ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr);
333   ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr);
334   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
335   PetscFunctionReturn(0);
336 }
337 
338 
339 #undef __FUNCT__
340 #define __FUNCT__ "DMPlexInsertBoundaryValues"
341 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec localX)
342 {
343   void        (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx);
344   void         **ctxs;
345   PetscFE       *fe;
346   PetscInt       numFields, f, numBd, b;
347   PetscErrorCode ierr;
348 
349   PetscFunctionBegin;
350   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
351   PetscValidHeaderSpecific(localX, VEC_CLASSID, 2);
352   ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr);
353   ierr = PetscMalloc3(numFields,&fe,numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr);
354   for (f = 0; f < numFields; ++f) {ierr = DMGetField(dm, f, (PetscObject *) &fe[f]);CHKERRQ(ierr);}
355   /* OPT: Could attempt to do multiple BCs at once */
356   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
357   for (b = 0; b < numBd; ++b) {
358     DMLabel         label;
359     const PetscInt *ids;
360     const char     *name;
361     PetscInt        numids, field;
362     PetscBool       isEssential;
363     void          (*func)();
364     void           *ctx;
365 
366     /* TODO: We need to set only the part indicated by the ids */
367     ierr = DMPlexGetBoundary(dm, b, &isEssential, &name, &field, &func, &numids, &ids, &ctx);CHKERRQ(ierr);
368     ierr = DMPlexGetLabel(dm, name, &label);CHKERRQ(ierr);
369     for (f = 0; f < numFields; ++f) {
370       funcs[f] = field == f ? (void (*)(const PetscReal[], PetscScalar *, void *)) func : NULL;
371       ctxs[f]  = field == f ? ctx : NULL;
372     }
373     ierr = DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
374   }
375   ierr = PetscFree3(fe,funcs,ctxs);CHKERRQ(ierr);
376   PetscFunctionReturn(0);
377 }
378 
379 #undef __FUNCT__
380 #define __FUNCT__ "DMPlexComputeL2Diff"
381 /*@C
382   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
383 
384   Input Parameters:
385 + dm    - The DM
386 . fe    - The PetscFE object for each field
387 . funcs - The functions to evaluate for each field component
388 . ctxs  - Optional array of contexts to pass to each function, or NULL.
389 - X     - The coefficient vector u_h
390 
391   Output Parameter:
392 . diff - The diff ||u - u_h||_2
393 
394   Level: developer
395 
396 .seealso: DMPlexProjectFunction(), DMPlexComputeL2GradientDiff()
397 @*/
398 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
399 {
400   const PetscInt  debug = 0;
401   PetscSection    section;
402   PetscQuadrature quad;
403   Vec             localX;
404   PetscScalar    *funcVal;
405   PetscReal      *coords, *v0, *J, *invJ, detJ;
406   PetscReal       localDiff = 0.0;
407   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
408   PetscErrorCode  ierr;
409 
410   PetscFunctionBegin;
411   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
412   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
413   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
414   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
415   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
416   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
417   for (field = 0; field < numFields; ++field) {
418     PetscInt Nc;
419 
420     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
421     numComponents += Nc;
422   }
423   ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
424   ierr = PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
425   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
426   ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr);
427   for (c = cStart; c < cEnd; ++c) {
428     PetscScalar *x = NULL;
429     PetscReal    elemDiff = 0.0;
430 
431     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
432     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
433     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
434 
435     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
436       void * const     ctx = ctxs ? ctxs[field] : NULL;
437       const PetscReal *quadPoints, *quadWeights;
438       PetscReal       *basis;
439       PetscInt         numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;
440 
441       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
442       ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr);
443       ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr);
444       ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr);
445       if (debug) {
446         char title[1024];
447         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
448         ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr);
449       }
450       for (q = 0; q < numQuadPoints; ++q) {
451         for (d = 0; d < dim; d++) {
452           coords[d] = v0[d];
453           for (e = 0; e < dim; e++) {
454             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
455           }
456         }
457         (*funcs[field])(coords, funcVal, ctx);
458         for (fc = 0; fc < numBasisComps; ++fc) {
459           PetscScalar interpolant = 0.0;
460 
461           for (f = 0; f < numBasisFuncs; ++f) {
462             const PetscInt fidx = f*numBasisComps+fc;
463             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
464           }
465           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);}
466           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
467         }
468       }
469       comp        += numBasisComps;
470       fieldOffset += numBasisFuncs*numBasisComps;
471     }
472     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
473     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
474     localDiff += elemDiff;
475   }
476   ierr  = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr);
477   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
478   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
479   *diff = PetscSqrtReal(*diff);
480   PetscFunctionReturn(0);
481 }
482 
483 #undef __FUNCT__
484 #define __FUNCT__ "DMPlexComputeL2GradientDiff"
485 /*@C
486   DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
487 
488   Input Parameters:
489 + dm    - The DM
490 . fe    - The PetscFE object for each field
491 . funcs - The gradient functions to evaluate for each field component
492 . ctxs  - Optional array of contexts to pass to each function, or NULL.
493 . X     - The coefficient vector u_h
494 - n     - The vector to project along
495 
496   Output Parameter:
497 . diff - The diff ||(grad u - grad u_h) . n||_2
498 
499   Level: developer
500 
501 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
502 @*/
503 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
504 {
505   const PetscInt  debug = 0;
506   PetscSection    section;
507   PetscQuadrature quad;
508   Vec             localX;
509   PetscScalar    *funcVal, *interpolantVec;
510   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
511   PetscReal       localDiff = 0.0;
512   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
513   PetscErrorCode  ierr;
514 
515   PetscFunctionBegin;
516   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
517   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
518   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
519   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
520   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
521   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
522   for (field = 0; field < numFields; ++field) {
523     PetscInt Nc;
524 
525     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
526     numComponents += Nc;
527   }
528   /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */
529   ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr);
530   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
531   ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr);
532   for (c = cStart; c < cEnd; ++c) {
533     PetscScalar *x = NULL;
534     PetscReal    elemDiff = 0.0;
535 
536     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
537     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
538     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
539 
540     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
541       void * const     ctx = ctxs ? ctxs[field] : NULL;
542       const PetscReal *quadPoints, *quadWeights;
543       PetscReal       *basisDer;
544       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;
545 
546       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
547       ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
548       ierr = PetscFEGetNumComponents(fe[field], &Ncomp);CHKERRQ(ierr);
549       ierr = PetscFEGetDefaultTabulation(fe[field], NULL, &basisDer, NULL);CHKERRQ(ierr);
550       if (debug) {
551         char title[1024];
552         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
553         ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr);
554       }
555       for (q = 0; q < numQuadPoints; ++q) {
556         for (d = 0; d < dim; d++) {
557           coords[d] = v0[d];
558           for (e = 0; e < dim; e++) {
559             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
560           }
561         }
562         (*funcs[field])(coords, n, funcVal, ctx);
563         for (fc = 0; fc < Ncomp; ++fc) {
564           PetscScalar interpolant = 0.0;
565 
566           for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
567           for (f = 0; f < Nb; ++f) {
568             const PetscInt fidx = f*Ncomp+fc;
569 
570             for (d = 0; d < dim; ++d) {
571               realSpaceDer[d] = 0.0;
572               for (g = 0; g < dim; ++g) {
573                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
574               }
575               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
576             }
577           }
578           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
579           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);}
580           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
581         }
582       }
583       comp        += Ncomp;
584       fieldOffset += Nb*Ncomp;
585     }
586     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
587     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
588     localDiff += elemDiff;
589   }
590   ierr  = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr);
591   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
592   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
593   *diff = PetscSqrtReal(*diff);
594   PetscFunctionReturn(0);
595 }
596 
597 #undef __FUNCT__
598 #define __FUNCT__ "DMPlexComputeResidualFEM"
599 /*@
600   DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
601 
602   Input Parameters:
603 + dm - The mesh
604 . X  - Local input vector
605 - user - The user context
606 
607   Output Parameter:
608 . F  - Local output vector
609 
610   Note:
611   The first member of the user context must be an FEMContext.
612 
613   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
614   like a GPU, or vectorize on a multicore machine.
615 
616   Level: developer
617 
618 .seealso: DMPlexComputeJacobianActionFEM()
619 @*/
620 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
621 {
622   DM_Plex          *mesh  = (DM_Plex *) dm->data;
623   PetscFEM         *fem   = (PetscFEM *) user;
624   PetscFE          *fe    = fem->fe;
625   PetscFE          *feAux = fem->feAux;
626   PetscFE          *feBd  = fem->feBd;
627   const char       *name  = "Residual";
628   DM                dmAux;
629   Vec               A;
630   PetscQuadrature   q;
631   PetscCellGeometry geom;
632   PetscSection      section, sectionAux;
633   PetscReal        *v0, *J, *invJ, *detJ;
634   PetscScalar      *elemVec, *u, *a = NULL;
635   PetscInt          dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c;
636   PetscInt          cellDof = 0, numComponents = 0;
637   PetscInt          cellDofAux = 0, numComponentsAux = 0;
638   PetscErrorCode    ierr;
639 
640   PetscFunctionBegin;
641   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
642   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
643   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
644   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
645   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
646   numCells = cEnd - cStart;
647   for (f = 0; f < Nf; ++f) {
648     PetscInt Nb, Nc;
649 
650     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
651     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
652     cellDof       += Nb*Nc;
653     numComponents += Nc;
654   }
655   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
656   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
657   if (dmAux) {
658     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
659     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
660   }
661   for (f = 0; f < NfAux; ++f) {
662     PetscInt Nb, Nc;
663 
664     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
665     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
666     cellDofAux       += Nb*Nc;
667     numComponentsAux += Nc;
668   }
669   ierr = DMPlexInsertBoundaryValues(dm, X);CHKERRQ(ierr);
670   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
671   ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr);
672   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
673   for (c = cStart; c < cEnd; ++c) {
674     PetscScalar *x = NULL;
675     PetscInt     i;
676 
677     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
678     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
679     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
680     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
681     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
682     if (dmAux) {
683       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
684       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
685       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
686     }
687   }
688   for (f = 0; f < Nf; ++f) {
689     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[f];
690     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[f];
691     PetscInt numQuadPoints, Nb;
692     /* Conforming batches */
693     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
694     /* Remainder */
695     PetscInt Nr, offset;
696 
697     ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr);
698     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
699     ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
700     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
701     blockSize = Nb*numQuadPoints;
702     batchSize = numBlocks * blockSize;
703     ierr =  PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
704     numChunks = numCells / (numBatches*batchSize);
705     Ne        = numChunks*numBatches*batchSize;
706     Nr        = numCells % (numBatches*batchSize);
707     offset    = numCells - Nr;
708     geom.v0   = v0;
709     geom.J    = J;
710     geom.invJ = invJ;
711     geom.detJ = detJ;
712     ierr = PetscFEIntegrateResidual(fe[f], Ne, Nf, fe, f, geom, u, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr);
713     geom.v0   = &v0[offset*dim];
714     geom.J    = &J[offset*dim*dim];
715     geom.invJ = &invJ[offset*dim*dim];
716     geom.detJ = &detJ[offset];
717     ierr = PetscFEIntegrateResidual(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
718   }
719   for (c = cStart; c < cEnd; ++c) {
720     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
721     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
722   }
723   ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
724   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
725   if (feBd) {
726     DMLabel         label, depth;
727     IS              pointIS;
728     const PetscInt *points;
729     PetscInt        dep, numPoints, p, numFaces;
730     PetscReal      *n;
731 
732     ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
733     ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
734     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
735     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
736     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
737     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
738       PetscInt Nb, Nc;
739 
740       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
741       ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr);
742       cellDof       += Nb*Nc;
743       numComponents += Nc;
744     }
745     for (p = 0, numFaces = 0; p < numPoints; ++p) {
746       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
747       if (dep == dim-1) ++numFaces;
748     }
749     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);
750     for (p = 0, f = 0; p < numPoints; ++p) {
751       const PetscInt point = points[p];
752       PetscScalar   *x     = NULL;
753       PetscInt       i;
754 
755       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
756       if (dep != dim-1) continue;
757       ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
758       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
759       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
760       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
761       for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
762       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
763       ++f;
764     }
765     for (f = 0; f < Nf; ++f) {
766       void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[f];
767       void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[f];
768       PetscInt numQuadPoints, Nb;
769       /* Conforming batches */
770       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
771       /* Remainder */
772       PetscInt Nr, offset;
773 
774       ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr);
775       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
776       ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
777       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
778       blockSize = Nb*numQuadPoints;
779       batchSize = numBlocks * blockSize;
780       ierr =  PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
781       numChunks = numFaces / (numBatches*batchSize);
782       Ne        = numChunks*numBatches*batchSize;
783       Nr        = numFaces % (numBatches*batchSize);
784       offset    = numFaces - Nr;
785       geom.v0   = v0;
786       geom.n    = n;
787       geom.J    = J;
788       geom.invJ = invJ;
789       geom.detJ = detJ;
790       ierr = PetscFEIntegrateBdResidual(feBd[f], Ne, Nf, feBd, f, geom, u, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr);
791       geom.v0   = &v0[offset*dim];
792       geom.n    = &n[offset*dim];
793       geom.J    = &J[offset*dim*dim];
794       geom.invJ = &invJ[offset*dim*dim];
795       geom.detJ = &detJ[offset];
796       ierr = PetscFEIntegrateBdResidual(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
797     }
798     for (p = 0, f = 0; p < numPoints; ++p) {
799       const PetscInt point = points[p];
800 
801       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
802       if (dep != dim-1) continue;
803       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);}
804       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr);
805       ++f;
806     }
807     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
808     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
809     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
810   }
811   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
812   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
813   PetscFunctionReturn(0);
814 }
815 
816 #undef __FUNCT__
817 #define __FUNCT__ "DMPlexComputeIFunctionFEM"
818 /*@
819   DMPlexComputeIFunctionFEM - Form the local implicit function F from the local input X, X_t using pointwise functions specified by the user
820 
821   Input Parameters:
822 + dm - The mesh
823 . time - The current time
824 . X  - Local input vector
825 . X_t  - Time derivative of the local input vector
826 - user - The user context
827 
828   Output Parameter:
829 . F  - Local output vector
830 
831   Note:
832   The first 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: DMPlexComputeResidualFEM()
840 @*/
841 PetscErrorCode DMPlexComputeIFunctionFEM(DM dm, PetscReal time, Vec X, Vec X_t, Vec F, 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   PetscFE          *feBd  = fem->feBd;
848   const char       *name  = "Residual";
849   DM                dmAux;
850   Vec               A;
851   PetscQuadrature   q;
852   PetscCellGeometry geom;
853   PetscSection      section, sectionAux;
854   PetscReal        *v0, *J, *invJ, *detJ;
855   PetscScalar      *elemVec, *u, *u_t, *a = NULL;
856   PetscInt          dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c;
857   PetscInt          cellDof = 0, numComponents = 0;
858   PetscInt          cellDofAux = 0, numComponentsAux = 0;
859   PetscErrorCode    ierr;
860 
861   PetscFunctionBegin;
862   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
863   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
864   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
865   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
866   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
867   numCells = cEnd - cStart;
868   for (f = 0; f < Nf; ++f) {
869     PetscInt Nb, Nc;
870 
871     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
872     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
873     cellDof       += Nb*Nc;
874     numComponents += Nc;
875   }
876   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
877   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
878   if (dmAux) {
879     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
880     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
881   }
882   for (f = 0; f < NfAux; ++f) {
883     PetscInt Nb, Nc;
884 
885     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
886     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
887     cellDofAux       += Nb*Nc;
888     numComponentsAux += Nc;
889   }
890   ierr = DMPlexInsertBoundaryValues(dm, X);CHKERRQ(ierr);
891   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
892   ierr = PetscMalloc7(numCells*cellDof,&u,numCells*cellDof,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr);
893   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
894   for (c = cStart; c < cEnd; ++c) {
895     PetscScalar *x = NULL, *x_t = NULL;
896     PetscInt     i;
897 
898     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
899     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
900     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
901     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
902     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
903     ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
904     for (i = 0; i < cellDof; ++i) u_t[c*cellDof+i] = x_t[i];
905     ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
906     if (dmAux) {
907       PetscScalar *x_a = NULL;
908       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr);
909       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x_a[i];
910       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr);
911     }
912   }
913   for (f = 0; f < Nf; ++f) {
914     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0IFuncs[f];
915     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1IFuncs[f];
916     PetscInt numQuadPoints, Nb;
917     /* Conforming batches */
918     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
919     /* Remainder */
920     PetscInt Nr, offset;
921 
922     ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr);
923     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
924     ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
925     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
926     blockSize = Nb*numQuadPoints;
927     batchSize = numBlocks * blockSize;
928     ierr =  PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
929     numChunks = numCells / (numBatches*batchSize);
930     Ne        = numChunks*numBatches*batchSize;
931     Nr        = numCells % (numBatches*batchSize);
932     offset    = numCells - Nr;
933     geom.v0   = v0;
934     geom.J    = J;
935     geom.invJ = invJ;
936     geom.detJ = detJ;
937     ierr = PetscFEIntegrateIFunction(fe[f], Ne, Nf, fe, f, geom, u, u_t, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr);
938     geom.v0   = &v0[offset*dim];
939     geom.J    = &J[offset*dim*dim];
940     geom.invJ = &invJ[offset*dim*dim];
941     geom.detJ = &detJ[offset];
942     ierr = PetscFEIntegrateIFunction(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], &u_t[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
943   }
944   for (c = cStart; c < cEnd; ++c) {
945     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
946     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
947   }
948   ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
949   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
950   if (feBd) {
951     DMLabel         label, depth;
952     IS              pointIS;
953     const PetscInt *points;
954     PetscInt        dep, numPoints, p, numFaces;
955     PetscReal      *n;
956 
957     ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
958     ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
959     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
960     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
961     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
962     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
963       PetscInt Nb, Nc;
964 
965       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
966       ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr);
967       cellDof       += Nb*Nc;
968       numComponents += Nc;
969     }
970     for (p = 0, numFaces = 0; p < numPoints; ++p) {
971       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
972       if (dep == dim-1) ++numFaces;
973     }
974     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);
975     ierr = PetscMalloc1(numFaces*cellDof,&u_t);CHKERRQ(ierr);
976     for (p = 0, f = 0; p < numPoints; ++p) {
977       const PetscInt point = points[p];
978       PetscScalar   *x     = NULL;
979       PetscInt       i;
980 
981       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
982       if (dep != dim-1) continue;
983       ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
984       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
985       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
986       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
987       for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
988       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
989       ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
990       for (i = 0; i < cellDof; ++i) u_t[f*cellDof+i] = x[i];
991       ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
992       ++f;
993     }
994     for (f = 0; f < Nf; ++f) {
995       void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdIFuncs[f];
996       void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdIFuncs[f];
997       PetscInt numQuadPoints, Nb;
998       /* Conforming batches */
999       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1000       /* Remainder */
1001       PetscInt Nr, offset;
1002 
1003       ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr);
1004       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
1005       ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1006       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1007       blockSize = Nb*numQuadPoints;
1008       batchSize = numBlocks * blockSize;
1009       ierr =  PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1010       numChunks = numFaces / (numBatches*batchSize);
1011       Ne        = numChunks*numBatches*batchSize;
1012       Nr        = numFaces % (numBatches*batchSize);
1013       offset    = numFaces - Nr;
1014       geom.v0   = v0;
1015       geom.n    = n;
1016       geom.J    = J;
1017       geom.invJ = invJ;
1018       geom.detJ = detJ;
1019       ierr = PetscFEIntegrateBdIFunction(feBd[f], Ne, Nf, feBd, f, geom, u, u_t, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr);
1020       geom.v0   = &v0[offset*dim];
1021       geom.n    = &n[offset*dim];
1022       geom.J    = &J[offset*dim*dim];
1023       geom.invJ = &invJ[offset*dim*dim];
1024       geom.detJ = &detJ[offset];
1025       ierr = PetscFEIntegrateBdIFunction(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], &u_t[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
1026     }
1027     for (p = 0, f = 0; p < numPoints; ++p) {
1028       const PetscInt point = points[p];
1029 
1030       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
1031       if (dep != dim-1) continue;
1032       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);}
1033       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr);
1034       ++f;
1035     }
1036     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
1037     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
1038     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
1039     ierr = PetscFree(u_t);CHKERRQ(ierr);
1040   }
1041   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
1042   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
1043   PetscFunctionReturn(0);
1044 }
1045 
1046 #undef __FUNCT__
1047 #define __FUNCT__ "DMPlexComputeJacobianActionFEM"
1048 /*@C
1049   DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user
1050 
1051   Input Parameters:
1052 + dm - The mesh
1053 . J  - The Jacobian shell matrix
1054 . X  - Local input vector
1055 - user - The user context
1056 
1057   Output Parameter:
1058 . F  - Local output vector
1059 
1060   Note:
1061   The first member of the user context must be an FEMContext.
1062 
1063   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1064   like a GPU, or vectorize on a multicore machine.
1065 
1066   Level: developer
1067 
1068 .seealso: DMPlexComputeResidualFEM()
1069 @*/
1070 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
1071 {
1072   DM_Plex          *mesh = (DM_Plex *) dm->data;
1073   PetscFEM         *fem  = (PetscFEM *) user;
1074   PetscFE          *fe   = fem->fe;
1075   PetscQuadrature   quad;
1076   PetscCellGeometry geom;
1077   PetscSection      section;
1078   JacActionCtx     *jctx;
1079   PetscReal        *v0, *J, *invJ, *detJ;
1080   PetscScalar      *elemVec, *u, *a;
1081   PetscInt          dim, numFields, field, numCells, cStart, cEnd, c;
1082   PetscInt          cellDof = 0;
1083   PetscErrorCode    ierr;
1084 
1085   PetscFunctionBegin;
1086   /* ierr = PetscLogEventBegin(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */
1087   ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
1088   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
1089   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1090   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
1091   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1092   numCells = cEnd - cStart;
1093   for (field = 0; field < numFields; ++field) {
1094     PetscInt Nb, Nc;
1095 
1096     ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
1097     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
1098     cellDof += Nb*Nc;
1099   }
1100   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
1101   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);
1102   for (c = cStart; c < cEnd; ++c) {
1103     PetscScalar *x = NULL;
1104     PetscInt     i;
1105 
1106     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1107     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1108     ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
1109     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
1110     ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
1111     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
1112     for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
1113     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
1114   }
1115   for (field = 0; field < numFields; ++field) {
1116     PetscInt numQuadPoints, Nb;
1117     /* Conforming batches */
1118     PetscInt numBlocks  = 1;
1119     PetscInt numBatches = 1;
1120     PetscInt numChunks, Ne, blockSize, batchSize;
1121     /* Remainder */
1122     PetscInt Nr, offset;
1123 
1124     ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr);
1125     ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
1126     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1127     blockSize = Nb*numQuadPoints;
1128     batchSize = numBlocks * blockSize;
1129     numChunks = numCells / (numBatches*batchSize);
1130     Ne        = numChunks*numBatches*batchSize;
1131     Nr        = numCells % (numBatches*batchSize);
1132     offset    = numCells - Nr;
1133     geom.v0   = v0;
1134     geom.J    = J;
1135     geom.invJ = invJ;
1136     geom.detJ = detJ;
1137     ierr = PetscFEIntegrateJacobianAction(fe[field], Ne, numFields, fe, field, geom, u, a, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr);
1138     geom.v0   = &v0[offset*dim];
1139     geom.J    = &J[offset*dim*dim];
1140     geom.invJ = &invJ[offset*dim*dim];
1141     geom.detJ = &detJ[offset];
1142     ierr = PetscFEIntegrateJacobianAction(fe[field], Nr, numFields, fe, field, geom, &u[offset*cellDof], &a[offset*cellDof],
1143                                           fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr);
1144   }
1145   for (c = cStart; c < cEnd; ++c) {
1146     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
1147     ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
1148   }
1149   ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
1150   if (mesh->printFEM) {
1151     PetscMPIInt rank, numProcs;
1152     PetscInt    p;
1153 
1154     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
1155     ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
1156     ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");CHKERRQ(ierr);
1157     for (p = 0; p < numProcs; ++p) {
1158       if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}
1159       ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
1160     }
1161   }
1162   /* ierr = PetscLogEventEnd(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */
1163   PetscFunctionReturn(0);
1164 }
1165 
1166 #undef __FUNCT__
1167 #define __FUNCT__ "DMPlexComputeJacobianFEM"
1168 /*@
1169   DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
1170 
1171   Input Parameters:
1172 + dm - The mesh
1173 . X  - Local input vector
1174 - user - The user context
1175 
1176   Output Parameter:
1177 . Jac  - Jacobian matrix
1178 
1179   Note:
1180   The first member of the user context must be an FEMContext.
1181 
1182   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1183   like a GPU, or vectorize on a multicore machine.
1184 
1185   Level: developer
1186 
1187 .seealso: FormFunctionLocal()
1188 @*/
1189 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1190 {
1191   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1192   PetscFEM         *fem   = (PetscFEM *) user;
1193   PetscFE          *fe    = fem->fe;
1194   PetscFE          *feAux = fem->feAux;
1195   const char       *name  = "Jacobian";
1196   DM                dmAux;
1197   Vec               A;
1198   PetscQuadrature   quad;
1199   PetscCellGeometry geom;
1200   PetscSection      section, globalSection, sectionAux;
1201   PetscReal        *v0, *J, *invJ, *detJ;
1202   PetscScalar      *elemMat, *u, *a;
1203   PetscInt          dim, Nf, NfAux = 0, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1204   PetscInt          cellDof = 0, numComponents = 0;
1205   PetscInt          cellDofAux = 0, numComponentsAux = 0;
1206   PetscBool         isShell;
1207   PetscErrorCode    ierr;
1208 
1209   PetscFunctionBegin;
1210   ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1211   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
1212   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1213   ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);
1214   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1215   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1216   numCells = cEnd - cStart;
1217   for (f = 0; f < Nf; ++f) {
1218     PetscInt Nb, Nc;
1219 
1220     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
1221     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1222     cellDof       += Nb*Nc;
1223     numComponents += Nc;
1224   }
1225   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1226   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1227   if (dmAux) {
1228     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1229     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
1230   }
1231   for (f = 0; f < NfAux; ++f) {
1232     PetscInt Nb, Nc;
1233 
1234     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
1235     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
1236     cellDofAux       += Nb*Nc;
1237     numComponentsAux += Nc;
1238   }
1239   ierr = DMPlexInsertBoundaryValues(dm, X);CHKERRQ(ierr);
1240   ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
1241   ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof*cellDof,&elemMat);CHKERRQ(ierr);
1242   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
1243   for (c = cStart; c < cEnd; ++c) {
1244     PetscScalar *x = NULL;
1245     PetscInt     i;
1246 
1247     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1248     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1249     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1250     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
1251     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1252     if (dmAux) {
1253       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1254       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
1255       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1256     }
1257   }
1258   ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1259   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1260     PetscInt numQuadPoints, Nb;
1261     /* Conforming batches */
1262     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1263     /* Remainder */
1264     PetscInt Nr, offset;
1265 
1266     ierr = PetscFEGetQuadrature(fe[fieldI], &quad);CHKERRQ(ierr);
1267     ierr = PetscFEGetDimension(fe[fieldI], &Nb);CHKERRQ(ierr);
1268     ierr = PetscFEGetTileSizes(fe[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1269     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1270     blockSize = Nb*numQuadPoints;
1271     batchSize = numBlocks * blockSize;
1272     ierr = PetscFESetTileSizes(fe[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1273     numChunks = numCells / (numBatches*batchSize);
1274     Ne        = numChunks*numBatches*batchSize;
1275     Nr        = numCells % (numBatches*batchSize);
1276     offset    = numCells - Nr;
1277     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1278       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ];
1279       void   (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*Nf+fieldJ];
1280       void   (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*Nf+fieldJ];
1281       void   (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*Nf+fieldJ];
1282 
1283       geom.v0   = v0;
1284       geom.J    = J;
1285       geom.invJ = invJ;
1286       geom.detJ = detJ;
1287       ierr = PetscFEIntegrateJacobian(fe[fieldI], Ne, Nf, fe, fieldI, fieldJ, geom, u, NfAux, feAux, a, g0, g1, g2, g3, elemMat);CHKERRQ(ierr);
1288       geom.v0   = &v0[offset*dim];
1289       geom.J    = &J[offset*dim*dim];
1290       geom.invJ = &invJ[offset*dim*dim];
1291       geom.detJ = &detJ[offset];
1292       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);
1293     }
1294   }
1295   for (c = cStart; c < cEnd; ++c) {
1296     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);}
1297     ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr);
1298   }
1299   ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1300   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1301   ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1302   ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1303   if (mesh->printFEM) {
1304     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1305     ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
1306     ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1307   }
1308   ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1309   ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr);
1310   if (isShell) {
1311     JacActionCtx *jctx;
1312 
1313     ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
1314     ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
1315   }
1316   PetscFunctionReturn(0);
1317 }
1318 
1319 #if 0
1320 
1321 static void g0_identity_1d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1322 {
1323   g3[0] = 1.0;
1324 }
1325 
1326 static void g0_identity_2d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1327 {
1328   g3[0] = g3[3] = 1.0;
1329 }
1330 
1331 static void g0_identity_3d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1332 {
1333   g3[0] = g3[4] = g3[8] = 1.0;
1334 }
1335 
1336 #undef __FUNCT__
1337 #define __FUNCT__ "DMPlexComputeInterpolatorFEMBroken"
1338 PetscErrorCode DMPlexComputeInterpolatorFEMBroken(DM dmc, DM dmf, Mat I, void *user)
1339 {
1340   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1341   PetscFEM         *fem   = (PetscFEM *) user;
1342   PetscFE          *fe    = fem->fe;
1343   const char       *name  = "Interpolator";
1344   PetscFE          *feRef;
1345   PetscQuadrature   quad, quadOld;
1346   PetscCellGeometry geom;
1347   PetscSection      fsection, fglobalSection;
1348   PetscSection      csection, cglobalSection;
1349   PetscReal        *v0, *J, *invJ, *detJ;
1350   PetscScalar      *elemMat;
1351   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1352   PetscInt          rCellDof = 0, cCellDof = 0, numComponents = 0;
1353   PetscErrorCode    ierr;
1354 
1355   PetscFunctionBegin;
1356 #if 0
1357   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1358 #endif
1359   ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr);
1360   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1361   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1362   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1363   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1364   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1365   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1366   numCells = cEnd - cStart;
1367   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1368   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1369     ierr = PetscFERefine(fe[fieldI], &feRef[fieldI]);CHKERRQ(ierr);
1370   }
1371   for (f = 0; f < Nf; ++f) {
1372     PetscInt rNb, cNb, Nc;
1373 
1374     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1375     ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr);
1376     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1377     numComponents += Nc;
1378     rCellDof += rNb*Nc;
1379     cCellDof += cNb*Nc;
1380   }
1381   ierr = MatZeroEntries(I);CHKERRQ(ierr);
1382   ierr = PetscMalloc5(numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*rCellDof*cCellDof,&elemMat);CHKERRQ(ierr);
1383   for (c = cStart; c < cEnd; ++c) {
1384     ierr = DMPlexComputeCellGeometry(dmc, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1385     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1386   }
1387   ierr = PetscMemzero(elemMat, numCells*rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1388   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1389     PetscInt numQuadPoints, Nb;
1390     /* Conforming batches */
1391     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1392     /* Remainder */
1393     PetscInt Nr, offset;
1394 
1395     /* Make new fine FE which refines the ref cell and the quadrature rule */
1396     ierr = PetscFEGetQuadrature(feRef[fieldI], &quad);CHKERRQ(ierr);
1397     ierr = PetscFEGetDimension(feRef[fieldI], &Nb);CHKERRQ(ierr);
1398     ierr = PetscFEGetTileSizes(feRef[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1399     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1400     blockSize = Nb*numQuadPoints;
1401     batchSize = numBlocks * blockSize;
1402     ierr = PetscFESetTileSizes(feRef[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1403     numChunks = numCells / (numBatches*batchSize);
1404     Ne        = numChunks*numBatches*batchSize;
1405     Nr        = numCells % (numBatches*batchSize);
1406     offset    = numCells - Nr;
1407 
1408     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1409       /* void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ]; */
1410       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = g0_identity_2d_static;
1411 
1412       /* Replace quadrature in coarse FE with refined quadrature */
1413       ierr = PetscFEGetQuadrature(fe[fieldJ], &quadOld);CHKERRQ(ierr);
1414       ierr = PetscObjectReference((PetscObject) quadOld);CHKERRQ(ierr);
1415       ierr = PetscFESetQuadrature(fe[fieldJ], quad);CHKERRQ(ierr);
1416       geom.v0   = v0;
1417       geom.J    = J;
1418       geom.invJ = invJ;
1419       geom.detJ = detJ;
1420       ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Ne, Nf, feRef, fieldI, fe, fieldJ, geom, g0, elemMat);CHKERRQ(ierr);
1421       geom.v0   = &v0[offset*dim];
1422       geom.J    = &J[offset*dim*dim];
1423       geom.invJ = &invJ[offset*dim*dim];
1424       geom.detJ = &detJ[offset];
1425       ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Nr, Nf, feRef, fieldI, fe, fieldJ, geom, g0, &elemMat[offset*rCellDof*cCellDof]);CHKERRQ(ierr);
1426       ierr = PetscFESetQuadrature(fe[fieldJ], quadOld);CHKERRQ(ierr);
1427     }
1428   }
1429   for (c = cStart; c < cEnd; ++c) {
1430     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, rCellDof, cCellDof, &elemMat[c*rCellDof*cCellDof]);CHKERRQ(ierr);}
1431     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, I, c, &elemMat[c*rCellDof*cCellDof], ADD_VALUES);CHKERRQ(ierr);
1432   }
1433   ierr = PetscFree5(v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1434   ierr = PetscFree(feRef);CHKERRQ(ierr);
1435   ierr = MatAssemblyBegin(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1436   ierr = MatAssemblyEnd(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1437   if (mesh->printFEM) {
1438     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1439     ierr = MatChop(I, 1.0e-10);CHKERRQ(ierr);
1440     ierr = MatView(I, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1441   }
1442 #if 0
1443   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1444 #endif
1445   PetscFunctionReturn(0);
1446 }
1447 #endif
1448 
1449 #undef __FUNCT__
1450 #define __FUNCT__ "DMPlexComputeInterpolatorFEM"
1451 /*@
1452   DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1453 
1454   Input Parameters:
1455 + dmf  - The fine mesh
1456 . dmc  - The coarse mesh
1457 - user - The user context
1458 
1459   Output Parameter:
1460 . In  - The interpolation matrix
1461 
1462   Note:
1463   The first member of the user context must be an FEMContext.
1464 
1465   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1466   like a GPU, or vectorize on a multicore machine.
1467 
1468   Level: developer
1469 
1470 .seealso: DMPlexComputeJacobianFEM()
1471 @*/
1472 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1473 {
1474   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1475   PetscFEM         *fem   = (PetscFEM *) user;
1476   PetscFE          *fe    = fem->fe;
1477   const char       *name  = "Interpolator";
1478   PetscFE          *feRef;
1479   PetscSection      fsection, fglobalSection;
1480   PetscSection      csection, cglobalSection;
1481   PetscScalar      *elemMat;
1482   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c;
1483   PetscInt          rCellDof = 0, cCellDof = 0;
1484   PetscErrorCode    ierr;
1485 
1486   PetscFunctionBegin;
1487 #if 0
1488   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1489 #endif
1490   ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr);
1491   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1492   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1493   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1494   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1495   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1496   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1497   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1498   for (f = 0; f < Nf; ++f) {
1499     PetscInt rNb, cNb, Nc;
1500 
1501     ierr = PetscFERefine(fe[f], &feRef[f]);CHKERRQ(ierr);
1502     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1503     ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr);
1504     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1505     rCellDof += rNb*Nc;
1506     cCellDof += cNb*Nc;
1507   }
1508   ierr = MatZeroEntries(In);CHKERRQ(ierr);
1509   ierr = PetscMalloc1(rCellDof*cCellDof,&elemMat);CHKERRQ(ierr);
1510   ierr = PetscMemzero(elemMat, rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1511   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1512     PetscDualSpace   Qref;
1513     PetscQuadrature  f;
1514     const PetscReal *qpoints, *qweights;
1515     PetscReal       *points;
1516     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;
1517 
1518     /* Compose points from all dual basis functionals */
1519     ierr = PetscFEGetNumComponents(fe[fieldI], &Nc);CHKERRQ(ierr);
1520     ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr);
1521     ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr);
1522     for (i = 0; i < fpdim; ++i) {
1523       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1524       ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr);
1525       npoints += Np;
1526     }
1527     ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr);
1528     for (i = 0, k = 0; i < fpdim; ++i) {
1529       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1530       ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr);
1531       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1532     }
1533 
1534     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1535       PetscReal *B;
1536       PetscInt   NcJ, cpdim, j;
1537 
1538       /* Evaluate basis at points */
1539       ierr = PetscFEGetNumComponents(fe[fieldJ], &NcJ);CHKERRQ(ierr);
1540       if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ);
1541       ierr = PetscFEGetDimension(fe[fieldJ], &cpdim);CHKERRQ(ierr);
1542       /* For now, fields only interpolate themselves */
1543       if (fieldI == fieldJ) {
1544         ierr = PetscFEGetTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);
1545         for (i = 0, k = 0; i < fpdim; ++i) {
1546           ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1547           ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr);
1548           for (p = 0; p < Np; ++p, ++k) {
1549             for (j = 0; j < cpdim; ++j) {
1550               for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cCellDof + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p];
1551             }
1552           }
1553         }
1554         ierr = PetscFERestoreTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr);
1555       }
1556       offsetJ += cpdim*NcJ;
1557     }
1558     offsetI += fpdim*Nc;
1559     ierr = PetscFree(points);CHKERRQ(ierr);
1560   }
1561   if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rCellDof, cCellDof, elemMat);CHKERRQ(ierr);}
1562   for (c = cStart; c < cEnd; ++c) {
1563     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr);
1564   }
1565   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);}
1566   ierr = PetscFree(feRef);CHKERRQ(ierr);
1567   ierr = PetscFree(elemMat);CHKERRQ(ierr);
1568   ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1569   ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1570   if (mesh->printFEM) {
1571     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1572     ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr);
1573     ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1574   }
1575 #if 0
1576   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1577 #endif
1578   PetscFunctionReturn(0);
1579 }
1580 
1581 #undef __FUNCT__
1582 #define __FUNCT__ "DMPlexAddBoundary"
1583 /* The ids can be overridden by the command line option -bc_<boundary name> */
1584 PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], PetscInt field, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx)
1585 {
1586   DM_Plex       *mesh = (DM_Plex *) dm->data;
1587   DMBoundary     b;
1588   PetscErrorCode ierr;
1589 
1590   PetscFunctionBegin;
1591   ierr = PetscNew(&b);CHKERRQ(ierr);
1592   ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr);
1593   ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr);
1594   ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr);
1595   b->essential   = isEssential;
1596   b->field       = field;
1597   b->func        = bcFunc;
1598   b->numids      = numids;
1599   b->ctx         = ctx;
1600   b->next        = mesh->boundary;
1601   mesh->boundary = b;
1602   PetscFunctionReturn(0);
1603 }
1604 
1605 #undef __FUNCT__
1606 #define __FUNCT__ "DMPlexGetNumBoundary"
1607 PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd)
1608 {
1609   DM_Plex   *mesh = (DM_Plex *) dm->data;
1610   DMBoundary b    = mesh->boundary;
1611 
1612   PetscFunctionBegin;
1613   *numBd = 0;
1614   while (b) {++(*numBd); b = b->next;}
1615   PetscFunctionReturn(0);
1616 }
1617 
1618 #undef __FUNCT__
1619 #define __FUNCT__ "DMPlexGetBoundary"
1620 PetscErrorCode DMPlexGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, PetscInt *field, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx)
1621 {
1622   DM_Plex   *mesh = (DM_Plex *) dm->data;
1623   DMBoundary b    = mesh->boundary;
1624   PetscInt   n    = 0;
1625 
1626   PetscFunctionBegin;
1627   while (b) {
1628     if (n == bd) break;
1629     b = b->next;
1630     ++n;
1631   }
1632   if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
1633   if (isEssential) {
1634     PetscValidPointer(isEssential, 3);
1635     *isEssential = b->essential;
1636   }
1637   if (name) {
1638     PetscValidPointer(name, 4);
1639     *name = b->name;
1640   }
1641   if (field) {
1642     PetscValidPointer(field, 5);
1643     *field = b->field;
1644   }
1645   if (func) {
1646     PetscValidPointer(func, 6);
1647     *func = b->func;
1648   }
1649   if (numids) {
1650     PetscValidPointer(numids, 7);
1651     *numids = b->numids;
1652   }
1653   if (ids) {
1654     PetscValidPointer(ids, 8);
1655     *ids = b->ids;
1656   }
1657   if (ctx) {
1658     PetscValidPointer(ctx, 9);
1659     *ctx = b->ctx;
1660   }
1661   PetscFunctionReturn(0);
1662 }
1663