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