xref: /petsc/src/dm/impls/plex/plexfem.c (revision 31383a9b1925af9a6e6adc8e5ad0d9cd5baf4d57) 
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__ "DMPlexSetMaxProjectionHeight"
186 /*@
187   DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs
188   are computed by associating the basis function with one of the mesh points in its transitively-closed support, and
189   evaluating the dual space basis of that point.  A basis function is associated with the point in its
190   transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum
191   projection height, which is set with this function.  By default, the maximum projection height is zero, which means
192   that only mesh cells are used to project basis functions.  A height of one, for example, evaluates a cell-interior
193   basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face.
194 
195   Input Parameters:
196 + dm - the DMPlex object
197 - height - the maximum projection height >= 0
198 
199   Level: advanced
200 
201 .seealso: DMPlexGetMaxProjectionHeight(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
202 @*/
203 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
204 {
205   DM_Plex *plex = (DM_Plex *) dm->data;
206 
207   PetscFunctionBegin;
208   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
209   plex->maxProjectionHeight = height;
210   PetscFunctionReturn(0);
211 }
212 
213 #undef __FUNCT__
214 #define __FUNCT__ "DMPlexGetMaxProjectionHeight"
215 /*@
216   DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in
217   DMPlexProjectXXXLocal() functions.
218 
219   Input Parameters:
220 . dm - the DMPlex object
221 
222   Output Parameters:
223 . height - the maximum projection height
224 
225   Level: intermediate
226 
227 .seealso: DMPlexSetMaxProjectionHeight(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
228 @*/
229 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
230 {
231   DM_Plex *plex = (DM_Plex *) dm->data;
232 
233   PetscFunctionBegin;
234   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
235   *height = plex->maxProjectionHeight;
236   PetscFunctionReturn(0);
237 }
238 
239 #undef __FUNCT__
240 #define __FUNCT__ "DMPlexProjectFunctionLabelLocal"
241 PetscErrorCode DMPlexProjectFunctionLabelLocal(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
242 {
243   PetscFE         fe;
244   PetscDualSpace *sp, *cellsp;
245   PetscSection    section;
246   PetscScalar    *values;
247   PetscBool      *fieldActive;
248   PetscInt        numFields, numComp, dim, dimEmbed, spDim, totDim, numValues, pStart, pEnd, f, d, v, i, comp, maxHeight, h;
249   PetscErrorCode  ierr;
250 
251   PetscFunctionBegin;
252   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
253   ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr);
254   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
255   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
256   ierr = PetscMalloc1(numFields,&sp);CHKERRQ(ierr);
257   ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr);
258   if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
259   if (maxHeight > 0) {
260     ierr = PetscMalloc1(numFields,&cellsp);CHKERRQ(ierr);
261   }
262   else {
263     cellsp = sp;
264   }
265   for (h = 0; h <= maxHeight; h++) {
266     ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr);
267     if (pEnd <= pStart) continue;
268     totDim = 0;
269     for (f = 0; f < numFields; ++f) {
270       ierr = DMGetField(dm, f, (PetscObject *) &fe);CHKERRQ(ierr);
271       if (!h) {
272         ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr);
273         sp[f] = cellsp[f];
274       }
275       else {
276         ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr);
277         if (!sp[f]) continue;
278       }
279       ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr);
280       ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
281       totDim += spDim*numComp;
282     }
283     ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr);
284     if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
285     if (!totDim) continue;
286     ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
287     ierr = DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr);
288     for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
289     for (i = 0; i < numIds; ++i) {
290       IS              pointIS;
291       const PetscInt *points;
292       PetscInt        n, p;
293 
294       ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr);
295       ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr);
296       ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
297       for (p = 0; p < n; ++p) {
298         const PetscInt    point = points[p];
299         PetscFECellGeom   geom;
300 
301         if ((point < pStart) || (point >= pEnd)) continue;
302         ierr          = DMPlexComputeCellGeometryFEM(dm, point, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr);
303         geom.dim      = dim - h;
304         geom.dimEmbed = dimEmbed;
305         for (f = 0, v = 0; f < numFields; ++f) {
306           void * const ctx = ctxs ? ctxs[f] : NULL;
307 
308           if (!sp[f]) continue;
309           ierr = DMGetField(dm, f, (PetscObject *) &fe);CHKERRQ(ierr);
310           ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr);
311           ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
312           for (d = 0; d < spDim; ++d) {
313             if (funcs[f]) {
314               ierr = PetscDualSpaceApply(sp[f], d, &geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr);
315             } else {
316               for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0;
317             }
318             v += numComp;
319           }
320         }
321         ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr);
322       }
323       ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
324       ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
325     }
326     if (h) {
327       for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);}
328     }
329   }
330   ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
331   ierr = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr);
332   ierr = PetscFree(sp);CHKERRQ(ierr);
333   if (maxHeight > 0) {
334     ierr = PetscFree(cellsp);CHKERRQ(ierr);
335   }
336   PetscFunctionReturn(0);
337 }
338 
339 #undef __FUNCT__
340 #define __FUNCT__ "DMPlexProjectFunctionLocal"
341 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
342 {
343   PetscDualSpace *sp, *cellsp;
344   PetscInt       *numComp;
345   PetscSection    section;
346   PetscScalar    *values;
347   PetscInt        numFields, dim, dimEmbed, spDim, totDim, numValues, pStart, pEnd, p, f, d, v, comp, h, maxHeight;
348   PetscErrorCode  ierr;
349 
350   PetscFunctionBegin;
351   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
352   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
353   ierr = PetscMalloc2(numFields, &sp, numFields, &numComp);CHKERRQ(ierr);
354   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
355   ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr);
356   ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr);
357   if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
358   if (maxHeight > 0) {
359     ierr = PetscMalloc1(numFields,&cellsp);CHKERRQ(ierr);
360   }
361   else {
362     cellsp = sp;
363   }
364   for (h = 0; h <= maxHeight; h++) {
365     ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr);
366     if (pEnd <= pStart) continue;
367     totDim = 0;
368     for (f = 0; f < numFields; ++f) {
369       PetscObject  obj;
370       PetscClassId id;
371 
372       ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr);
373       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
374       if (id == PETSCFE_CLASSID) {
375         PetscFE fe = (PetscFE) obj;
376 
377         ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr);
378         if (!h) {
379           ierr  = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr);
380           sp[f] = cellsp[f];
381           ierr  = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr);
382         }
383         else {
384           ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr);
385           if (!sp[f]) {
386             continue;
387           }
388         }
389       } else if (id == PETSCFV_CLASSID) {
390         PetscFV         fv = (PetscFV) obj;
391         PetscQuadrature q;
392 
393         if (h) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, "Projection height > 0 not supported for finite volume");
394         ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr);
395         ierr = PetscFVGetQuadrature(fv, &q);CHKERRQ(ierr);
396         ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) fv), &sp[f]);CHKERRQ(ierr);
397         ierr = PetscDualSpaceSetDM(sp[f], dm);CHKERRQ(ierr);
398         ierr = PetscDualSpaceSetType(sp[f], PETSCDUALSPACESIMPLE);CHKERRQ(ierr);
399         ierr = PetscDualSpaceSimpleSetDimension(sp[f], 1);CHKERRQ(ierr);
400         ierr = PetscDualSpaceSimpleSetFunctional(sp[f], 0, q);CHKERRQ(ierr);
401       } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
402       ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
403       totDim += spDim*numComp[f];
404     }
405     ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr);
406     if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
407     if (!totDim) continue;
408     ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
409     for (p = pStart; p < pEnd; ++p) {
410       PetscFECellGeom geom;
411 
412       ierr          = DMPlexComputeCellGeometryFEM(dm, p, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr);
413       geom.dim      = dim = h;
414       geom.dimEmbed = dimEmbed;
415       for (f = 0, v = 0; f < numFields; ++f) {
416         void * const ctx = ctxs ? ctxs[f] : NULL;
417 
418         if (!sp[f]) continue;
419         ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
420         for (d = 0; d < spDim; ++d) {
421           if (funcs[f]) {
422             ierr = PetscDualSpaceApply(sp[f], d, &geom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);
423           } else {
424             for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
425           }
426           v += numComp[f];
427         }
428       }
429       ierr = DMPlexVecSetClosure(dm, section, localX, p, values, mode);CHKERRQ(ierr);
430     }
431     ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
432     if (h || !maxHeight) {
433       for (f = 0; f < numFields; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);}
434     }
435   }
436   ierr = PetscFree2(sp, numComp);CHKERRQ(ierr);
437   if (maxHeight > 0) {
438     for (f = 0; f < numFields; f++) {ierr = PetscDualSpaceDestroy(&cellsp[f]);CHKERRQ(ierr);}
439     ierr = PetscFree(cellsp);CHKERRQ(ierr);
440   }
441   PetscFunctionReturn(0);
442 }
443 
444 #undef __FUNCT__
445 #define __FUNCT__ "DMPlexProjectFieldLocal"
446 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)
447 {
448   DM                dmAux;
449   PetscDS           prob, probAux;
450   Vec               A;
451   PetscSection      section, sectionAux;
452   PetscScalar      *values, *u, *u_x, *a, *a_x;
453   PetscReal        *x, *v0, *J, *invJ, detJ, **basisField, **basisFieldDer, **basisFieldAux, **basisFieldDerAux;
454   PetscInt          Nf, dim, spDim, totDim, numValues, cStart, cEnd, c, f, d, v, comp, maxHeight;
455   PetscErrorCode    ierr;
456 
457   PetscFunctionBegin;
458   ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr);
459   if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");}
460   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
461   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
462   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
463   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
464   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
465   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
466   ierr = PetscDSGetTabulation(prob, &basisField, &basisFieldDer);CHKERRQ(ierr);
467   ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr);
468   ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr);
469   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
470   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
471   if (dmAux) {
472     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
473     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
474     ierr = PetscDSGetTabulation(prob, &basisFieldAux, &basisFieldDerAux);CHKERRQ(ierr);
475     ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr);
476   }
477   ierr = DMPlexInsertBoundaryValues(dm, localU, 0.0, NULL, NULL, NULL);CHKERRQ(ierr);
478   ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr);
479   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
480   ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
481   ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
482   for (c = cStart; c < cEnd; ++c) {
483     PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
484 
485     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
486     ierr = DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr);
487     if (dmAux) {ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);}
488     for (f = 0, v = 0; f < Nf; ++f) {
489       PetscFE        fe;
490       PetscDualSpace sp;
491       PetscInt       Ncf;
492 
493       ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
494       ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr);
495       ierr = PetscFEGetNumComponents(fe, &Ncf);CHKERRQ(ierr);
496       ierr = PetscDualSpaceGetDimension(sp, &spDim);CHKERRQ(ierr);
497       for (d = 0; d < spDim; ++d) {
498         PetscQuadrature  quad;
499         const PetscReal *points, *weights;
500         PetscInt         numPoints, q;
501 
502         if (funcs[f]) {
503           ierr = PetscDualSpaceGetFunctional(sp, d, &quad);CHKERRQ(ierr);
504           ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr);
505           ierr = PetscFEGetTabulation(fe, numPoints, points, &basisField[f], &basisFieldDer[f], NULL);CHKERRQ(ierr);
506           for (q = 0; q < numPoints; ++q) {
507             CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x);
508             ierr = EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, coefficients,    NULL, u, u_x, NULL);CHKERRQ(ierr);
509             ierr = EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);CHKERRQ(ierr);
510             (*funcs[f])(u, NULL, u_x, a, NULL, a_x, x, &values[v]);
511           }
512           ierr = PetscFERestoreTabulation(fe, numPoints, points, &basisField[f], &basisFieldDer[f], NULL);CHKERRQ(ierr);
513         } else {
514           for (comp = 0; comp < Ncf; ++comp) values[v+comp] = 0.0;
515         }
516         v += Ncf;
517       }
518     }
519     ierr = DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr);
520     if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);}
521     ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr);
522   }
523   ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
524   ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr);
525   PetscFunctionReturn(0);
526 }
527 
528 #undef __FUNCT__
529 #define __FUNCT__ "DMPlexInsertBoundaryValues_FEM_Internal"
530 static PetscErrorCode DMPlexInsertBoundaryValues_FEM_Internal(DM dm, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], void (*func)(const PetscReal[], PetscScalar *, void *), void *ctx, Vec locX)
531 {
532   void        (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx);
533   void         **ctxs;
534   PetscInt       numFields;
535   PetscErrorCode ierr;
536 
537   PetscFunctionBegin;
538   ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr);
539   ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr);
540   funcs[field] = func;
541   ctxs[field]  = ctx;
542   ierr = DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr);
543   ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr);
544   PetscFunctionReturn(0);
545 }
546 
547 #undef __FUNCT__
548 #define __FUNCT__ "DMPlexInsertBoundaryValues_FVM_Internal"
549 static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad,
550                                                               PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
551 {
552   PetscFV            fv;
553   DM                 dmFace, dmCell, dmGrad;
554   const PetscScalar *facegeom, *cellgeom, *grad;
555   PetscScalar       *x, *fx;
556   PetscInt           dim, fStart, fEnd, pdim, i;
557   PetscErrorCode     ierr;
558 
559   PetscFunctionBegin;
560   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
561   ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
562   ierr = DMGetField(dm, field, (PetscObject *) &fv);CHKERRQ(ierr);
563   ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr);
564   ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr);
565   ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr);
566   ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);
567   if (Grad) {
568     ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr);
569     ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr);
570     ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr);
571     ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr);
572   }
573   ierr = VecGetArray(locX, &x);CHKERRQ(ierr);
574   for (i = 0; i < numids; ++i) {
575     IS              faceIS;
576     const PetscInt *faces;
577     PetscInt        numFaces, f;
578 
579     ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr);
580     if (!faceIS) continue; /* No points with that id on this process */
581     ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr);
582     ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr);
583     for (f = 0; f < numFaces; ++f) {
584       const PetscInt         face = faces[f], *cells;
585       const PetscFVFaceGeom *fg;
586 
587       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
588       ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr);
589       ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr);
590       if (Grad) {
591         const PetscFVCellGeom *cg;
592         const PetscScalar     *cx, *cgrad;
593         PetscScalar           *xG;
594         PetscReal              dx[3];
595         PetscInt               d;
596 
597         ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr);
598         ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr);
599         ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr);
600         ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr);
601         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
602         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
603         ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr);
604       } else {
605         const PetscScalar *xI;
606         PetscScalar       *xG;
607 
608         ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr);
609         ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr);
610         ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr);
611       }
612     }
613     ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr);
614     ierr = ISDestroy(&faceIS);CHKERRQ(ierr);
615   }
616   ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr);
617   if (Grad) {
618     ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr);
619     ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr);
620   }
621   ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);
622   ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr);
623   PetscFunctionReturn(0);
624 }
625 
626 #undef __FUNCT__
627 #define __FUNCT__ "DMPlexInsertBoundaryValues"
628 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
629 {
630   PetscInt       numBd, b;
631   PetscErrorCode ierr;
632 
633   PetscFunctionBegin;
634   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
635   PetscValidHeaderSpecific(locX, VEC_CLASSID, 2);
636   if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);}
637   if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);}
638   if (gradFVM)     {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);}
639   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
640   for (b = 0; b < numBd; ++b) {
641     PetscBool       isEssential;
642     const char     *labelname;
643     DMLabel         label;
644     PetscInt        field;
645     PetscObject     obj;
646     PetscClassId    id;
647     void          (*func)();
648     PetscInt        numids;
649     const PetscInt *ids;
650     void           *ctx;
651 
652     ierr = DMPlexGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, &func, &numids, &ids, &ctx);CHKERRQ(ierr);
653     if (!isEssential) continue;
654     ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr);
655     ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
656     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
657     if (id == PETSCFE_CLASSID) {
658       ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, field, label, numids, ids, (void (*)(const PetscReal[], PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr);
659     } else if (id == PETSCFV_CLASSID) {
660       ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM,
661                                                      field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr);
662     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
663   }
664   PetscFunctionReturn(0);
665 }
666 
667 #undef __FUNCT__
668 #define __FUNCT__ "DMPlexComputeL2Diff"
669 /*@C
670   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
671 
672   Input Parameters:
673 + dm    - The DM
674 . funcs - The functions to evaluate for each field component
675 . ctxs  - Optional array of contexts to pass to each function, or NULL.
676 - X     - The coefficient vector u_h
677 
678   Output Parameter:
679 . diff - The diff ||u - u_h||_2
680 
681   Level: developer
682 
683 .seealso: DMPlexProjectFunction(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff()
684 @*/
685 PetscErrorCode DMPlexComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
686 {
687   const PetscInt   debug = 0;
688   PetscSection     section;
689   PetscQuadrature  quad;
690   Vec              localX;
691   PetscScalar     *funcVal, *interpolant;
692   PetscReal       *coords, *v0, *J, *invJ, detJ;
693   PetscReal        localDiff = 0.0;
694   const PetscReal *quadPoints, *quadWeights;
695   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, c, field, fieldOffset;
696   PetscErrorCode   ierr;
697 
698   PetscFunctionBegin;
699   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
700   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
701   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
702   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
703   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
704   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
705   for (field = 0; field < numFields; ++field) {
706     PetscObject  obj;
707     PetscClassId id;
708     PetscInt     Nc;
709 
710     ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
711     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
712     if (id == PETSCFE_CLASSID) {
713       PetscFE fe = (PetscFE) obj;
714 
715       ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
716       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
717     } else if (id == PETSCFV_CLASSID) {
718       PetscFV fv = (PetscFV) obj;
719 
720       ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr);
721       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
722     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
723     numComponents += Nc;
724   }
725   ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
726   ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
727   ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
728   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
729   for (c = cStart; c < cEnd; ++c) {
730     PetscScalar *x = NULL;
731     PetscReal    elemDiff = 0.0;
732 
733     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
734     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
735     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
736 
737     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
738       PetscObject  obj;
739       PetscClassId id;
740       void * const ctx = ctxs ? ctxs[field] : NULL;
741       PetscInt     Nb, Nc, q, fc;
742 
743       ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
744       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
745       if (id == PETSCFE_CLASSID)      {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);}
746       else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;}
747       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
748       if (debug) {
749         char title[1024];
750         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
751         ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr);
752       }
753       for (q = 0; q < numQuadPoints; ++q) {
754         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
755         (*funcs[field])(coords, funcVal, ctx);
756         if (id == PETSCFE_CLASSID)      {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
757         else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
758         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
759         for (fc = 0; fc < Nc; ++fc) {
760           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);}
761           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
762         }
763       }
764       fieldOffset += Nb*Nc;
765     }
766     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
767     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
768     localDiff += elemDiff;
769   }
770   ierr  = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr);
771   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
772   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
773   *diff = PetscSqrtReal(*diff);
774   PetscFunctionReturn(0);
775 }
776 
777 #undef __FUNCT__
778 #define __FUNCT__ "DMPlexComputeL2GradientDiff"
779 /*@C
780   DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
781 
782   Input Parameters:
783 + dm    - The DM
784 . funcs - The gradient functions to evaluate for each field component
785 . ctxs  - Optional array of contexts to pass to each function, or NULL.
786 . X     - The coefficient vector u_h
787 - n     - The vector to project along
788 
789   Output Parameter:
790 . diff - The diff ||(grad u - grad u_h) . n||_2
791 
792   Level: developer
793 
794 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
795 @*/
796 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
797 {
798   const PetscInt  debug = 0;
799   PetscSection    section;
800   PetscQuadrature quad;
801   Vec             localX;
802   PetscScalar    *funcVal, *interpolantVec;
803   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
804   PetscReal       localDiff = 0.0;
805   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
806   PetscErrorCode  ierr;
807 
808   PetscFunctionBegin;
809   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
810   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
811   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
812   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
813   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
814   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
815   for (field = 0; field < numFields; ++field) {
816     PetscFE  fe;
817     PetscInt Nc;
818 
819     ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr);
820     ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
821     ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
822     numComponents += Nc;
823   }
824   /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */
825   ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr);
826   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
827   for (c = cStart; c < cEnd; ++c) {
828     PetscScalar *x = NULL;
829     PetscReal    elemDiff = 0.0;
830 
831     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
832     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
833     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
834 
835     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
836       PetscFE          fe;
837       void * const     ctx = ctxs ? ctxs[field] : NULL;
838       const PetscReal *quadPoints, *quadWeights;
839       PetscReal       *basisDer;
840       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;
841 
842       ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr);
843       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
844       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
845       ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr);
846       ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr);
847       if (debug) {
848         char title[1024];
849         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
850         ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr);
851       }
852       for (q = 0; q < numQuadPoints; ++q) {
853         for (d = 0; d < dim; d++) {
854           coords[d] = v0[d];
855           for (e = 0; e < dim; e++) {
856             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
857           }
858         }
859         (*funcs[field])(coords, n, funcVal, ctx);
860         for (fc = 0; fc < Ncomp; ++fc) {
861           PetscScalar interpolant = 0.0;
862 
863           for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
864           for (f = 0; f < Nb; ++f) {
865             const PetscInt fidx = f*Ncomp+fc;
866 
867             for (d = 0; d < dim; ++d) {
868               realSpaceDer[d] = 0.0;
869               for (g = 0; g < dim; ++g) {
870                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
871               }
872               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
873             }
874           }
875           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
876           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);}
877           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
878         }
879       }
880       comp        += Ncomp;
881       fieldOffset += Nb*Ncomp;
882     }
883     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
884     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
885     localDiff += elemDiff;
886   }
887   ierr  = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr);
888   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
889   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
890   *diff = PetscSqrtReal(*diff);
891   PetscFunctionReturn(0);
892 }
893 
894 #undef __FUNCT__
895 #define __FUNCT__ "DMPlexComputeL2FieldDiff"
896 /*@C
897   DMPlexComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
898 
899   Input Parameters:
900 + dm    - The DM
901 . funcs - The functions to evaluate for each field component
902 . ctxs  - Optional array of contexts to pass to each function, or NULL.
903 - X     - The coefficient vector u_h
904 
905   Output Parameter:
906 . diff - The array of differences, ||u^f - u^f_h||_2
907 
908   Level: developer
909 
910 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2GradientDiff()
911 @*/
912 PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
913 {
914   const PetscInt   debug = 0;
915   PetscSection     section;
916   PetscQuadrature  quad;
917   Vec              localX;
918   PetscScalar     *funcVal, *interpolant;
919   PetscReal       *coords, *v0, *J, *invJ, detJ;
920   PetscReal       *localDiff;
921   const PetscReal *quadPoints, *quadWeights;
922   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, c, field, fieldOffset;
923   PetscErrorCode   ierr;
924 
925   PetscFunctionBegin;
926   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
927   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
928   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
929   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
930   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
931   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
932   for (field = 0; field < numFields; ++field) {
933     PetscObject  obj;
934     PetscClassId id;
935     PetscInt     Nc;
936 
937     ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
938     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
939     if (id == PETSCFE_CLASSID) {
940       PetscFE fe = (PetscFE) obj;
941 
942       ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
943       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
944     } else if (id == PETSCFV_CLASSID) {
945       PetscFV fv = (PetscFV) obj;
946 
947       ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr);
948       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
949     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
950     numComponents += Nc;
951   }
952   ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
953   ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
954   ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
955   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
956   for (c = cStart; c < cEnd; ++c) {
957     PetscScalar *x = NULL;
958 
959     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
960     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
961     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
962 
963     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
964       PetscObject  obj;
965       PetscClassId id;
966       void * const ctx = ctxs ? ctxs[field] : NULL;
967       PetscInt     Nb, Nc, q, fc;
968 
969       PetscReal       elemDiff = 0.0;
970 
971       ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
972       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
973       if (id == PETSCFE_CLASSID)      {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);}
974       else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;}
975       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
976       if (debug) {
977         char title[1024];
978         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
979         ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr);
980       }
981       for (q = 0; q < numQuadPoints; ++q) {
982         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
983         (*funcs[field])(coords, funcVal, ctx);
984         if (id == PETSCFE_CLASSID)      {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
985         else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
986         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
987         for (fc = 0; fc < Nc; ++fc) {
988           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);}
989           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
990         }
991       }
992       fieldOffset += Nb*Nc;
993       localDiff[field] += elemDiff;
994     }
995     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
996   }
997   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
998   ierr = MPI_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
999   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
1000   ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr);
1001   PetscFunctionReturn(0);
1002 }
1003 
1004 #undef __FUNCT__
1005 #define __FUNCT__ "DMPlexComputeIntegralFEM"
1006 /*@
1007   DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user
1008 
1009   Input Parameters:
1010 + dm - The mesh
1011 . X  - Local input vector
1012 - user - The user context
1013 
1014   Output Parameter:
1015 . integral - Local integral for each field
1016 
1017   Level: developer
1018 
1019 .seealso: DMPlexComputeResidualFEM()
1020 @*/
1021 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
1022 {
1023   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1024   DM                dmAux;
1025   Vec               localX, A;
1026   PetscDS           prob, probAux = NULL;
1027   PetscQuadrature   q;
1028   PetscSection      section, sectionAux;
1029   PetscFECellGeom  *cgeom;
1030   PetscScalar      *u, *a = NULL;
1031   PetscInt          dim, Nf, f, numCells, cStart, cEnd, c;
1032   PetscInt          totDim, totDimAux;
1033   PetscErrorCode    ierr;
1034 
1035   PetscFunctionBegin;
1036   /*ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);*/
1037   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
1038   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
1039   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
1040   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1041   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1042   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1043   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1044   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1045   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1046   numCells = cEnd - cStart;
1047   for (f = 0; f < Nf; ++f) {integral[f]    = 0.0;}
1048   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1049   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1050   if (dmAux) {
1051     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1052     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
1053     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
1054   }
1055   ierr = DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr);
1056   ierr = PetscMalloc2(numCells*totDim,&u,numCells,&cgeom);CHKERRQ(ierr);
1057   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
1058   for (c = cStart; c < cEnd; ++c) {
1059     PetscScalar *x = NULL;
1060     PetscInt     i;
1061 
1062     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr);
1063     if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c);
1064     ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr);
1065     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1066     ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr);
1067     if (dmAux) {
1068       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1069       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1070       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1071     }
1072   }
1073   for (f = 0; f < Nf; ++f) {
1074     PetscFE  fe;
1075     PetscInt numQuadPoints, Nb;
1076     /* Conforming batches */
1077     PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1078     /* Remainder */
1079     PetscInt Nr, offset;
1080 
1081     ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
1082     ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
1083     ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1084     ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1085     ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr);
1086     blockSize = Nb*numQuadPoints;
1087     batchSize = numBlocks * blockSize;
1088     ierr =  PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1089     numChunks = numCells / (numBatches*batchSize);
1090     Ne        = numChunks*numBatches*batchSize;
1091     Nr        = numCells % (numBatches*batchSize);
1092     offset    = numCells - Nr;
1093     ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, integral);CHKERRQ(ierr);
1094     ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], integral);CHKERRQ(ierr);
1095   }
1096   ierr = PetscFree2(u,cgeom);CHKERRQ(ierr);
1097   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1098   if (mesh->printFEM) {
1099     ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr);
1100     for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", integral[f]);CHKERRQ(ierr);}
1101     ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr);
1102   }
1103   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
1104   /* TODO: Allreduce for integral */
1105   /*ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);*/
1106   PetscFunctionReturn(0);
1107 }
1108 
1109 #undef __FUNCT__
1110 #define __FUNCT__ "DMPlexComputeInterpolatorFEM"
1111 /*@
1112   DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1113 
1114   Input Parameters:
1115 + dmf  - The fine mesh
1116 . dmc  - The coarse mesh
1117 - user - The user context
1118 
1119   Output Parameter:
1120 . In  - The interpolation matrix
1121 
1122   Note:
1123   The first member of the user context must be an FEMContext.
1124 
1125   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1126   like a GPU, or vectorize on a multicore machine.
1127 
1128   Level: developer
1129 
1130 .seealso: DMPlexComputeJacobianFEM()
1131 @*/
1132 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1133 {
1134   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1135   const char       *name  = "Interpolator";
1136   PetscDS           prob;
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          cTotDim, rTotDim = 0;
1143   PetscErrorCode    ierr;
1144 
1145   PetscFunctionBegin;
1146   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1147   ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr);
1148   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1149   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1150   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1151   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1152   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1153   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1154   ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr);
1155   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1156   for (f = 0; f < Nf; ++f) {
1157     PetscFE  fe;
1158     PetscInt rNb, Nc;
1159 
1160     ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
1161     ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr);
1162     ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1163     ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
1164     rTotDim += rNb*Nc;
1165   }
1166   ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr);
1167   ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr);
1168   ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr);
1169   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1170     PetscDualSpace   Qref;
1171     PetscQuadrature  f;
1172     const PetscReal *qpoints, *qweights;
1173     PetscReal       *points;
1174     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;
1175 
1176     /* Compose points from all dual basis functionals */
1177     ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr);
1178     ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr);
1179     ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr);
1180     for (i = 0; i < fpdim; ++i) {
1181       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1182       ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr);
1183       npoints += Np;
1184     }
1185     ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr);
1186     for (i = 0, k = 0; i < fpdim; ++i) {
1187       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1188       ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr);
1189       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1190     }
1191 
1192     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1193       PetscFE    fe;
1194       PetscReal *B;
1195       PetscInt   NcJ, cpdim, j;
1196 
1197       /* Evaluate basis at points */
1198       ierr = PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);CHKERRQ(ierr);
1199       ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr);
1200       ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr);
1201       /* For now, fields only interpolate themselves */
1202       if (fieldI == fieldJ) {
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 = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);
1205         for (i = 0, k = 0; i < fpdim; ++i) {
1206           ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1207           ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr);
1208           for (p = 0; p < Np; ++p, ++k) {
1209             for (j = 0; j < cpdim; ++j) {
1210               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];
1211             }
1212           }
1213         }
1214         ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr);
1215       }
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, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);}
1222   /* Preallocate matrix */
1223   {
1224     PetscHashJK ht;
1225     PetscLayout rLayout;
1226     PetscInt   *dnz, *onz, *cellCIndices, *cellFIndices;
1227     PetscInt    locRows, rStart, rEnd, cell, r;
1228 
1229     ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr);
1230     ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr);
1231     ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
1232     ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
1233     ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
1234     ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
1235     ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
1236     ierr = PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr);
1237     ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr);
1238     for (cell = cStart; cell < cEnd; ++cell) {
1239       ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr);
1240       for (r = 0; r < rTotDim; ++r) {
1241         PetscHashJKKey  key;
1242         PetscHashJKIter missing, iter;
1243 
1244         key.j = cellFIndices[r];
1245         if (key.j < 0) continue;
1246         for (c = 0; c < cTotDim; ++c) {
1247           key.k = cellCIndices[c];
1248           if (key.k < 0) continue;
1249           ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr);
1250           if (missing) {
1251             ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr);
1252             if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1253             else                                     ++onz[key.j-rStart];
1254           }
1255         }
1256       }
1257     }
1258     ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr);
1259     ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr);
1260     ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
1261     ierr = PetscFree4(dnz,onz,cellCIndices,cellFIndices);CHKERRQ(ierr);
1262   }
1263   /* Fill matrix */
1264   ierr = MatZeroEntries(In);CHKERRQ(ierr);
1265   for (c = cStart; c < cEnd; ++c) {
1266     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr);
1267   }
1268   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);}
1269   ierr = PetscFree(feRef);CHKERRQ(ierr);
1270   ierr = PetscFree(elemMat);CHKERRQ(ierr);
1271   ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1272   ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1273   if (mesh->printFEM) {
1274     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1275     ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr);
1276     ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1277   }
1278   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1279   PetscFunctionReturn(0);
1280 }
1281 
1282 #undef __FUNCT__
1283 #define __FUNCT__ "DMPlexComputeInjectorFEM"
1284 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1285 {
1286   PetscDS        prob;
1287   PetscFE       *feRef;
1288   Vec            fv, cv;
1289   IS             fis, cis;
1290   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1291   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1292   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, c, dim, d, startC, offsetC, offsetF, m;
1293   PetscErrorCode ierr;
1294 
1295   PetscFunctionBegin;
1296   ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1297   ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr);
1298   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1299   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1300   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1301   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1302   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1303   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1304   ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr);
1305   ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr);
1306   for (f = 0; f < Nf; ++f) {
1307     PetscFE  fe;
1308     PetscInt fNb, Nc;
1309 
1310     ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr);
1311     ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr);
1312     ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr);
1313     ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
1314     fTotDim += fNb*Nc;
1315   }
1316   ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr);
1317   ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr);
1318   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1319     PetscFE        feC;
1320     PetscDualSpace QF, QC;
1321     PetscInt       NcF, NcC, fpdim, cpdim;
1322 
1323     ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr);
1324     ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr);
1325     ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr);
1326     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);
1327     ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr);
1328     ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr);
1329     ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr);
1330     ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr);
1331     for (c = 0; c < cpdim; ++c) {
1332       PetscQuadrature  cfunc;
1333       const PetscReal *cqpoints;
1334       PetscInt         NpC;
1335 
1336       ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr);
1337       ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr);
1338       if (NpC != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1339       for (f = 0; f < fpdim; ++f) {
1340         PetscQuadrature  ffunc;
1341         const PetscReal *fqpoints;
1342         PetscReal        sum = 0.0;
1343         PetscInt         NpF, comp;
1344 
1345         ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr);
1346         ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr);
1347         if (NpC != NpF) continue;
1348         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1349         if (sum > 1.0e-9) continue;
1350         for (comp = 0; comp < NcC; ++comp) {
1351           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1352         }
1353         break;
1354       }
1355     }
1356     offsetC += cpdim*NcC;
1357     offsetF += fpdim*NcF;
1358   }
1359   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);}
1360   ierr = PetscFree(feRef);CHKERRQ(ierr);
1361 
1362   ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr);
1363   ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr);
1364   ierr = VecGetOwnershipRange(cv, &startC, NULL);CHKERRQ(ierr);
1365   ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr);
1366   ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr);
1367   ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr);
1368   ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr);
1369   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1370   for (c = cStart; c < cEnd; ++c) {
1371     ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr);
1372     for (d = 0; d < cTotDim; ++d) {
1373       if (cellCIndices[d] < 0) continue;
1374       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]]);
1375       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1376       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1377     }
1378   }
1379   ierr = PetscFree(cmap);CHKERRQ(ierr);
1380   ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr);
1381 
1382   ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr);
1383   ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr);
1384   ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr);
1385   ierr = ISDestroy(&cis);CHKERRQ(ierr);
1386   ierr = ISDestroy(&fis);CHKERRQ(ierr);
1387   ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr);
1388   ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr);
1389   ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1390   PetscFunctionReturn(0);
1391 }
1392