xref: /petsc/src/dm/impls/plex/plexfem.c (revision bd5ce8098c519eb11f95642784d8bbf54bc60d52) 
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 . X  - Local input vector
824 . X_t  - Time derivative of the local input vector
825 - user - The user context
826 
827   Output Parameter:
828 . F  - Local output vector
829 
830   Note:
831   The first member of the user context must be an FEMContext.
832 
833   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
834   like a GPU, or vectorize on a multicore machine.
835 
836   Level: developer
837 
838 .seealso: DMPlexComputeResidualFEM()
839 @*/
840 PetscErrorCode DMPlexComputeIFunctionFEM(DM dm, Vec X, Vec X_t, Vec F, void *user)
841 {
842   DM_Plex          *mesh  = (DM_Plex *) dm->data;
843   PetscFEM         *fem   = (PetscFEM *) user;
844   PetscFE          *fe    = fem->fe;
845   PetscFE          *feAux = fem->feAux;
846   PetscFE          *feBd  = fem->feBd;
847   const char       *name  = "Residual";
848   DM                dmAux;
849   Vec               A;
850   PetscQuadrature   q;
851   PetscCellGeometry geom;
852   PetscSection      section, sectionAux;
853   PetscReal        *v0, *J, *invJ, *detJ;
854   PetscScalar      *elemVec, *u, *u_t, *a = NULL;
855   PetscInt          dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c;
856   PetscInt          cellDof = 0, numComponents = 0;
857   PetscInt          cellDofAux = 0, numComponentsAux = 0;
858   PetscErrorCode    ierr;
859 
860   PetscFunctionBegin;
861   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
862   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
863   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
864   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
865   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
866   numCells = cEnd - cStart;
867   for (f = 0; f < Nf; ++f) {
868     PetscInt Nb, Nc;
869 
870     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
871     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
872     cellDof       += Nb*Nc;
873     numComponents += Nc;
874   }
875   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
876   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
877   if (dmAux) {
878     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
879     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
880   }
881   for (f = 0; f < NfAux; ++f) {
882     PetscInt Nb, Nc;
883 
884     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
885     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
886     cellDofAux       += Nb*Nc;
887     numComponentsAux += Nc;
888   }
889   ierr = DMPlexInsertBoundaryValues(dm, X);CHKERRQ(ierr);
890   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
891   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);
892   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
893   for (c = cStart; c < cEnd; ++c) {
894     PetscScalar *x = NULL, *x_t = NULL;
895     PetscInt     i;
896 
897     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
898     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
899     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
900     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
901     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
902     ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
903     for (i = 0; i < cellDof; ++i) u_t[c*cellDof+i] = x_t[i];
904     ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr);
905     if (dmAux) {
906       PetscScalar *x_a = NULL;
907       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr);
908       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x_a[i];
909       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr);
910     }
911   }
912   for (f = 0; f < Nf; ++f) {
913     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0IFuncs[f];
914     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1IFuncs[f];
915     PetscInt numQuadPoints, Nb;
916     /* Conforming batches */
917     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
918     /* Remainder */
919     PetscInt Nr, offset;
920 
921     ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr);
922     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
923     ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
924     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
925     blockSize = Nb*numQuadPoints;
926     batchSize = numBlocks * blockSize;
927     ierr =  PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
928     numChunks = numCells / (numBatches*batchSize);
929     Ne        = numChunks*numBatches*batchSize;
930     Nr        = numCells % (numBatches*batchSize);
931     offset    = numCells - Nr;
932     geom.v0   = v0;
933     geom.J    = J;
934     geom.invJ = invJ;
935     geom.detJ = detJ;
936     ierr = PetscFEIntegrateIFunction(fe[f], Ne, Nf, fe, f, geom, u, u_t, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr);
937     geom.v0   = &v0[offset*dim];
938     geom.J    = &J[offset*dim*dim];
939     geom.invJ = &invJ[offset*dim*dim];
940     geom.detJ = &detJ[offset];
941     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);
942   }
943   for (c = cStart; c < cEnd; ++c) {
944     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
945     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
946   }
947   ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
948   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
949   if (feBd) {
950     DMLabel         label, depth;
951     IS              pointIS;
952     const PetscInt *points;
953     PetscInt        dep, numPoints, p, numFaces;
954     PetscReal      *n;
955 
956     ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
957     ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
958     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
959     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
960     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
961     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
962       PetscInt Nb, Nc;
963 
964       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
965       ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr);
966       cellDof       += Nb*Nc;
967       numComponents += Nc;
968     }
969     for (p = 0, numFaces = 0; p < numPoints; ++p) {
970       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
971       if (dep == dim-1) ++numFaces;
972     }
973     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);
974     ierr = PetscMalloc1(numFaces*cellDof,&u_t);CHKERRQ(ierr);
975     for (p = 0, f = 0; p < numPoints; ++p) {
976       const PetscInt point = points[p];
977       PetscScalar   *x     = NULL;
978       PetscInt       i;
979 
980       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
981       if (dep != dim-1) continue;
982       ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
983       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
984       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
985       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
986       for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
987       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
988       ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
989       for (i = 0; i < cellDof; ++i) u_t[f*cellDof+i] = x[i];
990       ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr);
991       ++f;
992     }
993     for (f = 0; f < Nf; ++f) {
994       void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdIFuncs[f];
995       void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdIFuncs[f];
996       PetscInt numQuadPoints, Nb;
997       /* Conforming batches */
998       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
999       /* Remainder */
1000       PetscInt Nr, offset;
1001 
1002       ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr);
1003       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
1004       ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1005       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1006       blockSize = Nb*numQuadPoints;
1007       batchSize = numBlocks * blockSize;
1008       ierr =  PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1009       numChunks = numFaces / (numBatches*batchSize);
1010       Ne        = numChunks*numBatches*batchSize;
1011       Nr        = numFaces % (numBatches*batchSize);
1012       offset    = numFaces - Nr;
1013       geom.v0   = v0;
1014       geom.n    = n;
1015       geom.J    = J;
1016       geom.invJ = invJ;
1017       geom.detJ = detJ;
1018       ierr = PetscFEIntegrateBdIFunction(feBd[f], Ne, Nf, feBd, f, geom, u, u_t, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr);
1019       geom.v0   = &v0[offset*dim];
1020       geom.n    = &n[offset*dim];
1021       geom.J    = &J[offset*dim*dim];
1022       geom.invJ = &invJ[offset*dim*dim];
1023       geom.detJ = &detJ[offset];
1024       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);
1025     }
1026     for (p = 0, f = 0; p < numPoints; ++p) {
1027       const PetscInt point = points[p];
1028 
1029       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
1030       if (dep != dim-1) continue;
1031       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);}
1032       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr);
1033       ++f;
1034     }
1035     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
1036     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
1037     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
1038     ierr = PetscFree(u_t);CHKERRQ(ierr);
1039   }
1040   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
1041   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
1042   PetscFunctionReturn(0);
1043 }
1044 
1045 #undef __FUNCT__
1046 #define __FUNCT__ "DMPlexComputeJacobianActionFEM"
1047 /*@C
1048   DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user
1049 
1050   Input Parameters:
1051 + dm - The mesh
1052 . J  - The Jacobian shell matrix
1053 . X  - Local input vector
1054 - user - The user context
1055 
1056   Output Parameter:
1057 . F  - Local output vector
1058 
1059   Note:
1060   The first member of the user context must be an FEMContext.
1061 
1062   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1063   like a GPU, or vectorize on a multicore machine.
1064 
1065   Level: developer
1066 
1067 .seealso: DMPlexComputeResidualFEM()
1068 @*/
1069 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
1070 {
1071   DM_Plex          *mesh = (DM_Plex *) dm->data;
1072   PetscFEM         *fem  = (PetscFEM *) user;
1073   PetscFE          *fe   = fem->fe;
1074   PetscQuadrature   quad;
1075   PetscCellGeometry geom;
1076   PetscSection      section;
1077   JacActionCtx     *jctx;
1078   PetscReal        *v0, *J, *invJ, *detJ;
1079   PetscScalar      *elemVec, *u, *a;
1080   PetscInt          dim, numFields, field, numCells, cStart, cEnd, c;
1081   PetscInt          cellDof = 0;
1082   PetscErrorCode    ierr;
1083 
1084   PetscFunctionBegin;
1085   /* ierr = PetscLogEventBegin(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */
1086   ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
1087   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
1088   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1089   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
1090   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1091   numCells = cEnd - cStart;
1092   for (field = 0; field < numFields; ++field) {
1093     PetscInt Nb, Nc;
1094 
1095     ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
1096     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
1097     cellDof += Nb*Nc;
1098   }
1099   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
1100   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);
1101   for (c = cStart; c < cEnd; ++c) {
1102     PetscScalar *x = NULL;
1103     PetscInt     i;
1104 
1105     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1106     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1107     ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
1108     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
1109     ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
1110     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
1111     for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
1112     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
1113   }
1114   for (field = 0; field < numFields; ++field) {
1115     PetscInt numQuadPoints, Nb;
1116     /* Conforming batches */
1117     PetscInt numBlocks  = 1;
1118     PetscInt numBatches = 1;
1119     PetscInt numChunks, Ne, blockSize, batchSize;
1120     /* Remainder */
1121     PetscInt Nr, offset;
1122 
1123     ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr);
1124     ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
1125     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1126     blockSize = Nb*numQuadPoints;
1127     batchSize = numBlocks * blockSize;
1128     numChunks = numCells / (numBatches*batchSize);
1129     Ne        = numChunks*numBatches*batchSize;
1130     Nr        = numCells % (numBatches*batchSize);
1131     offset    = numCells - Nr;
1132     geom.v0   = v0;
1133     geom.J    = J;
1134     geom.invJ = invJ;
1135     geom.detJ = detJ;
1136     ierr = PetscFEIntegrateJacobianAction(fe[field], Ne, numFields, fe, field, geom, u, a, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr);
1137     geom.v0   = &v0[offset*dim];
1138     geom.J    = &J[offset*dim*dim];
1139     geom.invJ = &invJ[offset*dim*dim];
1140     geom.detJ = &detJ[offset];
1141     ierr = PetscFEIntegrateJacobianAction(fe[field], Nr, numFields, fe, field, geom, &u[offset*cellDof], &a[offset*cellDof],
1142                                           fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr);
1143   }
1144   for (c = cStart; c < cEnd; ++c) {
1145     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
1146     ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
1147   }
1148   ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
1149   if (mesh->printFEM) {
1150     PetscMPIInt rank, numProcs;
1151     PetscInt    p;
1152 
1153     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
1154     ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
1155     ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");CHKERRQ(ierr);
1156     for (p = 0; p < numProcs; ++p) {
1157       if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}
1158       ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
1159     }
1160   }
1161   /* ierr = PetscLogEventEnd(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */
1162   PetscFunctionReturn(0);
1163 }
1164 
1165 #undef __FUNCT__
1166 #define __FUNCT__ "DMPlexComputeJacobianFEM"
1167 /*@
1168   DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
1169 
1170   Input Parameters:
1171 + dm - The mesh
1172 . X  - Local input vector
1173 - user - The user context
1174 
1175   Output Parameter:
1176 . Jac  - Jacobian matrix
1177 
1178   Note:
1179   The first member of the user context must be an FEMContext.
1180 
1181   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1182   like a GPU, or vectorize on a multicore machine.
1183 
1184   Level: developer
1185 
1186 .seealso: FormFunctionLocal()
1187 @*/
1188 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1189 {
1190   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1191   PetscFEM         *fem   = (PetscFEM *) user;
1192   PetscFE          *fe    = fem->fe;
1193   PetscFE          *feAux = fem->feAux;
1194   const char       *name  = "Jacobian";
1195   DM                dmAux;
1196   Vec               A;
1197   PetscQuadrature   quad;
1198   PetscCellGeometry geom;
1199   PetscSection      section, globalSection, sectionAux;
1200   PetscReal        *v0, *J, *invJ, *detJ;
1201   PetscScalar      *elemMat, *u, *a;
1202   PetscInt          dim, Nf, NfAux = 0, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1203   PetscInt          cellDof = 0, numComponents = 0;
1204   PetscInt          cellDofAux = 0, numComponentsAux = 0;
1205   PetscBool         isShell;
1206   PetscErrorCode    ierr;
1207 
1208   PetscFunctionBegin;
1209   ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1210   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
1211   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1212   ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);
1213   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1214   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1215   numCells = cEnd - cStart;
1216   for (f = 0; f < Nf; ++f) {
1217     PetscInt Nb, Nc;
1218 
1219     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
1220     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1221     cellDof       += Nb*Nc;
1222     numComponents += Nc;
1223   }
1224   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1225   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1226   if (dmAux) {
1227     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1228     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
1229   }
1230   for (f = 0; f < NfAux; ++f) {
1231     PetscInt Nb, Nc;
1232 
1233     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
1234     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
1235     cellDofAux       += Nb*Nc;
1236     numComponentsAux += Nc;
1237   }
1238   ierr = DMPlexInsertBoundaryValues(dm, X);CHKERRQ(ierr);
1239   ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
1240   ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof*cellDof,&elemMat);CHKERRQ(ierr);
1241   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
1242   for (c = cStart; c < cEnd; ++c) {
1243     PetscScalar *x = NULL;
1244     PetscInt     i;
1245 
1246     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1247     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1248     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1249     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
1250     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
1251     if (dmAux) {
1252       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1253       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
1254       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1255     }
1256   }
1257   ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1258   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1259     PetscInt numQuadPoints, Nb;
1260     /* Conforming batches */
1261     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1262     /* Remainder */
1263     PetscInt Nr, offset;
1264 
1265     ierr = PetscFEGetQuadrature(fe[fieldI], &quad);CHKERRQ(ierr);
1266     ierr = PetscFEGetDimension(fe[fieldI], &Nb);CHKERRQ(ierr);
1267     ierr = PetscFEGetTileSizes(fe[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1268     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1269     blockSize = Nb*numQuadPoints;
1270     batchSize = numBlocks * blockSize;
1271     ierr = PetscFESetTileSizes(fe[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1272     numChunks = numCells / (numBatches*batchSize);
1273     Ne        = numChunks*numBatches*batchSize;
1274     Nr        = numCells % (numBatches*batchSize);
1275     offset    = numCells - Nr;
1276     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1277       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ];
1278       void   (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*Nf+fieldJ];
1279       void   (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*Nf+fieldJ];
1280       void   (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*Nf+fieldJ];
1281 
1282       geom.v0   = v0;
1283       geom.J    = J;
1284       geom.invJ = invJ;
1285       geom.detJ = detJ;
1286       ierr = PetscFEIntegrateJacobian(fe[fieldI], Ne, Nf, fe, fieldI, fieldJ, geom, u, NfAux, feAux, a, g0, g1, g2, g3, elemMat);CHKERRQ(ierr);
1287       geom.v0   = &v0[offset*dim];
1288       geom.J    = &J[offset*dim*dim];
1289       geom.invJ = &invJ[offset*dim*dim];
1290       geom.detJ = &detJ[offset];
1291       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);
1292     }
1293   }
1294   for (c = cStart; c < cEnd; ++c) {
1295     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);}
1296     ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr);
1297   }
1298   ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1299   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1300   ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1301   ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1302   if (mesh->printFEM) {
1303     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1304     ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
1305     ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1306   }
1307   ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
1308   ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr);
1309   if (isShell) {
1310     JacActionCtx *jctx;
1311 
1312     ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
1313     ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
1314   }
1315   PetscFunctionReturn(0);
1316 }
1317 
1318 #if 0
1319 
1320 static void g0_identity_1d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1321 {
1322   g3[0] = 1.0;
1323 }
1324 
1325 static void g0_identity_2d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1326 {
1327   g3[0] = g3[3] = 1.0;
1328 }
1329 
1330 static void g0_identity_3d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
1331 {
1332   g3[0] = g3[4] = g3[8] = 1.0;
1333 }
1334 
1335 #undef __FUNCT__
1336 #define __FUNCT__ "DMPlexComputeInterpolatorFEMBroken"
1337 PetscErrorCode DMPlexComputeInterpolatorFEMBroken(DM dmc, DM dmf, Mat I, void *user)
1338 {
1339   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1340   PetscFEM         *fem   = (PetscFEM *) user;
1341   PetscFE          *fe    = fem->fe;
1342   const char       *name  = "Interpolator";
1343   PetscFE          *feRef;
1344   PetscQuadrature   quad, quadOld;
1345   PetscCellGeometry geom;
1346   PetscSection      fsection, fglobalSection;
1347   PetscSection      csection, cglobalSection;
1348   PetscReal        *v0, *J, *invJ, *detJ;
1349   PetscScalar      *elemMat;
1350   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1351   PetscInt          rCellDof = 0, cCellDof = 0, numComponents = 0;
1352   PetscErrorCode    ierr;
1353 
1354   PetscFunctionBegin;
1355 #if 0
1356   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1357 #endif
1358   ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr);
1359   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1360   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1361   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1362   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1363   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1364   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1365   numCells = cEnd - cStart;
1366   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1367   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1368     ierr = PetscFERefine(fe[fieldI], &feRef[fieldI]);CHKERRQ(ierr);
1369   }
1370   for (f = 0; f < Nf; ++f) {
1371     PetscInt rNb, cNb, Nc;
1372 
1373     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1374     ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr);
1375     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1376     numComponents += Nc;
1377     rCellDof += rNb*Nc;
1378     cCellDof += cNb*Nc;
1379   }
1380   ierr = MatZeroEntries(I);CHKERRQ(ierr);
1381   ierr = PetscMalloc5(numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*rCellDof*cCellDof,&elemMat);CHKERRQ(ierr);
1382   for (c = cStart; c < cEnd; ++c) {
1383     ierr = DMPlexComputeCellGeometry(dmc, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1384     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1385   }
1386   ierr = PetscMemzero(elemMat, numCells*rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1387   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1388     PetscInt numQuadPoints, Nb;
1389     /* Conforming batches */
1390     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1391     /* Remainder */
1392     PetscInt Nr, offset;
1393 
1394     /* Make new fine FE which refines the ref cell and the quadrature rule */
1395     ierr = PetscFEGetQuadrature(feRef[fieldI], &quad);CHKERRQ(ierr);
1396     ierr = PetscFEGetDimension(feRef[fieldI], &Nb);CHKERRQ(ierr);
1397     ierr = PetscFEGetTileSizes(feRef[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1398     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1399     blockSize = Nb*numQuadPoints;
1400     batchSize = numBlocks * blockSize;
1401     ierr = PetscFESetTileSizes(feRef[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1402     numChunks = numCells / (numBatches*batchSize);
1403     Ne        = numChunks*numBatches*batchSize;
1404     Nr        = numCells % (numBatches*batchSize);
1405     offset    = numCells - Nr;
1406 
1407     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1408       /* void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ]; */
1409       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = g0_identity_2d_static;
1410 
1411       /* Replace quadrature in coarse FE with refined quadrature */
1412       ierr = PetscFEGetQuadrature(fe[fieldJ], &quadOld);CHKERRQ(ierr);
1413       ierr = PetscObjectReference((PetscObject) quadOld);CHKERRQ(ierr);
1414       ierr = PetscFESetQuadrature(fe[fieldJ], quad);CHKERRQ(ierr);
1415       geom.v0   = v0;
1416       geom.J    = J;
1417       geom.invJ = invJ;
1418       geom.detJ = detJ;
1419       ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Ne, Nf, feRef, fieldI, fe, fieldJ, geom, g0, elemMat);CHKERRQ(ierr);
1420       geom.v0   = &v0[offset*dim];
1421       geom.J    = &J[offset*dim*dim];
1422       geom.invJ = &invJ[offset*dim*dim];
1423       geom.detJ = &detJ[offset];
1424       ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Nr, Nf, feRef, fieldI, fe, fieldJ, geom, g0, &elemMat[offset*rCellDof*cCellDof]);CHKERRQ(ierr);
1425       ierr = PetscFESetQuadrature(fe[fieldJ], quadOld);CHKERRQ(ierr);
1426     }
1427   }
1428   for (c = cStart; c < cEnd; ++c) {
1429     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, rCellDof, cCellDof, &elemMat[c*rCellDof*cCellDof]);CHKERRQ(ierr);}
1430     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, I, c, &elemMat[c*rCellDof*cCellDof], ADD_VALUES);CHKERRQ(ierr);
1431   }
1432   ierr = PetscFree5(v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1433   ierr = PetscFree(feRef);CHKERRQ(ierr);
1434   ierr = MatAssemblyBegin(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1435   ierr = MatAssemblyEnd(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1436   if (mesh->printFEM) {
1437     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1438     ierr = MatChop(I, 1.0e-10);CHKERRQ(ierr);
1439     ierr = MatView(I, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1440   }
1441 #if 0
1442   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1443 #endif
1444   PetscFunctionReturn(0);
1445 }
1446 #endif
1447 
1448 #undef __FUNCT__
1449 #define __FUNCT__ "DMPlexComputeInterpolatorFEM"
1450 /*@
1451   DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1452 
1453   Input Parameters:
1454 + dmf  - The fine mesh
1455 . dmc  - The coarse mesh
1456 - user - The user context
1457 
1458   Output Parameter:
1459 . In  - The interpolation matrix
1460 
1461   Note:
1462   The first member of the user context must be an FEMContext.
1463 
1464   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1465   like a GPU, or vectorize on a multicore machine.
1466 
1467   Level: developer
1468 
1469 .seealso: DMPlexComputeJacobianFEM()
1470 @*/
1471 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1472 {
1473   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1474   PetscFEM         *fem   = (PetscFEM *) user;
1475   PetscFE          *fe    = fem->fe;
1476   const char       *name  = "Interpolator";
1477   PetscFE          *feRef;
1478   PetscSection      fsection, fglobalSection;
1479   PetscSection      csection, cglobalSection;
1480   PetscScalar      *elemMat;
1481   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c;
1482   PetscInt          rCellDof = 0, cCellDof = 0;
1483   PetscErrorCode    ierr;
1484 
1485   PetscFunctionBegin;
1486 #if 0
1487   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1488 #endif
1489   ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr);
1490   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1491   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1492   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1493   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1494   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1495   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1496   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1497   for (f = 0; f < Nf; ++f) {
1498     PetscInt rNb, cNb, Nc;
1499 
1500     ierr = PetscFERefine(fe[f], &feRef[f]);CHKERRQ(ierr);
1501     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1502     ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr);
1503     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1504     rCellDof += rNb*Nc;
1505     cCellDof += cNb*Nc;
1506   }
1507   ierr = MatZeroEntries(In);CHKERRQ(ierr);
1508   ierr = PetscMalloc1(rCellDof*cCellDof,&elemMat);CHKERRQ(ierr);
1509   ierr = PetscMemzero(elemMat, rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1510   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1511     PetscDualSpace   Qref;
1512     PetscQuadrature  f;
1513     const PetscReal *qpoints, *qweights;
1514     PetscReal       *points;
1515     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;
1516 
1517     /* Compose points from all dual basis functionals */
1518     ierr = PetscFEGetNumComponents(fe[fieldI], &Nc);CHKERRQ(ierr);
1519     ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr);
1520     ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr);
1521     for (i = 0; i < fpdim; ++i) {
1522       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1523       ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr);
1524       npoints += Np;
1525     }
1526     ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr);
1527     for (i = 0, k = 0; i < fpdim; ++i) {
1528       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1529       ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr);
1530       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1531     }
1532 
1533     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1534       PetscReal *B;
1535       PetscInt   NcJ, cpdim, j;
1536 
1537       /* Evaluate basis at points */
1538       ierr = PetscFEGetNumComponents(fe[fieldJ], &NcJ);CHKERRQ(ierr);
1539       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);
1540       ierr = PetscFEGetDimension(fe[fieldJ], &cpdim);CHKERRQ(ierr);
1541       /* For now, fields only interpolate themselves */
1542       if (fieldI == fieldJ) {
1543         ierr = PetscFEGetTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);
1544         for (i = 0, k = 0; i < fpdim; ++i) {
1545           ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1546           ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr);
1547           for (p = 0; p < Np; ++p, ++k) {
1548             for (j = 0; j < cpdim; ++j) {
1549               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];
1550             }
1551           }
1552         }
1553         ierr = PetscFERestoreTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr);
1554       }
1555       offsetJ += cpdim*NcJ;
1556     }
1557     offsetI += fpdim*Nc;
1558     ierr = PetscFree(points);CHKERRQ(ierr);
1559   }
1560   if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rCellDof, cCellDof, elemMat);CHKERRQ(ierr);}
1561   for (c = cStart; c < cEnd; ++c) {
1562     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr);
1563   }
1564   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);}
1565   ierr = PetscFree(feRef);CHKERRQ(ierr);
1566   ierr = PetscFree(elemMat);CHKERRQ(ierr);
1567   ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1568   ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1569   if (mesh->printFEM) {
1570     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1571     ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr);
1572     ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1573   }
1574 #if 0
1575   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1576 #endif
1577   PetscFunctionReturn(0);
1578 }
1579 
1580 #undef __FUNCT__
1581 #define __FUNCT__ "DMPlexAddBoundary"
1582 /* The ids can be overridden by the command line option -bc_<boundary name> */
1583 PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], PetscInt field, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx)
1584 {
1585   DM_Plex       *mesh = (DM_Plex *) dm->data;
1586   DMBoundary     b;
1587   PetscErrorCode ierr;
1588 
1589   PetscFunctionBegin;
1590   ierr = PetscNew(&b);CHKERRQ(ierr);
1591   ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr);
1592   ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr);
1593   ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr);
1594   b->essential   = isEssential;
1595   b->field       = field;
1596   b->func        = bcFunc;
1597   b->numids      = numids;
1598   b->ctx         = ctx;
1599   b->next        = mesh->boundary;
1600   mesh->boundary = b;
1601   PetscFunctionReturn(0);
1602 }
1603 
1604 #undef __FUNCT__
1605 #define __FUNCT__ "DMPlexGetNumBoundary"
1606 PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd)
1607 {
1608   DM_Plex   *mesh = (DM_Plex *) dm->data;
1609   DMBoundary b    = mesh->boundary;
1610 
1611   PetscFunctionBegin;
1612   *numBd = 0;
1613   while (b) {++(*numBd); b = b->next;}
1614   PetscFunctionReturn(0);
1615 }
1616 
1617 #undef __FUNCT__
1618 #define __FUNCT__ "DMPlexGetBoundary"
1619 PetscErrorCode DMPlexGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, PetscInt *field, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx)
1620 {
1621   DM_Plex   *mesh = (DM_Plex *) dm->data;
1622   DMBoundary b    = mesh->boundary;
1623   PetscInt   n    = 0;
1624 
1625   PetscFunctionBegin;
1626   while (b) {
1627     if (n == bd) break;
1628     b = b->next;
1629     ++n;
1630   }
1631   if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
1632   if (isEssential) {
1633     PetscValidPointer(isEssential, 3);
1634     *isEssential = b->essential;
1635   }
1636   if (name) {
1637     PetscValidPointer(name, 4);
1638     *name = b->name;
1639   }
1640   if (field) {
1641     PetscValidPointer(field, 5);
1642     *field = b->field;
1643   }
1644   if (func) {
1645     PetscValidPointer(func, 6);
1646     *func = b->func;
1647   }
1648   if (numids) {
1649     PetscValidPointer(numids, 7);
1650     *numids = b->numids;
1651   }
1652   if (ids) {
1653     PetscValidPointer(ids, 8);
1654     *ids = b->ids;
1655   }
1656   if (ctx) {
1657     PetscValidPointer(ctx, 9);
1658     *ctx = b->ctx;
1659   }
1660   PetscFunctionReturn(0);
1661 }
1662