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