xref: /petsc/src/dm/impls/plex/plexfem.c (revision 7fdeb8b9be692273ce99bf256278a44482a052ca) 
1 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
3 #include <petscfe.h>
4 
5 #undef __FUNCT__
6 #define __FUNCT__ "DMPlexGetScale"
7 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
8 {
9   DM_Plex *mesh = (DM_Plex*) dm->data;
10 
11   PetscFunctionBegin;
12   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
13   PetscValidPointer(scale, 3);
14   *scale = mesh->scale[unit];
15   PetscFunctionReturn(0);
16 }
17 
18 #undef __FUNCT__
19 #define __FUNCT__ "DMPlexSetScale"
20 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
21 {
22   DM_Plex *mesh = (DM_Plex*) dm->data;
23 
24   PetscFunctionBegin;
25   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
26   mesh->scale[unit] = scale;
27   PetscFunctionReturn(0);
28 }
29 
30 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
31 {
32   switch (i) {
33   case 0:
34     switch (j) {
35     case 0: return 0;
36     case 1:
37       switch (k) {
38       case 0: return 0;
39       case 1: return 0;
40       case 2: return 1;
41       }
42     case 2:
43       switch (k) {
44       case 0: return 0;
45       case 1: return -1;
46       case 2: return 0;
47       }
48     }
49   case 1:
50     switch (j) {
51     case 0:
52       switch (k) {
53       case 0: return 0;
54       case 1: return 0;
55       case 2: return -1;
56       }
57     case 1: return 0;
58     case 2:
59       switch (k) {
60       case 0: return 1;
61       case 1: return 0;
62       case 2: return 0;
63       }
64     }
65   case 2:
66     switch (j) {
67     case 0:
68       switch (k) {
69       case 0: return 0;
70       case 1: return 1;
71       case 2: return 0;
72       }
73     case 1:
74       switch (k) {
75       case 0: return -1;
76       case 1: return 0;
77       case 2: return 0;
78       }
79     case 2: return 0;
80     }
81   }
82   return 0;
83 }
84 
85 #undef __FUNCT__
86 #define __FUNCT__ "DMPlexCreateRigidBody"
87 /*@C
88   DMPlexCreateRigidBody - create rigid body modes from coordinates
89 
90   Collective on DM
91 
92   Input Arguments:
93 + dm - the DM
94 . section - the local section associated with the rigid field, or NULL for the default section
95 - globalSection - the global section associated with the rigid field, or NULL for the default section
96 
97   Output Argument:
98 . sp - the null space
99 
100   Note: This is necessary to take account of Dirichlet conditions on the displacements
101 
102   Level: advanced
103 
104 .seealso: MatNullSpaceCreate()
105 @*/
106 PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp)
107 {
108   MPI_Comm       comm;
109   Vec            coordinates, localMode, mode[6];
110   PetscSection   coordSection;
111   PetscScalar   *coords;
112   PetscInt       dim, vStart, vEnd, v, n, m, d, i, j;
113   PetscErrorCode ierr;
114 
115   PetscFunctionBegin;
116   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
117   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
118   if (dim == 1) {
119     ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr);
120     PetscFunctionReturn(0);
121   }
122   if (!section)       {ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);}
123   if (!globalSection) {ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);}
124   ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr);
125   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
126   ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
127   ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
128   m    = (dim*(dim+1))/2;
129   ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr);
130   ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr);
131   ierr = VecSetUp(mode[0]);CHKERRQ(ierr);
132   for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);}
133   /* Assume P1 */
134   ierr = DMGetLocalVector(dm, &localMode);CHKERRQ(ierr);
135   for (d = 0; d < dim; ++d) {
136     PetscScalar values[3] = {0.0, 0.0, 0.0};
137 
138     values[d] = 1.0;
139     ierr      = VecSet(localMode, 0.0);CHKERRQ(ierr);
140     for (v = vStart; v < vEnd; ++v) {
141       ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
142     }
143     ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
144     ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
145   }
146   ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
147   for (d = dim; d < dim*(dim+1)/2; ++d) {
148     PetscInt i, j, k = dim > 2 ? d - dim : d;
149 
150     ierr = VecSet(localMode, 0.0);CHKERRQ(ierr);
151     for (v = vStart; v < vEnd; ++v) {
152       PetscScalar values[3] = {0.0, 0.0, 0.0};
153       PetscInt    off;
154 
155       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
156       for (i = 0; i < dim; ++i) {
157         for (j = 0; j < dim; ++j) {
158           values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]);
159         }
160       }
161       ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr);
162     }
163     ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
164     ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr);
165   }
166   ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
167   ierr = DMRestoreLocalVector(dm, &localMode);CHKERRQ(ierr);
168   for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);}
169   /* Orthonormalize system */
170   for (i = dim; i < m; ++i) {
171     PetscScalar dots[6];
172 
173     ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr);
174     for (j = 0; j < i; ++j) dots[j] *= -1.0;
175     ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr);
176     ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);
177   }
178   ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr);
179   for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);}
180   PetscFunctionReturn(0);
181 }
182 
183 #undef __FUNCT__
184 #define __FUNCT__ "DMPlexProjectFunctionLocal"
185 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, 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, c, f, d, v;
192   PetscErrorCode  ierr;
193 
194   PetscFunctionBegin;
195   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
196   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
197   ierr = PetscMalloc1(numFields, &sp);CHKERRQ(ierr);
198   for (f = 0; f < numFields; ++f) {
199     ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr);
200     ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr);
201     ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
202     totDim += spDim*numComp;
203   }
204   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
205   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
206   ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr);
207   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
208   ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
209   ierr = PetscMalloc2(dim,&v0,dim*dim,&J);CHKERRQ(ierr);
210   for (c = cStart; c < cEnd; ++c) {
211     PetscCellGeometry geom;
212 
213     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr);
214     geom.v0   = v0;
215     geom.J    = J;
216     geom.detJ = &detJ;
217     for (f = 0, v = 0; f < numFields; ++f) {
218       void * const ctx = ctxs ? ctxs[f] : NULL;
219       ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr);
220       ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
221       for (d = 0; d < spDim; ++d) {
222         ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr);
223         v += numComp;
224       }
225     }
226     ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr);
227   }
228   ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
229   ierr = PetscFree2(v0,J);CHKERRQ(ierr);
230   ierr = PetscFree(sp);CHKERRQ(ierr);
231   PetscFunctionReturn(0);
232 }
233 
234 #undef __FUNCT__
235 #define __FUNCT__ "DMPlexProjectFunction"
236 /*@C
237   DMPlexProjectFunction - This projects the given function into the function space provided.
238 
239   Input Parameters:
240 + dm      - The DM
241 . fe      - The PetscFE associated with the field
242 . funcs   - The coordinate functions to evaluate, one per field
243 . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
244 - mode    - The insertion mode for values
245 
246   Output Parameter:
247 . X - vector
248 
249   Level: developer
250 
251 .seealso: DMPlexComputeL2Diff()
252 @*/
253 PetscErrorCode DMPlexProjectFunction(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
254 {
255   Vec            localX;
256   PetscErrorCode ierr;
257 
258   PetscFunctionBegin;
259   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
260   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
261   ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, mode, localX);CHKERRQ(ierr);
262   ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr);
263   ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr);
264   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
265   PetscFunctionReturn(0);
266 }
267 
268 #undef __FUNCT__
269 #define __FUNCT__ "DMPlexComputeL2Diff"
270 /*@C
271   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
272 
273   Input Parameters:
274 + dm    - The DM
275 . fe    - The PetscFE object for each field
276 . funcs - The functions to evaluate for each field component
277 . ctxs  - Optional array of contexts to pass to each function, or NULL.
278 - X     - The coefficient vector u_h
279 
280   Output Parameter:
281 . diff - The diff ||u - u_h||_2
282 
283   Level: developer
284 
285 .seealso: DMPlexProjectFunction()
286 @*/
287 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
288 {
289   const PetscInt  debug = 0;
290   PetscSection    section;
291   PetscQuadrature quad;
292   Vec             localX;
293   PetscScalar    *funcVal;
294   PetscReal      *coords, *v0, *J, *invJ, detJ;
295   PetscReal       localDiff = 0.0;
296   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
297   PetscErrorCode  ierr;
298 
299   PetscFunctionBegin;
300   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
301   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
302   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
303   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
304   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
305   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
306   for (field = 0; field < numFields; ++field) {
307     PetscInt Nc;
308 
309     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
310     numComponents += Nc;
311   }
312   ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
313   ierr = PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
314   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
315   ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr);
316   for (c = cStart; c < cEnd; ++c) {
317     PetscScalar *x = NULL;
318     PetscReal    elemDiff = 0.0;
319 
320     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
321     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
322     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
323 
324     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
325       void * const     ctx = ctxs ? ctxs[field] : NULL;
326       const PetscReal *quadPoints, *quadWeights;
327       PetscReal       *basis;
328       PetscInt         numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f;
329 
330       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
331       ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr);
332       ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr);
333       ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr);
334       if (debug) {
335         char title[1024];
336         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
337         ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr);
338       }
339       for (q = 0; q < numQuadPoints; ++q) {
340         for (d = 0; d < dim; d++) {
341           coords[d] = v0[d];
342           for (e = 0; e < dim; e++) {
343             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
344           }
345         }
346         (*funcs[field])(coords, funcVal, ctx);
347         for (fc = 0; fc < numBasisComps; ++fc) {
348           PetscScalar interpolant = 0.0;
349 
350           for (f = 0; f < numBasisFuncs; ++f) {
351             const PetscInt fidx = f*numBasisComps+fc;
352             interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
353           }
354           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);}
355           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
356         }
357       }
358       comp        += numBasisComps;
359       fieldOffset += numBasisFuncs*numBasisComps;
360     }
361     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
362     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
363     localDiff += elemDiff;
364   }
365   ierr  = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr);
366   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
367   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
368   *diff = PetscSqrtReal(*diff);
369   PetscFunctionReturn(0);
370 }
371 
372 #undef __FUNCT__
373 #define __FUNCT__ "DMPlexComputeL2GradientDiff"
374 /*@C
375   DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
376 
377   Input Parameters:
378 + dm    - The DM
379 . fe    - The PetscFE object for each field
380 . funcs - The gradient functions to evaluate for each field component
381 . ctxs  - Optional array of contexts to pass to each function, or NULL.
382 . X     - The coefficient vector u_h
383 - n     - The vector to project along
384 
385   Output Parameter:
386 . diff - The diff ||(grad u - grad u_h) . n||_2
387 
388   Level: developer
389 
390 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
391 @*/
392 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
393 {
394   const PetscInt  debug = 0;
395   PetscSection    section;
396   PetscQuadrature quad;
397   Vec             localX;
398   PetscScalar    *funcVal, *interpolantVec;
399   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
400   PetscReal       localDiff = 0.0;
401   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
402   PetscErrorCode  ierr;
403 
404   PetscFunctionBegin;
405   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
406   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
407   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
408   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
409   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
410   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
411   for (field = 0; field < numFields; ++field) {
412     PetscInt Nc;
413 
414     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
415     numComponents += Nc;
416   }
417   /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */
418   ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr);
419   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
420   ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr);
421   for (c = cStart; c < cEnd; ++c) {
422     PetscScalar *x = NULL;
423     PetscReal    elemDiff = 0.0;
424 
425     ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr);
426     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
427     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
428 
429     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
430       void * const     ctx = ctxs ? ctxs[field] : NULL;
431       const PetscReal *quadPoints, *quadWeights;
432       PetscReal       *basisDer;
433       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;
434 
435       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
436       ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
437       ierr = PetscFEGetNumComponents(fe[field], &Ncomp);CHKERRQ(ierr);
438       ierr = PetscFEGetDefaultTabulation(fe[field], NULL, &basisDer, NULL);CHKERRQ(ierr);
439       if (debug) {
440         char title[1024];
441         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
442         ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr);
443       }
444       for (q = 0; q < numQuadPoints; ++q) {
445         for (d = 0; d < dim; d++) {
446           coords[d] = v0[d];
447           for (e = 0; e < dim; e++) {
448             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
449           }
450         }
451         (*funcs[field])(coords, n, funcVal, ctx);
452         for (fc = 0; fc < Ncomp; ++fc) {
453           PetscScalar interpolant = 0.0;
454 
455           for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
456           for (f = 0; f < Nb; ++f) {
457             const PetscInt fidx = f*Ncomp+fc;
458 
459             for (d = 0; d < dim; ++d) {
460               realSpaceDer[d] = 0.0;
461               for (g = 0; g < dim; ++g) {
462                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
463               }
464               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
465             }
466           }
467           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
468           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);}
469           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
470         }
471       }
472       comp        += Ncomp;
473       fieldOffset += Nb*Ncomp;
474     }
475     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
476     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
477     localDiff += elemDiff;
478   }
479     ierr  = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr);
480   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
481   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
482   *diff = PetscSqrtReal(*diff);
483   PetscFunctionReturn(0);
484 }
485 
486 #undef __FUNCT__
487 #define __FUNCT__ "DMPlexComputeResidualFEM"
488 /*@
489   DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
490 
491   Input Parameters:
492 + dm - The mesh
493 . X  - Local input vector
494 - user - The user context
495 
496   Output Parameter:
497 . F  - Local output vector
498 
499   Note:
500   The first member of the user context must be an FEMContext.
501 
502   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
503   like a GPU, or vectorize on a multicore machine.
504 
505   Level: developer
506 
507 .seealso: DMPlexComputeJacobianActionFEM()
508 @*/
509 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
510 {
511   DM_Plex          *mesh  = (DM_Plex *) dm->data;
512   PetscFEM         *fem   = (PetscFEM *) user;
513   PetscFE          *fe    = fem->fe;
514   PetscFE          *feAux = fem->feAux;
515   PetscFE          *feBd  = fem->feBd;
516   const char       *name  = "Residual";
517   DM                dmAux;
518   Vec               A;
519   PetscQuadrature   q;
520   PetscCellGeometry geom;
521   PetscSection      section, sectionAux;
522   PetscReal        *v0, *J, *invJ, *detJ;
523   PetscScalar      *elemVec, *u, *a = NULL;
524   PetscInt          dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c;
525   PetscInt          cellDof = 0, numComponents = 0;
526   PetscInt          cellDofAux = 0, numComponentsAux = 0;
527   PetscErrorCode    ierr;
528 
529   PetscFunctionBegin;
530   ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
531   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
532   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
533   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
534   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
535   numCells = cEnd - cStart;
536   for (f = 0; f < Nf; ++f) {
537     PetscInt Nb, Nc;
538 
539     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
540     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
541     cellDof       += Nb*Nc;
542     numComponents += Nc;
543   }
544   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
545   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
546   if (dmAux) {
547     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
548     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
549   }
550   for (f = 0; f < NfAux; ++f) {
551     PetscInt Nb, Nc;
552 
553     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
554     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
555     cellDofAux       += Nb*Nc;
556     numComponentsAux += Nc;
557   }
558   ierr = DMPlexProjectFunctionLocal(dm, fe, fem->bcFuncs, fem->bcCtxs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
559   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
560   ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr);
561   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
562   for (c = cStart; c < cEnd; ++c) {
563     PetscScalar *x = NULL;
564     PetscInt     i;
565 
566     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
567     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
568     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
569     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
570     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
571     if (dmAux) {
572       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
573       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
574       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
575     }
576   }
577   for (f = 0; f < Nf; ++f) {
578     void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[f];
579     void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[f];
580     PetscInt numQuadPoints, Nb;
581     /* Conforming batches */
582     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
583     /* Remainder */
584     PetscInt Nr, offset;
585 
586     ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr);
587     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
588     ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
589     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
590     blockSize = Nb*numQuadPoints;
591     batchSize = numBlocks * blockSize;
592     ierr =  PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
593     numChunks = numCells / (numBatches*batchSize);
594     Ne        = numChunks*numBatches*batchSize;
595     Nr        = numCells % (numBatches*batchSize);
596     offset    = numCells - Nr;
597     geom.v0   = v0;
598     geom.J    = J;
599     geom.invJ = invJ;
600     geom.detJ = detJ;
601     ierr = PetscFEIntegrateResidual(fe[f], Ne, Nf, fe, f, geom, u, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr);
602     geom.v0   = &v0[offset*dim];
603     geom.J    = &J[offset*dim*dim];
604     geom.invJ = &invJ[offset*dim*dim];
605     geom.detJ = &detJ[offset];
606     ierr = PetscFEIntegrateResidual(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
607   }
608   for (c = cStart; c < cEnd; ++c) {
609     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
610     ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
611   }
612   ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
613   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
614   if (feBd) {
615     DMLabel         label, depth;
616     IS              pointIS;
617     const PetscInt *points;
618     PetscInt        dep, numPoints, p, numFaces;
619     PetscReal      *n;
620 
621     ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr);
622     ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr);
623     ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr);
624     ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr);
625     ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
626     for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) {
627       PetscInt Nb, Nc;
628 
629       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
630       ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr);
631       cellDof       += Nb*Nc;
632       numComponents += Nc;
633     }
634     for (p = 0, numFaces = 0; p < numPoints; ++p) {
635       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
636       if (dep == dim-1) ++numFaces;
637     }
638     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);
639     for (p = 0, f = 0; p < numPoints; ++p) {
640       const PetscInt point = points[p];
641       PetscScalar   *x     = NULL;
642       PetscInt       i;
643 
644       ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr);
645       if (dep != dim-1) continue;
646       ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr);
647       ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]);
648       if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point);
649       ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
650       for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i];
651       ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr);
652       ++f;
653     }
654     for (f = 0; f < Nf; ++f) {
655       void   (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[f];
656       void   (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[f];
657       PetscInt numQuadPoints, Nb;
658       /* Conforming batches */
659       PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
660       /* Remainder */
661       PetscInt Nr, offset;
662 
663       ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr);
664       ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr);
665       ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
666       ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
667       blockSize = Nb*numQuadPoints;
668       batchSize = numBlocks * blockSize;
669       ierr =  PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
670       numChunks = numFaces / (numBatches*batchSize);
671       Ne        = numChunks*numBatches*batchSize;
672       Nr        = numFaces % (numBatches*batchSize);
673       offset    = numFaces - Nr;
674       geom.v0   = v0;
675       geom.n    = n;
676       geom.J    = J;
677       geom.invJ = invJ;
678       geom.detJ = detJ;
679       ierr = PetscFEIntegrateBdResidual(feBd[f], Ne, Nf, feBd, f, geom, u, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr);
680       geom.v0   = &v0[offset*dim];
681       geom.n    = &n[offset*dim];
682       geom.J    = &J[offset*dim*dim];
683       geom.invJ = &invJ[offset*dim*dim];
684       geom.detJ = &detJ[offset];
685       ierr = PetscFEIntegrateBdResidual(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr);
686     }
687     for (p = 0, f = 0; p < numPoints; ++p) {
688       const PetscInt point = points[p];
689 
690       ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr);
691       if (dep != dim-1) continue;
692       if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);}
693       ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr);
694       ++f;
695     }
696     ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
697     ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
698     ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr);
699   }
700   if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);}
701   ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr);
702   PetscFunctionReturn(0);
703 }
704 
705 #undef __FUNCT__
706 #define __FUNCT__ "DMPlexComputeJacobianActionFEM"
707 /*@C
708   DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user
709 
710   Input Parameters:
711 + dm - The mesh
712 . J  - The Jacobian shell matrix
713 . X  - Local input vector
714 - user - The user context
715 
716   Output Parameter:
717 . F  - Local output vector
718 
719   Note:
720   The first member of the user context must be an FEMContext.
721 
722   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
723   like a GPU, or vectorize on a multicore machine.
724 
725   Level: developer
726 
727 .seealso: DMPlexComputeResidualFEM()
728 @*/
729 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
730 {
731   DM_Plex          *mesh = (DM_Plex *) dm->data;
732   PetscFEM         *fem  = (PetscFEM *) user;
733   PetscFE          *fe   = fem->fe;
734   PetscQuadrature   quad;
735   PetscCellGeometry geom;
736   PetscSection      section;
737   JacActionCtx     *jctx;
738   PetscReal        *v0, *J, *invJ, *detJ;
739   PetscScalar      *elemVec, *u, *a;
740   PetscInt          dim, numFields, field, numCells, cStart, cEnd, c;
741   PetscInt          cellDof = 0;
742   PetscErrorCode    ierr;
743 
744   PetscFunctionBegin;
745   /* ierr = PetscLogEventBegin(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */
746   ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
747   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
748   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
749   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
750   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
751   numCells = cEnd - cStart;
752   for (field = 0; field < numFields; ++field) {
753     PetscInt Nb, Nc;
754 
755     ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
756     ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr);
757     cellDof += Nb*Nc;
758   }
759   ierr = VecSet(F, 0.0);CHKERRQ(ierr);
760   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);
761   for (c = cStart; c < cEnd; ++c) {
762     PetscScalar *x = NULL;
763     PetscInt     i;
764 
765     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
766     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
767     ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
768     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
769     ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr);
770     ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
771     for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
772     ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr);
773   }
774   for (field = 0; field < numFields; ++field) {
775     PetscInt numQuadPoints, Nb;
776     /* Conforming batches */
777     PetscInt numBlocks  = 1;
778     PetscInt numBatches = 1;
779     PetscInt numChunks, Ne, blockSize, batchSize;
780     /* Remainder */
781     PetscInt Nr, offset;
782 
783     ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr);
784     ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr);
785     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
786     blockSize = Nb*numQuadPoints;
787     batchSize = numBlocks * blockSize;
788     numChunks = numCells / (numBatches*batchSize);
789     Ne        = numChunks*numBatches*batchSize;
790     Nr        = numCells % (numBatches*batchSize);
791     offset    = numCells - Nr;
792     geom.v0   = v0;
793     geom.J    = J;
794     geom.invJ = invJ;
795     geom.detJ = detJ;
796     ierr = PetscFEIntegrateJacobianAction(fe[field], Ne, numFields, fe, field, geom, u, a, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr);
797     geom.v0   = &v0[offset*dim];
798     geom.J    = &J[offset*dim*dim];
799     geom.invJ = &invJ[offset*dim*dim];
800     geom.detJ = &detJ[offset];
801     ierr = PetscFEIntegrateJacobianAction(fe[field], Nr, numFields, fe, field, geom, &u[offset*cellDof], &a[offset*cellDof],
802                                           fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr);
803   }
804   for (c = cStart; c < cEnd; ++c) {
805     if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
806     ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
807   }
808   ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr);
809   if (mesh->printFEM) {
810     PetscMPIInt rank, numProcs;
811     PetscInt    p;
812 
813     ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
814     ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr);
815     ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");CHKERRQ(ierr);
816     for (p = 0; p < numProcs; ++p) {
817       if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);}
818       ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
819     }
820   }
821   /* ierr = PetscLogEventEnd(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */
822   PetscFunctionReturn(0);
823 }
824 
825 #undef __FUNCT__
826 #define __FUNCT__ "DMPlexComputeJacobianFEM"
827 /*@
828   DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
829 
830   Input Parameters:
831 + dm - The mesh
832 . X  - Local input vector
833 - user - The user context
834 
835   Output Parameter:
836 . Jac  - Jacobian matrix
837 
838   Note:
839   The first member of the user context must be an FEMContext.
840 
841   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
842   like a GPU, or vectorize on a multicore machine.
843 
844   Level: developer
845 
846 .seealso: FormFunctionLocal()
847 @*/
848 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
849 {
850   DM_Plex          *mesh  = (DM_Plex *) dm->data;
851   PetscFEM         *fem   = (PetscFEM *) user;
852   PetscFE          *fe    = fem->fe;
853   PetscFE          *feAux = fem->feAux;
854   const char       *name  = "Jacobian";
855   DM                dmAux;
856   Vec               A;
857   PetscQuadrature   quad;
858   PetscCellGeometry geom;
859   PetscSection      section, globalSection, sectionAux;
860   PetscReal        *v0, *J, *invJ, *detJ;
861   PetscScalar      *elemMat, *u, *a;
862   PetscInt          dim, Nf, NfAux = 0, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
863   PetscInt          cellDof = 0, numComponents = 0;
864   PetscInt          cellDofAux = 0, numComponentsAux = 0;
865   PetscBool         isShell;
866   PetscErrorCode    ierr;
867 
868   PetscFunctionBegin;
869   ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
870   ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr);
871   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
872   ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);
873   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
874   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
875   numCells = cEnd - cStart;
876   for (f = 0; f < Nf; ++f) {
877     PetscInt Nb, Nc;
878 
879     ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr);
880     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
881     cellDof       += Nb*Nc;
882     numComponents += Nc;
883   }
884   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
885   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
886   if (dmAux) {
887     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
888     ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr);
889   }
890   for (f = 0; f < NfAux; ++f) {
891     PetscInt Nb, Nc;
892 
893     ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr);
894     ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr);
895     cellDofAux       += Nb*Nc;
896     numComponentsAux += Nc;
897   }
898   ierr = DMPlexProjectFunctionLocal(dm, fe, fem->bcFuncs, fem->bcCtxs, INSERT_BC_VALUES, X);CHKERRQ(ierr);
899   ierr = MatZeroEntries(JacP);CHKERRQ(ierr);
900   ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof*cellDof,&elemMat);CHKERRQ(ierr);
901   if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);}
902   for (c = cStart; c < cEnd; ++c) {
903     PetscScalar *x = NULL;
904     PetscInt     i;
905 
906     ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
907     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
908     ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
909     for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
910     ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr);
911     if (dmAux) {
912       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
913       for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i];
914       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
915     }
916   }
917   ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr);
918   for (fieldI = 0; fieldI < Nf; ++fieldI) {
919     PetscInt numQuadPoints, Nb;
920     /* Conforming batches */
921     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
922     /* Remainder */
923     PetscInt Nr, offset;
924 
925     ierr = PetscFEGetQuadrature(fe[fieldI], &quad);CHKERRQ(ierr);
926     ierr = PetscFEGetDimension(fe[fieldI], &Nb);CHKERRQ(ierr);
927     ierr = PetscFEGetTileSizes(fe[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
928     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
929     blockSize = Nb*numQuadPoints;
930     batchSize = numBlocks * blockSize;
931     ierr = PetscFESetTileSizes(fe[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
932     numChunks = numCells / (numBatches*batchSize);
933     Ne        = numChunks*numBatches*batchSize;
934     Nr        = numCells % (numBatches*batchSize);
935     offset    = numCells - Nr;
936     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
937       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ];
938       void   (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*Nf+fieldJ];
939       void   (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*Nf+fieldJ];
940       void   (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*Nf+fieldJ];
941 
942       geom.v0   = v0;
943       geom.J    = J;
944       geom.invJ = invJ;
945       geom.detJ = detJ;
946       ierr = PetscFEIntegrateJacobian(fe[fieldI], Ne, Nf, fe, fieldI, fieldJ, geom, u, NfAux, feAux, a, g0, g1, g2, g3, elemMat);CHKERRQ(ierr);
947       geom.v0   = &v0[offset*dim];
948       geom.J    = &J[offset*dim*dim];
949       geom.invJ = &invJ[offset*dim*dim];
950       geom.detJ = &detJ[offset];
951       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);
952     }
953   }
954   for (c = cStart; c < cEnd; ++c) {
955     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);}
956     ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr);
957   }
958   ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
959   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
960   ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
961   ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
962   if (mesh->printFEM) {
963     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
964     ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr);
965     ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
966   }
967   ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr);
968   ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr);
969   if (isShell) {
970     JacActionCtx *jctx;
971 
972     ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr);
973     ierr = VecCopy(X, jctx->u);CHKERRQ(ierr);
974   }
975   PetscFunctionReturn(0);
976 }
977 
978 #if 0
979 
980 static void g0_identity_1d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
981 {
982   g3[0] = 1.0;
983 }
984 
985 static void g0_identity_2d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
986 {
987   g3[0] = g3[3] = 1.0;
988 }
989 
990 static void g0_identity_3d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[])
991 {
992   g3[0] = g3[4] = g3[8] = 1.0;
993 }
994 
995 #undef __FUNCT__
996 #define __FUNCT__ "DMPlexComputeInterpolatorFEMBroken"
997 PetscErrorCode DMPlexComputeInterpolatorFEMBroken(DM dmc, DM dmf, Mat I, void *user)
998 {
999   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1000   PetscFEM         *fem   = (PetscFEM *) user;
1001   PetscFE          *fe    = fem->fe;
1002   const char       *name  = "Interpolator";
1003   PetscFE          *feRef;
1004   PetscQuadrature   quad, quadOld;
1005   PetscCellGeometry geom;
1006   PetscSection      fsection, fglobalSection;
1007   PetscSection      csection, cglobalSection;
1008   PetscReal        *v0, *J, *invJ, *detJ;
1009   PetscScalar      *elemMat;
1010   PetscInt          dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c;
1011   PetscInt          rCellDof = 0, cCellDof = 0, numComponents = 0;
1012   PetscErrorCode    ierr;
1013 
1014   PetscFunctionBegin;
1015 #if 0
1016   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1017 #endif
1018   ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr);
1019   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1020   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1021   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1022   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1023   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1024   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1025   numCells = cEnd - cStart;
1026   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1027   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1028     ierr = PetscFERefine(fe[fieldI], &feRef[fieldI]);CHKERRQ(ierr);
1029   }
1030   for (f = 0; f < Nf; ++f) {
1031     PetscInt rNb, cNb, Nc;
1032 
1033     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1034     ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr);
1035     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1036     numComponents += Nc;
1037     rCellDof += rNb*Nc;
1038     cCellDof += cNb*Nc;
1039   }
1040   ierr = MatZeroEntries(I);CHKERRQ(ierr);
1041   ierr = PetscMalloc5(numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*rCellDof*cCellDof,&elemMat);CHKERRQ(ierr);
1042   for (c = cStart; c < cEnd; ++c) {
1043     ierr = DMPlexComputeCellGeometry(dmc, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr);
1044     if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
1045   }
1046   ierr = PetscMemzero(elemMat, numCells*rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1047   for (fieldI = 0; fieldI < Nf; ++fieldI) {
1048     PetscInt numQuadPoints, Nb;
1049     /* Conforming batches */
1050     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1051     /* Remainder */
1052     PetscInt Nr, offset;
1053 
1054     /* Make new fine FE which refines the ref cell and the quadrature rule */
1055     ierr = PetscFEGetQuadrature(feRef[fieldI], &quad);CHKERRQ(ierr);
1056     ierr = PetscFEGetDimension(feRef[fieldI], &Nb);CHKERRQ(ierr);
1057     ierr = PetscFEGetTileSizes(feRef[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1058     ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1059     blockSize = Nb*numQuadPoints;
1060     batchSize = numBlocks * blockSize;
1061     ierr = PetscFESetTileSizes(feRef[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1062     numChunks = numCells / (numBatches*batchSize);
1063     Ne        = numChunks*numBatches*batchSize;
1064     Nr        = numCells % (numBatches*batchSize);
1065     offset    = numCells - Nr;
1066 
1067     for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1068       /* void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ]; */
1069       void   (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = g0_identity_2d_static;
1070 
1071       /* Replace quadrature in coarse FE with refined quadrature */
1072       ierr = PetscFEGetQuadrature(fe[fieldJ], &quadOld);CHKERRQ(ierr);
1073       ierr = PetscObjectReference((PetscObject) quadOld);CHKERRQ(ierr);
1074       ierr = PetscFESetQuadrature(fe[fieldJ], quad);CHKERRQ(ierr);
1075       geom.v0   = v0;
1076       geom.J    = J;
1077       geom.invJ = invJ;
1078       geom.detJ = detJ;
1079       ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Ne, Nf, feRef, fieldI, fe, fieldJ, geom, g0, elemMat);CHKERRQ(ierr);
1080       geom.v0   = &v0[offset*dim];
1081       geom.J    = &J[offset*dim*dim];
1082       geom.invJ = &invJ[offset*dim*dim];
1083       geom.detJ = &detJ[offset];
1084       ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Nr, Nf, feRef, fieldI, fe, fieldJ, geom, g0, &elemMat[offset*rCellDof*cCellDof]);CHKERRQ(ierr);
1085       ierr = PetscFESetQuadrature(fe[fieldJ], quadOld);CHKERRQ(ierr);
1086     }
1087   }
1088   for (c = cStart; c < cEnd; ++c) {
1089     if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, rCellDof, cCellDof, &elemMat[c*rCellDof*cCellDof]);CHKERRQ(ierr);}
1090     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, I, c, &elemMat[c*rCellDof*cCellDof], ADD_VALUES);CHKERRQ(ierr);
1091   }
1092   ierr = PetscFree5(v0,J,invJ,detJ,elemMat);CHKERRQ(ierr);
1093   ierr = PetscFree(feRef);CHKERRQ(ierr);
1094   ierr = MatAssemblyBegin(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1095   ierr = MatAssemblyEnd(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1096   if (mesh->printFEM) {
1097     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1098     ierr = MatChop(I, 1.0e-10);CHKERRQ(ierr);
1099     ierr = MatView(I, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1100   }
1101 #if 0
1102   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1103 #endif
1104   PetscFunctionReturn(0);
1105 }
1106 #endif
1107 
1108 #undef __FUNCT__
1109 #define __FUNCT__ "DMPlexComputeInterpolatorFEM"
1110 /*@
1111   DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1112 
1113   Input Parameters:
1114 + dmf  - The fine mesh
1115 . dmc  - The coarse mesh
1116 - user - The user context
1117 
1118   Output Parameter:
1119 . In  - The interpolation matrix
1120 
1121   Note:
1122   The first member of the user context must be an FEMContext.
1123 
1124   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1125   like a GPU, or vectorize on a multicore machine.
1126 
1127   Level: developer
1128 
1129 .seealso: DMPlexComputeJacobianFEM()
1130 @*/
1131 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1132 {
1133   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1134   PetscFEM         *fem   = (PetscFEM *) user;
1135   PetscFE          *fe    = fem->fe;
1136   const char       *name  = "Interpolator";
1137   PetscFE          *feRef;
1138   PetscSection      fsection, fglobalSection;
1139   PetscSection      csection, cglobalSection;
1140   PetscScalar      *elemMat;
1141   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c;
1142   PetscInt          rCellDof = 0, cCellDof = 0;
1143   PetscErrorCode    ierr;
1144 
1145   PetscFunctionBegin;
1146 #if 0
1147   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1148 #endif
1149   ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr);
1150   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1151   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1152   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1153   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1154   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1155   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1156   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1157   for (f = 0; f < Nf; ++f) {
1158     PetscInt rNb, cNb, Nc;
1159 
1160     ierr = PetscFERefine(fe[f], &feRef[f]);CHKERRQ(ierr);
1161     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1162     ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr);
1163     ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr);
1164     rCellDof += rNb*Nc;
1165     cCellDof += cNb*Nc;
1166   }
1167   ierr = MatZeroEntries(In);CHKERRQ(ierr);
1168   ierr = PetscMalloc1(rCellDof*cCellDof,&elemMat);CHKERRQ(ierr);
1169   ierr = PetscMemzero(elemMat, rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr);
1170   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1171     PetscDualSpace   Qref;
1172     PetscQuadrature  f;
1173     const PetscReal *qpoints, *qweights;
1174     PetscReal       *points;
1175     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;
1176 
1177     /* Compose points from all dual basis functionals */
1178     ierr = PetscFEGetNumComponents(fe[fieldI], &Nc);CHKERRQ(ierr);
1179     ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr);
1180     ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr);
1181     for (i = 0; i < fpdim; ++i) {
1182       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1183       ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr);
1184       npoints += Np;
1185     }
1186     ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr);
1187     for (i = 0, k = 0; i < fpdim; ++i) {
1188       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1189       ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr);
1190       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1191     }
1192 
1193     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1194       PetscSpace P;
1195       PetscReal *B;
1196       PetscInt   NcJ, cpdim, j;
1197 
1198       /* For now, fields only interpolate themselves */
1199       if (fieldI != fieldJ) {offsetJ += cpdim; continue;}
1200       /* Evaluate basis at points */
1201       ierr = PetscFEGetBasisSpace(fe[fieldJ], &P);CHKERRQ(ierr);
1202       ierr = PetscFEGetNumComponents(fe[fieldJ], &NcJ);CHKERRQ(ierr);
1203       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);
1204       ierr = PetscFEGetDimension(fe[fieldJ], &cpdim);CHKERRQ(ierr);
1205       ierr = PetscFEGetTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);
1206       for (i = 0, k = 0; i < fpdim; ++i) {
1207         ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1208         ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr);
1209         for (p = 0; p < Np; ++p, ++k) {
1210           for (j = 0; j < cpdim; ++j) {
1211             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];
1212           }
1213         }
1214       }
1215       ierr = PetscFERestoreTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr);
1216       offsetJ += cpdim*NcJ;
1217     }
1218     offsetI += fpdim*Nc;
1219     ierr = PetscFree(points);CHKERRQ(ierr);
1220   }
1221   if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rCellDof, cCellDof, elemMat);CHKERRQ(ierr);}
1222   for (c = cStart; c < cEnd; ++c) {
1223     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr);
1224   }
1225   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);}
1226   ierr = PetscFree(feRef);CHKERRQ(ierr);
1227   ierr = PetscFree(elemMat);CHKERRQ(ierr);
1228   ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1229   ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1230   if (mesh->printFEM) {
1231     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1232     ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr);
1233     ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1234   }
1235 #if 0
1236   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1237 #endif
1238   PetscFunctionReturn(0);
1239 }
1240 
1241 #undef __FUNCT__
1242 #define __FUNCT__ "DMPlexAddBoundary"
1243 /* The ids can be overridden by the command line option -bc_<boundary name> */
1244 PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], PetscInt field, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx)
1245 {
1246   DM_Plex       *mesh = (DM_Plex *) dm->data;
1247   DMBoundary     b;
1248   PetscErrorCode ierr;
1249 
1250   PetscFunctionBegin;
1251   ierr = PetscNew(&b);CHKERRQ(ierr);
1252   ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr);
1253   ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr);
1254   ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr);
1255   b->essential   = isEssential;
1256   b->field       = field;
1257   b->func        = bcFunc;
1258   b->numids      = numids;
1259   b->ctx         = ctx;
1260   b->next        = mesh->boundary;
1261   mesh->boundary = b;
1262   PetscFunctionReturn(0);
1263 }
1264 
1265 #undef __FUNCT__
1266 #define __FUNCT__ "DMPlexGetNumBoundary"
1267 PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd)
1268 {
1269   DM_Plex   *mesh = (DM_Plex *) dm->data;
1270   DMBoundary b    = mesh->boundary;
1271 
1272   PetscFunctionBegin;
1273   *numBd = 0;
1274   while (b) {++(*numBd); b = b->next;}
1275   PetscFunctionReturn(0);
1276 }
1277 
1278 #undef __FUNCT__
1279 #define __FUNCT__ "DMPlexGetBoundary"
1280 PetscErrorCode DMPlexGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, PetscInt *field, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx)
1281 {
1282   DM_Plex   *mesh = (DM_Plex *) dm->data;
1283   DMBoundary b    = mesh->boundary;
1284   PetscInt   n    = 0;
1285 
1286   PetscFunctionBegin;
1287   while (b) {
1288     if (n == bd) break;
1289     b = b->next;
1290     ++n;
1291   }
1292   if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
1293   if (isEssential) {
1294     PetscValidPointer(isEssential, 3);
1295     *isEssential = b->essential;
1296   }
1297   if (name) {
1298     PetscValidPointer(name, 4);
1299     *name = b->name;
1300   }
1301   if (field) {
1302     PetscValidPointer(field, 5);
1303     *field = b->field;
1304   }
1305   if (func) {
1306     PetscValidPointer(func, 6);
1307     *func = b->func;
1308   }
1309   if (numids) {
1310     PetscValidPointer(numids, 7);
1311     *numids = b->numids;
1312   }
1313   if (ids) {
1314     PetscValidPointer(ids, 8);
1315     *ids = b->ids;
1316   }
1317   if (ctx) {
1318     PetscValidPointer(ctx, 9);
1319     *ctx = b->ctx;
1320   }
1321   PetscFunctionReturn(0);
1322 }
1323