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