xref: /petsc/src/dm/impls/plex/plexfem.c (revision 0e0f25f23fb62c92bf5bd72e43fa8dee8cfd8e10) 
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]);CHKERRQ(ierr);
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]);CHKERRQ(ierr);
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]);CHKERRQ(ierr);
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 /* This ignores numcomps/comps */
602 static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad,
603                                                               PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX)
604 {
605   PetscDS            prob;
606   PetscSF            sf;
607   DM                 dmFace, dmCell, dmGrad;
608   const PetscScalar *facegeom, *cellgeom, *grad;
609   const PetscInt    *leaves;
610   PetscScalar       *x, *fx;
611   PetscInt           dim, nleaves, loc, fStart, fEnd, pdim, i;
612   PetscErrorCode     ierr;
613 
614   PetscFunctionBegin;
615   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
616   ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr);
617   nleaves = PetscMax(0, nleaves);
618   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
619   ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
620   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
621   ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr);
622   ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr);
623   ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr);
624   ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);
625   if (Grad) {
626     PetscFV fv;
627 
628     ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr);
629     ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr);
630     ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr);
631     ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr);
632     ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr);
633   }
634   ierr = VecGetArray(locX, &x);CHKERRQ(ierr);
635   for (i = 0; i < numids; ++i) {
636     IS              faceIS;
637     const PetscInt *faces;
638     PetscInt        numFaces, f;
639 
640     ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr);
641     if (!faceIS) continue; /* No points with that id on this process */
642     ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr);
643     ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr);
644     for (f = 0; f < numFaces; ++f) {
645       const PetscInt         face = faces[f], *cells;
646       const PetscFVFaceGeom *fg;
647 
648       if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */
649       ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr);
650       if (loc >= 0) continue;
651       ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr);
652       ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr);
653       if (Grad) {
654         const PetscFVCellGeom *cg;
655         const PetscScalar     *cx, *cgrad;
656         PetscScalar           *xG;
657         PetscReal              dx[3];
658         PetscInt               d;
659 
660         ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr);
661         ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr);
662         ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr);
663         ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr);
664         DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx);
665         for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx);
666         ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr);
667       } else {
668         const PetscScalar *xI;
669         PetscScalar       *xG;
670 
671         ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr);
672         ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr);
673         ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr);
674       }
675     }
676     ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr);
677     ierr = ISDestroy(&faceIS);CHKERRQ(ierr);
678   }
679   ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr);
680   if (Grad) {
681     ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr);
682     ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr);
683   }
684   ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);
685   ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr);
686   PetscFunctionReturn(0);
687 }
688 
689 #undef __FUNCT__
690 #define __FUNCT__ "DMPlexInsertBoundaryValues"
691 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM)
692 {
693   PetscInt       numBd, b;
694   PetscErrorCode ierr;
695 
696   PetscFunctionBegin;
697   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
698   PetscValidHeaderSpecific(locX, VEC_CLASSID, 2);
699   if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);}
700   if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);}
701   if (gradFVM)     {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);}
702   ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr);
703   for (b = 0; b < numBd; ++b) {
704     PetscBool       isEssential;
705     const char     *labelname;
706     DMLabel         label;
707     PetscInt        field;
708     PetscObject     obj;
709     PetscClassId    id;
710     void          (*func)();
711     PetscInt        numids;
712     const PetscInt *ids;
713     void           *ctx;
714 
715     ierr = DMPlexGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);CHKERRQ(ierr);
716     if (!isEssential) continue;
717     ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr);
718     ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
719     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
720     if (id == PETSCFE_CLASSID) {
721       ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, field, label, numids, ids, (void (*)(const PetscReal[], PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr);
722     } else if (id == PETSCFV_CLASSID) {
723       if (!faceGeomFVM) continue;
724       ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM,
725                                                      field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr);
726     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
727   }
728   PetscFunctionReturn(0);
729 }
730 
731 #undef __FUNCT__
732 #define __FUNCT__ "DMPlexComputeL2Diff"
733 /*@C
734   DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
735 
736   Input Parameters:
737 + dm    - The DM
738 . funcs - The functions to evaluate for each field component
739 . ctxs  - Optional array of contexts to pass to each function, or NULL.
740 - X     - The coefficient vector u_h
741 
742   Output Parameter:
743 . diff - The diff ||u - u_h||_2
744 
745   Level: developer
746 
747 .seealso: DMPlexProjectFunction(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff()
748 @*/
749 PetscErrorCode DMPlexComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
750 {
751   const PetscInt   debug = 0;
752   PetscSection     section;
753   PetscQuadrature  quad;
754   Vec              localX;
755   PetscScalar     *funcVal, *interpolant;
756   PetscReal       *coords, *v0, *J, *invJ, detJ;
757   PetscReal        localDiff = 0.0;
758   const PetscReal *quadPoints, *quadWeights;
759   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
760   PetscErrorCode   ierr;
761 
762   PetscFunctionBegin;
763   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
764   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
765   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
766   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
767   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
768   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
769   for (field = 0; field < numFields; ++field) {
770     PetscObject  obj;
771     PetscClassId id;
772     PetscInt     Nc;
773 
774     ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
775     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
776     if (id == PETSCFE_CLASSID) {
777       PetscFE fe = (PetscFE) obj;
778 
779       ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
780       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
781     } else if (id == PETSCFV_CLASSID) {
782       PetscFV fv = (PetscFV) obj;
783 
784       ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr);
785       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
786     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
787     numComponents += Nc;
788   }
789   ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
790   ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
791   ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
792   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
793   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
794   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
795   for (c = cStart; c < cEnd; ++c) {
796     PetscScalar *x = NULL;
797     PetscReal    elemDiff = 0.0;
798 
799     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
800     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
801     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
802 
803     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
804       PetscObject  obj;
805       PetscClassId id;
806       void * const ctx = ctxs ? ctxs[field] : NULL;
807       PetscInt     Nb, Nc, q, fc;
808 
809       ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
810       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
811       if (id == PETSCFE_CLASSID)      {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);}
812       else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;}
813       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
814       if (debug) {
815         char title[1024];
816         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
817         ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr);
818       }
819       for (q = 0; q < numQuadPoints; ++q) {
820         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
821         (*funcs[field])(coords, funcVal, ctx);
822         if (id == PETSCFE_CLASSID)      {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
823         else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
824         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
825         for (fc = 0; fc < Nc; ++fc) {
826           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);}
827           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
828         }
829       }
830       fieldOffset += Nb*Nc;
831     }
832     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
833     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
834     localDiff += elemDiff;
835   }
836   ierr  = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr);
837   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
838   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
839   *diff = PetscSqrtReal(*diff);
840   PetscFunctionReturn(0);
841 }
842 
843 #undef __FUNCT__
844 #define __FUNCT__ "DMPlexComputeL2GradientDiff"
845 /*@C
846   DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.
847 
848   Input Parameters:
849 + dm    - The DM
850 . funcs - The gradient functions to evaluate for each field component
851 . ctxs  - Optional array of contexts to pass to each function, or NULL.
852 . X     - The coefficient vector u_h
853 - n     - The vector to project along
854 
855   Output Parameter:
856 . diff - The diff ||(grad u - grad u_h) . n||_2
857 
858   Level: developer
859 
860 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff()
861 @*/
862 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
863 {
864   const PetscInt  debug = 0;
865   PetscSection    section;
866   PetscQuadrature quad;
867   Vec             localX;
868   PetscScalar    *funcVal, *interpolantVec;
869   PetscReal      *coords, *realSpaceDer, *v0, *J, *invJ, detJ;
870   PetscReal       localDiff = 0.0;
871   PetscInt        dim, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp;
872   PetscErrorCode  ierr;
873 
874   PetscFunctionBegin;
875   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
876   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
877   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
878   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
879   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
880   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
881   for (field = 0; field < numFields; ++field) {
882     PetscFE  fe;
883     PetscInt Nc;
884 
885     ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr);
886     ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
887     ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
888     numComponents += Nc;
889   }
890   /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */
891   ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr);
892   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
893   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
894   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
895   for (c = cStart; c < cEnd; ++c) {
896     PetscScalar *x = NULL;
897     PetscReal    elemDiff = 0.0;
898 
899     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
900     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
901     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
902 
903     for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
904       PetscFE          fe;
905       void * const     ctx = ctxs ? ctxs[field] : NULL;
906       const PetscReal *quadPoints, *quadWeights;
907       PetscReal       *basisDer;
908       PetscInt         numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g;
909 
910       ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr);
911       ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
912       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
913       ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr);
914       ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr);
915       if (debug) {
916         char title[1024];
917         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
918         ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr);
919       }
920       for (q = 0; q < numQuadPoints; ++q) {
921         for (d = 0; d < dim; d++) {
922           coords[d] = v0[d];
923           for (e = 0; e < dim; e++) {
924             coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
925           }
926         }
927         (*funcs[field])(coords, n, funcVal, ctx);
928         for (fc = 0; fc < Ncomp; ++fc) {
929           PetscScalar interpolant = 0.0;
930 
931           for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0;
932           for (f = 0; f < Nb; ++f) {
933             const PetscInt fidx = f*Ncomp+fc;
934 
935             for (d = 0; d < dim; ++d) {
936               realSpaceDer[d] = 0.0;
937               for (g = 0; g < dim; ++g) {
938                 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g];
939               }
940               interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d];
941             }
942           }
943           for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d];
944           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);}
945           elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ;
946         }
947       }
948       comp        += Ncomp;
949       fieldOffset += Nb*Ncomp;
950     }
951     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
952     if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);}
953     localDiff += elemDiff;
954   }
955   ierr  = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr);
956   ierr  = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
957   ierr  = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
958   *diff = PetscSqrtReal(*diff);
959   PetscFunctionReturn(0);
960 }
961 
962 #undef __FUNCT__
963 #define __FUNCT__ "DMPlexComputeL2FieldDiff"
964 /*@C
965   DMPlexComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.
966 
967   Input Parameters:
968 + dm    - The DM
969 . funcs - The functions to evaluate for each field component
970 . ctxs  - Optional array of contexts to pass to each function, or NULL.
971 - X     - The coefficient vector u_h
972 
973   Output Parameter:
974 . diff - The array of differences, ||u^f - u^f_h||_2
975 
976   Level: developer
977 
978 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2GradientDiff()
979 @*/
980 PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
981 {
982   const PetscInt   debug = 0;
983   PetscSection     section;
984   PetscQuadrature  quad;
985   Vec              localX;
986   PetscScalar     *funcVal, *interpolant;
987   PetscReal       *coords, *v0, *J, *invJ, detJ;
988   PetscReal       *localDiff;
989   const PetscReal *quadPoints, *quadWeights;
990   PetscInt         dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, cEndInterior, c, field, fieldOffset;
991   PetscErrorCode   ierr;
992 
993   PetscFunctionBegin;
994   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
995   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
996   ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr);
997   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
998   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
999   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
1000   for (field = 0; field < numFields; ++field) {
1001     PetscObject  obj;
1002     PetscClassId id;
1003     PetscInt     Nc;
1004 
1005     ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
1006     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1007     if (id == PETSCFE_CLASSID) {
1008       PetscFE fe = (PetscFE) obj;
1009 
1010       ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
1011       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
1012     } else if (id == PETSCFV_CLASSID) {
1013       PetscFV fv = (PetscFV) obj;
1014 
1015       ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr);
1016       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
1017     } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1018     numComponents += Nc;
1019   }
1020   ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr);
1021   ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr);
1022   ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr);
1023   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1024   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
1025   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1026   for (c = cStart; c < cEnd; ++c) {
1027     PetscScalar *x = NULL;
1028 
1029     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr);
1030     if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
1031     ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
1032 
1033     for (field = 0, fieldOffset = 0; field < numFields; ++field) {
1034       PetscObject  obj;
1035       PetscClassId id;
1036       void * const ctx = ctxs ? ctxs[field] : NULL;
1037       PetscInt     Nb, Nc, q, fc;
1038 
1039       PetscReal       elemDiff = 0.0;
1040 
1041       ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr);
1042       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1043       if (id == PETSCFE_CLASSID)      {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);}
1044       else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;}
1045       else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1046       if (debug) {
1047         char title[1024];
1048         ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr);
1049         ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr);
1050       }
1051       for (q = 0; q < numQuadPoints; ++q) {
1052         CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords);
1053         (*funcs[field])(coords, funcVal, ctx);
1054         if (id == PETSCFE_CLASSID)      {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
1055         else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);}
1056         else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field);
1057         for (fc = 0; fc < Nc; ++fc) {
1058           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);}
1059           elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ;
1060         }
1061       }
1062       fieldOffset += Nb*Nc;
1063       localDiff[field] += elemDiff;
1064     }
1065     ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr);
1066   }
1067   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
1068   ierr = MPI_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
1069   for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]);
1070   ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr);
1071   PetscFunctionReturn(0);
1072 }
1073 
1074 #undef __FUNCT__
1075 #define __FUNCT__ "DMPlexComputeIntegralFEM"
1076 /*@
1077   DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user
1078 
1079   Input Parameters:
1080 + dm - The mesh
1081 . X  - Local input vector
1082 - user - The user context
1083 
1084   Output Parameter:
1085 . integral - Local integral for each field
1086 
1087   Level: developer
1088 
1089 .seealso: DMPlexComputeResidualFEM()
1090 @*/
1091 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user)
1092 {
1093   DM_Plex          *mesh  = (DM_Plex *) dm->data;
1094   DM                dmAux;
1095   Vec               localX, A;
1096   PetscDS           prob, probAux = NULL;
1097   PetscSection      section, sectionAux;
1098   PetscFECellGeom  *cgeom;
1099   PetscScalar      *u, *a = NULL;
1100   PetscReal        *lintegral, *vol;
1101   PetscInt          dim, Nf, f, numCells, cStart, cEnd, cEndInterior, c;
1102   PetscInt          totDim, totDimAux;
1103   PetscErrorCode    ierr;
1104 
1105   PetscFunctionBegin;
1106   ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);
1107   ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr);
1108   ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
1109   ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr);
1110   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1111   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
1112   ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1113   ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr);
1114   ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr);
1115   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1116   ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
1117   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1118   numCells = cEnd - cStart;
1119   ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr);
1120   ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr);
1121   if (dmAux) {
1122     ierr = DMGetDefaultSection(dmAux, &sectionAux);CHKERRQ(ierr);
1123     ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
1124     ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr);
1125   }
1126   ierr = DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr);
1127   ierr = PetscMalloc4(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeom,numCells,&vol);CHKERRQ(ierr);
1128   if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);}
1129   for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;}
1130   for (c = cStart; c < cEnd; ++c) {
1131     PetscScalar *x = NULL;
1132     PetscInt     i;
1133 
1134     ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr);
1135     ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol[c], NULL, NULL);CHKERRQ(ierr);
1136     if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c);
1137     ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr);
1138     for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1139     ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr);
1140     if (dmAux) {
1141       ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1142       for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1143       ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr);
1144     }
1145   }
1146   for (f = 0; f < Nf; ++f) {
1147     PetscObject  obj;
1148     PetscClassId id;
1149     PetscInt     numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1150 
1151     ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);
1152     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1153     if (id == PETSCFE_CLASSID) {
1154       PetscFE         fe = (PetscFE) obj;
1155       PetscQuadrature q;
1156       PetscInt        Nq, Nb;
1157 
1158       ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr);
1159       ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
1160       ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr);
1161       ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr);
1162       blockSize = Nb*Nq;
1163       batchSize = numBlocks * blockSize;
1164       ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr);
1165       numChunks = numCells / (numBatches*batchSize);
1166       Ne        = numChunks*numBatches*batchSize;
1167       Nr        = numCells % (numBatches*batchSize);
1168       offset    = numCells - Nr;
1169       ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, lintegral);CHKERRQ(ierr);
1170       ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr);
1171     } else if (id == PETSCFV_CLASSID) {
1172       /* PetscFV  fv = (PetscFV) obj; */
1173       PetscInt       foff;
1174       PetscPointFunc obj_func;
1175       PetscScalar    lint;
1176 
1177       ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr);
1178       ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr);
1179       if (obj_func) {
1180         for (c = 0; c < numCells; ++c) {
1181           /* TODO: Need full pointwise interpolation and get centroid for x argument */
1182           obj_func(dim, Nf, 0, NULL, NULL, &u[totDim*c+foff], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, &lint);
1183           lintegral[f] = PetscRealPart(lint)*vol[c];
1184         }
1185       }
1186     } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1187   }
1188   if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);}
1189   if (mesh->printFEM) {
1190     ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr);
1191     for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);}
1192     ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr);
1193   }
1194   ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr);
1195   ierr = MPI_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr);
1196   ierr = PetscFree4(lintegral,u,cgeom,vol);CHKERRQ(ierr);
1197   ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);
1198   PetscFunctionReturn(0);
1199 }
1200 
1201 #undef __FUNCT__
1202 #define __FUNCT__ "DMPlexComputeInterpolatorFEM"
1203 /*@
1204   DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM.
1205 
1206   Input Parameters:
1207 + dmf  - The fine mesh
1208 . dmc  - The coarse mesh
1209 - user - The user context
1210 
1211   Output Parameter:
1212 . In  - The interpolation matrix
1213 
1214   Note:
1215   The first member of the user context must be an FEMContext.
1216 
1217   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1218   like a GPU, or vectorize on a multicore machine.
1219 
1220   Level: developer
1221 
1222 .seealso: DMPlexComputeJacobianFEM()
1223 @*/
1224 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user)
1225 {
1226   DM_Plex          *mesh  = (DM_Plex *) dmc->data;
1227   const char       *name  = "Interpolator";
1228   PetscDS           prob;
1229   PetscFE          *feRef;
1230   PetscFV          *fvRef;
1231   PetscSection      fsection, fglobalSection;
1232   PetscSection      csection, cglobalSection;
1233   PetscScalar      *elemMat;
1234   PetscInt          dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c;
1235   PetscInt          cTotDim, rTotDim = 0;
1236   PetscErrorCode    ierr;
1237 
1238   PetscFunctionBegin;
1239   ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1240   ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr);
1241   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1242   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1243   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1244   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1245   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1246   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1247   ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
1248   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1249   ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr);
1250   ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr);
1251   for (f = 0; f < Nf; ++f) {
1252     PetscObject  obj;
1253     PetscClassId id;
1254     PetscInt     rNb = 0, Nc = 0;
1255 
1256     ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);
1257     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1258     if (id == PETSCFE_CLASSID) {
1259       PetscFE fe = (PetscFE) obj;
1260 
1261       ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr);
1262       ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr);
1263       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
1264     } else if (id == PETSCFV_CLASSID) {
1265       PetscFV        fv = (PetscFV) obj;
1266       PetscDualSpace Q;
1267 
1268       ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr);
1269       ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr);
1270       ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr);
1271       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
1272     }
1273     rTotDim += rNb*Nc;
1274   }
1275   ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr);
1276   ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr);
1277   ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr);
1278   for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) {
1279     PetscDualSpace   Qref;
1280     PetscQuadrature  f;
1281     const PetscReal *qpoints, *qweights;
1282     PetscReal       *points;
1283     PetscInt         npoints = 0, Nc, Np, fpdim, i, k, p, d;
1284 
1285     /* Compose points from all dual basis functionals */
1286     if (feRef[fieldI]) {
1287       ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr);
1288       ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr);
1289     } else {
1290       ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr);
1291       ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr);
1292     }
1293     ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr);
1294     for (i = 0; i < fpdim; ++i) {
1295       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1296       ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr);
1297       npoints += Np;
1298     }
1299     ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr);
1300     for (i = 0, k = 0; i < fpdim; ++i) {
1301       ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1302       ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr);
1303       for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d];
1304     }
1305 
1306     for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) {
1307       PetscObject  obj;
1308       PetscClassId id;
1309       PetscReal   *B;
1310       PetscInt     NcJ = 0, cpdim = 0, j;
1311 
1312       ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr);
1313       ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1314       if (id == PETSCFE_CLASSID) {
1315         PetscFE fe = (PetscFE) obj;
1316 
1317         /* Evaluate basis at points */
1318         ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr);
1319         ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr);
1320         /* For now, fields only interpolate themselves */
1321         if (fieldI == fieldJ) {
1322           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);
1323           ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);
1324           for (i = 0, k = 0; i < fpdim; ++i) {
1325             ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1326             ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr);
1327             for (p = 0; p < Np; ++p, ++k) {
1328               for (j = 0; j < cpdim; ++j) {
1329                 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];
1330               }
1331             }
1332           }
1333           ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr);
1334         }
1335       } else if (id == PETSCFV_CLASSID) {
1336         PetscFV        fv = (PetscFV) obj;
1337 
1338         /* Evaluate constant function at points */
1339         ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr);
1340         cpdim = 1;
1341         /* For now, fields only interpolate themselves */
1342         if (fieldI == fieldJ) {
1343           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);
1344           for (i = 0, k = 0; i < fpdim; ++i) {
1345             ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr);
1346             ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr);
1347             for (p = 0; p < Np; ++p, ++k) {
1348               for (j = 0; j < cpdim; ++j) {
1349                 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p];
1350               }
1351             }
1352           }
1353         }
1354       }
1355       offsetJ += cpdim*NcJ;
1356     }
1357     offsetI += fpdim*Nc;
1358     ierr = PetscFree(points);CHKERRQ(ierr);
1359   }
1360   if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);}
1361   /* Preallocate matrix */
1362   {
1363     PetscHashJK ht;
1364     PetscLayout rLayout;
1365     PetscInt   *dnz, *onz, *cellCIndices, *cellFIndices;
1366     PetscInt    locRows, rStart, rEnd, cell, r;
1367 
1368     ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr);
1369     ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr);
1370     ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr);
1371     ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr);
1372     ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr);
1373     ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr);
1374     ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr);
1375     ierr = PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr);
1376     ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr);
1377     for (cell = cStart; cell < cEnd; ++cell) {
1378       ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr);
1379       for (r = 0; r < rTotDim; ++r) {
1380         PetscHashJKKey  key;
1381         PetscHashJKIter missing, iter;
1382 
1383         key.j = cellFIndices[r];
1384         if (key.j < 0) continue;
1385         for (c = 0; c < cTotDim; ++c) {
1386           key.k = cellCIndices[c];
1387           if (key.k < 0) continue;
1388           ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr);
1389           if (missing) {
1390             ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr);
1391             if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart];
1392             else                                     ++onz[key.j-rStart];
1393           }
1394         }
1395       }
1396     }
1397     ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr);
1398     ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr);
1399     ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
1400     ierr = PetscFree4(dnz,onz,cellCIndices,cellFIndices);CHKERRQ(ierr);
1401   }
1402   /* Fill matrix */
1403   ierr = MatZeroEntries(In);CHKERRQ(ierr);
1404   for (c = cStart; c < cEnd; ++c) {
1405     ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr);
1406   }
1407   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);}
1408   ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr);
1409   ierr = PetscFree(elemMat);CHKERRQ(ierr);
1410   ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1411   ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1412   if (mesh->printFEM) {
1413     ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr);
1414     ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr);
1415     ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1416   }
1417   ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1418   PetscFunctionReturn(0);
1419 }
1420 
1421 #undef __FUNCT__
1422 #define __FUNCT__ "DMPlexComputeInjectorFEM"
1423 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user)
1424 {
1425   PetscDS        prob;
1426   PetscFE       *feRef;
1427   PetscFV       *fvRef;
1428   Vec            fv, cv;
1429   IS             fis, cis;
1430   PetscSection   fsection, fglobalSection, csection, cglobalSection;
1431   PetscInt      *cmap, *cellCIndices, *cellFIndices, *cindices, *findices;
1432   PetscInt       cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, offsetC, offsetF, m;
1433   PetscErrorCode ierr;
1434 
1435   PetscFunctionBegin;
1436   ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1437   ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr);
1438   ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr);
1439   ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr);
1440   ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr);
1441   ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr);
1442   ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr);
1443   ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr);
1444   ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr);
1445   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1446   ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr);
1447   ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr);
1448   for (f = 0; f < Nf; ++f) {
1449     PetscObject  obj;
1450     PetscClassId id;
1451     PetscInt     fNb = 0, Nc = 0;
1452 
1453     ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr);
1454     ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr);
1455     if (id == PETSCFE_CLASSID) {
1456       PetscFE fe = (PetscFE) obj;
1457 
1458       ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr);
1459       ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr);
1460       ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
1461     } else if (id == PETSCFV_CLASSID) {
1462       PetscFV        fv = (PetscFV) obj;
1463       PetscDualSpace Q;
1464 
1465       ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr);
1466       ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr);
1467       ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr);
1468       ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr);
1469     }
1470     fTotDim += fNb*Nc;
1471   }
1472   ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr);
1473   ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr);
1474   for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) {
1475     PetscFE        feC;
1476     PetscFV        fvC;
1477     PetscDualSpace QF, QC;
1478     PetscInt       NcF, NcC, fpdim, cpdim;
1479 
1480     if (feRef[field]) {
1481       ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr);
1482       ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr);
1483       ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr);
1484       ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr);
1485       ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr);
1486       ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr);
1487       ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr);
1488     } else {
1489       ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr);
1490       ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr);
1491       ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr);
1492       ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr);
1493       ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr);
1494       ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr);
1495       ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr);
1496     }
1497     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);
1498     for (c = 0; c < cpdim; ++c) {
1499       PetscQuadrature  cfunc;
1500       const PetscReal *cqpoints;
1501       PetscInt         NpC;
1502       PetscBool        found = PETSC_FALSE;
1503 
1504       ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr);
1505       ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr);
1506       if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments");
1507       for (f = 0; f < fpdim; ++f) {
1508         PetscQuadrature  ffunc;
1509         const PetscReal *fqpoints;
1510         PetscReal        sum = 0.0;
1511         PetscInt         NpF, comp;
1512 
1513         ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr);
1514         ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr);
1515         if (NpC != NpF) continue;
1516         for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]);
1517         if (sum > 1.0e-9) continue;
1518         for (comp = 0; comp < NcC; ++comp) {
1519           cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp;
1520         }
1521         found = PETSC_TRUE;
1522         break;
1523       }
1524       if (!found) {
1525         /* TODO We really want the average here, but some asshole put VecScatter in the interface */
1526         if (fvRef[field]) {
1527           PetscInt comp;
1528           for (comp = 0; comp < NcC; ++comp) {
1529             cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp;
1530           }
1531         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection");
1532       }
1533     }
1534     offsetC += cpdim*NcC;
1535     offsetF += fpdim*NcF;
1536   }
1537   for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);}
1538   ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr);
1539 
1540   ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr);
1541   ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr);
1542   ierr = VecGetOwnershipRange(cv, &startC, NULL);CHKERRQ(ierr);
1543   ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr);
1544   ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr);
1545   ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr);
1546   ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr);
1547   for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1;
1548   for (c = cStart; c < cEnd; ++c) {
1549     ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr);
1550     for (d = 0; d < cTotDim; ++d) {
1551       if (cellCIndices[d] < 0) continue;
1552       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]]);
1553       cindices[cellCIndices[d]-startC] = cellCIndices[d];
1554       findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]];
1555     }
1556   }
1557   ierr = PetscFree(cmap);CHKERRQ(ierr);
1558   ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr);
1559 
1560   ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr);
1561   ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr);
1562   ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr);
1563   ierr = ISDestroy(&cis);CHKERRQ(ierr);
1564   ierr = ISDestroy(&fis);CHKERRQ(ierr);
1565   ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr);
1566   ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr);
1567   ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr);
1568   PetscFunctionReturn(0);
1569 }
1570