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