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