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