xref: /petsc/src/dm/impls/plex/plexfem.c (revision a9b180a608ab5b0c8e2e0dfcd432ec21ea26b82a) 
1 #include <petsc-private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 #include <petscsf.h>
3 
4 #include <petsc-private/petscfeimpl.h>
5 #include <petsc-private/petscfvimpl.h>
6 
7 #undef __FUNCT__
8 #define __FUNCT__ "DMPlexGetScale"
9 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
10 {
11   DM_Plex *mesh = (DM_Plex*) dm->data;
12 
13   PetscFunctionBegin;
14   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
15   PetscValidPointer(scale, 3);
16   *scale = mesh->scale[unit];
17   PetscFunctionReturn(0);
18 }
19 
20 #undef __FUNCT__
21 #define __FUNCT__ "DMPlexSetScale"
22 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
23 {
24   DM_Plex *mesh = (DM_Plex*) dm->data;
25 
26   PetscFunctionBegin;
27   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
28   mesh->scale[unit] = scale;
29   PetscFunctionReturn(0);
30 }
31 
32 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
33 {
34   switch (i) {
35   case 0:
36     switch (j) {
37     case 0: return 0;
38     case 1:
39       switch (k) {
40       case 0: return 0;
41       case 1: return 0;
42       case 2: return 1;
43       }
44     case 2:
45       switch (k) {
46       case 0: return 0;
47       case 1: return -1;
48       case 2: return 0;
49       }
50     }
51   case 1:
52     switch (j) {
53     case 0:
54       switch (k) {
55       case 0: return 0;
56       case 1: return 0;
57       case 2: return -1;
58       }
59     case 1: return 0;
60     case 2:
61       switch (k) {
62       case 0: return 1;
63       case 1: return 0;
64       case 2: return 0;
65       }
66     }
67   case 2:
68     switch (j) {
69     case 0:
70       switch (k) {
71       case 0: return 0;
72       case 1: return 1;
73       case 2: return 0;
74       }
75     case 1:
76       switch (k) {
77       case 0: return -1;
78       case 1: return 0;
79       case 2: return 0;
80       }
81     case 2: return 0;
82     }
83   }
84   return 0;
85 }
86 
87 #undef __FUNCT__
88 #define __FUNCT__ "DMPlexProjectRigidBody"
89 static void DMPlexProjectRigidBody(const PetscReal X[], PetscScalar *mode, void *ctx)
90 {
91   PetscInt       *ctxInt = (PetscInt *)ctx;
92   PetscInt       dim     = ctxInt[0];
93   PetscInt       d       = ctxInt[1];
94   PetscInt       i, j, k = dim > 2 ? d - dim : d;
95 
96   for (i = 0; i < dim; i++) mode[i] = 0.;
97   if (d < dim) {
98     mode[d] = 1.;
99   }
100   else {
101     for (i = 0; i < dim; i++) {
102       for (j = 0; j < dim; j++) {
103         mode[j] += epsilon(i, j, k)*X[i];
104       }
105     }
106   }
107 }
108 
109 #undef __FUNCT__
110 #define __FUNCT__ "DMPlexCreateRigidBody"
111 /*@C
112   DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates
113 
114   Collective on DM
115 
116   Input Arguments:
117 . dm - the DM
118 
119   Output Argument:
120 . sp - the null space
121 
122   Note: This is necessary to take account of Dirichlet conditions on the displacements
123 
124   Level: advanced
125 
126 .seealso: MatNullSpaceCreate()
127 @*/
128 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp)
129 {
130   MPI_Comm       comm;
131   Vec            mode[6];
132   PetscSection   section, globalSection;
133   PetscInt       dim, n, m, d, i, j;
134   PetscErrorCode ierr;
135 
136   PetscFunctionBegin;
137   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
138   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
139   if (dim == 1) {
140     ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr);
141     PetscFunctionReturn(0);
142   }
143   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
144   ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);
145   ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr);
146   m    = (dim*(dim+1))/2;
147   ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr);
148   ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr);
149   ierr = VecSetUp(mode[0]);CHKERRQ(ierr);
150   for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);}
151   for (d = 0; d < m; d++) {
152     PetscInt ctx[2];
153     void     *voidctx = (void *)(&ctx[0]);
154     void     (*func)(const PetscReal *,PetscScalar *,void *) = DMPlexProjectRigidBody;
155 
156     ctx[0] = dim;
157     ctx[1] = d;
158     ierr = DMPlexProjectFunction(dm,&func,&voidctx,INSERT_VALUES,mode[d]);CHKERRQ(ierr);
159   }
160   for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);}
161   /* Orthonormalize system */
162   for (i = dim; i < m; ++i) {
163     PetscScalar dots[6];
164 
165     ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr);
166     for (j = 0; j < i; ++j) dots[j] *= -1.0;
167     ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr);
168     ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);
169   }
170   ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr);
171   for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);}
172   PetscFunctionReturn(0);
173 }
174 
175 #undef __FUNCT__
176 #define __FUNCT__ "DMPlexSetMaxProjectionHeight"
177 /*@
178   DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs
179   are computed by associating the basis function with one of the mesh points in its transitively-closed support, and
180   evaluating the dual space basis of that point.  A basis function is associated with the point in its
181   transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum
182   projection height, which is set with this function.  By default, the maximum projection height is zero, which means
183   that only mesh cells are used to project basis functions.  A height of one, for example, evaluates a cell-interior
184   basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face.
185 
186   Input Parameters:
187 + dm - the DMPlex object
188 - height - the maximum projection height >= 0
189 
190   Level: advanced
191 
192 .seealso: DMPlexGetMaxProjectionHeight(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
193 @*/
194 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height)
195 {
196   DM_Plex *plex = (DM_Plex *) dm->data;
197 
198   PetscFunctionBegin;
199   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
200   plex->maxProjectionHeight = height;
201   PetscFunctionReturn(0);
202 }
203 
204 #undef __FUNCT__
205 #define __FUNCT__ "DMPlexGetMaxProjectionHeight"
206 /*@
207   DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in
208   DMPlexProjectXXXLocal() functions.
209 
210   Input Parameters:
211 . dm - the DMPlex object
212 
213   Output Parameters:
214 . height - the maximum projection height
215 
216   Level: intermediate
217 
218 .seealso: DMPlexSetMaxProjectionHeight(), DMPlexProjectFunctionLocal(), DMPlexProjectFunctionLabelLocal()
219 @*/
220 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height)
221 {
222   DM_Plex *plex = (DM_Plex *) dm->data;
223 
224   PetscFunctionBegin;
225   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
226   *height = plex->maxProjectionHeight;
227   PetscFunctionReturn(0);
228 }
229 
230 #undef __FUNCT__
231 #define __FUNCT__ "DMPlexProjectFunctionLabelLocal"
232 PetscErrorCode DMPlexProjectFunctionLabelLocal(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
233 {
234   PetscDualSpace *sp, *cellsp;
235   PetscInt       *numComp;
236   PetscSection    section;
237   PetscScalar    *values;
238   PetscBool      *fieldActive;
239   PetscInt        numFields, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, cStart, cEnd, cEndInterior, f, d, v, i, comp, maxHeight, h;
240   PetscErrorCode  ierr;
241 
242   PetscFunctionBegin;
243   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
244   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
245   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
246   if (cEnd <= cStart) PetscFunctionReturn(0);
247   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
248   ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr);
249   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
250   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
251   ierr = PetscMalloc2(numFields,&sp,numFields,&numComp);CHKERRQ(ierr);
252   ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr);
253   if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
254   if (maxHeight > 0) {ierr = PetscMalloc1(numFields,&cellsp);CHKERRQ(ierr);}
255   else               {cellsp = sp;}
256   for (h = 0; h <= maxHeight; h++) {
257     ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr);
258     if (!h) {pStart = cStart; pEnd = cEnd;}
259     if (pEnd <= pStart) continue;
260     totDim = 0;
261     for (f = 0; f < numFields; ++f) {
262       PetscObject  obj;
263       PetscClassId id;
264 
265       ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr);
266       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
267       if (id == PETSCFE_CLASSID) {
268         PetscFE fe = (PetscFE) obj;
269 
270         ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr);
271         if (!h) {
272           ierr = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr);
273           sp[f] = cellsp[f];
274           ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr);
275         } else {
276           ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr);
277           if (!sp[f]) continue;
278         }
279       } else if (id == PETSCFV_CLASSID) {
280         PetscFV fv = (PetscFV) obj;
281 
282         ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr);
283         ierr = PetscFVGetDualSpace(fv, &sp[f]);
284         ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr);
285       } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
286       ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
287       totDim += spDim*numComp[f];
288     }
289     ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr);
290     if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
291     if (!totDim) continue;
292     ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
293     ierr = DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr);
294     for (f = 0; f < numFields; ++f) fieldActive[f] = (funcs[f] && sp[f]) ? PETSC_TRUE : PETSC_FALSE;
295     for (i = 0; i < numIds; ++i) {
296       IS              pointIS;
297       const PetscInt *points;
298       PetscInt        n, p;
299 
300       ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr);
301       ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr);
302       ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr);
303       for (p = 0; p < n; ++p) {
304         const PetscInt    point = points[p];
305         PetscFECellGeom   geom;
306 
307         if ((point < pStart) || (point >= pEnd)) continue;
308         ierr          = DMPlexComputeCellGeometryFEM(dm, point, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr);
309         geom.dim      = dim - h;
310         geom.dimEmbed = dimEmbed;
311         for (f = 0, v = 0; f < numFields; ++f) {
312           void * const ctx = ctxs ? ctxs[f] : NULL;
313 
314           if (!sp[f]) continue;
315           ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
316           for (d = 0; d < spDim; ++d) {
317             if (funcs[f]) {
318               ierr = PetscDualSpaceApply(sp[f], d, &geom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);
319             } else {
320               for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
321             }
322             v += numComp[f];
323           }
324         }
325         ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr);
326       }
327       ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr);
328       ierr = ISDestroy(&pointIS);CHKERRQ(ierr);
329     }
330     if (h) {
331       for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);}
332     }
333   }
334   ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
335   ierr = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr);
336   for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);}
337   ierr = PetscFree2(sp, numComp);CHKERRQ(ierr);
338   if (maxHeight > 0) {
339     ierr = PetscFree(cellsp);CHKERRQ(ierr);
340   }
341   PetscFunctionReturn(0);
342 }
343 
344 #undef __FUNCT__
345 #define __FUNCT__ "DMPlexProjectFunctionLocal"
346 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
347 {
348   PetscDualSpace *sp, *cellsp;
349   PetscInt       *numComp;
350   PetscSection    section;
351   PetscScalar    *values;
352   PetscInt        Nf, dim, dimEmbed, spDim, totDim = 0, numValues, pStart, pEnd, p, cStart, cEnd, cEndInterior, f, d, v, comp, h, maxHeight;
353   PetscErrorCode  ierr;
354 
355   PetscFunctionBegin;
356   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
357   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
358   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
359   if (cEnd <= cStart) PetscFunctionReturn(0);
360   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
361   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
362   ierr = PetscMalloc2(Nf, &sp, Nf, &numComp);CHKERRQ(ierr);
363   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
364   ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr);
365   ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr);
366   if (maxHeight < 0 || maxHeight > dim) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"maximum projection height %d not in [0, %d)\n", maxHeight,dim);}
367   if (maxHeight > 0) {
368     ierr = PetscMalloc1(Nf,&cellsp);CHKERRQ(ierr);
369   }
370   else {
371     cellsp = sp;
372   }
373   for (h = 0; h <= maxHeight; h++) {
374     ierr = DMPlexGetHeightStratum(dm, h, &pStart, &pEnd);CHKERRQ(ierr);
375     if (!h) {pStart = cStart; pEnd = cEnd;}
376     if (pEnd <= pStart) continue;
377     totDim = 0;
378     for (f = 0; f < Nf; ++f) {
379       PetscObject  obj;
380       PetscClassId id;
381 
382       ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr);
383       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
384       if (id == PETSCFE_CLASSID) {
385         PetscFE fe = (PetscFE) obj;
386 
387         ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr);
388         if (!h) {
389           ierr  = PetscFEGetDualSpace(fe, &cellsp[f]);CHKERRQ(ierr);
390           sp[f] = cellsp[f];
391           ierr  = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr);
392         }
393         else {
394           ierr = PetscDualSpaceGetHeightSubspace(cellsp[f], h, &sp[f]);CHKERRQ(ierr);
395           if (!sp[f]) {
396             continue;
397           }
398         }
399       } else if (id == PETSCFV_CLASSID) {
400         PetscFV fv = (PetscFV) obj;
401 
402         if (h) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP, "Projection height > 0 not supported for finite volume");
403         ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr);
404         ierr = PetscFVGetDualSpace(fv, &sp[f]);
405         ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr);
406       } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
407       ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
408       totDim += spDim*numComp[f];
409     }
410     ierr = DMPlexVecGetClosure(dm, section, localX, pStart, &numValues, NULL);CHKERRQ(ierr);
411     if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section point closure size %d != dual space dimension %d", numValues, totDim);
412     if (!totDim) continue;
413     ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
414     for (p = pStart; p < pEnd; ++p) {
415       PetscFECellGeom geom;
416 
417       ierr          = DMPlexComputeCellGeometryFEM(dm, p, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr);
418       geom.dim      = dim - h;
419       geom.dimEmbed = dimEmbed;
420       for (f = 0, v = 0; f < Nf; ++f) {
421         void * const ctx = ctxs ? ctxs[f] : NULL;
422 
423         if (!sp[f]) continue;
424         ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
425         for (d = 0; d < spDim; ++d) {
426           if (funcs[f]) {
427             ierr = PetscDualSpaceApply(sp[f], d, &geom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr);
428           } else {
429             for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0;
430           }
431           v += numComp[f];
432         }
433       }
434       ierr = DMPlexVecSetClosure(dm, section, localX, p, values, mode);CHKERRQ(ierr);
435     }
436     ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
437     if (h || !maxHeight) {
438       for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);}
439     }
440   }
441   ierr = PetscFree2(sp, numComp);CHKERRQ(ierr);
442   if (maxHeight > 0) {
443     for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&cellsp[f]);CHKERRQ(ierr);}
444     ierr = PetscFree(cellsp);CHKERRQ(ierr);
445   }
446   PetscFunctionReturn(0);
447 }
448 
449 #undef __FUNCT__
450 #define __FUNCT__ "DMPlexProjectFunction"
451 /*@C
452   DMPlexProjectFunction - This projects the given function into the function space provided.
453 
454   Input Parameters:
455 + dm      - The DM
456 . funcs   - The coordinate functions to evaluate, one per field
457 . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
458 - mode    - The insertion mode for values
459 
460   Output Parameter:
461 . X - vector
462 
463   Level: developer
464 
465 .seealso: DMPlexComputeL2Diff()
466 @*/
467 PetscErrorCode DMPlexProjectFunction(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
468 {
469   Vec            localX;
470   PetscErrorCode ierr;
471 
472   PetscFunctionBegin;
473   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
474   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
475   ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, mode, localX);CHKERRQ(ierr);
476   ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr);
477   ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr);
478   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
479   PetscFunctionReturn(0);
480 }
481 
482 #undef __FUNCT__
483 #define __FUNCT__ "DMPlexProjectFieldLocal"
484 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)
485 {
486   DM                dmAux;
487   PetscDS           prob, probAux = NULL;
488   Vec               A;
489   PetscSection      section, sectionAux = NULL;
490   PetscDualSpace   *sp;
491   PetscInt         *Ncf;
492   PetscScalar      *values, *u, *u_x, *a, *a_x;
493   PetscReal        *x, *v0, *J, *invJ, detJ;
494   PetscInt          Nf, dim, spDim, totDim, numValues, cStart, cEnd, cEndInterior, c, f, d, v, comp, maxHeight;
495   PetscErrorCode    ierr;
496 
497   PetscFunctionBegin;
498   ierr = DMPlexGetMaxProjectionHeight(dm,&maxHeight);CHKERRQ(ierr);
499   if (maxHeight > 0) {SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Field projection for height > 0 not supported yet");}
500   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
501   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
502   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
503   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
504   ierr = PetscMalloc2(Nf, &sp, Nf, &Ncf);CHKERRQ(ierr);
505   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
506   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
507   ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr);
508   ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr);
509   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
510   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
511   if (dmAux) {
512     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
513     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
514     ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr);
515   }
516   ierr = DMPlexInsertBoundaryValues(dm, localU, 0.0, NULL, NULL, NULL);CHKERRQ(ierr);
517   ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr);
518   if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim);
519   ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
520   ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
521   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
522   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
523   for (c = cStart; c < cEnd; ++c) {
524     PetscScalar *coefficients = NULL, *coefficientsAux = NULL;
525 
526     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
527     ierr = DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr);
528     if (dmAux) {ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);}
529     for (f = 0, v = 0; f < Nf; ++f) {
530       PetscObject  obj;
531       PetscClassId id;
532 
533       ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);
534       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
535       if (id == PETSCFE_CLASSID) {
536         PetscFE fe = (PetscFE) obj;
537 
538         ierr = PetscFEGetDualSpace(fe, &sp[f]);CHKERRQ(ierr);
539         ierr  = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr);
540         ierr = PetscFEGetNumComponents(fe, &Ncf[f]);CHKERRQ(ierr);
541       } else if (id == PETSCFV_CLASSID) {
542         PetscFV fv = (PetscFV) obj;
543 
544         ierr = PetscFVGetNumComponents(fv, &Ncf[f]);CHKERRQ(ierr);
545         ierr = PetscFVGetDualSpace(fv, &sp[f]);
546         ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr);
547       }
548       ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr);
549       for (d = 0; d < spDim; ++d) {
550         PetscQuadrature  quad;
551         const PetscReal *points, *weights;
552         PetscInt         numPoints, q;
553 
554         if (funcs[f]) {
555           ierr = PetscDualSpaceGetFunctional(sp[f], d, &quad);CHKERRQ(ierr);
556           ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr);
557           for (q = 0; q < numPoints; ++q) {
558             CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x);
559             ierr = EvaluateFieldJets(prob,    PETSC_FALSE, q, invJ, coefficients,    NULL, u, u_x, NULL);CHKERRQ(ierr);
560             ierr = EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);CHKERRQ(ierr);
561             (*funcs[f])(u, NULL, u_x, a, NULL, a_x, x, &values[v]);
562           }
563         } else {
564           for (comp = 0; comp < Ncf[f]; ++comp) values[v+comp] = 0.0;
565         }
566         v += Ncf[f];
567       }
568     }
569     ierr = DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr);
570     if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);}
571     ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr);
572   }
573   ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr);
574   ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr);
575   for (f = 0; f < Nf; f++) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);}
576   ierr = PetscFree2(sp, Ncf);CHKERRQ(ierr);
577   PetscFunctionReturn(0);
578 }
579 
580 #undef __FUNCT__
581 #define __FUNCT__ "DMPlexInsertBoundaryValues_FEM_Internal"
582 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)
583 {
584   void        (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx);
585   void         **ctxs;
586   PetscInt       numFields;
587   PetscErrorCode ierr;
588 
589   PetscFunctionBegin;
590   ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr);
591   ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr);
592   funcs[field] = func;
593   ctxs[field]  = ctx;
594   ierr = DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr);
595   ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr);
596   PetscFunctionReturn(0);
597 }
598 
599 #undef __FUNCT__
600 #define __FUNCT__ "DMPlexInsertBoundaryValues_FVM_Internal"
601 static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad,
602                                                               PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
603 {
604   PetscDS            prob;
605   PetscSF            sf;
606   DM                 dmFace, dmCell, dmGrad;
607   const PetscScalar *facegeom, *cellgeom, *grad;
608   const PetscInt    *leaves;
609   PetscScalar       *x, *fx;
610   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
611   PetscErrorCode     ierr;
612 
613   PetscFunctionBegin;
614   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
615   ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr);
616   nleaves = PetscMax(0, nleaves);
617   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
618   ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
619   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
620   ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr);
621   ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr);
622   ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr);
623   ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);
624   if (Grad) {
625     PetscFV fv;
626 
627     ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr);
628     ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr);
629     ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr);
630     ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr);
631     ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr);
632   }
633   ierr = VecGetArray(locX, &x);CHKERRQ(ierr);
634   for (i = 0; i < numids; ++i) {
635     IS              faceIS;
636     const PetscInt *faces;
637     PetscInt        numFaces, f;
638 
639     ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr);
640     if (!faceIS) continue; /* No points with that id on this process */
641     ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr);
642     ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr);
643     for (f = 0; f < numFaces; ++f) {
644       const PetscInt         face = faces[f], *cells;
645       const PetscFVFaceGeom *fg;
646 
647       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
648       ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr);
649       if (loc >= 0) continue;
650       ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr);
651       ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr);
652       if (Grad) {
653         const PetscFVCellGeom *cg;
654         const PetscScalar     *cx, *cgrad;
655         PetscScalar           *xG;
656         PetscReal              dx[3];
657         PetscInt               d;
658 
659         ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr);
660         ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr);
661         ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr);
662         ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr);
663         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
664         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
665         ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr);
666       } else {
667         const PetscScalar *xI;
668         PetscScalar       *xG;
669 
670         ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr);
671         ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr);
672         ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr);
673       }
674     }
675     ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr);
676     ierr = ISDestroy(&faceIS);CHKERRQ(ierr);
677   }
678   ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr);
679   if (Grad) {
680     ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr);
681     ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr);
682   }
683   ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);
684   ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr);
685   PetscFunctionReturn(0);
686 }
687 
688 #undef __FUNCT__
689 #define __FUNCT__ "DMPlexInsertBoundaryValues"
690 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
691 {
692   PetscInt       numBd, b;
693   PetscErrorCode ierr;
694 
695   PetscFunctionBegin;
696   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
697   PetscValidHeaderSpecific(locX, VEC_CLASSID, 2);
698   if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);}
699   if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);}
700   if (gradFVM)     {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);}
701   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
702   for (b = 0; b < numBd; ++b) {
703     PetscBool       isEssential;
704     const char     *labelname;
705     DMLabel         label;
706     PetscInt        field;
707     PetscObject     obj;
708     PetscClassId    id;
709     void          (*func)();
710     PetscInt        numids;
711     const PetscInt *ids;
712     void           *ctx;
713 
714     ierr = DMPlexGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, &func, &numids, &ids, &ctx);CHKERRQ(ierr);
715     if (!isEssential) continue;
716     ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr);
717     ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
718     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
719     if (id == PETSCFE_CLASSID) {
720       ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, field, label, numids, ids, (void (*)(const PetscReal[], PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr);
721     } else if (id == PETSCFV_CLASSID) {
722       if (!faceGeomFVM) continue;
723       ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM,
724                                                      field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr);
725     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
726   }
727   PetscFunctionReturn(0);
728 }
729 
730 #undef __FUNCT__
731 #define __FUNCT__ "DMPlexComputeL2Diff"
732 /*@C
733   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
734 
735   Input Parameters:
736 + dm    - The DM
737 . funcs - The functions to evaluate for each field component
738 . ctxs  - Optional array of contexts to pass to each function, or NULL.
739 - X     - The coefficient vector u_h
740 
741   Output Parameter:
742 . diff - The diff ||u - u_h||_2
743 
744   Level: developer
745 
746 .seealso: DMPlexProjectFunction(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff()
747 @*/
748 PetscErrorCode DMPlexComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
749 {
750   const PetscInt   debug = 0;
751   PetscSection     section;
752   PetscQuadrature  quad;
753   Vec              localX;
754   PetscScalar     *funcVal, *interpolant;
755   PetscReal       *coords, *v0, *J, *invJ, detJ;
756   PetscReal        localDiff = 0.0;
757   const PetscReal *quadPoints, *quadWeights;
758   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
759   PetscErrorCode   ierr;
760 
761   PetscFunctionBegin;
762   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
763   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
764   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
765   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
766   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
767   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
768   for (field = 0; field < numFields; ++field) {
769     PetscObject  obj;
770     PetscClassId id;
771     PetscInt     Nc;
772 
773     ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
774     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
775     if (id == PETSCFE_CLASSID) {
776       PetscFE fe = (PetscFE) obj;
777 
778       ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
779       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
780     } else if (id == PETSCFV_CLASSID) {
781       PetscFV fv = (PetscFV) obj;
782 
783       ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr);
784       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
785     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
786     numComponents += Nc;
787   }
788   ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
789   ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
790   ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
791   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
792   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
793   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
794   for (c = cStart; c < cEnd; ++c) {
795     PetscScalar *x = NULL;
796     PetscReal    elemDiff = 0.0;
797 
798     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
799     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
800     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
801 
802     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
803       PetscObject  obj;
804       PetscClassId id;
805       void * const ctx = ctxs ? ctxs[field] : NULL;
806       PetscInt     Nb, Nc, q, fc;
807 
808       ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
809       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
810       if (id == PETSCFE_CLASSID)      {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);}
811       else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;}
812       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
813       if (debug) {
814         char title[1024];
815         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
816         ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr);
817       }
818       for (q = 0; q < numQuadPoints; ++q) {
819         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
820         (*funcs[field])(coords, funcVal, ctx);
821         if (id == PETSCFE_CLASSID)      {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
822         else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
823         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
824         for (fc = 0; fc < Nc; ++fc) {
825           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);}
826           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
827         }
828       }
829       fieldOffset += Nb*Nc;
830     }
831     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
832     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
833     localDiff += elemDiff;
834   }
835   ierr  = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr);
836   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
837   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
838   *diff = PetscSqrtReal(*diff);
839   PetscFunctionReturn(0);
840 }
841 
842 #undef __FUNCT__
843 #define __FUNCT__ "DMPlexComputeL2GradientDiff"
844 /*@C
845   DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
846 
847   Input Parameters:
848 + dm    - The DM
849 . funcs - The gradient functions to evaluate for each field component
850 . ctxs  - Optional array of contexts to pass to each function, or NULL.
851 . X     - The coefficient vector u_h
852 - n     - The vector to project along
853 
854   Output Parameter:
855 . diff - The diff ||(grad u - grad u_h) . n||_2
856 
857   Level: developer
858 
859 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
860 @*/
861 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
862 {
863   const PetscInt  debug = 0;
864   PetscSection    section;
865   PetscQuadrature quad;
866   Vec             localX;
867   PetscScalar    *funcVal, *interpolantVec;
868   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
869   PetscReal       localDiff = 0.0;
870   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp;
871   PetscErrorCode  ierr;
872 
873   PetscFunctionBegin;
874   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
875   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
876   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
877   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
878   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
879   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
880   for (field = 0; field < numFields; ++field) {
881     PetscFE  fe;
882     PetscInt Nc;
883 
884     ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr);
885     ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
886     ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
887     numComponents += Nc;
888   }
889   /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */
890   ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr);
891   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
892   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
893   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
894   for (c = cStart; c < cEnd; ++c) {
895     PetscScalar *x = NULL;
896     PetscReal    elemDiff = 0.0;
897 
898     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
899     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
900     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
901 
902     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
903       PetscFE          fe;
904       void * const     ctx = ctxs ? ctxs[field] : NULL;
905       const PetscReal *quadPoints, *quadWeights;
906       PetscReal       *basisDer;
907       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;
908 
909       ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr);
910       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
911       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
912       ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr);
913       ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr);
914       if (debug) {
915         char title[1024];
916         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
917         ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr);
918       }
919       for (q = 0; q < numQuadPoints; ++q) {
920         for (d = 0; d < dim; d++) {
921           coords[d] = v0[d];
922           for (e = 0; e < dim; e++) {
923             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
924           }
925         }
926         (*funcs[field])(coords, n, funcVal, ctx);
927         for (fc = 0; fc < Ncomp; ++fc) {
928           PetscScalar interpolant = 0.0;
929 
930           for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
931           for (f = 0; f < Nb; ++f) {
932             const PetscInt fidx = f*Ncomp+fc;
933 
934             for (d = 0; d < dim; ++d) {
935               realSpaceDer[d] = 0.0;
936               for (g = 0; g < dim; ++g) {
937                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
938               }
939               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
940             }
941           }
942           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
943           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);}
944           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
945         }
946       }
947       comp        += Ncomp;
948       fieldOffset += Nb*Ncomp;
949     }
950     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
951     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
952     localDiff += elemDiff;
953   }
954   ierr  = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr);
955   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
956   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
957   *diff = PetscSqrtReal(*diff);
958   PetscFunctionReturn(0);
959 }
960 
961 #undef __FUNCT__
962 #define __FUNCT__ "DMPlexComputeL2FieldDiff"
963 /*@C
964   DMPlexComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
965 
966   Input Parameters:
967 + dm    - The DM
968 . funcs - The functions to evaluate for each field component
969 . ctxs  - Optional array of contexts to pass to each function, or NULL.
970 - X     - The coefficient vector u_h
971 
972   Output Parameter:
973 . diff - The array of differences, ||u^f - u^f_h||_2
974 
975   Level: developer
976 
977 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2GradientDiff()
978 @*/
979 PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
980 {
981   const PetscInt   debug = 0;
982   PetscSection     section;
983   PetscQuadrature  quad;
984   Vec              localX;
985   PetscScalar     *funcVal, *interpolant;
986   PetscReal       *coords, *v0, *J, *invJ, detJ;
987   PetscReal       *localDiff;
988   const PetscReal *quadPoints, *quadWeights;
989   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
990   PetscErrorCode   ierr;
991 
992   PetscFunctionBegin;
993   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
994   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
995   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
996   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
997   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
998   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
999   for (field = 0; field < numFields; ++field) {
1000     PetscObject  obj;
1001     PetscClassId id;
1002     PetscInt     Nc;
1003 
1004     ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
1005     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1006     if (id == PETSCFE_CLASSID) {
1007       PetscFE fe = (PetscFE) obj;
1008 
1009       ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
1010       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
1011     } else if (id == PETSCFV_CLASSID) {
1012       PetscFV fv = (PetscFV) obj;
1013 
1014       ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr);
1015       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
1016     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1017     numComponents += Nc;
1018   }
1019   ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
1020   ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
1021   ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
1022   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1023   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
1024   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1025   for (c = cStart; c < cEnd; ++c) {
1026     PetscScalar *x = NULL;
1027 
1028     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
1029     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1030     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
1031 
1032     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1033       PetscObject  obj;
1034       PetscClassId id;
1035       void * const ctx = ctxs ? ctxs[field] : NULL;
1036       PetscInt     Nb, Nc, q, fc;
1037 
1038       PetscReal       elemDiff = 0.0;
1039 
1040       ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
1041       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1042       if (id == PETSCFE_CLASSID)      {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);}
1043       else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;}
1044       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1045       if (debug) {
1046         char title[1024];
1047         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
1048         ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr);
1049       }
1050       for (q = 0; q < numQuadPoints; ++q) {
1051         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
1052         (*funcs[field])(coords, funcVal, ctx);
1053         if (id == PETSCFE_CLASSID)      {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
1054         else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
1055         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1056         for (fc = 0; fc < Nc; ++fc) {
1057           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);}
1058           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
1059         }
1060       }
1061       fieldOffset += Nb*Nc;
1062       localDiff[field] += elemDiff;
1063     }
1064     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
1065   }
1066   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
1067   ierr = MPI_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
1068   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
1069   ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr);
1070   PetscFunctionReturn(0);
1071 }
1072 
1073 #undef __FUNCT__
1074 #define __FUNCT__ "DMPlexComputeIntegralFEM"
1075 /*@
1076   DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user
1077 
1078   Input Parameters:
1079 + dm - The mesh
1080 . X  - Local input vector
1081 - user - The user context
1082 
1083   Output Parameter:
1084 . integral - Local integral for each field
1085 
1086   Level: developer
1087 
1088 .seealso: DMPlexComputeResidualFEM()
1089 @*/
1090 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
1091 {
1092   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1093   DM                dmAux;
1094   Vec               localX, A;
1095   PetscDS           prob, probAux = NULL;
1096   PetscSection      section, sectionAux;
1097   PetscFECellGeom  *cgeom;
1098   PetscScalar      *u, *a = NULL;
1099   PetscReal        *lintegral, *vol;
1100   PetscInt          dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c;
1101   PetscInt          totDim, totDimAux;
1102   PetscErrorCode    ierr;
1103 
1104   PetscFunctionBegin;
1105   ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);
1106   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
1107   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
1108   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
1109   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1110   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1111   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1112   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1113   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1114   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1115   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
1116   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1117   numCells = cEnd - cStart;
1118   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1119   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1120   if (dmAux) {
1121     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1122     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
1123     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
1124   }
1125   ierr = DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr);
1126   ierr = PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);CHKERRQ(ierr);
1127   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
1128   for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
1129   for (c = cStart; c < cEnd; ++c) {
1130     PetscScalar *x = NULL;
1131     PetscInt     i;
1132 
1133     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr);
1134     ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);CHKERRQ(ierr);
1135     if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c);
1136     ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr);
1137     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1138     ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr);
1139     if (dmAux) {
1140       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1141       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1142       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1143     }
1144   }
1145   for (f = 0; f < Nf; ++f) {
1146     PetscObject  obj;
1147     PetscClassId id;
1148     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1149 
1150     ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);
1151     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1152     if (id == PETSCFE_CLASSID) {
1153       PetscFE         fe = (PetscFE) obj;
1154       PetscQuadrature q;
1155       PetscInt        Nq, Nb;
1156 
1157       ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1158       ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
1159       ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr);
1160       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1161       blockSize = Nb*Nq;
1162       batchSize = numBlocks * blockSize;
1163       ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1164       numChunks = numCells / (numBatches*batchSize);
1165       Ne        = numChunks*numBatches*batchSize;
1166       Nr        = numCells % (numBatches*batchSize);
1167       offset    = numCells - Nr;
1168       ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);CHKERRQ(ierr);
1169       ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr);
1170     } else if (id == PETSCFV_CLASSID) {
1171       /* PetscFV  fv = (PetscFV) obj; */
1172       PetscInt    foff;
1173       void      (*obj_func)(const PetscScalar u[], const PetscScalar u_x[], const PetscScalar u_t[], const PetscScalar a[], const PetscScalar a_x[], const PetscScalar a_t[], const PetscReal x[], PetscScalar obj[]);
1174       PetscScalar lint;
1175 
1176       ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr);
1177       ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr);
1178       if (obj_func) {
1179         for (c = 0; c < numCells; ++c) {
1180           /* TODO: Need full pointwise interpolation and get centroid for x argument */
1181           obj_func(&u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, &lint);
1182           lintegral[f] = PetscRealPart(lint)*vol[c];
1183         }
1184       }
1185     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1186   }
1187   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1188   if (mesh->printFEM) {
1189     ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr);
1190     for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);}
1191     ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr);
1192   }
1193   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
1194   ierr = MPI_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr);
1195   ierr = PetscFree4(lintegral,u,cgeom,vol);CHKERRQ(ierr);
1196   ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);
1197   PetscFunctionReturn(0);
1198 }
1199 
1200 #undef __FUNCT__
1201 #define __FUNCT__ "DMPlexComputeInterpolatorFEM"
1202 /*@
1203   DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1204 
1205   Input Parameters:
1206 + dmf  - The fine mesh
1207 . dmc  - The coarse mesh
1208 - user - The user context
1209 
1210   Output Parameter:
1211 . In  - The interpolation matrix
1212 
1213   Note:
1214   The first member of the user context must be an FEMContext.
1215 
1216   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1217   like a GPU, or vectorize on a multicore machine.
1218 
1219   Level: developer
1220 
1221 .seealso: DMPlexComputeJacobianFEM()
1222 @*/
1223 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1224 {
1225   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1226   const char       *name  = "Interpolator";
1227   PetscDS           prob;
1228   PetscFE          *feRef;
1229   PetscFV          *fvRef;
1230   PetscSection      fsection, fglobalSection;
1231   PetscSection      csection, cglobalSection;
1232   PetscScalar      *elemMat;
1233   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1234   PetscInt          cTotDim, rTotDim = 0;
1235   PetscErrorCode    ierr;
1236 
1237   PetscFunctionBegin;
1238   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1239   ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr);
1240   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1241   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1242   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1243   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1244   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1245   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1246   ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
1247   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1248   ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr);
1249   ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr);
1250   for (f = 0; f < Nf; ++f) {
1251     PetscObject  obj;
1252     PetscClassId id;
1253     PetscInt     rNb = 0, Nc = 0;
1254 
1255     ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);
1256     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1257     if (id == PETSCFE_CLASSID) {
1258       PetscFE fe = (PetscFE) obj;
1259 
1260       ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr);
1261       ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1262       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
1263     } else if (id == PETSCFV_CLASSID) {
1264       PetscFV        fv = (PetscFV) obj;
1265       PetscDualSpace Q;
1266 
1267       ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr);
1268       ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr);
1269       ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr);
1270       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
1271     }
1272     rTotDim += rNb*Nc;
1273   }
1274   ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr);
1275   ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr);
1276   ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr);
1277   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1278     PetscDualSpace   Qref;
1279     PetscQuadrature  f;
1280     const PetscReal *qpoints, *qweights;
1281     PetscReal       *points;
1282     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;
1283 
1284     /* Compose points from all dual basis functionals */
1285     if (feRef[fieldI]) {
1286       ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr);
1287       ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr);
1288     } else {
1289       ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr);
1290       ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr);
1291     }
1292     ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr);
1293     for (i = 0; i < fpdim; ++i) {
1294       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1295       ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr);
1296       npoints += Np;
1297     }
1298     ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr);
1299     for (i = 0, k = 0; i < fpdim; ++i) {
1300       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1301       ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr);
1302       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1303     }
1304 
1305     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1306       PetscObject  obj;
1307       PetscClassId id;
1308       PetscReal   *B;
1309       PetscInt     NcJ = 0, cpdim = 0, j;
1310 
1311       ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr);
1312       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1313       if (id == PETSCFE_CLASSID) {
1314         PetscFE fe = (PetscFE) obj;
1315 
1316         /* Evaluate basis at points */
1317         ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr);
1318         ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr);
1319         /* For now, fields only interpolate themselves */
1320         if (fieldI == fieldJ) {
1321           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);
1322           ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);
1323           for (i = 0, k = 0; i < fpdim; ++i) {
1324             ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1325             ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr);
1326             for (p = 0; p < Np; ++p, ++k) {
1327               for (j = 0; j < cpdim; ++j) {
1328                 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];
1329               }
1330             }
1331           }
1332           ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr);
1333         }
1334       } else if (id == PETSCFV_CLASSID) {
1335         PetscFV        fv = (PetscFV) obj;
1336 
1337         /* Evaluate constant function at points */
1338         ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr);
1339         cpdim = 1;
1340         /* For now, fields only interpolate themselves */
1341         if (fieldI == fieldJ) {
1342           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);
1343           for (i = 0, k = 0; i < fpdim; ++i) {
1344             ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1345             ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr);
1346             for (p = 0; p < Np; ++p, ++k) {
1347               for (j = 0; j < cpdim; ++j) {
1348                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p];
1349               }
1350             }
1351           }
1352         }
1353       }
1354       offsetJ += cpdim*NcJ;
1355     }
1356     offsetI += fpdim*Nc;
1357     ierr = PetscFree(points);CHKERRQ(ierr);
1358   }
1359   if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);}
1360   /* Preallocate matrix */
1361   {
1362     PetscHashJK ht;
1363     PetscLayout rLayout;
1364     PetscInt   *dnz, *onz, *cellCIndices, *cellFIndices;
1365     PetscInt    locRows, rStart, rEnd, cell, r;
1366 
1367     ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr);
1368     ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr);
1369     ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
1370     ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
1371     ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
1372     ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
1373     ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
1374     ierr = PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr);
1375     ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr);
1376     for (cell = cStart; cell < cEnd; ++cell) {
1377       ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr);
1378       for (r = 0; r < rTotDim; ++r) {
1379         PetscHashJKKey  key;
1380         PetscHashJKIter missing, iter;
1381 
1382         key.j = cellFIndices[r];
1383         if (key.j < 0) continue;
1384         for (c = 0; c < cTotDim; ++c) {
1385           key.k = cellCIndices[c];
1386           if (key.k < 0) continue;
1387           ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr);
1388           if (missing) {
1389             ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr);
1390             if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1391             else                                     ++onz[key.j-rStart];
1392           }
1393         }
1394       }
1395     }
1396     ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr);
1397     ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr);
1398     ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
1399     ierr = PetscFree4(dnz,onz,cellCIndices,cellFIndices);CHKERRQ(ierr);
1400   }
1401   /* Fill matrix */
1402   ierr = MatZeroEntries(In);CHKERRQ(ierr);
1403   for (c = cStart; c < cEnd; ++c) {
1404     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr);
1405   }
1406   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);}
1407   ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr);
1408   ierr = PetscFree(elemMat);CHKERRQ(ierr);
1409   ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1410   ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1411   if (mesh->printFEM) {
1412     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1413     ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr);
1414     ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1415   }
1416   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1417   PetscFunctionReturn(0);
1418 }
1419 
1420 #undef __FUNCT__
1421 #define __FUNCT__ "DMPlexComputeInjectorFEM"
1422 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1423 {
1424   PetscDS        prob;
1425   PetscFE       *feRef;
1426   PetscFV       *fvRef;
1427   Vec            fv, cv;
1428   IS             fis, cis;
1429   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1430   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1431   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m;
1432   PetscErrorCode ierr;
1433 
1434   PetscFunctionBegin;
1435   ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1436   ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr);
1437   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1438   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1439   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1440   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1441   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1442   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1443   ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
1444   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1445   ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr);
1446   ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr);
1447   for (f = 0; f < Nf; ++f) {
1448     PetscObject  obj;
1449     PetscClassId id;
1450     PetscInt     fNb = 0, Nc = 0;
1451 
1452     ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);
1453     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1454     if (id == PETSCFE_CLASSID) {
1455       PetscFE fe = (PetscFE) obj;
1456 
1457       ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr);
1458       ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr);
1459       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
1460     } else if (id == PETSCFV_CLASSID) {
1461       PetscFV        fv = (PetscFV) obj;
1462       PetscDualSpace Q;
1463 
1464       ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr);
1465       ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr);
1466       ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr);
1467       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
1468     }
1469     fTotDim += fNb*Nc;
1470   }
1471   ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr);
1472   ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr);
1473   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1474     PetscFE        feC;
1475     PetscFV        fvC;
1476     PetscDualSpace QF, QC;
1477     PetscInt       NcF, NcC, fpdim, cpdim;
1478 
1479     if (feRef[field]) {
1480       ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr);
1481       ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr);
1482       ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr);
1483       ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr);
1484       ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr);
1485       ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr);
1486       ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr);
1487     } else {
1488       ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr);
1489       ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr);
1490       ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr);
1491       ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr);
1492       ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr);
1493       ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr);
1494       ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr);
1495     }
1496     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);
1497     for (c = 0; c < cpdim; ++c) {
1498       PetscQuadrature  cfunc;
1499       const PetscReal *cqpoints;
1500       PetscInt         NpC;
1501       PetscBool        found = PETSC_FALSE;
1502 
1503       ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr);
1504       ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr);
1505       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1506       for (f = 0; f < fpdim; ++f) {
1507         PetscQuadrature  ffunc;
1508         const PetscReal *fqpoints;
1509         PetscReal        sum = 0.0;
1510         PetscInt         NpF, comp;
1511 
1512         ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr);
1513         ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr);
1514         if (NpC != NpF) continue;
1515         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1516         if (sum > 1.0e-9) continue;
1517         for (comp = 0; comp < NcC; ++comp) {
1518           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1519         }
1520         found = PETSC_TRUE;
1521         break;
1522       }
1523       if (!found) {
1524         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1525         if (fvRef[field]) {
1526           PetscInt comp;
1527           for (comp = 0; comp < NcC; ++comp) {
1528             cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1529           }
1530         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1531       }
1532     }
1533     offsetC += cpdim*NcC;
1534     offsetF += fpdim*NcF;
1535   }
1536   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);}
1537   ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr);
1538 
1539   ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr);
1540   ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr);
1541   ierr = VecGetOwnershipRange(cv, &startC, NULL);CHKERRQ(ierr);
1542   ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr);
1543   ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr);
1544   ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr);
1545   ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr);
1546   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1547   for (c = cStart; c < cEnd; ++c) {
1548     ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr);
1549     for (d = 0; d < cTotDim; ++d) {
1550       if (cellCIndices[d] < 0) continue;
1551       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]]);
1552       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1553       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1554     }
1555   }
1556   ierr = PetscFree(cmap);CHKERRQ(ierr);
1557   ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr);
1558 
1559   ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr);
1560   ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr);
1561   ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr);
1562   ierr = ISDestroy(&cis);CHKERRQ(ierr);
1563   ierr = ISDestroy(&fis);CHKERRQ(ierr);
1564   ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr);
1565   ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr);
1566   ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1567   PetscFunctionReturn(0);
1568 }
1569