Actual source code: dspacelagrange.c

  1: #include <petsc/private/petscfeimpl.h>
  2: #include <petscdmplex.h>
  3: #include <petscblaslapack.h>

  5: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]);

  7: struct _n_Petsc1DNodeFamily {
  8:   PetscInt        refct;
  9:   PetscDTNodeType nodeFamily;
 10:   PetscReal       gaussJacobiExp;
 11:   PetscInt        nComputed;
 12:   PetscReal     **nodesets;
 13:   PetscBool       endpoints;
 14: };

 16: /* users set node families for PETSCDUALSPACELAGRANGE with just the inputs to this function, but internally we create
 17:  * an object that can cache the computations across multiple dual spaces */
 18: static PetscErrorCode Petsc1DNodeFamilyCreate(PetscDTNodeType family, PetscReal gaussJacobiExp, PetscBool endpoints, Petsc1DNodeFamily *nf)
 19: {
 20:   Petsc1DNodeFamily f;

 22:   PetscFunctionBegin;
 23:   PetscCall(PetscNew(&f));
 24:   switch (family) {
 25:   case PETSCDTNODES_GAUSSJACOBI:
 26:   case PETSCDTNODES_EQUISPACED:
 27:     f->nodeFamily = family;
 28:     break;
 29:   default:
 30:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family");
 31:   }
 32:   f->endpoints      = endpoints;
 33:   f->gaussJacobiExp = 0.;
 34:   if (family == PETSCDTNODES_GAUSSJACOBI) {
 35:     PetscCheck(gaussJacobiExp > -1., PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Gauss-Jacobi exponent must be > -1.");
 36:     f->gaussJacobiExp = gaussJacobiExp;
 37:   }
 38:   f->refct = 1;
 39:   *nf      = f;
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: static PetscErrorCode Petsc1DNodeFamilyReference(Petsc1DNodeFamily nf)
 44: {
 45:   PetscFunctionBegin;
 46:   if (nf) nf->refct++;
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

 50: static PetscErrorCode Petsc1DNodeFamilyDestroy(Petsc1DNodeFamily *nf)
 51: {
 52:   PetscInt i, nc;

 54:   PetscFunctionBegin;
 55:   if (!(*nf)) PetscFunctionReturn(PETSC_SUCCESS);
 56:   if (--(*nf)->refct > 0) {
 57:     *nf = NULL;
 58:     PetscFunctionReturn(PETSC_SUCCESS);
 59:   }
 60:   nc = (*nf)->nComputed;
 61:   for (i = 0; i < nc; i++) PetscCall(PetscFree((*nf)->nodesets[i]));
 62:   PetscCall(PetscFree((*nf)->nodesets));
 63:   PetscCall(PetscFree(*nf));
 64:   *nf = NULL;
 65:   PetscFunctionReturn(PETSC_SUCCESS);
 66: }

 68: static PetscErrorCode Petsc1DNodeFamilyGetNodeSets(Petsc1DNodeFamily f, PetscInt degree, PetscReal ***nodesets)
 69: {
 70:   PetscInt nc;

 72:   PetscFunctionBegin;
 73:   nc = f->nComputed;
 74:   if (degree >= nc) {
 75:     PetscInt    i, j;
 76:     PetscReal **new_nodesets;
 77:     PetscReal  *w;

 79:     PetscCall(PetscMalloc1(degree + 1, &new_nodesets));
 80:     PetscCall(PetscArraycpy(new_nodesets, f->nodesets, nc));
 81:     PetscCall(PetscFree(f->nodesets));
 82:     f->nodesets = new_nodesets;
 83:     PetscCall(PetscMalloc1(degree + 1, &w));
 84:     for (i = nc; i < degree + 1; i++) {
 85:       PetscCall(PetscMalloc1(i + 1, &(f->nodesets[i])));
 86:       if (!i) {
 87:         f->nodesets[i][0] = 0.5;
 88:       } else {
 89:         switch (f->nodeFamily) {
 90:         case PETSCDTNODES_EQUISPACED:
 91:           if (f->endpoints) {
 92:             for (j = 0; j <= i; j++) f->nodesets[i][j] = (PetscReal)j / (PetscReal)i;
 93:           } else {
 94:             /* these nodes are at the centroids of the small simplices created by the equispaced nodes that include
 95:              * the endpoints */
 96:             for (j = 0; j <= i; j++) f->nodesets[i][j] = ((PetscReal)j + 0.5) / ((PetscReal)i + 1.);
 97:           }
 98:           break;
 99:         case PETSCDTNODES_GAUSSJACOBI:
100:           if (f->endpoints) {
101:             PetscCall(PetscDTGaussLobattoJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w));
102:           } else {
103:             PetscCall(PetscDTGaussJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w));
104:           }
105:           break;
106:         default:
107:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family");
108:         }
109:       }
110:     }
111:     PetscCall(PetscFree(w));
112:     f->nComputed = degree + 1;
113:   }
114:   *nodesets = f->nodesets;
115:   PetscFunctionReturn(PETSC_SUCCESS);
116: }

118: /* http://arxiv.org/abs/2002.09421 for details */
119: static PetscErrorCode PetscNodeRecursive_Internal(PetscInt dim, PetscInt degree, PetscReal **nodesets, PetscInt tup[], PetscReal node[])
120: {
121:   PetscReal w;
122:   PetscInt  i, j;

124:   PetscFunctionBeginHot;
125:   w = 0.;
126:   if (dim == 1) {
127:     node[0] = nodesets[degree][tup[0]];
128:     node[1] = nodesets[degree][tup[1]];
129:   } else {
130:     for (i = 0; i < dim + 1; i++) node[i] = 0.;
131:     for (i = 0; i < dim + 1; i++) {
132:       PetscReal wi = nodesets[degree][degree - tup[i]];

134:       for (j = 0; j < dim + 1; j++) tup[dim + 1 + j] = tup[j + (j >= i)];
135:       PetscCall(PetscNodeRecursive_Internal(dim - 1, degree - tup[i], nodesets, &tup[dim + 1], &node[dim + 1]));
136:       for (j = 0; j < dim + 1; j++) node[j + (j >= i)] += wi * node[dim + 1 + j];
137:       w += wi;
138:     }
139:     for (i = 0; i < dim + 1; i++) node[i] /= w;
140:   }
141:   PetscFunctionReturn(PETSC_SUCCESS);
142: }

144: /* compute simplex nodes for the biunit simplex from the 1D node family */
145: static PetscErrorCode Petsc1DNodeFamilyComputeSimplexNodes(Petsc1DNodeFamily f, PetscInt dim, PetscInt degree, PetscReal points[])
146: {
147:   PetscInt   *tup;
148:   PetscInt    k;
149:   PetscInt    npoints;
150:   PetscReal **nodesets = NULL;
151:   PetscInt    worksize;
152:   PetscReal  *nodework;
153:   PetscInt   *tupwork;

155:   PetscFunctionBegin;
156:   PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have non-negative dimension");
157:   PetscCheck(degree >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have non-negative degree");
158:   if (!dim) PetscFunctionReturn(PETSC_SUCCESS);
159:   PetscCall(PetscCalloc1(dim + 2, &tup));
160:   k = 0;
161:   PetscCall(PetscDTBinomialInt(degree + dim, dim, &npoints));
162:   PetscCall(Petsc1DNodeFamilyGetNodeSets(f, degree, &nodesets));
163:   worksize = ((dim + 2) * (dim + 3)) / 2;
164:   PetscCall(PetscCalloc2(worksize, &nodework, worksize, &tupwork));
165:   /* loop over the tuples of length dim with sum at most degree */
166:   for (k = 0; k < npoints; k++) {
167:     PetscInt i;

169:     /* turn thm into tuples of length dim + 1 with sum equal to degree (barycentric indice) */
170:     tup[0] = degree;
171:     for (i = 0; i < dim; i++) tup[0] -= tup[i + 1];
172:     switch (f->nodeFamily) {
173:     case PETSCDTNODES_EQUISPACED:
174:       /* compute equispaces nodes on the unit reference triangle */
175:       if (f->endpoints) {
176:         PetscCheck(degree > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have positive degree");
177:         for (i = 0; i < dim; i++) points[dim * k + i] = (PetscReal)tup[i + 1] / (PetscReal)degree;
178:       } else {
179:         for (i = 0; i < dim; i++) {
180:           /* these nodes are at the centroids of the small simplices created by the equispaced nodes that include
181:            * the endpoints */
182:           points[dim * k + i] = ((PetscReal)tup[i + 1] + 1. / (dim + 1.)) / (PetscReal)(degree + 1.);
183:         }
184:       }
185:       break;
186:     default:
187:       /* compute equispaced nodes on the barycentric reference triangle (the trace on the first dim dimensions are the
188:        * unit reference triangle nodes */
189:       for (i = 0; i < dim + 1; i++) tupwork[i] = tup[i];
190:       PetscCall(PetscNodeRecursive_Internal(dim, degree, nodesets, tupwork, nodework));
191:       for (i = 0; i < dim; i++) points[dim * k + i] = nodework[i + 1];
192:       break;
193:     }
194:     PetscCall(PetscDualSpaceLatticePointLexicographic_Internal(dim, degree, &tup[1]));
195:   }
196:   /* map from unit simplex to biunit simplex */
197:   for (k = 0; k < npoints * dim; k++) points[k] = points[k] * 2. - 1.;
198:   PetscCall(PetscFree2(nodework, tupwork));
199:   PetscCall(PetscFree(tup));
200:   PetscFunctionReturn(PETSC_SUCCESS);
201: }

203: /* If we need to get the dofs from a mesh point, or add values into dofs at a mesh point, and there is more than one dof
204:  * on that mesh point, we have to be careful about getting/adding everything in the right place.
205:  *
206:  * With nodal dofs like PETSCDUALSPACELAGRANGE makes, the general approach to calculate the value of dofs associate
207:  * with a node A is
208:  * - transform the node locations x(A) by the map that takes the mesh point to its reorientation, x' = phi(x(A))
209:  * - figure out which node was originally at the location of the transformed point, A' = idx(x')
210:  * - if the dofs are not scalars, figure out how to represent the transformed dofs in terms of the basis
211:  *   of dofs at A' (using pushforward/pullback rules)
212:  *
213:  * The one sticky point with this approach is the "A' = idx(x')" step: trying to go from real valued coordinates
214:  * back to indices.  I don't want to rely on floating point tolerances.  Additionally, PETSCDUALSPACELAGRANGE may
215:  * eventually support quasi-Lagrangian dofs, which could involve quadrature at multiple points, so the location "x(A)"
216:  * would be ambiguous.
217:  *
218:  * So each dof gets an integer value coordinate (nodeIdx in the structure below).  The choice of integer coordinates
219:  * is somewhat arbitrary, as long as all of the relevant symmetries of the mesh point correspond to *permutations* of
220:  * the integer coordinates, which do not depend on numerical precision.
221:  *
222:  * So
223:  *
224:  * - DMPlexGetTransitiveClosure_Internal() tells me how an orientation turns into a permutation of the vertices of a
225:  *   mesh point
226:  * - The permutation of the vertices, and the nodeIdx values assigned to them, tells what permutation in index space
227:  *   is associated with the orientation
228:  * - I uses that permutation to get xi' = phi(xi(A)), the integer coordinate of the transformed dof
229:  * - I can without numerical issues compute A' = idx(xi')
230:  *
231:  * Here are some examples of how the process works
232:  *
233:  * - With a triangle:
234:  *
235:  *   The triangle has the following integer coordinates for vertices, taken from the barycentric triangle
236:  *
237:  *     closure order 2
238:  *     nodeIdx (0,0,1)
239:  *      \
240:  *       +
241:  *       |\
242:  *       | \
243:  *       |  \
244:  *       |   \    closure order 1
245:  *       |    \ / nodeIdx (0,1,0)
246:  *       +-----+
247:  *        \
248:  *      closure order 0
249:  *      nodeIdx (1,0,0)
250:  *
251:  *   If I do DMPlexGetTransitiveClosure_Internal() with orientation 1, the vertices would appear
252:  *   in the order (1, 2, 0)
253:  *
254:  *   If I list the nodeIdx of each vertex in closure order for orientation 0 (0, 1, 2) and orientation 1 (1, 2, 0), I
255:  *   see
256:  *
257:  *   orientation 0  | orientation 1
258:  *
259:  *   [0] (1,0,0)      [1] (0,1,0)
260:  *   [1] (0,1,0)      [2] (0,0,1)
261:  *   [2] (0,0,1)      [0] (1,0,0)
262:  *          A                B
263:  *
264:  *   In other words, B is the result of a row permutation of A.  But, there is also
265:  *   a column permutation that accomplishes the same result, (2,0,1).
266:  *
267:  *   So if a dof has nodeIdx coordinate (a,b,c), after the transformation its nodeIdx coordinate
268:  *   is (c,a,b), and the transformed degree of freedom will be a linear combination of dofs
269:  *   that originally had coordinate (c,a,b).
270:  *
271:  * - With a quadrilateral:
272:  *
273:  *   The quadrilateral has the following integer coordinates for vertices, taken from concatenating barycentric
274:  *   coordinates for two segments:
275:  *
276:  *     closure order 3      closure order 2
277:  *     nodeIdx (1,0,0,1)    nodeIdx (0,1,0,1)
278:  *                   \      /
279:  *                    +----+
280:  *                    |    |
281:  *                    |    |
282:  *                    +----+
283:  *                   /      \
284:  *     closure order 0      closure order 1
285:  *     nodeIdx (1,0,1,0)    nodeIdx (0,1,1,0)
286:  *
287:  *   If I do DMPlexGetTransitiveClosure_Internal() with orientation 1, the vertices would appear
288:  *   in the order (1, 2, 3, 0)
289:  *
290:  *   If I list the nodeIdx of each vertex in closure order for orientation 0 (0, 1, 2, 3) and
291:  *   orientation 1 (1, 2, 3, 0), I see
292:  *
293:  *   orientation 0  | orientation 1
294:  *
295:  *   [0] (1,0,1,0)    [1] (0,1,1,0)
296:  *   [1] (0,1,1,0)    [2] (0,1,0,1)
297:  *   [2] (0,1,0,1)    [3] (1,0,0,1)
298:  *   [3] (1,0,0,1)    [0] (1,0,1,0)
299:  *          A                B
300:  *
301:  *   The column permutation that accomplishes the same result is (3,2,0,1).
302:  *
303:  *   So if a dof has nodeIdx coordinate (a,b,c,d), after the transformation its nodeIdx coordinate
304:  *   is (d,c,a,b), and the transformed degree of freedom will be a linear combination of dofs
305:  *   that originally had coordinate (d,c,a,b).
306:  *
307:  * Previously PETSCDUALSPACELAGRANGE had hardcoded symmetries for the triangle and quadrilateral,
308:  * but this approach will work for any polytope, such as the wedge (triangular prism).
309:  */
310: struct _n_PetscLagNodeIndices {
311:   PetscInt   refct;
312:   PetscInt   nodeIdxDim;
313:   PetscInt   nodeVecDim;
314:   PetscInt   nNodes;
315:   PetscInt  *nodeIdx; /* for each node an index of size nodeIdxDim */
316:   PetscReal *nodeVec; /* for each node a vector of size nodeVecDim */
317:   PetscInt  *perm;    /* if these are vertices, perm takes DMPlex point index to closure order;
318:                               if these are nodes, perm lists nodes in index revlex order */
319: };

321: /* this is just here so I can access the values in tests/ex1.c outside the library */
322: PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices ni, PetscInt *nodeIdxDim, PetscInt *nodeVecDim, PetscInt *nNodes, const PetscInt *nodeIdx[], const PetscReal *nodeVec[])
323: {
324:   PetscFunctionBegin;
325:   *nodeIdxDim = ni->nodeIdxDim;
326:   *nodeVecDim = ni->nodeVecDim;
327:   *nNodes     = ni->nNodes;
328:   *nodeIdx    = ni->nodeIdx;
329:   *nodeVec    = ni->nodeVec;
330:   PetscFunctionReturn(PETSC_SUCCESS);
331: }

333: static PetscErrorCode PetscLagNodeIndicesReference(PetscLagNodeIndices ni)
334: {
335:   PetscFunctionBegin;
336:   if (ni) ni->refct++;
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

340: static PetscErrorCode PetscLagNodeIndicesDuplicate(PetscLagNodeIndices ni, PetscLagNodeIndices *niNew)
341: {
342:   PetscFunctionBegin;
343:   PetscCall(PetscNew(niNew));
344:   (*niNew)->refct      = 1;
345:   (*niNew)->nodeIdxDim = ni->nodeIdxDim;
346:   (*niNew)->nodeVecDim = ni->nodeVecDim;
347:   (*niNew)->nNodes     = ni->nNodes;
348:   PetscCall(PetscMalloc1(ni->nNodes * ni->nodeIdxDim, &((*niNew)->nodeIdx)));
349:   PetscCall(PetscArraycpy((*niNew)->nodeIdx, ni->nodeIdx, ni->nNodes * ni->nodeIdxDim));
350:   PetscCall(PetscMalloc1(ni->nNodes * ni->nodeVecDim, &((*niNew)->nodeVec)));
351:   PetscCall(PetscArraycpy((*niNew)->nodeVec, ni->nodeVec, ni->nNodes * ni->nodeVecDim));
352:   (*niNew)->perm = NULL;
353:   PetscFunctionReturn(PETSC_SUCCESS);
354: }

356: static PetscErrorCode PetscLagNodeIndicesDestroy(PetscLagNodeIndices *ni)
357: {
358:   PetscFunctionBegin;
359:   if (!(*ni)) PetscFunctionReturn(PETSC_SUCCESS);
360:   if (--(*ni)->refct > 0) {
361:     *ni = NULL;
362:     PetscFunctionReturn(PETSC_SUCCESS);
363:   }
364:   PetscCall(PetscFree((*ni)->nodeIdx));
365:   PetscCall(PetscFree((*ni)->nodeVec));
366:   PetscCall(PetscFree((*ni)->perm));
367:   PetscCall(PetscFree(*ni));
368:   *ni = NULL;
369:   PetscFunctionReturn(PETSC_SUCCESS);
370: }

372: /* The vertices are given nodeIdx coordinates (e.g. the corners of the barycentric triangle).  Those coordinates are
373:  * in some other order, and to understand the effect of different symmetries, we need them to be in closure order.
374:  *
375:  * If sortIdx is PETSC_FALSE, the coordinates are already in revlex order, otherwise we must sort them
376:  * to that order before we do the real work of this function, which is
377:  *
378:  * - mark the vertices in closure order
379:  * - sort them in revlex order
380:  * - use the resulting permutation to list the vertex coordinates in closure order
381:  */
382: static PetscErrorCode PetscLagNodeIndicesComputeVertexOrder(DM dm, PetscLagNodeIndices ni, PetscBool sortIdx)
383: {
384:   PetscInt           v, w, vStart, vEnd, c, d;
385:   PetscInt           nVerts;
386:   PetscInt           closureSize = 0;
387:   PetscInt          *closure     = NULL;
388:   PetscInt          *closureOrder;
389:   PetscInt          *invClosureOrder;
390:   PetscInt          *revlexOrder;
391:   PetscInt          *newNodeIdx;
392:   PetscInt           dim;
393:   Vec                coordVec;
394:   const PetscScalar *coords;

396:   PetscFunctionBegin;
397:   PetscCall(DMGetDimension(dm, &dim));
398:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
399:   nVerts = vEnd - vStart;
400:   PetscCall(PetscMalloc1(nVerts, &closureOrder));
401:   PetscCall(PetscMalloc1(nVerts, &invClosureOrder));
402:   PetscCall(PetscMalloc1(nVerts, &revlexOrder));
403:   if (sortIdx) { /* bubble sort nodeIdx into revlex order */
404:     PetscInt  nodeIdxDim = ni->nodeIdxDim;
405:     PetscInt *idxOrder;

407:     PetscCall(PetscMalloc1(nVerts * nodeIdxDim, &newNodeIdx));
408:     PetscCall(PetscMalloc1(nVerts, &idxOrder));
409:     for (v = 0; v < nVerts; v++) idxOrder[v] = v;
410:     for (v = 0; v < nVerts; v++) {
411:       for (w = v + 1; w < nVerts; w++) {
412:         const PetscInt *iv   = &(ni->nodeIdx[idxOrder[v] * nodeIdxDim]);
413:         const PetscInt *iw   = &(ni->nodeIdx[idxOrder[w] * nodeIdxDim]);
414:         PetscInt        diff = 0;

416:         for (d = nodeIdxDim - 1; d >= 0; d--)
417:           if ((diff = (iv[d] - iw[d]))) break;
418:         if (diff > 0) {
419:           PetscInt swap = idxOrder[v];

421:           idxOrder[v] = idxOrder[w];
422:           idxOrder[w] = swap;
423:         }
424:       }
425:     }
426:     for (v = 0; v < nVerts; v++) {
427:       for (d = 0; d < nodeIdxDim; d++) newNodeIdx[v * ni->nodeIdxDim + d] = ni->nodeIdx[idxOrder[v] * nodeIdxDim + d];
428:     }
429:     PetscCall(PetscFree(ni->nodeIdx));
430:     ni->nodeIdx = newNodeIdx;
431:     newNodeIdx  = NULL;
432:     PetscCall(PetscFree(idxOrder));
433:   }
434:   PetscCall(DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure));
435:   c = closureSize - nVerts;
436:   for (v = 0; v < nVerts; v++) closureOrder[v] = closure[2 * (c + v)] - vStart;
437:   for (v = 0; v < nVerts; v++) invClosureOrder[closureOrder[v]] = v;
438:   PetscCall(DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure));
439:   PetscCall(DMGetCoordinatesLocal(dm, &coordVec));
440:   PetscCall(VecGetArrayRead(coordVec, &coords));
441:   /* bubble sort closure vertices by coordinates in revlex order */
442:   for (v = 0; v < nVerts; v++) revlexOrder[v] = v;
443:   for (v = 0; v < nVerts; v++) {
444:     for (w = v + 1; w < nVerts; w++) {
445:       const PetscScalar *cv   = &coords[closureOrder[revlexOrder[v]] * dim];
446:       const PetscScalar *cw   = &coords[closureOrder[revlexOrder[w]] * dim];
447:       PetscReal          diff = 0;

449:       for (d = dim - 1; d >= 0; d--)
450:         if ((diff = PetscRealPart(cv[d] - cw[d])) != 0.) break;
451:       if (diff > 0.) {
452:         PetscInt swap = revlexOrder[v];

454:         revlexOrder[v] = revlexOrder[w];
455:         revlexOrder[w] = swap;
456:       }
457:     }
458:   }
459:   PetscCall(VecRestoreArrayRead(coordVec, &coords));
460:   PetscCall(PetscMalloc1(ni->nodeIdxDim * nVerts, &newNodeIdx));
461:   /* reorder nodeIdx to be in closure order */
462:   for (v = 0; v < nVerts; v++) {
463:     for (d = 0; d < ni->nodeIdxDim; d++) newNodeIdx[revlexOrder[v] * ni->nodeIdxDim + d] = ni->nodeIdx[v * ni->nodeIdxDim + d];
464:   }
465:   PetscCall(PetscFree(ni->nodeIdx));
466:   ni->nodeIdx = newNodeIdx;
467:   ni->perm    = invClosureOrder;
468:   PetscCall(PetscFree(revlexOrder));
469:   PetscCall(PetscFree(closureOrder));
470:   PetscFunctionReturn(PETSC_SUCCESS);
471: }

473: /* the coordinates of the simplex vertices are the corners of the barycentric simplex.
474:  * When we stack them on top of each other in revlex order, they look like the identity matrix */
475: static PetscErrorCode PetscLagNodeIndicesCreateSimplexVertices(DM dm, PetscLagNodeIndices *nodeIndices)
476: {
477:   PetscLagNodeIndices ni;
478:   PetscInt            dim, d;

480:   PetscFunctionBegin;
481:   PetscCall(PetscNew(&ni));
482:   PetscCall(DMGetDimension(dm, &dim));
483:   ni->nodeIdxDim = dim + 1;
484:   ni->nodeVecDim = 0;
485:   ni->nNodes     = dim + 1;
486:   ni->refct      = 1;
487:   PetscCall(PetscCalloc1((dim + 1) * (dim + 1), &(ni->nodeIdx)));
488:   for (d = 0; d < dim + 1; d++) ni->nodeIdx[d * (dim + 2)] = 1;
489:   PetscCall(PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_FALSE));
490:   *nodeIndices = ni;
491:   PetscFunctionReturn(PETSC_SUCCESS);
492: }

494: /* A polytope that is a tensor product of a facet and a segment.
495:  * We take whatever coordinate system was being used for the facet
496:  * and we concatenate the barycentric coordinates for the vertices
497:  * at the end of the segment, (1,0) and (0,1), to get a coordinate
498:  * system for the tensor product element */
499: static PetscErrorCode PetscLagNodeIndicesCreateTensorVertices(DM dm, PetscLagNodeIndices facetni, PetscLagNodeIndices *nodeIndices)
500: {
501:   PetscLagNodeIndices ni;
502:   PetscInt            nodeIdxDim, subNodeIdxDim = facetni->nodeIdxDim;
503:   PetscInt            nVerts, nSubVerts         = facetni->nNodes;
504:   PetscInt            dim, d, e, f, g;

506:   PetscFunctionBegin;
507:   PetscCall(PetscNew(&ni));
508:   PetscCall(DMGetDimension(dm, &dim));
509:   ni->nodeIdxDim = nodeIdxDim = subNodeIdxDim + 2;
510:   ni->nodeVecDim              = 0;
511:   ni->nNodes = nVerts = 2 * nSubVerts;
512:   ni->refct           = 1;
513:   PetscCall(PetscCalloc1(nodeIdxDim * nVerts, &(ni->nodeIdx)));
514:   for (f = 0, d = 0; d < 2; d++) {
515:     for (e = 0; e < nSubVerts; e++, f++) {
516:       for (g = 0; g < subNodeIdxDim; g++) ni->nodeIdx[f * nodeIdxDim + g] = facetni->nodeIdx[e * subNodeIdxDim + g];
517:       ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim]     = (1 - d);
518:       ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim + 1] = d;
519:     }
520:   }
521:   PetscCall(PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_TRUE));
522:   *nodeIndices = ni;
523:   PetscFunctionReturn(PETSC_SUCCESS);
524: }

526: /* This helps us compute symmetries, and it also helps us compute coordinates for dofs that are being pushed
527:  * forward from a boundary mesh point.
528:  *
529:  * Input:
530:  *
531:  * dm - the target reference cell where we want new coordinates and dof directions to be valid
532:  * vert - the vertex coordinate system for the target reference cell
533:  * p - the point in the target reference cell that the dofs are coming from
534:  * vertp - the vertex coordinate system for p's reference cell
535:  * ornt - the resulting coordinates and dof vectors will be for p under this orientation
536:  * nodep - the node coordinates and dof vectors in p's reference cell
537:  * formDegree - the form degree that the dofs transform as
538:  *
539:  * Output:
540:  *
541:  * pfNodeIdx - the node coordinates for p's dofs, in the dm reference cell, from the ornt perspective
542:  * pfNodeVec - the node dof vectors for p's dofs, in the dm reference cell, from the ornt perspective
543:  */
544: static PetscErrorCode PetscLagNodeIndicesPushForward(DM dm, PetscLagNodeIndices vert, PetscInt p, PetscLagNodeIndices vertp, PetscLagNodeIndices nodep, PetscInt ornt, PetscInt formDegree, PetscInt pfNodeIdx[], PetscReal pfNodeVec[])
545: {
546:   PetscInt          *closureVerts;
547:   PetscInt           closureSize = 0;
548:   PetscInt          *closure     = NULL;
549:   PetscInt           dim, pdim, c, i, j, k, n, v, vStart, vEnd;
550:   PetscInt           nSubVert      = vertp->nNodes;
551:   PetscInt           nodeIdxDim    = vert->nodeIdxDim;
552:   PetscInt           subNodeIdxDim = vertp->nodeIdxDim;
553:   PetscInt           nNodes        = nodep->nNodes;
554:   const PetscInt    *vertIdx       = vert->nodeIdx;
555:   const PetscInt    *subVertIdx    = vertp->nodeIdx;
556:   const PetscInt    *nodeIdx       = nodep->nodeIdx;
557:   const PetscReal   *nodeVec       = nodep->nodeVec;
558:   PetscReal         *J, *Jstar;
559:   PetscReal          detJ;
560:   PetscInt           depth, pdepth, Nk, pNk;
561:   Vec                coordVec;
562:   PetscScalar       *newCoords = NULL;
563:   const PetscScalar *oldCoords = NULL;

565:   PetscFunctionBegin;
566:   PetscCall(DMGetDimension(dm, &dim));
567:   PetscCall(DMPlexGetDepth(dm, &depth));
568:   PetscCall(DMGetCoordinatesLocal(dm, &coordVec));
569:   PetscCall(DMPlexGetPointDepth(dm, p, &pdepth));
570:   pdim = pdepth != depth ? pdepth != 0 ? pdepth : 0 : dim;
571:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
572:   PetscCall(DMGetWorkArray(dm, nSubVert, MPIU_INT, &closureVerts));
573:   PetscCall(DMPlexGetTransitiveClosure_Internal(dm, p, ornt, PETSC_TRUE, &closureSize, &closure));
574:   c = closureSize - nSubVert;
575:   /* we want which cell closure indices the closure of this point corresponds to */
576:   for (v = 0; v < nSubVert; v++) closureVerts[v] = vert->perm[closure[2 * (c + v)] - vStart];
577:   PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure));
578:   /* push forward indices */
579:   for (i = 0; i < nodeIdxDim; i++) { /* for every component of the target index space */
580:     /* check if this is a component that all vertices around this point have in common */
581:     for (j = 1; j < nSubVert; j++) {
582:       if (vertIdx[closureVerts[j] * nodeIdxDim + i] != vertIdx[closureVerts[0] * nodeIdxDim + i]) break;
583:     }
584:     if (j == nSubVert) { /* all vertices have this component in common, directly copy to output */
585:       PetscInt val = vertIdx[closureVerts[0] * nodeIdxDim + i];
586:       for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = val;
587:     } else {
588:       PetscInt subi = -1;
589:       /* there must be a component in vertp that looks the same */
590:       for (k = 0; k < subNodeIdxDim; k++) {
591:         for (j = 0; j < nSubVert; j++) {
592:           if (vertIdx[closureVerts[j] * nodeIdxDim + i] != subVertIdx[j * subNodeIdxDim + k]) break;
593:         }
594:         if (j == nSubVert) {
595:           subi = k;
596:           break;
597:         }
598:       }
599:       PetscCheck(subi >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not find matching coordinate");
600:       /* that component in the vertp system becomes component i in the vert system for each dof */
601:       for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = nodeIdx[n * subNodeIdxDim + subi];
602:     }
603:   }
604:   /* push forward vectors */
605:   PetscCall(DMGetWorkArray(dm, dim * dim, MPIU_REAL, &J));
606:   if (ornt != 0) { /* temporarily change the coordinate vector so
607:                       DMPlexComputeCellGeometryAffineFEM gives us the Jacobian we want */
608:     PetscInt  closureSize2 = 0;
609:     PetscInt *closure2     = NULL;

611:     PetscCall(DMPlexGetTransitiveClosure_Internal(dm, p, 0, PETSC_TRUE, &closureSize2, &closure2));
612:     PetscCall(PetscMalloc1(dim * nSubVert, &newCoords));
613:     PetscCall(VecGetArrayRead(coordVec, &oldCoords));
614:     for (v = 0; v < nSubVert; v++) {
615:       PetscInt d;
616:       for (d = 0; d < dim; d++) newCoords[(closure2[2 * (c + v)] - vStart) * dim + d] = oldCoords[closureVerts[v] * dim + d];
617:     }
618:     PetscCall(VecRestoreArrayRead(coordVec, &oldCoords));
619:     PetscCall(DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize2, &closure2));
620:     PetscCall(VecPlaceArray(coordVec, newCoords));
621:   }
622:   PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, p, NULL, J, NULL, &detJ));
623:   if (ornt != 0) {
624:     PetscCall(VecResetArray(coordVec));
625:     PetscCall(PetscFree(newCoords));
626:   }
627:   PetscCall(DMRestoreWorkArray(dm, nSubVert, MPIU_INT, &closureVerts));
628:   /* compactify */
629:   for (i = 0; i < dim; i++)
630:     for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
631:   /* We have the Jacobian mapping the point's reference cell to this reference cell:
632:    * pulling back a function to the point and applying the dof is what we want,
633:    * so we get the pullback matrix and multiply the dof by that matrix on the right */
634:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk));
635:   PetscCall(PetscDTBinomialInt(pdim, PetscAbsInt(formDegree), &pNk));
636:   PetscCall(DMGetWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar));
637:   PetscCall(PetscDTAltVPullbackMatrix(pdim, dim, J, formDegree, Jstar));
638:   for (n = 0; n < nNodes; n++) {
639:     for (i = 0; i < Nk; i++) {
640:       PetscReal val = 0.;
641:       for (j = 0; j < pNk; j++) val += nodeVec[n * pNk + j] * Jstar[j * Nk + i];
642:       pfNodeVec[n * Nk + i] = val;
643:     }
644:   }
645:   PetscCall(DMRestoreWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar));
646:   PetscCall(DMRestoreWorkArray(dm, dim * dim, MPIU_REAL, &J));
647:   PetscFunctionReturn(PETSC_SUCCESS);
648: }

650: /* given to sets of nodes, take the tensor product, where the product of the dof indices is concatenation and the
651:  * product of the dof vectors is the wedge product */
652: static PetscErrorCode PetscLagNodeIndicesTensor(PetscLagNodeIndices tracei, PetscInt dimT, PetscInt kT, PetscLagNodeIndices fiberi, PetscInt dimF, PetscInt kF, PetscLagNodeIndices *nodeIndices)
653: {
654:   PetscInt            dim = dimT + dimF;
655:   PetscInt            nodeIdxDim, nNodes;
656:   PetscInt            formDegree = kT + kF;
657:   PetscInt            Nk, NkT, NkF;
658:   PetscInt            MkT, MkF;
659:   PetscLagNodeIndices ni;
660:   PetscInt            i, j, l;
661:   PetscReal          *projF, *projT;
662:   PetscReal          *projFstar, *projTstar;
663:   PetscReal          *workF, *workF2, *workT, *workT2, *work, *work2;
664:   PetscReal          *wedgeMat;
665:   PetscReal           sign;

667:   PetscFunctionBegin;
668:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk));
669:   PetscCall(PetscDTBinomialInt(dimT, PetscAbsInt(kT), &NkT));
670:   PetscCall(PetscDTBinomialInt(dimF, PetscAbsInt(kF), &NkF));
671:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(kT), &MkT));
672:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(kF), &MkF));
673:   PetscCall(PetscNew(&ni));
674:   ni->nodeIdxDim = nodeIdxDim = tracei->nodeIdxDim + fiberi->nodeIdxDim;
675:   ni->nodeVecDim              = Nk;
676:   ni->nNodes = nNodes = tracei->nNodes * fiberi->nNodes;
677:   ni->refct           = 1;
678:   PetscCall(PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx)));
679:   /* first concatenate the indices */
680:   for (l = 0, j = 0; j < fiberi->nNodes; j++) {
681:     for (i = 0; i < tracei->nNodes; i++, l++) {
682:       PetscInt m, n = 0;

684:       for (m = 0; m < tracei->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = tracei->nodeIdx[i * tracei->nodeIdxDim + m];
685:       for (m = 0; m < fiberi->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = fiberi->nodeIdx[j * fiberi->nodeIdxDim + m];
686:     }
687:   }

689:   /* now wedge together the push-forward vectors */
690:   PetscCall(PetscMalloc1(nNodes * Nk, &(ni->nodeVec)));
691:   PetscCall(PetscCalloc2(dimT * dim, &projT, dimF * dim, &projF));
692:   for (i = 0; i < dimT; i++) projT[i * (dim + 1)] = 1.;
693:   for (i = 0; i < dimF; i++) projF[i * (dim + dimT + 1) + dimT] = 1.;
694:   PetscCall(PetscMalloc2(MkT * NkT, &projTstar, MkF * NkF, &projFstar));
695:   PetscCall(PetscDTAltVPullbackMatrix(dim, dimT, projT, kT, projTstar));
696:   PetscCall(PetscDTAltVPullbackMatrix(dim, dimF, projF, kF, projFstar));
697:   PetscCall(PetscMalloc6(MkT, &workT, MkT, &workT2, MkF, &workF, MkF, &workF2, Nk, &work, Nk, &work2));
698:   PetscCall(PetscMalloc1(Nk * MkT, &wedgeMat));
699:   sign = (PetscAbsInt(kT * kF) & 1) ? -1. : 1.;
700:   for (l = 0, j = 0; j < fiberi->nNodes; j++) {
701:     PetscInt d, e;

703:     /* push forward fiber k-form */
704:     for (d = 0; d < MkF; d++) {
705:       PetscReal val = 0.;
706:       for (e = 0; e < NkF; e++) val += projFstar[d * NkF + e] * fiberi->nodeVec[j * NkF + e];
707:       workF[d] = val;
708:     }
709:     /* Hodge star to proper form if necessary */
710:     if (kF < 0) {
711:       for (d = 0; d < MkF; d++) workF2[d] = workF[d];
712:       PetscCall(PetscDTAltVStar(dim, PetscAbsInt(kF), 1, workF2, workF));
713:     }
714:     /* Compute the matrix that wedges this form with one of the trace k-form */
715:     PetscCall(PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kF), PetscAbsInt(kT), workF, wedgeMat));
716:     for (i = 0; i < tracei->nNodes; i++, l++) {
717:       /* push forward trace k-form */
718:       for (d = 0; d < MkT; d++) {
719:         PetscReal val = 0.;
720:         for (e = 0; e < NkT; e++) val += projTstar[d * NkT + e] * tracei->nodeVec[i * NkT + e];
721:         workT[d] = val;
722:       }
723:       /* Hodge star to proper form if necessary */
724:       if (kT < 0) {
725:         for (d = 0; d < MkT; d++) workT2[d] = workT[d];
726:         PetscCall(PetscDTAltVStar(dim, PetscAbsInt(kT), 1, workT2, workT));
727:       }
728:       /* compute the wedge product of the push-forward trace form and firer forms */
729:       for (d = 0; d < Nk; d++) {
730:         PetscReal val = 0.;
731:         for (e = 0; e < MkT; e++) val += wedgeMat[d * MkT + e] * workT[e];
732:         work[d] = val;
733:       }
734:       /* inverse Hodge star from proper form if necessary */
735:       if (formDegree < 0) {
736:         for (d = 0; d < Nk; d++) work2[d] = work[d];
737:         PetscCall(PetscDTAltVStar(dim, PetscAbsInt(formDegree), -1, work2, work));
738:       }
739:       /* insert into the array (adjusting for sign) */
740:       for (d = 0; d < Nk; d++) ni->nodeVec[l * Nk + d] = sign * work[d];
741:     }
742:   }
743:   PetscCall(PetscFree(wedgeMat));
744:   PetscCall(PetscFree6(workT, workT2, workF, workF2, work, work2));
745:   PetscCall(PetscFree2(projTstar, projFstar));
746:   PetscCall(PetscFree2(projT, projF));
747:   *nodeIndices = ni;
748:   PetscFunctionReturn(PETSC_SUCCESS);
749: }

751: /* simple union of two sets of nodes */
752: static PetscErrorCode PetscLagNodeIndicesMerge(PetscLagNodeIndices niA, PetscLagNodeIndices niB, PetscLagNodeIndices *nodeIndices)
753: {
754:   PetscLagNodeIndices ni;
755:   PetscInt            nodeIdxDim, nodeVecDim, nNodes;

757:   PetscFunctionBegin;
758:   PetscCall(PetscNew(&ni));
759:   ni->nodeIdxDim = nodeIdxDim = niA->nodeIdxDim;
760:   PetscCheck(niB->nodeIdxDim == nodeIdxDim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeIdxDim");
761:   ni->nodeVecDim = nodeVecDim = niA->nodeVecDim;
762:   PetscCheck(niB->nodeVecDim == nodeVecDim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeVecDim");
763:   ni->nNodes = nNodes = niA->nNodes + niB->nNodes;
764:   ni->refct           = 1;
765:   PetscCall(PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx)));
766:   PetscCall(PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec)));
767:   PetscCall(PetscArraycpy(ni->nodeIdx, niA->nodeIdx, niA->nNodes * nodeIdxDim));
768:   PetscCall(PetscArraycpy(ni->nodeVec, niA->nodeVec, niA->nNodes * nodeVecDim));
769:   PetscCall(PetscArraycpy(&(ni->nodeIdx[niA->nNodes * nodeIdxDim]), niB->nodeIdx, niB->nNodes * nodeIdxDim));
770:   PetscCall(PetscArraycpy(&(ni->nodeVec[niA->nNodes * nodeVecDim]), niB->nodeVec, niB->nNodes * nodeVecDim));
771:   *nodeIndices = ni;
772:   PetscFunctionReturn(PETSC_SUCCESS);
773: }

775: #define PETSCTUPINTCOMPREVLEX(N) \
776:   static int PetscConcat_(PetscTupIntCompRevlex_, N)(const void *a, const void *b) \
777:   { \
778:     const PetscInt *A = (const PetscInt *)a; \
779:     const PetscInt *B = (const PetscInt *)b; \
780:     int             i; \
781:     PetscInt        diff = 0; \
782:     for (i = 0; i < N; i++) { \
783:       diff = A[N - i] - B[N - i]; \
784:       if (diff) break; \
785:     } \
786:     return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1; \
787:   }

789: PETSCTUPINTCOMPREVLEX(3)
790: PETSCTUPINTCOMPREVLEX(4)
791: PETSCTUPINTCOMPREVLEX(5)
792: PETSCTUPINTCOMPREVLEX(6)
793: PETSCTUPINTCOMPREVLEX(7)

795: static int PetscTupIntCompRevlex_N(const void *a, const void *b)
796: {
797:   const PetscInt *A = (const PetscInt *)a;
798:   const PetscInt *B = (const PetscInt *)b;
799:   int             i;
800:   int             N    = A[0];
801:   PetscInt        diff = 0;
802:   for (i = 0; i < N; i++) {
803:     diff = A[N - i] - B[N - i];
804:     if (diff) break;
805:   }
806:   return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1;
807: }

809: /* The nodes are not necessarily in revlex order wrt nodeIdx: get the permutation
810:  * that puts them in that order */
811: static PetscErrorCode PetscLagNodeIndicesGetPermutation(PetscLagNodeIndices ni, PetscInt *perm[])
812: {
813:   PetscFunctionBegin;
814:   if (!(ni->perm)) {
815:     PetscInt *sorter;
816:     PetscInt  m          = ni->nNodes;
817:     PetscInt  nodeIdxDim = ni->nodeIdxDim;
818:     PetscInt  i, j, k, l;
819:     PetscInt *prm;
820:     int (*comp)(const void *, const void *);

822:     PetscCall(PetscMalloc1((nodeIdxDim + 2) * m, &sorter));
823:     for (k = 0, l = 0, i = 0; i < m; i++) {
824:       sorter[k++] = nodeIdxDim + 1;
825:       sorter[k++] = i;
826:       for (j = 0; j < nodeIdxDim; j++) sorter[k++] = ni->nodeIdx[l++];
827:     }
828:     switch (nodeIdxDim) {
829:     case 2:
830:       comp = PetscTupIntCompRevlex_3;
831:       break;
832:     case 3:
833:       comp = PetscTupIntCompRevlex_4;
834:       break;
835:     case 4:
836:       comp = PetscTupIntCompRevlex_5;
837:       break;
838:     case 5:
839:       comp = PetscTupIntCompRevlex_6;
840:       break;
841:     case 6:
842:       comp = PetscTupIntCompRevlex_7;
843:       break;
844:     default:
845:       comp = PetscTupIntCompRevlex_N;
846:       break;
847:     }
848:     qsort(sorter, m, (nodeIdxDim + 2) * sizeof(PetscInt), comp);
849:     PetscCall(PetscMalloc1(m, &prm));
850:     for (i = 0; i < m; i++) prm[i] = sorter[(nodeIdxDim + 2) * i + 1];
851:     ni->perm = prm;
852:     PetscCall(PetscFree(sorter));
853:   }
854:   *perm = ni->perm;
855:   PetscFunctionReturn(PETSC_SUCCESS);
856: }

858: static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
859: {
860:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

862:   PetscFunctionBegin;
863:   if (lag->symperms) {
864:     PetscInt **selfSyms = lag->symperms[0];

866:     if (selfSyms) {
867:       PetscInt i, **allocated = &selfSyms[-lag->selfSymOff];

869:       for (i = 0; i < lag->numSelfSym; i++) PetscCall(PetscFree(allocated[i]));
870:       PetscCall(PetscFree(allocated));
871:     }
872:     PetscCall(PetscFree(lag->symperms));
873:   }
874:   if (lag->symflips) {
875:     PetscScalar **selfSyms = lag->symflips[0];

877:     if (selfSyms) {
878:       PetscInt      i;
879:       PetscScalar **allocated = &selfSyms[-lag->selfSymOff];

881:       for (i = 0; i < lag->numSelfSym; i++) PetscCall(PetscFree(allocated[i]));
882:       PetscCall(PetscFree(allocated));
883:     }
884:     PetscCall(PetscFree(lag->symflips));
885:   }
886:   PetscCall(Petsc1DNodeFamilyDestroy(&(lag->nodeFamily)));
887:   PetscCall(PetscLagNodeIndicesDestroy(&(lag->vertIndices)));
888:   PetscCall(PetscLagNodeIndicesDestroy(&(lag->intNodeIndices)));
889:   PetscCall(PetscLagNodeIndicesDestroy(&(lag->allNodeIndices)));
890:   PetscCall(PetscFree(lag));
891:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL));
892:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL));
893:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", NULL));
894:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", NULL));
895:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL));
896:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL));
897:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL));
898:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL));
899:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", NULL));
900:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", NULL));
901:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", NULL));
902:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", NULL));
903:   PetscFunctionReturn(PETSC_SUCCESS);
904: }

906: static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
907: {
908:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

910:   PetscFunctionBegin;
911:   PetscCall(PetscViewerASCIIPrintf(viewer, "%s %s%sLagrange dual space\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "", lag->trimmed ? "trimmed " : ""));
912:   PetscFunctionReturn(PETSC_SUCCESS);
913: }

915: static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
916: {
917:   PetscBool iascii;

919:   PetscFunctionBegin;
922:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
923:   if (iascii) PetscCall(PetscDualSpaceLagrangeView_Ascii(sp, viewer));
924:   PetscFunctionReturn(PETSC_SUCCESS);
925: }

927: static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscDualSpace sp, PetscOptionItems *PetscOptionsObject)
928: {
929:   PetscBool       continuous, tensor, trimmed, flg, flg2, flg3;
930:   PetscDTNodeType nodeType;
931:   PetscReal       nodeExponent;
932:   PetscInt        momentOrder;
933:   PetscBool       nodeEndpoints, useMoments;

935:   PetscFunctionBegin;
936:   PetscCall(PetscDualSpaceLagrangeGetContinuity(sp, &continuous));
937:   PetscCall(PetscDualSpaceLagrangeGetTensor(sp, &tensor));
938:   PetscCall(PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed));
939:   PetscCall(PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &nodeEndpoints, &nodeExponent));
940:   if (nodeType == PETSCDTNODES_DEFAULT) nodeType = PETSCDTNODES_GAUSSJACOBI;
941:   PetscCall(PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments));
942:   PetscCall(PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder));
943:   PetscOptionsHeadBegin(PetscOptionsObject, "PetscDualSpace Lagrange Options");
944:   PetscCall(PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg));
945:   if (flg) PetscCall(PetscDualSpaceLagrangeSetContinuity(sp, continuous));
946:   PetscCall(PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetTensor", tensor, &tensor, &flg));
947:   if (flg) PetscCall(PetscDualSpaceLagrangeSetTensor(sp, tensor));
948:   PetscCall(PetscOptionsBool("-petscdualspace_lagrange_trimmed", "Flag for trimmed dual space", "PetscDualSpaceLagrangeSetTrimmed", trimmed, &trimmed, &flg));
949:   if (flg) PetscCall(PetscDualSpaceLagrangeSetTrimmed(sp, trimmed));
950:   PetscCall(PetscOptionsEnum("-petscdualspace_lagrange_node_type", "Lagrange node location type", "PetscDualSpaceLagrangeSetNodeType", PetscDTNodeTypes, (PetscEnum)nodeType, (PetscEnum *)&nodeType, &flg));
951:   PetscCall(PetscOptionsBool("-petscdualspace_lagrange_node_endpoints", "Flag for nodes that include endpoints", "PetscDualSpaceLagrangeSetNodeType", nodeEndpoints, &nodeEndpoints, &flg2));
952:   flg3 = PETSC_FALSE;
953:   if (nodeType == PETSCDTNODES_GAUSSJACOBI) PetscCall(PetscOptionsReal("-petscdualspace_lagrange_node_exponent", "Gauss-Jacobi weight function exponent", "PetscDualSpaceLagrangeSetNodeType", nodeExponent, &nodeExponent, &flg3));
954:   if (flg || flg2 || flg3) PetscCall(PetscDualSpaceLagrangeSetNodeType(sp, nodeType, nodeEndpoints, nodeExponent));
955:   PetscCall(PetscOptionsBool("-petscdualspace_lagrange_use_moments", "Use moments (where appropriate) for functionals", "PetscDualSpaceLagrangeSetUseMoments", useMoments, &useMoments, &flg));
956:   if (flg) PetscCall(PetscDualSpaceLagrangeSetUseMoments(sp, useMoments));
957:   PetscCall(PetscOptionsInt("-petscdualspace_lagrange_moment_order", "Quadrature order for moment functionals", "PetscDualSpaceLagrangeSetMomentOrder", momentOrder, &momentOrder, &flg));
958:   if (flg) PetscCall(PetscDualSpaceLagrangeSetMomentOrder(sp, momentOrder));
959:   PetscOptionsHeadEnd();
960:   PetscFunctionReturn(PETSC_SUCCESS);
961: }

963: static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace spNew)
964: {
965:   PetscBool           cont, tensor, trimmed, boundary;
966:   PetscDTNodeType     nodeType;
967:   PetscReal           exponent;
968:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

970:   PetscFunctionBegin;
971:   PetscCall(PetscDualSpaceLagrangeGetContinuity(sp, &cont));
972:   PetscCall(PetscDualSpaceLagrangeSetContinuity(spNew, cont));
973:   PetscCall(PetscDualSpaceLagrangeGetTensor(sp, &tensor));
974:   PetscCall(PetscDualSpaceLagrangeSetTensor(spNew, tensor));
975:   PetscCall(PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed));
976:   PetscCall(PetscDualSpaceLagrangeSetTrimmed(spNew, trimmed));
977:   PetscCall(PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &boundary, &exponent));
978:   PetscCall(PetscDualSpaceLagrangeSetNodeType(spNew, nodeType, boundary, exponent));
979:   if (lag->nodeFamily) {
980:     PetscDualSpace_Lag *lagnew = (PetscDualSpace_Lag *)spNew->data;

982:     PetscCall(Petsc1DNodeFamilyReference(lag->nodeFamily));
983:     lagnew->nodeFamily = lag->nodeFamily;
984:   }
985:   PetscFunctionReturn(PETSC_SUCCESS);
986: }

988: /* for making tensor product spaces: take a dual space and product a segment space that has all the same
989:  * specifications (trimmed, continuous, order, node set), except for the form degree */
990: static PetscErrorCode PetscDualSpaceCreateEdgeSubspace_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt k, PetscInt Nc, PetscBool interiorOnly, PetscDualSpace *bdsp)
991: {
992:   DM                  K;
993:   PetscDualSpace_Lag *newlag;

995:   PetscFunctionBegin;
996:   PetscCall(PetscDualSpaceDuplicate(sp, bdsp));
997:   PetscCall(PetscDualSpaceSetFormDegree(*bdsp, k));
998:   PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(1, PETSC_TRUE), &K));
999:   PetscCall(PetscDualSpaceSetDM(*bdsp, K));
1000:   PetscCall(DMDestroy(&K));
1001:   PetscCall(PetscDualSpaceSetOrder(*bdsp, order));
1002:   PetscCall(PetscDualSpaceSetNumComponents(*bdsp, Nc));
1003:   newlag               = (PetscDualSpace_Lag *)(*bdsp)->data;
1004:   newlag->interiorOnly = interiorOnly;
1005:   PetscCall(PetscDualSpaceSetUp(*bdsp));
1006:   PetscFunctionReturn(PETSC_SUCCESS);
1007: }

1009: /* just the points, weights aren't handled */
1010: static PetscErrorCode PetscQuadratureCreateTensor(PetscQuadrature trace, PetscQuadrature fiber, PetscQuadrature *product)
1011: {
1012:   PetscInt         dimTrace, dimFiber;
1013:   PetscInt         numPointsTrace, numPointsFiber;
1014:   PetscInt         dim, numPoints;
1015:   const PetscReal *pointsTrace;
1016:   const PetscReal *pointsFiber;
1017:   PetscReal       *points;
1018:   PetscInt         i, j, k, p;

1020:   PetscFunctionBegin;
1021:   PetscCall(PetscQuadratureGetData(trace, &dimTrace, NULL, &numPointsTrace, &pointsTrace, NULL));
1022:   PetscCall(PetscQuadratureGetData(fiber, &dimFiber, NULL, &numPointsFiber, &pointsFiber, NULL));
1023:   dim       = dimTrace + dimFiber;
1024:   numPoints = numPointsFiber * numPointsTrace;
1025:   PetscCall(PetscMalloc1(numPoints * dim, &points));
1026:   for (p = 0, j = 0; j < numPointsFiber; j++) {
1027:     for (i = 0; i < numPointsTrace; i++, p++) {
1028:       for (k = 0; k < dimTrace; k++) points[p * dim + k] = pointsTrace[i * dimTrace + k];
1029:       for (k = 0; k < dimFiber; k++) points[p * dim + dimTrace + k] = pointsFiber[j * dimFiber + k];
1030:     }
1031:   }
1032:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, product));
1033:   PetscCall(PetscQuadratureSetData(*product, dim, 0, numPoints, points, NULL));
1034:   PetscFunctionReturn(PETSC_SUCCESS);
1035: }

1037: /* Kronecker tensor product where matrix is considered a matrix of k-forms, so that
1038:  * the entries in the product matrix are wedge products of the entries in the original matrices */
1039: static PetscErrorCode MatTensorAltV(Mat trace, Mat fiber, PetscInt dimTrace, PetscInt kTrace, PetscInt dimFiber, PetscInt kFiber, Mat *product)
1040: {
1041:   PetscInt     mTrace, nTrace, mFiber, nFiber, m, n, k, i, j, l;
1042:   PetscInt     dim, NkTrace, NkFiber, Nk;
1043:   PetscInt     dT, dF;
1044:   PetscInt    *nnzTrace, *nnzFiber, *nnz;
1045:   PetscInt     iT, iF, jT, jF, il, jl;
1046:   PetscReal   *workT, *workT2, *workF, *workF2, *work, *workstar;
1047:   PetscReal   *projT, *projF;
1048:   PetscReal   *projTstar, *projFstar;
1049:   PetscReal   *wedgeMat;
1050:   PetscReal    sign;
1051:   PetscScalar *workS;
1052:   Mat          prod;
1053:   /* this produces dof groups that look like the identity */

1055:   PetscFunctionBegin;
1056:   PetscCall(MatGetSize(trace, &mTrace, &nTrace));
1057:   PetscCall(PetscDTBinomialInt(dimTrace, PetscAbsInt(kTrace), &NkTrace));
1058:   PetscCheck(nTrace % NkTrace == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "point value space of trace matrix is not a multiple of k-form size");
1059:   PetscCall(MatGetSize(fiber, &mFiber, &nFiber));
1060:   PetscCall(PetscDTBinomialInt(dimFiber, PetscAbsInt(kFiber), &NkFiber));
1061:   PetscCheck(nFiber % NkFiber == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "point value space of fiber matrix is not a multiple of k-form size");
1062:   PetscCall(PetscMalloc2(mTrace, &nnzTrace, mFiber, &nnzFiber));
1063:   for (i = 0; i < mTrace; i++) {
1064:     PetscCall(MatGetRow(trace, i, &(nnzTrace[i]), NULL, NULL));
1065:     PetscCheck(nnzTrace[i] % NkTrace == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in trace matrix are not in k-form size blocks");
1066:   }
1067:   for (i = 0; i < mFiber; i++) {
1068:     PetscCall(MatGetRow(fiber, i, &(nnzFiber[i]), NULL, NULL));
1069:     PetscCheck(nnzFiber[i] % NkFiber == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in fiber matrix are not in k-form size blocks");
1070:   }
1071:   dim = dimTrace + dimFiber;
1072:   k   = kFiber + kTrace;
1073:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk));
1074:   m = mTrace * mFiber;
1075:   PetscCall(PetscMalloc1(m, &nnz));
1076:   for (l = 0, j = 0; j < mFiber; j++)
1077:     for (i = 0; i < mTrace; i++, l++) nnz[l] = (nnzTrace[i] / NkTrace) * (nnzFiber[j] / NkFiber) * Nk;
1078:   n = (nTrace / NkTrace) * (nFiber / NkFiber) * Nk;
1079:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &prod));
1080:   PetscCall(PetscFree(nnz));
1081:   PetscCall(PetscFree2(nnzTrace, nnzFiber));
1082:   /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1083:   PetscCall(MatSetOption(prod, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
1084:   /* compute pullbacks */
1085:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(kTrace), &dT));
1086:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(kFiber), &dF));
1087:   PetscCall(PetscMalloc4(dimTrace * dim, &projT, dimFiber * dim, &projF, dT * NkTrace, &projTstar, dF * NkFiber, &projFstar));
1088:   PetscCall(PetscArrayzero(projT, dimTrace * dim));
1089:   for (i = 0; i < dimTrace; i++) projT[i * (dim + 1)] = 1.;
1090:   PetscCall(PetscArrayzero(projF, dimFiber * dim));
1091:   for (i = 0; i < dimFiber; i++) projF[i * (dim + 1) + dimTrace] = 1.;
1092:   PetscCall(PetscDTAltVPullbackMatrix(dim, dimTrace, projT, kTrace, projTstar));
1093:   PetscCall(PetscDTAltVPullbackMatrix(dim, dimFiber, projF, kFiber, projFstar));
1094:   PetscCall(PetscMalloc5(dT, &workT, dF, &workF, Nk, &work, Nk, &workstar, Nk, &workS));
1095:   PetscCall(PetscMalloc2(dT, &workT2, dF, &workF2));
1096:   PetscCall(PetscMalloc1(Nk * dT, &wedgeMat));
1097:   sign = (PetscAbsInt(kTrace * kFiber) & 1) ? -1. : 1.;
1098:   for (i = 0, iF = 0; iF < mFiber; iF++) {
1099:     PetscInt           ncolsF, nformsF;
1100:     const PetscInt    *colsF;
1101:     const PetscScalar *valsF;

1103:     PetscCall(MatGetRow(fiber, iF, &ncolsF, &colsF, &valsF));
1104:     nformsF = ncolsF / NkFiber;
1105:     for (iT = 0; iT < mTrace; iT++, i++) {
1106:       PetscInt           ncolsT, nformsT;
1107:       const PetscInt    *colsT;
1108:       const PetscScalar *valsT;

1110:       PetscCall(MatGetRow(trace, iT, &ncolsT, &colsT, &valsT));
1111:       nformsT = ncolsT / NkTrace;
1112:       for (j = 0, jF = 0; jF < nformsF; jF++) {
1113:         PetscInt colF = colsF[jF * NkFiber] / NkFiber;

1115:         for (il = 0; il < dF; il++) {
1116:           PetscReal val = 0.;
1117:           for (jl = 0; jl < NkFiber; jl++) val += projFstar[il * NkFiber + jl] * PetscRealPart(valsF[jF * NkFiber + jl]);
1118:           workF[il] = val;
1119:         }
1120:         if (kFiber < 0) {
1121:           for (il = 0; il < dF; il++) workF2[il] = workF[il];
1122:           PetscCall(PetscDTAltVStar(dim, PetscAbsInt(kFiber), 1, workF2, workF));
1123:         }
1124:         PetscCall(PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kFiber), PetscAbsInt(kTrace), workF, wedgeMat));
1125:         for (jT = 0; jT < nformsT; jT++, j++) {
1126:           PetscInt           colT = colsT[jT * NkTrace] / NkTrace;
1127:           PetscInt           col  = colF * (nTrace / NkTrace) + colT;
1128:           const PetscScalar *vals;

1130:           for (il = 0; il < dT; il++) {
1131:             PetscReal val = 0.;
1132:             for (jl = 0; jl < NkTrace; jl++) val += projTstar[il * NkTrace + jl] * PetscRealPart(valsT[jT * NkTrace + jl]);
1133:             workT[il] = val;
1134:           }
1135:           if (kTrace < 0) {
1136:             for (il = 0; il < dT; il++) workT2[il] = workT[il];
1137:             PetscCall(PetscDTAltVStar(dim, PetscAbsInt(kTrace), 1, workT2, workT));
1138:           }

1140:           for (il = 0; il < Nk; il++) {
1141:             PetscReal val = 0.;
1142:             for (jl = 0; jl < dT; jl++) val += sign * wedgeMat[il * dT + jl] * workT[jl];
1143:             work[il] = val;
1144:           }
1145:           if (k < 0) {
1146:             PetscCall(PetscDTAltVStar(dim, PetscAbsInt(k), -1, work, workstar));
1147: #if defined(PETSC_USE_COMPLEX)
1148:             for (l = 0; l < Nk; l++) workS[l] = workstar[l];
1149:             vals = &workS[0];
1150: #else
1151:             vals = &workstar[0];
1152: #endif
1153:           } else {
1154: #if defined(PETSC_USE_COMPLEX)
1155:             for (l = 0; l < Nk; l++) workS[l] = work[l];
1156:             vals = &workS[0];
1157: #else
1158:             vals = &work[0];
1159: #endif
1160:           }
1161:           for (l = 0; l < Nk; l++) PetscCall(MatSetValue(prod, i, col * Nk + l, vals[l], INSERT_VALUES)); /* Nk */
1162:         }                                                                                                 /* jT */
1163:       }                                                                                                   /* jF */
1164:       PetscCall(MatRestoreRow(trace, iT, &ncolsT, &colsT, &valsT));
1165:     } /* iT */
1166:     PetscCall(MatRestoreRow(fiber, iF, &ncolsF, &colsF, &valsF));
1167:   } /* iF */
1168:   PetscCall(PetscFree(wedgeMat));
1169:   PetscCall(PetscFree4(projT, projF, projTstar, projFstar));
1170:   PetscCall(PetscFree2(workT2, workF2));
1171:   PetscCall(PetscFree5(workT, workF, work, workstar, workS));
1172:   PetscCall(MatAssemblyBegin(prod, MAT_FINAL_ASSEMBLY));
1173:   PetscCall(MatAssemblyEnd(prod, MAT_FINAL_ASSEMBLY));
1174:   *product = prod;
1175:   PetscFunctionReturn(PETSC_SUCCESS);
1176: }

1178: /* Union of quadrature points, with an attempt to identify common points in the two sets */
1179: static PetscErrorCode PetscQuadraturePointsMerge(PetscQuadrature quadA, PetscQuadrature quadB, PetscQuadrature *quadJoint, PetscInt *aToJoint[], PetscInt *bToJoint[])
1180: {
1181:   PetscInt         dimA, dimB;
1182:   PetscInt         nA, nB, nJoint, i, j, d;
1183:   const PetscReal *pointsA;
1184:   const PetscReal *pointsB;
1185:   PetscReal       *pointsJoint;
1186:   PetscInt        *aToJ, *bToJ;
1187:   PetscQuadrature  qJ;

1189:   PetscFunctionBegin;
1190:   PetscCall(PetscQuadratureGetData(quadA, &dimA, NULL, &nA, &pointsA, NULL));
1191:   PetscCall(PetscQuadratureGetData(quadB, &dimB, NULL, &nB, &pointsB, NULL));
1192:   PetscCheck(dimA == dimB, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Quadrature points must be in the same dimension");
1193:   nJoint = nA;
1194:   PetscCall(PetscMalloc1(nA, &aToJ));
1195:   for (i = 0; i < nA; i++) aToJ[i] = i;
1196:   PetscCall(PetscMalloc1(nB, &bToJ));
1197:   for (i = 0; i < nB; i++) {
1198:     for (j = 0; j < nA; j++) {
1199:       bToJ[i] = -1;
1200:       for (d = 0; d < dimA; d++)
1201:         if (PetscAbsReal(pointsB[i * dimA + d] - pointsA[j * dimA + d]) > PETSC_SMALL) break;
1202:       if (d == dimA) {
1203:         bToJ[i] = j;
1204:         break;
1205:       }
1206:     }
1207:     if (bToJ[i] == -1) bToJ[i] = nJoint++;
1208:   }
1209:   *aToJoint = aToJ;
1210:   *bToJoint = bToJ;
1211:   PetscCall(PetscMalloc1(nJoint * dimA, &pointsJoint));
1212:   PetscCall(PetscArraycpy(pointsJoint, pointsA, nA * dimA));
1213:   for (i = 0; i < nB; i++) {
1214:     if (bToJ[i] >= nA) {
1215:       for (d = 0; d < dimA; d++) pointsJoint[bToJ[i] * dimA + d] = pointsB[i * dimA + d];
1216:     }
1217:   }
1218:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &qJ));
1219:   PetscCall(PetscQuadratureSetData(qJ, dimA, 0, nJoint, pointsJoint, NULL));
1220:   *quadJoint = qJ;
1221:   PetscFunctionReturn(PETSC_SUCCESS);
1222: }

1224: /* Matrices matA and matB are both quadrature -> dof matrices: produce a matrix that is joint quadrature -> union of
1225:  * dofs, where the joint quadrature was produced by PetscQuadraturePointsMerge */
1226: static PetscErrorCode MatricesMerge(Mat matA, Mat matB, PetscInt dim, PetscInt k, PetscInt numMerged, const PetscInt aToMerged[], const PetscInt bToMerged[], Mat *matMerged)
1227: {
1228:   PetscInt  m, n, mA, nA, mB, nB, Nk, i, j, l;
1229:   Mat       M;
1230:   PetscInt *nnz;
1231:   PetscInt  maxnnz;
1232:   PetscInt *work;

1234:   PetscFunctionBegin;
1235:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk));
1236:   PetscCall(MatGetSize(matA, &mA, &nA));
1237:   PetscCheck(nA % Nk == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matA column space not a multiple of k-form size");
1238:   PetscCall(MatGetSize(matB, &mB, &nB));
1239:   PetscCheck(nB % Nk == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matB column space not a multiple of k-form size");
1240:   m = mA + mB;
1241:   n = numMerged * Nk;
1242:   PetscCall(PetscMalloc1(m, &nnz));
1243:   maxnnz = 0;
1244:   for (i = 0; i < mA; i++) {
1245:     PetscCall(MatGetRow(matA, i, &(nnz[i]), NULL, NULL));
1246:     PetscCheck(nnz[i] % Nk == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matA are not in k-form size blocks");
1247:     maxnnz = PetscMax(maxnnz, nnz[i]);
1248:   }
1249:   for (i = 0; i < mB; i++) {
1250:     PetscCall(MatGetRow(matB, i, &(nnz[i + mA]), NULL, NULL));
1251:     PetscCheck(nnz[i + mA] % Nk == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matB are not in k-form size blocks");
1252:     maxnnz = PetscMax(maxnnz, nnz[i + mA]);
1253:   }
1254:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &M));
1255:   PetscCall(PetscFree(nnz));
1256:   /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1257:   PetscCall(MatSetOption(M, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
1258:   PetscCall(PetscMalloc1(maxnnz, &work));
1259:   for (i = 0; i < mA; i++) {
1260:     const PetscInt    *cols;
1261:     const PetscScalar *vals;
1262:     PetscInt           nCols;
1263:     PetscCall(MatGetRow(matA, i, &nCols, &cols, &vals));
1264:     for (j = 0; j < nCols / Nk; j++) {
1265:       PetscInt newCol = aToMerged[cols[j * Nk] / Nk];
1266:       for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1267:     }
1268:     PetscCall(MatSetValuesBlocked(M, 1, &i, nCols, work, vals, INSERT_VALUES));
1269:     PetscCall(MatRestoreRow(matA, i, &nCols, &cols, &vals));
1270:   }
1271:   for (i = 0; i < mB; i++) {
1272:     const PetscInt    *cols;
1273:     const PetscScalar *vals;

1275:     PetscInt row = i + mA;
1276:     PetscInt nCols;
1277:     PetscCall(MatGetRow(matB, i, &nCols, &cols, &vals));
1278:     for (j = 0; j < nCols / Nk; j++) {
1279:       PetscInt newCol = bToMerged[cols[j * Nk] / Nk];
1280:       for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1281:     }
1282:     PetscCall(MatSetValuesBlocked(M, 1, &row, nCols, work, vals, INSERT_VALUES));
1283:     PetscCall(MatRestoreRow(matB, i, &nCols, &cols, &vals));
1284:   }
1285:   PetscCall(PetscFree(work));
1286:   PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
1287:   PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
1288:   *matMerged = M;
1289:   PetscFunctionReturn(PETSC_SUCCESS);
1290: }

1292: /* Take a dual space and product a segment space that has all the same specifications (trimmed, continuous, order,
1293:  * node set), except for the form degree.  For computing boundary dofs and for making tensor product spaces */
1294: static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Lagrange(PetscDualSpace sp, DM K, PetscInt f, PetscInt k, PetscInt Ncopies, PetscBool interiorOnly, PetscDualSpace *bdsp)
1295: {
1296:   PetscInt            Nknew, Ncnew;
1297:   PetscInt            dim, pointDim = -1;
1298:   PetscInt            depth;
1299:   DM                  dm;
1300:   PetscDualSpace_Lag *newlag;

1302:   PetscFunctionBegin;
1303:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
1304:   PetscCall(DMGetDimension(dm, &dim));
1305:   PetscCall(DMPlexGetDepth(dm, &depth));
1306:   PetscCall(PetscDualSpaceDuplicate(sp, bdsp));
1307:   PetscCall(PetscDualSpaceSetFormDegree(*bdsp, k));
1308:   if (!K) {
1309:     if (depth == dim) {
1310:       DMPolytopeType ct;

1312:       pointDim = dim - 1;
1313:       PetscCall(DMPlexGetCellType(dm, f, &ct));
1314:       PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
1315:     } else if (depth == 1) {
1316:       pointDim = 0;
1317:       PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DM_POLYTOPE_POINT, &K));
1318:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
1319:   } else {
1320:     PetscCall(PetscObjectReference((PetscObject)K));
1321:     PetscCall(DMGetDimension(K, &pointDim));
1322:   }
1323:   PetscCall(PetscDualSpaceSetDM(*bdsp, K));
1324:   PetscCall(DMDestroy(&K));
1325:   PetscCall(PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew));
1326:   Ncnew = Nknew * Ncopies;
1327:   PetscCall(PetscDualSpaceSetNumComponents(*bdsp, Ncnew));
1328:   newlag               = (PetscDualSpace_Lag *)(*bdsp)->data;
1329:   newlag->interiorOnly = interiorOnly;
1330:   PetscCall(PetscDualSpaceSetUp(*bdsp));
1331:   PetscFunctionReturn(PETSC_SUCCESS);
1332: }

1334: /* Construct simplex nodes from a nodefamily, add Nk dof vectors of length Nk at each node.
1335:  * Return the (quadrature, matrix) form of the dofs and the nodeIndices form as well.
1336:  *
1337:  * Sometimes we want a set of nodes to be contained in the interior of the element,
1338:  * even when the node scheme puts nodes on the boundaries.  numNodeSkip tells
1339:  * the routine how many "layers" of nodes need to be skipped.
1340:  * */
1341: static PetscErrorCode PetscDualSpaceLagrangeCreateSimplexNodeMat(Petsc1DNodeFamily nodeFamily, PetscInt dim, PetscInt sum, PetscInt Nk, PetscInt numNodeSkip, PetscQuadrature *iNodes, Mat *iMat, PetscLagNodeIndices *nodeIndices)
1342: {
1343:   PetscReal          *extraNodeCoords, *nodeCoords;
1344:   PetscInt            nNodes, nExtraNodes;
1345:   PetscInt            i, j, k, extraSum = sum + numNodeSkip * (1 + dim);
1346:   PetscQuadrature     intNodes;
1347:   Mat                 intMat;
1348:   PetscLagNodeIndices ni;

1350:   PetscFunctionBegin;
1351:   PetscCall(PetscDTBinomialInt(dim + sum, dim, &nNodes));
1352:   PetscCall(PetscDTBinomialInt(dim + extraSum, dim, &nExtraNodes));

1354:   PetscCall(PetscMalloc1(dim * nExtraNodes, &extraNodeCoords));
1355:   PetscCall(PetscNew(&ni));
1356:   ni->nodeIdxDim = dim + 1;
1357:   ni->nodeVecDim = Nk;
1358:   ni->nNodes     = nNodes * Nk;
1359:   ni->refct      = 1;
1360:   PetscCall(PetscMalloc1(nNodes * Nk * (dim + 1), &(ni->nodeIdx)));
1361:   PetscCall(PetscMalloc1(nNodes * Nk * Nk, &(ni->nodeVec)));
1362:   for (i = 0; i < nNodes; i++)
1363:     for (j = 0; j < Nk; j++)
1364:       for (k = 0; k < Nk; k++) ni->nodeVec[(i * Nk + j) * Nk + k] = (j == k) ? 1. : 0.;
1365:   PetscCall(Petsc1DNodeFamilyComputeSimplexNodes(nodeFamily, dim, extraSum, extraNodeCoords));
1366:   if (numNodeSkip) {
1367:     PetscInt  k;
1368:     PetscInt *tup;

1370:     PetscCall(PetscMalloc1(dim * nNodes, &nodeCoords));
1371:     PetscCall(PetscMalloc1(dim + 1, &tup));
1372:     for (k = 0; k < nNodes; k++) {
1373:       PetscInt j, c;
1374:       PetscInt index;

1376:       PetscCall(PetscDTIndexToBary(dim + 1, sum, k, tup));
1377:       for (j = 0; j < dim + 1; j++) tup[j] += numNodeSkip;
1378:       for (c = 0; c < Nk; c++) {
1379:         for (j = 0; j < dim + 1; j++) ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1380:       }
1381:       PetscCall(PetscDTBaryToIndex(dim + 1, extraSum, tup, &index));
1382:       for (j = 0; j < dim; j++) nodeCoords[k * dim + j] = extraNodeCoords[index * dim + j];
1383:     }
1384:     PetscCall(PetscFree(tup));
1385:     PetscCall(PetscFree(extraNodeCoords));
1386:   } else {
1387:     PetscInt  k;
1388:     PetscInt *tup;

1390:     nodeCoords = extraNodeCoords;
1391:     PetscCall(PetscMalloc1(dim + 1, &tup));
1392:     for (k = 0; k < nNodes; k++) {
1393:       PetscInt j, c;

1395:       PetscCall(PetscDTIndexToBary(dim + 1, sum, k, tup));
1396:       for (c = 0; c < Nk; c++) {
1397:         for (j = 0; j < dim + 1; j++) {
1398:           /* barycentric indices can have zeros, but we don't want to push forward zeros because it makes it harder to
1399:            * determine which nodes correspond to which under symmetries, so we increase by 1.  This is fine
1400:            * because the nodeIdx coordinates don't have any meaning other than helping to identify symmetries */
1401:           ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1402:         }
1403:       }
1404:     }
1405:     PetscCall(PetscFree(tup));
1406:   }
1407:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &intNodes));
1408:   PetscCall(PetscQuadratureSetData(intNodes, dim, 0, nNodes, nodeCoords, NULL));
1409:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes * Nk, nNodes * Nk, Nk, NULL, &intMat));
1410:   PetscCall(MatSetOption(intMat, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
1411:   for (j = 0; j < nNodes * Nk; j++) {
1412:     PetscInt rem = j % Nk;
1413:     PetscInt a, aprev = j - rem;
1414:     PetscInt anext = aprev + Nk;

1416:     for (a = aprev; a < anext; a++) PetscCall(MatSetValue(intMat, j, a, (a == j) ? 1. : 0., INSERT_VALUES));
1417:   }
1418:   PetscCall(MatAssemblyBegin(intMat, MAT_FINAL_ASSEMBLY));
1419:   PetscCall(MatAssemblyEnd(intMat, MAT_FINAL_ASSEMBLY));
1420:   *iNodes      = intNodes;
1421:   *iMat        = intMat;
1422:   *nodeIndices = ni;
1423:   PetscFunctionReturn(PETSC_SUCCESS);
1424: }

1426: /* once the nodeIndices have been created for the interior of the reference cell, and for all of the boundary cells,
1427:  * push forward the boundary dofs and concatenate them into the full node indices for the dual space */
1428: static PetscErrorCode PetscDualSpaceLagrangeCreateAllNodeIdx(PetscDualSpace sp)
1429: {
1430:   DM                  dm;
1431:   PetscInt            dim, nDofs;
1432:   PetscSection        section;
1433:   PetscInt            pStart, pEnd, p;
1434:   PetscInt            formDegree, Nk;
1435:   PetscInt            nodeIdxDim, spintdim;
1436:   PetscDualSpace_Lag *lag;
1437:   PetscLagNodeIndices ni, verti;

1439:   PetscFunctionBegin;
1440:   lag   = (PetscDualSpace_Lag *)sp->data;
1441:   verti = lag->vertIndices;
1442:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
1443:   PetscCall(DMGetDimension(dm, &dim));
1444:   PetscCall(PetscDualSpaceGetFormDegree(sp, &formDegree));
1445:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk));
1446:   PetscCall(PetscDualSpaceGetSection(sp, &section));
1447:   PetscCall(PetscSectionGetStorageSize(section, &nDofs));
1448:   PetscCall(PetscNew(&ni));
1449:   ni->nodeIdxDim = nodeIdxDim = verti->nodeIdxDim;
1450:   ni->nodeVecDim              = Nk;
1451:   ni->nNodes                  = nDofs;
1452:   ni->refct                   = 1;
1453:   PetscCall(PetscMalloc1(nodeIdxDim * nDofs, &(ni->nodeIdx)));
1454:   PetscCall(PetscMalloc1(Nk * nDofs, &(ni->nodeVec)));
1455:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1456:   PetscCall(PetscSectionGetDof(section, 0, &spintdim));
1457:   if (spintdim) {
1458:     PetscCall(PetscArraycpy(ni->nodeIdx, lag->intNodeIndices->nodeIdx, spintdim * nodeIdxDim));
1459:     PetscCall(PetscArraycpy(ni->nodeVec, lag->intNodeIndices->nodeVec, spintdim * Nk));
1460:   }
1461:   for (p = pStart + 1; p < pEnd; p++) {
1462:     PetscDualSpace      psp = sp->pointSpaces[p];
1463:     PetscDualSpace_Lag *plag;
1464:     PetscInt            dof, off;

1466:     PetscCall(PetscSectionGetDof(section, p, &dof));
1467:     if (!dof) continue;
1468:     plag = (PetscDualSpace_Lag *)psp->data;
1469:     PetscCall(PetscSectionGetOffset(section, p, &off));
1470:     PetscCall(PetscLagNodeIndicesPushForward(dm, verti, p, plag->vertIndices, plag->intNodeIndices, 0, formDegree, &(ni->nodeIdx[off * nodeIdxDim]), &(ni->nodeVec[off * Nk])));
1471:   }
1472:   lag->allNodeIndices = ni;
1473:   PetscFunctionReturn(PETSC_SUCCESS);
1474: }

1476: /* once the (quadrature, Matrix) forms of the dofs have been created for the interior of the
1477:  * reference cell and for the boundary cells, jk
1478:  * push forward the boundary data and concatenate them into the full (quadrature, matrix) data
1479:  * for the dual space */
1480: static PetscErrorCode PetscDualSpaceCreateAllDataFromInteriorData(PetscDualSpace sp)
1481: {
1482:   DM              dm;
1483:   PetscSection    section;
1484:   PetscInt        pStart, pEnd, p, k, Nk, dim, Nc;
1485:   PetscInt        nNodes;
1486:   PetscInt        countNodes;
1487:   Mat             allMat;
1488:   PetscQuadrature allNodes;
1489:   PetscInt        nDofs;
1490:   PetscInt        maxNzforms, j;
1491:   PetscScalar    *work;
1492:   PetscReal      *L, *J, *Jinv, *v0, *pv0;
1493:   PetscInt       *iwork;
1494:   PetscReal      *nodes;

1496:   PetscFunctionBegin;
1497:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
1498:   PetscCall(DMGetDimension(dm, &dim));
1499:   PetscCall(PetscDualSpaceGetSection(sp, &section));
1500:   PetscCall(PetscSectionGetStorageSize(section, &nDofs));
1501:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
1502:   PetscCall(PetscDualSpaceGetFormDegree(sp, &k));
1503:   PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
1504:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk));
1505:   for (p = pStart, nNodes = 0, maxNzforms = 0; p < pEnd; p++) {
1506:     PetscDualSpace  psp;
1507:     DM              pdm;
1508:     PetscInt        pdim, pNk;
1509:     PetscQuadrature intNodes;
1510:     Mat             intMat;

1512:     PetscCall(PetscDualSpaceGetPointSubspace(sp, p, &psp));
1513:     if (!psp) continue;
1514:     PetscCall(PetscDualSpaceGetDM(psp, &pdm));
1515:     PetscCall(DMGetDimension(pdm, &pdim));
1516:     if (pdim < PetscAbsInt(k)) continue;
1517:     PetscCall(PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk));
1518:     PetscCall(PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat));
1519:     if (intNodes) {
1520:       PetscInt nNodesp;

1522:       PetscCall(PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, NULL, NULL));
1523:       nNodes += nNodesp;
1524:     }
1525:     if (intMat) {
1526:       PetscInt maxNzsp;
1527:       PetscInt maxNzformsp;

1529:       PetscCall(MatSeqAIJGetMaxRowNonzeros(intMat, &maxNzsp));
1530:       PetscCheck(maxNzsp % pNk == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms");
1531:       maxNzformsp = maxNzsp / pNk;
1532:       maxNzforms  = PetscMax(maxNzforms, maxNzformsp);
1533:     }
1534:   }
1535:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, nDofs, nNodes * Nc, maxNzforms * Nk, NULL, &allMat));
1536:   PetscCall(MatSetOption(allMat, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
1537:   PetscCall(PetscMalloc7(dim, &v0, dim, &pv0, dim * dim, &J, dim * dim, &Jinv, Nk * Nk, &L, maxNzforms * Nk, &work, maxNzforms * Nk, &iwork));
1538:   for (j = 0; j < dim; j++) pv0[j] = -1.;
1539:   PetscCall(PetscMalloc1(dim * nNodes, &nodes));
1540:   for (p = pStart, countNodes = 0; p < pEnd; p++) {
1541:     PetscDualSpace  psp;
1542:     PetscQuadrature intNodes;
1543:     DM              pdm;
1544:     PetscInt        pdim, pNk;
1545:     PetscInt        countNodesIn = countNodes;
1546:     PetscReal       detJ;
1547:     Mat             intMat;

1549:     PetscCall(PetscDualSpaceGetPointSubspace(sp, p, &psp));
1550:     if (!psp) continue;
1551:     PetscCall(PetscDualSpaceGetDM(psp, &pdm));
1552:     PetscCall(DMGetDimension(pdm, &pdim));
1553:     if (pdim < PetscAbsInt(k)) continue;
1554:     PetscCall(PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat));
1555:     if (intNodes == NULL && intMat == NULL) continue;
1556:     PetscCall(PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk));
1557:     if (p) {
1558:       PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, Jinv, &detJ));
1559:     } else { /* identity */
1560:       PetscInt i, j;

1562:       for (i = 0; i < dim; i++)
1563:         for (j = 0; j < dim; j++) J[i * dim + j] = Jinv[i * dim + j] = 0.;
1564:       for (i = 0; i < dim; i++) J[i * dim + i] = Jinv[i * dim + i] = 1.;
1565:       for (i = 0; i < dim; i++) v0[i] = -1.;
1566:     }
1567:     if (pdim != dim) { /* compactify Jacobian */
1568:       PetscInt i, j;

1570:       for (i = 0; i < dim; i++)
1571:         for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
1572:     }
1573:     PetscCall(PetscDTAltVPullbackMatrix(pdim, dim, J, k, L));
1574:     if (intNodes) { /* push forward quadrature locations by the affine transformation */
1575:       PetscInt         nNodesp;
1576:       const PetscReal *nodesp;
1577:       PetscInt         j;

1579:       PetscCall(PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, &nodesp, NULL));
1580:       for (j = 0; j < nNodesp; j++, countNodes++) {
1581:         PetscInt d, e;

1583:         for (d = 0; d < dim; d++) {
1584:           nodes[countNodes * dim + d] = v0[d];
1585:           for (e = 0; e < pdim; e++) nodes[countNodes * dim + d] += J[d * pdim + e] * (nodesp[j * pdim + e] - pv0[e]);
1586:         }
1587:       }
1588:     }
1589:     if (intMat) {
1590:       PetscInt nrows;
1591:       PetscInt off;

1593:       PetscCall(PetscSectionGetDof(section, p, &nrows));
1594:       PetscCall(PetscSectionGetOffset(section, p, &off));
1595:       for (j = 0; j < nrows; j++) {
1596:         PetscInt           ncols;
1597:         const PetscInt    *cols;
1598:         const PetscScalar *vals;
1599:         PetscInt           l, d, e;
1600:         PetscInt           row = j + off;

1602:         PetscCall(MatGetRow(intMat, j, &ncols, &cols, &vals));
1603:         PetscCheck(ncols % pNk == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms");
1604:         for (l = 0; l < ncols / pNk; l++) {
1605:           PetscInt blockcol;

1607:           for (d = 0; d < pNk; d++) PetscCheck((cols[l * pNk + d] % pNk) == d, PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms");
1608:           blockcol = cols[l * pNk] / pNk;
1609:           for (d = 0; d < Nk; d++) iwork[l * Nk + d] = (blockcol + countNodesIn) * Nk + d;
1610:           for (d = 0; d < Nk; d++) work[l * Nk + d] = 0.;
1611:           for (d = 0; d < Nk; d++) {
1612:             for (e = 0; e < pNk; e++) {
1613:               /* "push forward" dof by pulling back a k-form to be evaluated on the point: multiply on the right by L */
1614:               work[l * Nk + d] += vals[l * pNk + e] * L[e * Nk + d];
1615:             }
1616:           }
1617:         }
1618:         PetscCall(MatSetValues(allMat, 1, &row, (ncols / pNk) * Nk, iwork, work, INSERT_VALUES));
1619:         PetscCall(MatRestoreRow(intMat, j, &ncols, &cols, &vals));
1620:       }
1621:     }
1622:   }
1623:   PetscCall(MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY));
1624:   PetscCall(MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY));
1625:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &allNodes));
1626:   PetscCall(PetscQuadratureSetData(allNodes, dim, 0, nNodes, nodes, NULL));
1627:   PetscCall(PetscFree7(v0, pv0, J, Jinv, L, work, iwork));
1628:   PetscCall(MatDestroy(&(sp->allMat)));
1629:   sp->allMat = allMat;
1630:   PetscCall(PetscQuadratureDestroy(&(sp->allNodes)));
1631:   sp->allNodes = allNodes;
1632:   PetscFunctionReturn(PETSC_SUCCESS);
1633: }

1635: /* rather than trying to get all data from the functionals, we create
1636:  * the functionals from rows of the quadrature -> dof matrix.
1637:  *
1638:  * Ideally most of the uses of PetscDualSpace in PetscFE will switch
1639:  * to using intMat and allMat, so that the individual functionals
1640:  * don't need to be constructed at all */
1641: static PetscErrorCode PetscDualSpaceComputeFunctionalsFromAllData(PetscDualSpace sp)
1642: {
1643:   PetscQuadrature  allNodes;
1644:   Mat              allMat;
1645:   PetscInt         nDofs;
1646:   PetscInt         dim, k, Nk, Nc, f;
1647:   DM               dm;
1648:   PetscInt         nNodes, spdim;
1649:   const PetscReal *nodes = NULL;
1650:   PetscSection     section;
1651:   PetscBool        useMoments;

1653:   PetscFunctionBegin;
1654:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
1655:   PetscCall(DMGetDimension(dm, &dim));
1656:   PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
1657:   PetscCall(PetscDualSpaceGetFormDegree(sp, &k));
1658:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk));
1659:   PetscCall(PetscDualSpaceGetAllData(sp, &allNodes, &allMat));
1660:   nNodes = 0;
1661:   if (allNodes) PetscCall(PetscQuadratureGetData(allNodes, NULL, NULL, &nNodes, &nodes, NULL));
1662:   PetscCall(MatGetSize(allMat, &nDofs, NULL));
1663:   PetscCall(PetscDualSpaceGetSection(sp, &section));
1664:   PetscCall(PetscSectionGetStorageSize(section, &spdim));
1665:   PetscCheck(spdim == nDofs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "incompatible all matrix size");
1666:   PetscCall(PetscMalloc1(nDofs, &(sp->functional)));
1667:   PetscCall(PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments));
1668:   if (useMoments) {
1669:     Mat              allMat;
1670:     PetscInt         momentOrder, i;
1671:     PetscBool        tensor = PETSC_FALSE;
1672:     const PetscReal *weights;
1673:     PetscScalar     *array;

1675:     PetscCheck(nDofs == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "We do not yet support moments beyond P0, nDofs == %" PetscInt_FMT, nDofs);
1676:     PetscCall(PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder));
1677:     PetscCall(PetscDualSpaceLagrangeGetTensor(sp, &tensor));
1678:     if (!tensor) PetscCall(PetscDTStroudConicalQuadrature(dim, Nc, PetscMax(momentOrder + 1, 1), -1.0, 1.0, &(sp->functional[0])));
1679:     else PetscCall(PetscDTGaussTensorQuadrature(dim, Nc, PetscMax(momentOrder + 1, 1), -1.0, 1.0, &(sp->functional[0])));
1680:     /* Need to replace allNodes and allMat */
1681:     PetscCall(PetscObjectReference((PetscObject)sp->functional[0]));
1682:     PetscCall(PetscQuadratureDestroy(&(sp->allNodes)));
1683:     sp->allNodes = sp->functional[0];
1684:     PetscCall(PetscQuadratureGetData(sp->allNodes, NULL, NULL, &nNodes, NULL, &weights));
1685:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, nDofs, nNodes * Nc, NULL, &allMat));
1686:     PetscCall(MatDenseGetArrayWrite(allMat, &array));
1687:     for (i = 0; i < nNodes * Nc; ++i) array[i] = weights[i];
1688:     PetscCall(MatDenseRestoreArrayWrite(allMat, &array));
1689:     PetscCall(MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY));
1690:     PetscCall(MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY));
1691:     PetscCall(MatDestroy(&(sp->allMat)));
1692:     sp->allMat = allMat;
1693:     PetscFunctionReturn(PETSC_SUCCESS);
1694:   }
1695:   for (f = 0; f < nDofs; f++) {
1696:     PetscInt           ncols, c;
1697:     const PetscInt    *cols;
1698:     const PetscScalar *vals;
1699:     PetscReal         *nodesf;
1700:     PetscReal         *weightsf;
1701:     PetscInt           nNodesf;
1702:     PetscInt           countNodes;

1704:     PetscCall(MatGetRow(allMat, f, &ncols, &cols, &vals));
1705:     PetscCheck(ncols % Nk == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "all matrix is not laid out as blocks of k-forms");
1706:     for (c = 1, nNodesf = 1; c < ncols; c++) {
1707:       if ((cols[c] / Nc) != (cols[c - 1] / Nc)) nNodesf++;
1708:     }
1709:     PetscCall(PetscMalloc1(dim * nNodesf, &nodesf));
1710:     PetscCall(PetscMalloc1(Nc * nNodesf, &weightsf));
1711:     for (c = 0, countNodes = 0; c < ncols; c++) {
1712:       if (!c || ((cols[c] / Nc) != (cols[c - 1] / Nc))) {
1713:         PetscInt d;

1715:         for (d = 0; d < Nc; d++) weightsf[countNodes * Nc + d] = 0.;
1716:         for (d = 0; d < dim; d++) nodesf[countNodes * dim + d] = nodes[(cols[c] / Nc) * dim + d];
1717:         countNodes++;
1718:       }
1719:       weightsf[(countNodes - 1) * Nc + (cols[c] % Nc)] = PetscRealPart(vals[c]);
1720:     }
1721:     PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &(sp->functional[f])));
1722:     PetscCall(PetscQuadratureSetData(sp->functional[f], dim, Nc, nNodesf, nodesf, weightsf));
1723:     PetscCall(MatRestoreRow(allMat, f, &ncols, &cols, &vals));
1724:   }
1725:   PetscFunctionReturn(PETSC_SUCCESS);
1726: }

1728: /* take a matrix meant for k-forms and expand it to one for Ncopies */
1729: static PetscErrorCode PetscDualSpaceLagrangeMatrixCreateCopies(Mat A, PetscInt Nk, PetscInt Ncopies, Mat *Abs)
1730: {
1731:   PetscInt m, n, i, j, k;
1732:   PetscInt maxnnz, *nnz, *iwork;
1733:   Mat      Ac;

1735:   PetscFunctionBegin;
1736:   PetscCall(MatGetSize(A, &m, &n));
1737:   PetscCheck(n % Nk == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of columns in A %" PetscInt_FMT " is not a multiple of Nk %" PetscInt_FMT, n, Nk);
1738:   PetscCall(PetscMalloc1(m * Ncopies, &nnz));
1739:   for (i = 0, maxnnz = 0; i < m; i++) {
1740:     PetscInt innz;
1741:     PetscCall(MatGetRow(A, i, &innz, NULL, NULL));
1742:     PetscCheck(innz % Nk == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A row %" PetscInt_FMT " nnzs is not a multiple of Nk %" PetscInt_FMT, innz, Nk);
1743:     for (j = 0; j < Ncopies; j++) nnz[i * Ncopies + j] = innz;
1744:     maxnnz = PetscMax(maxnnz, innz);
1745:   }
1746:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, m * Ncopies, n * Ncopies, 0, nnz, &Ac));
1747:   PetscCall(MatSetOption(Ac, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE));
1748:   PetscCall(PetscFree(nnz));
1749:   PetscCall(PetscMalloc1(maxnnz, &iwork));
1750:   for (i = 0; i < m; i++) {
1751:     PetscInt           innz;
1752:     const PetscInt    *cols;
1753:     const PetscScalar *vals;

1755:     PetscCall(MatGetRow(A, i, &innz, &cols, &vals));
1756:     for (j = 0; j < innz; j++) iwork[j] = (cols[j] / Nk) * (Nk * Ncopies) + (cols[j] % Nk);
1757:     for (j = 0; j < Ncopies; j++) {
1758:       PetscInt row = i * Ncopies + j;

1760:       PetscCall(MatSetValues(Ac, 1, &row, innz, iwork, vals, INSERT_VALUES));
1761:       for (k = 0; k < innz; k++) iwork[k] += Nk;
1762:     }
1763:     PetscCall(MatRestoreRow(A, i, &innz, &cols, &vals));
1764:   }
1765:   PetscCall(PetscFree(iwork));
1766:   PetscCall(MatAssemblyBegin(Ac, MAT_FINAL_ASSEMBLY));
1767:   PetscCall(MatAssemblyEnd(Ac, MAT_FINAL_ASSEMBLY));
1768:   *Abs = Ac;
1769:   PetscFunctionReturn(PETSC_SUCCESS);
1770: }

1772: /* check if a cell is a tensor product of the segment with a facet,
1773:  * specifically checking if f and f2 can be the "endpoints" (like the triangles
1774:  * at either end of a wedge) */
1775: static PetscErrorCode DMPlexPointIsTensor_Internal_Given(DM dm, PetscInt p, PetscInt f, PetscInt f2, PetscBool *isTensor)
1776: {
1777:   PetscInt        coneSize, c;
1778:   const PetscInt *cone;
1779:   const PetscInt *fCone;
1780:   const PetscInt *f2Cone;
1781:   PetscInt        fs[2];
1782:   PetscInt        meetSize, nmeet;
1783:   const PetscInt *meet;

1785:   PetscFunctionBegin;
1786:   fs[0] = f;
1787:   fs[1] = f2;
1788:   PetscCall(DMPlexGetMeet(dm, 2, fs, &meetSize, &meet));
1789:   nmeet = meetSize;
1790:   PetscCall(DMPlexRestoreMeet(dm, 2, fs, &meetSize, &meet));
1791:   /* two points that have a non-empty meet cannot be at opposite ends of a cell */
1792:   if (nmeet) {
1793:     *isTensor = PETSC_FALSE;
1794:     PetscFunctionReturn(PETSC_SUCCESS);
1795:   }
1796:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1797:   PetscCall(DMPlexGetCone(dm, p, &cone));
1798:   PetscCall(DMPlexGetCone(dm, f, &fCone));
1799:   PetscCall(DMPlexGetCone(dm, f2, &f2Cone));
1800:   for (c = 0; c < coneSize; c++) {
1801:     PetscInt        e, ef;
1802:     PetscInt        d = -1, d2 = -1;
1803:     PetscInt        dcount, d2count;
1804:     PetscInt        t = cone[c];
1805:     PetscInt        tConeSize;
1806:     PetscBool       tIsTensor;
1807:     const PetscInt *tCone;

1809:     if (t == f || t == f2) continue;
1810:     /* for every other facet in the cone, check that is has
1811:      * one ridge in common with each end */
1812:     PetscCall(DMPlexGetConeSize(dm, t, &tConeSize));
1813:     PetscCall(DMPlexGetCone(dm, t, &tCone));

1815:     dcount  = 0;
1816:     d2count = 0;
1817:     for (e = 0; e < tConeSize; e++) {
1818:       PetscInt q = tCone[e];
1819:       for (ef = 0; ef < coneSize - 2; ef++) {
1820:         if (fCone[ef] == q) {
1821:           if (dcount) {
1822:             *isTensor = PETSC_FALSE;
1823:             PetscFunctionReturn(PETSC_SUCCESS);
1824:           }
1825:           d = q;
1826:           dcount++;
1827:         } else if (f2Cone[ef] == q) {
1828:           if (d2count) {
1829:             *isTensor = PETSC_FALSE;
1830:             PetscFunctionReturn(PETSC_SUCCESS);
1831:           }
1832:           d2 = q;
1833:           d2count++;
1834:         }
1835:       }
1836:     }
1837:     /* if the whole cell is a tensor with the segment, then this
1838:      * facet should be a tensor with the segment */
1839:     PetscCall(DMPlexPointIsTensor_Internal_Given(dm, t, d, d2, &tIsTensor));
1840:     if (!tIsTensor) {
1841:       *isTensor = PETSC_FALSE;
1842:       PetscFunctionReturn(PETSC_SUCCESS);
1843:     }
1844:   }
1845:   *isTensor = PETSC_TRUE;
1846:   PetscFunctionReturn(PETSC_SUCCESS);
1847: }

1849: /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair
1850:  * that could be the opposite ends */
1851: static PetscErrorCode DMPlexPointIsTensor_Internal(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB)
1852: {
1853:   PetscInt        coneSize, c, c2;
1854:   const PetscInt *cone;

1856:   PetscFunctionBegin;
1857:   PetscCall(DMPlexGetConeSize(dm, p, &coneSize));
1858:   if (!coneSize) {
1859:     if (isTensor) *isTensor = PETSC_FALSE;
1860:     if (endA) *endA = -1;
1861:     if (endB) *endB = -1;
1862:   }
1863:   PetscCall(DMPlexGetCone(dm, p, &cone));
1864:   for (c = 0; c < coneSize; c++) {
1865:     PetscInt f = cone[c];
1866:     PetscInt fConeSize;

1868:     PetscCall(DMPlexGetConeSize(dm, f, &fConeSize));
1869:     if (fConeSize != coneSize - 2) continue;

1871:     for (c2 = c + 1; c2 < coneSize; c2++) {
1872:       PetscInt  f2 = cone[c2];
1873:       PetscBool isTensorff2;
1874:       PetscInt  f2ConeSize;

1876:       PetscCall(DMPlexGetConeSize(dm, f2, &f2ConeSize));
1877:       if (f2ConeSize != coneSize - 2) continue;

1879:       PetscCall(DMPlexPointIsTensor_Internal_Given(dm, p, f, f2, &isTensorff2));
1880:       if (isTensorff2) {
1881:         if (isTensor) *isTensor = PETSC_TRUE;
1882:         if (endA) *endA = f;
1883:         if (endB) *endB = f2;
1884:         PetscFunctionReturn(PETSC_SUCCESS);
1885:       }
1886:     }
1887:   }
1888:   if (isTensor) *isTensor = PETSC_FALSE;
1889:   if (endA) *endA = -1;
1890:   if (endB) *endB = -1;
1891:   PetscFunctionReturn(PETSC_SUCCESS);
1892: }

1894: /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair
1895:  * that could be the opposite ends */
1896: static PetscErrorCode DMPlexPointIsTensor(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB)
1897: {
1898:   DMPlexInterpolatedFlag interpolated;

1900:   PetscFunctionBegin;
1901:   PetscCall(DMPlexIsInterpolated(dm, &interpolated));
1902:   PetscCheck(interpolated == DMPLEX_INTERPOLATED_FULL, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Only for interpolated DMPlex's");
1903:   PetscCall(DMPlexPointIsTensor_Internal(dm, p, isTensor, endA, endB));
1904:   PetscFunctionReturn(PETSC_SUCCESS);
1905: }

1907: /* Let k = formDegree and k' = -sign(k) * dim + k.  Transform a symmetric frame for k-forms on the biunit simplex into
1908:  * a symmetric frame for k'-forms on the biunit simplex.
1909:  *
1910:  * A frame is "symmetric" if the pullback of every symmetry of the biunit simplex is a permutation of the frame.
1911:  *
1912:  * forms in the symmetric frame are used as dofs in the untrimmed simplex spaces.  This way, symmetries of the
1913:  * reference cell result in permutations of dofs grouped by node.
1914:  *
1915:  * Use T to transform dof matrices for k'-forms into dof matrices for k-forms as a block diagonal transformation on
1916:  * the right.
1917:  */
1918: static PetscErrorCode BiunitSimplexSymmetricFormTransformation(PetscInt dim, PetscInt formDegree, PetscReal T[])
1919: {
1920:   PetscInt   k  = formDegree;
1921:   PetscInt   kd = k < 0 ? dim + k : k - dim;
1922:   PetscInt   Nk;
1923:   PetscReal *biToEq, *eqToBi, *biToEqStar, *eqToBiStar;
1924:   PetscInt   fact;

1926:   PetscFunctionBegin;
1927:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk));
1928:   PetscCall(PetscCalloc4(dim * dim, &biToEq, dim * dim, &eqToBi, Nk * Nk, &biToEqStar, Nk * Nk, &eqToBiStar));
1929:   /* fill in biToEq: Jacobian of the transformation from the biunit simplex to the equilateral simplex */
1930:   fact = 0;
1931:   for (PetscInt i = 0; i < dim; i++) {
1932:     biToEq[i * dim + i] = PetscSqrtReal(((PetscReal)i + 2.) / (2. * ((PetscReal)i + 1.)));
1933:     fact += 4 * (i + 1);
1934:     for (PetscInt j = i + 1; j < dim; j++) biToEq[i * dim + j] = PetscSqrtReal(1. / (PetscReal)fact);
1935:   }
1936:   /* fill in eqToBi: Jacobian of the transformation from the equilateral simplex to the biunit simplex */
1937:   fact = 0;
1938:   for (PetscInt j = 0; j < dim; j++) {
1939:     eqToBi[j * dim + j] = PetscSqrtReal(2. * ((PetscReal)j + 1.) / ((PetscReal)j + 2));
1940:     fact += j + 1;
1941:     for (PetscInt i = 0; i < j; i++) eqToBi[i * dim + j] = -PetscSqrtReal(1. / (PetscReal)fact);
1942:   }
1943:   PetscCall(PetscDTAltVPullbackMatrix(dim, dim, biToEq, kd, biToEqStar));
1944:   PetscCall(PetscDTAltVPullbackMatrix(dim, dim, eqToBi, k, eqToBiStar));
1945:   /* product of pullbacks simulates the following steps
1946:    *
1947:    * 1. start with frame W = [w_1, w_2, ..., w_m] of k forms that is symmetric on the biunit simplex:
1948:           if J is the Jacobian of a symmetry of the biunit simplex, then J_k* W = [J_k*w_1, ..., J_k*w_m]
1949:           is a permutation of W.
1950:           Even though a k' form --- a (dim - k) form represented by its Hodge star --- has the same geometric
1951:           content as a k form, W is not a symmetric frame of k' forms on the biunit simplex.  That's because,
1952:           for general Jacobian J, J_k* != J_k'*.
1953:    * 2. pullback W to the equilateral triangle using the k pullback, W_eq = eqToBi_k* W.  All symmetries of the
1954:           equilateral simplex have orthonormal Jacobians.  For an orthonormal Jacobian O, J_k* = J_k'*, so W_eq is
1955:           also a symmetric frame for k' forms on the equilateral simplex.
1956:      3. pullback W_eq back to the biunit simplex using the k' pulback, V = biToEq_k'* W_eq = biToEq_k'* eqToBi_k* W.
1957:           V is a symmetric frame for k' forms on the biunit simplex.
1958:    */
1959:   for (PetscInt i = 0; i < Nk; i++) {
1960:     for (PetscInt j = 0; j < Nk; j++) {
1961:       PetscReal val = 0.;
1962:       for (PetscInt k = 0; k < Nk; k++) val += biToEqStar[i * Nk + k] * eqToBiStar[k * Nk + j];
1963:       T[i * Nk + j] = val;
1964:     }
1965:   }
1966:   PetscCall(PetscFree4(biToEq, eqToBi, biToEqStar, eqToBiStar));
1967:   PetscFunctionReturn(PETSC_SUCCESS);
1968: }

1970: /* permute a quadrature -> dof matrix so that its rows are in revlex order by nodeIdx */
1971: static PetscErrorCode MatPermuteByNodeIdx(Mat A, PetscLagNodeIndices ni, Mat *Aperm)
1972: {
1973:   PetscInt   m, n, i, j;
1974:   PetscInt   nodeIdxDim = ni->nodeIdxDim;
1975:   PetscInt   nodeVecDim = ni->nodeVecDim;
1976:   PetscInt  *perm;
1977:   IS         permIS;
1978:   IS         id;
1979:   PetscInt  *nIdxPerm;
1980:   PetscReal *nVecPerm;

1982:   PetscFunctionBegin;
1983:   PetscCall(PetscLagNodeIndicesGetPermutation(ni, &perm));
1984:   PetscCall(MatGetSize(A, &m, &n));
1985:   PetscCall(PetscMalloc1(nodeIdxDim * m, &nIdxPerm));
1986:   PetscCall(PetscMalloc1(nodeVecDim * m, &nVecPerm));
1987:   for (i = 0; i < m; i++)
1988:     for (j = 0; j < nodeIdxDim; j++) nIdxPerm[i * nodeIdxDim + j] = ni->nodeIdx[perm[i] * nodeIdxDim + j];
1989:   for (i = 0; i < m; i++)
1990:     for (j = 0; j < nodeVecDim; j++) nVecPerm[i * nodeVecDim + j] = ni->nodeVec[perm[i] * nodeVecDim + j];
1991:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, perm, PETSC_USE_POINTER, &permIS));
1992:   PetscCall(ISSetPermutation(permIS));
1993:   PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &id));
1994:   PetscCall(ISSetPermutation(id));
1995:   PetscCall(MatPermute(A, permIS, id, Aperm));
1996:   PetscCall(ISDestroy(&permIS));
1997:   PetscCall(ISDestroy(&id));
1998:   for (i = 0; i < m; i++) perm[i] = i;
1999:   PetscCall(PetscFree(ni->nodeIdx));
2000:   PetscCall(PetscFree(ni->nodeVec));
2001:   ni->nodeIdx = nIdxPerm;
2002:   ni->nodeVec = nVecPerm;
2003:   PetscFunctionReturn(PETSC_SUCCESS);
2004: }

2006: static PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
2007: {
2008:   PetscDualSpace_Lag    *lag   = (PetscDualSpace_Lag *)sp->data;
2009:   DM                     dm    = sp->dm;
2010:   DM                     dmint = NULL;
2011:   PetscInt               order;
2012:   PetscInt               Nc = sp->Nc;
2013:   MPI_Comm               comm;
2014:   PetscBool              continuous;
2015:   PetscSection           section;
2016:   PetscInt               depth, dim, pStart, pEnd, cStart, cEnd, p, *pStratStart, *pStratEnd, d;
2017:   PetscInt               formDegree, Nk, Ncopies;
2018:   PetscInt               tensorf = -1, tensorf2 = -1;
2019:   PetscBool              tensorCell, tensorSpace;
2020:   PetscBool              uniform, trimmed;
2021:   Petsc1DNodeFamily      nodeFamily;
2022:   PetscInt               numNodeSkip;
2023:   DMPlexInterpolatedFlag interpolated;
2024:   PetscBool              isbdm;

2026:   PetscFunctionBegin;
2027:   /* step 1: sanitize input */
2028:   PetscCall(PetscObjectGetComm((PetscObject)sp, &comm));
2029:   PetscCall(DMGetDimension(dm, &dim));
2030:   PetscCall(PetscObjectTypeCompare((PetscObject)sp, PETSCDUALSPACEBDM, &isbdm));
2031:   if (isbdm) {
2032:     sp->k = -(dim - 1); /* form degree of H-div */
2033:     PetscCall(PetscObjectChangeTypeName((PetscObject)sp, PETSCDUALSPACELAGRANGE));
2034:   }
2035:   PetscCall(PetscDualSpaceGetFormDegree(sp, &formDegree));
2036:   PetscCheck(PetscAbsInt(formDegree) <= dim, comm, PETSC_ERR_ARG_OUTOFRANGE, "Form degree must be bounded by dimension");
2037:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk));
2038:   if (sp->Nc <= 0 && lag->numCopies > 0) sp->Nc = Nk * lag->numCopies;
2039:   Nc = sp->Nc;
2040:   PetscCheck(Nc % Nk == 0, comm, PETSC_ERR_ARG_INCOMP, "Number of components is not a multiple of form degree size");
2041:   if (lag->numCopies <= 0) lag->numCopies = Nc / Nk;
2042:   Ncopies = lag->numCopies;
2043:   PetscCheck(Nc / Nk == Ncopies, comm, PETSC_ERR_ARG_INCOMP, "Number of copies * (dim choose k) != Nc");
2044:   if (!dim) sp->order = 0;
2045:   order   = sp->order;
2046:   uniform = sp->uniform;
2047:   PetscCheck(uniform, PETSC_COMM_SELF, PETSC_ERR_SUP, "Variable order not supported yet");
2048:   if (lag->trimmed && !formDegree) lag->trimmed = PETSC_FALSE; /* trimmed spaces are the same as full spaces for 0-forms */
2049:   if (lag->nodeType == PETSCDTNODES_DEFAULT) {
2050:     lag->nodeType     = PETSCDTNODES_GAUSSJACOBI;
2051:     lag->nodeExponent = 0.;
2052:     /* trimmed spaces don't include corner vertices, so don't use end nodes by default */
2053:     lag->endNodes = lag->trimmed ? PETSC_FALSE : PETSC_TRUE;
2054:   }
2055:   /* If a trimmed space and the user did choose nodes with endpoints, skip them by default */
2056:   if (lag->numNodeSkip < 0) lag->numNodeSkip = (lag->trimmed && lag->endNodes) ? 1 : 0;
2057:   numNodeSkip = lag->numNodeSkip;
2058:   PetscCheck(!lag->trimmed || order, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zeroth order trimmed elements");
2059:   if (lag->trimmed && PetscAbsInt(formDegree) == dim) { /* convert trimmed n-forms to untrimmed of one polynomial order less */
2060:     sp->order--;
2061:     order--;
2062:     lag->trimmed = PETSC_FALSE;
2063:   }
2064:   trimmed = lag->trimmed;
2065:   if (!order || PetscAbsInt(formDegree) == dim) lag->continuous = PETSC_FALSE;
2066:   continuous = lag->continuous;
2067:   PetscCall(DMPlexGetDepth(dm, &depth));
2068:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
2069:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2070:   PetscCheck(pStart == 0 && cStart == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Expect DM with chart starting at zero and cells first");
2071:   PetscCheck(cEnd == 1, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Use PETSCDUALSPACEREFINED for multi-cell reference meshes");
2072:   PetscCall(DMPlexIsInterpolated(dm, &interpolated));
2073:   if (interpolated != DMPLEX_INTERPOLATED_FULL) {
2074:     PetscCall(DMPlexInterpolate(dm, &dmint));
2075:   } else {
2076:     PetscCall(PetscObjectReference((PetscObject)dm));
2077:     dmint = dm;
2078:   }
2079:   tensorCell = PETSC_FALSE;
2080:   if (dim > 1) PetscCall(DMPlexPointIsTensor(dmint, 0, &tensorCell, &tensorf, &tensorf2));
2081:   lag->tensorCell = tensorCell;
2082:   if (dim < 2 || !lag->tensorCell) lag->tensorSpace = PETSC_FALSE;
2083:   tensorSpace = lag->tensorSpace;
2084:   if (!lag->nodeFamily) PetscCall(Petsc1DNodeFamilyCreate(lag->nodeType, lag->nodeExponent, lag->endNodes, &lag->nodeFamily));
2085:   nodeFamily = lag->nodeFamily;
2086:   PetscCheck(interpolated == DMPLEX_INTERPOLATED_FULL || !continuous || (PetscAbsInt(formDegree) <= 0 && order <= 1), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Reference element won't support all boundary nodes");

2088:   /* step 2: construct the boundary spaces */
2089:   PetscCall(PetscMalloc2(depth + 1, &pStratStart, depth + 1, &pStratEnd));
2090:   PetscCall(PetscCalloc1(pEnd, &(sp->pointSpaces)));
2091:   for (d = 0; d <= depth; ++d) PetscCall(DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]));
2092:   PetscCall(PetscDualSpaceSectionCreate_Internal(sp, &section));
2093:   sp->pointSection = section;
2094:   if (continuous && !(lag->interiorOnly)) {
2095:     PetscInt h;

2097:     for (p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */
2098:       PetscReal      v0[3];
2099:       DMPolytopeType ptype;
2100:       PetscReal      J[9], detJ;
2101:       PetscInt       q;

2103:       PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, NULL, &detJ));
2104:       PetscCall(DMPlexGetCellType(dm, p, &ptype));

2106:       /* compare to previous facets: if computed, reference that dualspace */
2107:       for (q = pStratStart[depth - 1]; q < p; q++) {
2108:         DMPolytopeType qtype;

2110:         PetscCall(DMPlexGetCellType(dm, q, &qtype));
2111:         if (qtype == ptype) break;
2112:       }
2113:       if (q < p) { /* this facet has the same dual space as that one */
2114:         PetscCall(PetscObjectReference((PetscObject)sp->pointSpaces[q]));
2115:         sp->pointSpaces[p] = sp->pointSpaces[q];
2116:         continue;
2117:       }
2118:       /* if not, recursively compute this dual space */
2119:       PetscCall(PetscDualSpaceCreateFacetSubspace_Lagrange(sp, NULL, p, formDegree, Ncopies, PETSC_FALSE, &sp->pointSpaces[p]));
2120:     }
2121:     for (h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */
2122:       PetscInt hd   = depth - h;
2123:       PetscInt hdim = dim - h;

2125:       if (hdim < PetscAbsInt(formDegree)) break;
2126:       for (p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
2127:         PetscInt        suppSize, s;
2128:         const PetscInt *supp;

2130:         PetscCall(DMPlexGetSupportSize(dm, p, &suppSize));
2131:         PetscCall(DMPlexGetSupport(dm, p, &supp));
2132:         for (s = 0; s < suppSize; s++) {
2133:           DM              qdm;
2134:           PetscDualSpace  qsp, psp;
2135:           PetscInt        c, coneSize, q;
2136:           const PetscInt *cone;
2137:           const PetscInt *refCone;

2139:           q   = supp[0];
2140:           qsp = sp->pointSpaces[q];
2141:           PetscCall(DMPlexGetConeSize(dm, q, &coneSize));
2142:           PetscCall(DMPlexGetCone(dm, q, &cone));
2143:           for (c = 0; c < coneSize; c++)
2144:             if (cone[c] == p) break;
2145:           PetscCheck(c != coneSize, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "cone/support mismatch");
2146:           PetscCall(PetscDualSpaceGetDM(qsp, &qdm));
2147:           PetscCall(DMPlexGetCone(qdm, 0, &refCone));
2148:           /* get the equivalent dual space from the support dual space */
2149:           PetscCall(PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp));
2150:           if (!s) {
2151:             PetscCall(PetscObjectReference((PetscObject)psp));
2152:             sp->pointSpaces[p] = psp;
2153:           }
2154:         }
2155:       }
2156:     }
2157:     for (p = 1; p < pEnd; p++) {
2158:       PetscInt pspdim;
2159:       if (!sp->pointSpaces[p]) continue;
2160:       PetscCall(PetscDualSpaceGetInteriorDimension(sp->pointSpaces[p], &pspdim));
2161:       PetscCall(PetscSectionSetDof(section, p, pspdim));
2162:     }
2163:   }

2165:   if (Ncopies > 1) {
2166:     Mat                 intMatScalar, allMatScalar;
2167:     PetscDualSpace      scalarsp;
2168:     PetscDualSpace_Lag *scalarlag;

2170:     PetscCall(PetscDualSpaceDuplicate(sp, &scalarsp));
2171:     /* Setting the number of components to Nk is a space with 1 copy of each k-form */
2172:     PetscCall(PetscDualSpaceSetNumComponents(scalarsp, Nk));
2173:     PetscCall(PetscDualSpaceSetUp(scalarsp));
2174:     PetscCall(PetscDualSpaceGetInteriorData(scalarsp, &(sp->intNodes), &intMatScalar));
2175:     PetscCall(PetscObjectReference((PetscObject)(sp->intNodes)));
2176:     if (intMatScalar) PetscCall(PetscDualSpaceLagrangeMatrixCreateCopies(intMatScalar, Nk, Ncopies, &(sp->intMat)));
2177:     PetscCall(PetscDualSpaceGetAllData(scalarsp, &(sp->allNodes), &allMatScalar));
2178:     PetscCall(PetscObjectReference((PetscObject)(sp->allNodes)));
2179:     PetscCall(PetscDualSpaceLagrangeMatrixCreateCopies(allMatScalar, Nk, Ncopies, &(sp->allMat)));
2180:     sp->spdim    = scalarsp->spdim * Ncopies;
2181:     sp->spintdim = scalarsp->spintdim * Ncopies;
2182:     scalarlag    = (PetscDualSpace_Lag *)scalarsp->data;
2183:     PetscCall(PetscLagNodeIndicesReference(scalarlag->vertIndices));
2184:     lag->vertIndices = scalarlag->vertIndices;
2185:     PetscCall(PetscLagNodeIndicesReference(scalarlag->intNodeIndices));
2186:     lag->intNodeIndices = scalarlag->intNodeIndices;
2187:     PetscCall(PetscLagNodeIndicesReference(scalarlag->allNodeIndices));
2188:     lag->allNodeIndices = scalarlag->allNodeIndices;
2189:     PetscCall(PetscDualSpaceDestroy(&scalarsp));
2190:     PetscCall(PetscSectionSetDof(section, 0, sp->spintdim));
2191:     PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, section));
2192:     PetscCall(PetscDualSpaceComputeFunctionalsFromAllData(sp));
2193:     PetscCall(PetscFree2(pStratStart, pStratEnd));
2194:     PetscCall(DMDestroy(&dmint));
2195:     PetscFunctionReturn(PETSC_SUCCESS);
2196:   }

2198:   if (trimmed && !continuous) {
2199:     /* the dofs of a trimmed space don't have a nice tensor/lattice structure:
2200:      * just construct the continuous dual space and copy all of the data over,
2201:      * allocating it all to the cell instead of splitting it up between the boundaries */
2202:     PetscDualSpace      spcont;
2203:     PetscInt            spdim, f;
2204:     PetscQuadrature     allNodes;
2205:     PetscDualSpace_Lag *lagc;
2206:     Mat                 allMat;

2208:     PetscCall(PetscDualSpaceDuplicate(sp, &spcont));
2209:     PetscCall(PetscDualSpaceLagrangeSetContinuity(spcont, PETSC_TRUE));
2210:     PetscCall(PetscDualSpaceSetUp(spcont));
2211:     PetscCall(PetscDualSpaceGetDimension(spcont, &spdim));
2212:     sp->spdim = sp->spintdim = spdim;
2213:     PetscCall(PetscSectionSetDof(section, 0, spdim));
2214:     PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, section));
2215:     PetscCall(PetscMalloc1(spdim, &(sp->functional)));
2216:     for (f = 0; f < spdim; f++) {
2217:       PetscQuadrature fn;

2219:       PetscCall(PetscDualSpaceGetFunctional(spcont, f, &fn));
2220:       PetscCall(PetscObjectReference((PetscObject)fn));
2221:       sp->functional[f] = fn;
2222:     }
2223:     PetscCall(PetscDualSpaceGetAllData(spcont, &allNodes, &allMat));
2224:     PetscCall(PetscObjectReference((PetscObject)allNodes));
2225:     PetscCall(PetscObjectReference((PetscObject)allNodes));
2226:     sp->allNodes = sp->intNodes = allNodes;
2227:     PetscCall(PetscObjectReference((PetscObject)allMat));
2228:     PetscCall(PetscObjectReference((PetscObject)allMat));
2229:     sp->allMat = sp->intMat = allMat;
2230:     lagc                    = (PetscDualSpace_Lag *)spcont->data;
2231:     PetscCall(PetscLagNodeIndicesReference(lagc->vertIndices));
2232:     lag->vertIndices = lagc->vertIndices;
2233:     PetscCall(PetscLagNodeIndicesReference(lagc->allNodeIndices));
2234:     PetscCall(PetscLagNodeIndicesReference(lagc->allNodeIndices));
2235:     lag->intNodeIndices = lagc->allNodeIndices;
2236:     lag->allNodeIndices = lagc->allNodeIndices;
2237:     PetscCall(PetscDualSpaceDestroy(&spcont));
2238:     PetscCall(PetscFree2(pStratStart, pStratEnd));
2239:     PetscCall(DMDestroy(&dmint));
2240:     PetscFunctionReturn(PETSC_SUCCESS);
2241:   }

2243:   /* step 3: construct intNodes, and intMat, and combine it with boundray data to make allNodes and allMat */
2244:   if (!tensorSpace) {
2245:     if (!tensorCell) PetscCall(PetscLagNodeIndicesCreateSimplexVertices(dm, &(lag->vertIndices)));

2247:     if (trimmed) {
2248:       /* there is one dof in the interior of the a trimmed element for each full polynomial of with degree at most
2249:        * order + k - dim - 1 */
2250:       if (order + PetscAbsInt(formDegree) > dim) {
2251:         PetscInt sum = order + PetscAbsInt(formDegree) - dim - 1;
2252:         PetscInt nDofs;

2254:         PetscCall(PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices)));
2255:         PetscCall(MatGetSize(sp->intMat, &nDofs, NULL));
2256:         PetscCall(PetscSectionSetDof(section, 0, nDofs));
2257:       }
2258:       PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, section));
2259:       PetscCall(PetscDualSpaceCreateAllDataFromInteriorData(sp));
2260:       PetscCall(PetscDualSpaceLagrangeCreateAllNodeIdx(sp));
2261:     } else {
2262:       if (!continuous) {
2263:         /* if discontinuous just construct one node for each set of dofs (a set of dofs is a basis for the k-form
2264:          * space) */
2265:         PetscInt sum = order;
2266:         PetscInt nDofs;

2268:         PetscCall(PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices)));
2269:         PetscCall(MatGetSize(sp->intMat, &nDofs, NULL));
2270:         PetscCall(PetscSectionSetDof(section, 0, nDofs));
2271:         PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, section));
2272:         PetscCall(PetscObjectReference((PetscObject)(sp->intNodes)));
2273:         sp->allNodes = sp->intNodes;
2274:         PetscCall(PetscObjectReference((PetscObject)(sp->intMat)));
2275:         sp->allMat = sp->intMat;
2276:         PetscCall(PetscLagNodeIndicesReference(lag->intNodeIndices));
2277:         lag->allNodeIndices = lag->intNodeIndices;
2278:       } else {
2279:         /* there is one dof in the interior of the a full element for each trimmed polynomial of with degree at most
2280:          * order + k - dim, but with complementary form degree */
2281:         if (order + PetscAbsInt(formDegree) > dim) {
2282:           PetscDualSpace      trimmedsp;
2283:           PetscDualSpace_Lag *trimmedlag;
2284:           PetscQuadrature     intNodes;
2285:           PetscInt            trFormDegree = formDegree >= 0 ? formDegree - dim : dim - PetscAbsInt(formDegree);
2286:           PetscInt            nDofs;
2287:           Mat                 intMat;

2289:           PetscCall(PetscDualSpaceDuplicate(sp, &trimmedsp));
2290:           PetscCall(PetscDualSpaceLagrangeSetTrimmed(trimmedsp, PETSC_TRUE));
2291:           PetscCall(PetscDualSpaceSetOrder(trimmedsp, order + PetscAbsInt(formDegree) - dim));
2292:           PetscCall(PetscDualSpaceSetFormDegree(trimmedsp, trFormDegree));
2293:           trimmedlag              = (PetscDualSpace_Lag *)trimmedsp->data;
2294:           trimmedlag->numNodeSkip = numNodeSkip + 1;
2295:           PetscCall(PetscDualSpaceSetUp(trimmedsp));
2296:           PetscCall(PetscDualSpaceGetAllData(trimmedsp, &intNodes, &intMat));
2297:           PetscCall(PetscObjectReference((PetscObject)intNodes));
2298:           sp->intNodes = intNodes;
2299:           PetscCall(PetscLagNodeIndicesReference(trimmedlag->allNodeIndices));
2300:           lag->intNodeIndices = trimmedlag->allNodeIndices;
2301:           PetscCall(PetscObjectReference((PetscObject)intMat));
2302:           if (PetscAbsInt(formDegree) > 0 && PetscAbsInt(formDegree) < dim) {
2303:             PetscReal   *T;
2304:             PetscScalar *work;
2305:             PetscInt     nCols, nRows;
2306:             Mat          intMatT;

2308:             PetscCall(MatDuplicate(intMat, MAT_COPY_VALUES, &intMatT));
2309:             PetscCall(MatGetSize(intMat, &nRows, &nCols));
2310:             PetscCall(PetscMalloc2(Nk * Nk, &T, nCols, &work));
2311:             PetscCall(BiunitSimplexSymmetricFormTransformation(dim, formDegree, T));
2312:             for (PetscInt row = 0; row < nRows; row++) {
2313:               PetscInt           nrCols;
2314:               const PetscInt    *rCols;
2315:               const PetscScalar *rVals;

2317:               PetscCall(MatGetRow(intMat, row, &nrCols, &rCols, &rVals));
2318:               PetscCheck(nrCols % Nk == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in intMat matrix are not in k-form size blocks");
2319:               for (PetscInt b = 0; b < nrCols; b += Nk) {
2320:                 const PetscScalar *v = &rVals[b];
2321:                 PetscScalar       *w = &work[b];
2322:                 for (PetscInt j = 0; j < Nk; j++) {
2323:                   w[j] = 0.;
2324:                   for (PetscInt i = 0; i < Nk; i++) w[j] += v[i] * T[i * Nk + j];
2325:                 }
2326:               }
2327:               PetscCall(MatSetValuesBlocked(intMatT, 1, &row, nrCols, rCols, work, INSERT_VALUES));
2328:               PetscCall(MatRestoreRow(intMat, row, &nrCols, &rCols, &rVals));
2329:             }
2330:             PetscCall(MatAssemblyBegin(intMatT, MAT_FINAL_ASSEMBLY));
2331:             PetscCall(MatAssemblyEnd(intMatT, MAT_FINAL_ASSEMBLY));
2332:             PetscCall(MatDestroy(&intMat));
2333:             intMat = intMatT;
2334:             PetscCall(PetscLagNodeIndicesDestroy(&(lag->intNodeIndices)));
2335:             PetscCall(PetscLagNodeIndicesDuplicate(trimmedlag->allNodeIndices, &(lag->intNodeIndices)));
2336:             {
2337:               PetscInt         nNodes     = lag->intNodeIndices->nNodes;
2338:               PetscReal       *newNodeVec = lag->intNodeIndices->nodeVec;
2339:               const PetscReal *oldNodeVec = trimmedlag->allNodeIndices->nodeVec;

2341:               for (PetscInt n = 0; n < nNodes; n++) {
2342:                 PetscReal       *w = &newNodeVec[n * Nk];
2343:                 const PetscReal *v = &oldNodeVec[n * Nk];

2345:                 for (PetscInt j = 0; j < Nk; j++) {
2346:                   w[j] = 0.;
2347:                   for (PetscInt i = 0; i < Nk; i++) w[j] += v[i] * T[i * Nk + j];
2348:                 }
2349:               }
2350:             }
2351:             PetscCall(PetscFree2(T, work));
2352:           }
2353:           sp->intMat = intMat;
2354:           PetscCall(MatGetSize(sp->intMat, &nDofs, NULL));
2355:           PetscCall(PetscDualSpaceDestroy(&trimmedsp));
2356:           PetscCall(PetscSectionSetDof(section, 0, nDofs));
2357:         }
2358:         PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, section));
2359:         PetscCall(PetscDualSpaceCreateAllDataFromInteriorData(sp));
2360:         PetscCall(PetscDualSpaceLagrangeCreateAllNodeIdx(sp));
2361:       }
2362:     }
2363:   } else {
2364:     PetscQuadrature     intNodesTrace  = NULL;
2365:     PetscQuadrature     intNodesFiber  = NULL;
2366:     PetscQuadrature     intNodes       = NULL;
2367:     PetscLagNodeIndices intNodeIndices = NULL;
2368:     Mat                 intMat         = NULL;

2370:     if (PetscAbsInt(formDegree) < dim) { /* get the trace k-forms on the first facet, and the 0-forms on the edge,
2371:                                             and wedge them together to create some of the k-form dofs */
2372:       PetscDualSpace      trace, fiber;
2373:       PetscDualSpace_Lag *tracel, *fiberl;
2374:       Mat                 intMatTrace, intMatFiber;

2376:       if (sp->pointSpaces[tensorf]) {
2377:         PetscCall(PetscObjectReference((PetscObject)(sp->pointSpaces[tensorf])));
2378:         trace = sp->pointSpaces[tensorf];
2379:       } else {
2380:         PetscCall(PetscDualSpaceCreateFacetSubspace_Lagrange(sp, NULL, tensorf, formDegree, Ncopies, PETSC_TRUE, &trace));
2381:       }
2382:       PetscCall(PetscDualSpaceCreateEdgeSubspace_Lagrange(sp, order, 0, 1, PETSC_TRUE, &fiber));
2383:       tracel = (PetscDualSpace_Lag *)trace->data;
2384:       fiberl = (PetscDualSpace_Lag *)fiber->data;
2385:       PetscCall(PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices)));
2386:       PetscCall(PetscDualSpaceGetInteriorData(trace, &intNodesTrace, &intMatTrace));
2387:       PetscCall(PetscDualSpaceGetInteriorData(fiber, &intNodesFiber, &intMatFiber));
2388:       if (intNodesTrace && intNodesFiber) {
2389:         PetscCall(PetscQuadratureCreateTensor(intNodesTrace, intNodesFiber, &intNodes));
2390:         PetscCall(MatTensorAltV(intMatTrace, intMatFiber, dim - 1, formDegree, 1, 0, &intMat));
2391:         PetscCall(PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, formDegree, fiberl->intNodeIndices, 1, 0, &intNodeIndices));
2392:       }
2393:       PetscCall(PetscObjectReference((PetscObject)intNodesTrace));
2394:       PetscCall(PetscObjectReference((PetscObject)intNodesFiber));
2395:       PetscCall(PetscDualSpaceDestroy(&fiber));
2396:       PetscCall(PetscDualSpaceDestroy(&trace));
2397:     }
2398:     if (PetscAbsInt(formDegree) > 0) { /* get the trace (k-1)-forms on the first facet, and the 1-forms on the edge,
2399:                                           and wedge them together to create the remaining k-form dofs */
2400:       PetscDualSpace      trace, fiber;
2401:       PetscDualSpace_Lag *tracel, *fiberl;
2402:       PetscQuadrature     intNodesTrace2, intNodesFiber2, intNodes2;
2403:       PetscLagNodeIndices intNodeIndices2;
2404:       Mat                 intMatTrace, intMatFiber, intMat2;
2405:       PetscInt            traceDegree = formDegree > 0 ? formDegree - 1 : formDegree + 1;
2406:       PetscInt            fiberDegree = formDegree > 0 ? 1 : -1;

2408:       PetscCall(PetscDualSpaceCreateFacetSubspace_Lagrange(sp, NULL, tensorf, traceDegree, Ncopies, PETSC_TRUE, &trace));
2409:       PetscCall(PetscDualSpaceCreateEdgeSubspace_Lagrange(sp, order, fiberDegree, 1, PETSC_TRUE, &fiber));
2410:       tracel = (PetscDualSpace_Lag *)trace->data;
2411:       fiberl = (PetscDualSpace_Lag *)fiber->data;
2412:       if (!lag->vertIndices) PetscCall(PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices)));
2413:       PetscCall(PetscDualSpaceGetInteriorData(trace, &intNodesTrace2, &intMatTrace));
2414:       PetscCall(PetscDualSpaceGetInteriorData(fiber, &intNodesFiber2, &intMatFiber));
2415:       if (intNodesTrace2 && intNodesFiber2) {
2416:         PetscCall(PetscQuadratureCreateTensor(intNodesTrace2, intNodesFiber2, &intNodes2));
2417:         PetscCall(MatTensorAltV(intMatTrace, intMatFiber, dim - 1, traceDegree, 1, fiberDegree, &intMat2));
2418:         PetscCall(PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, traceDegree, fiberl->intNodeIndices, 1, fiberDegree, &intNodeIndices2));
2419:         if (!intMat) {
2420:           intMat         = intMat2;
2421:           intNodes       = intNodes2;
2422:           intNodeIndices = intNodeIndices2;
2423:         } else {
2424:           /* merge the matrices, quadrature points, and nodes */
2425:           PetscInt            nM;
2426:           PetscInt            nDof, nDof2;
2427:           PetscInt           *toMerged = NULL, *toMerged2 = NULL;
2428:           PetscQuadrature     merged               = NULL;
2429:           PetscLagNodeIndices intNodeIndicesMerged = NULL;
2430:           Mat                 matMerged            = NULL;

2432:           PetscCall(MatGetSize(intMat, &nDof, NULL));
2433:           PetscCall(MatGetSize(intMat2, &nDof2, NULL));
2434:           PetscCall(PetscQuadraturePointsMerge(intNodes, intNodes2, &merged, &toMerged, &toMerged2));
2435:           PetscCall(PetscQuadratureGetData(merged, NULL, NULL, &nM, NULL, NULL));
2436:           PetscCall(MatricesMerge(intMat, intMat2, dim, formDegree, nM, toMerged, toMerged2, &matMerged));
2437:           PetscCall(PetscLagNodeIndicesMerge(intNodeIndices, intNodeIndices2, &intNodeIndicesMerged));
2438:           PetscCall(PetscFree(toMerged));
2439:           PetscCall(PetscFree(toMerged2));
2440:           PetscCall(MatDestroy(&intMat));
2441:           PetscCall(MatDestroy(&intMat2));
2442:           PetscCall(PetscQuadratureDestroy(&intNodes));
2443:           PetscCall(PetscQuadratureDestroy(&intNodes2));
2444:           PetscCall(PetscLagNodeIndicesDestroy(&intNodeIndices));
2445:           PetscCall(PetscLagNodeIndicesDestroy(&intNodeIndices2));
2446:           intNodes       = merged;
2447:           intMat         = matMerged;
2448:           intNodeIndices = intNodeIndicesMerged;
2449:           if (!trimmed) {
2450:             /* I think users expect that, when a node has a full basis for the k-forms,
2451:              * they should be consecutive dofs.  That isn't the case for trimmed spaces,
2452:              * but is for some of the nodes in untrimmed spaces, so in that case we
2453:              * sort them to group them by node */
2454:             Mat intMatPerm;

2456:             PetscCall(MatPermuteByNodeIdx(intMat, intNodeIndices, &intMatPerm));
2457:             PetscCall(MatDestroy(&intMat));
2458:             intMat = intMatPerm;
2459:           }
2460:         }
2461:       }
2462:       PetscCall(PetscDualSpaceDestroy(&fiber));
2463:       PetscCall(PetscDualSpaceDestroy(&trace));
2464:     }
2465:     PetscCall(PetscQuadratureDestroy(&intNodesTrace));
2466:     PetscCall(PetscQuadratureDestroy(&intNodesFiber));
2467:     sp->intNodes        = intNodes;
2468:     sp->intMat          = intMat;
2469:     lag->intNodeIndices = intNodeIndices;
2470:     {
2471:       PetscInt nDofs = 0;

2473:       if (intMat) PetscCall(MatGetSize(intMat, &nDofs, NULL));
2474:       PetscCall(PetscSectionSetDof(section, 0, nDofs));
2475:     }
2476:     PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, section));
2477:     if (continuous) {
2478:       PetscCall(PetscDualSpaceCreateAllDataFromInteriorData(sp));
2479:       PetscCall(PetscDualSpaceLagrangeCreateAllNodeIdx(sp));
2480:     } else {
2481:       PetscCall(PetscObjectReference((PetscObject)intNodes));
2482:       sp->allNodes = intNodes;
2483:       PetscCall(PetscObjectReference((PetscObject)intMat));
2484:       sp->allMat = intMat;
2485:       PetscCall(PetscLagNodeIndicesReference(intNodeIndices));
2486:       lag->allNodeIndices = intNodeIndices;
2487:     }
2488:   }
2489:   PetscCall(PetscSectionGetStorageSize(section, &sp->spdim));
2490:   PetscCall(PetscSectionGetConstrainedStorageSize(section, &sp->spintdim));
2491:   PetscCall(PetscDualSpaceComputeFunctionalsFromAllData(sp));
2492:   PetscCall(PetscFree2(pStratStart, pStratEnd));
2493:   PetscCall(DMDestroy(&dmint));
2494:   PetscFunctionReturn(PETSC_SUCCESS);
2495: }

2497: /* Create a matrix that represents the transformation that DMPlexVecGetClosure() would need
2498:  * to get the representation of the dofs for a mesh point if the mesh point had this orientation
2499:  * relative to the cell */
2500: PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat)
2501: {
2502:   PetscDualSpace_Lag *lag;
2503:   DM                  dm;
2504:   PetscLagNodeIndices vertIndices, intNodeIndices;
2505:   PetscLagNodeIndices ni;
2506:   PetscInt            nodeIdxDim, nodeVecDim, nNodes;
2507:   PetscInt            formDegree;
2508:   PetscInt           *perm, *permOrnt;
2509:   PetscInt           *nnz;
2510:   PetscInt            n;
2511:   PetscInt            maxGroupSize;
2512:   PetscScalar        *V, *W, *work;
2513:   Mat                 A;

2515:   PetscFunctionBegin;
2516:   if (!sp->spintdim) {
2517:     *symMat = NULL;
2518:     PetscFunctionReturn(PETSC_SUCCESS);
2519:   }
2520:   lag            = (PetscDualSpace_Lag *)sp->data;
2521:   vertIndices    = lag->vertIndices;
2522:   intNodeIndices = lag->intNodeIndices;
2523:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
2524:   PetscCall(PetscDualSpaceGetFormDegree(sp, &formDegree));
2525:   PetscCall(PetscNew(&ni));
2526:   ni->refct      = 1;
2527:   ni->nodeIdxDim = nodeIdxDim = intNodeIndices->nodeIdxDim;
2528:   ni->nodeVecDim = nodeVecDim = intNodeIndices->nodeVecDim;
2529:   ni->nNodes = nNodes = intNodeIndices->nNodes;
2530:   PetscCall(PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx)));
2531:   PetscCall(PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec)));
2532:   /* push forward the dofs by the symmetry of the reference element induced by ornt */
2533:   PetscCall(PetscLagNodeIndicesPushForward(dm, vertIndices, 0, vertIndices, intNodeIndices, ornt, formDegree, ni->nodeIdx, ni->nodeVec));
2534:   /* get the revlex order for both the original and transformed dofs */
2535:   PetscCall(PetscLagNodeIndicesGetPermutation(intNodeIndices, &perm));
2536:   PetscCall(PetscLagNodeIndicesGetPermutation(ni, &permOrnt));
2537:   PetscCall(PetscMalloc1(nNodes, &nnz));
2538:   for (n = 0, maxGroupSize = 0; n < nNodes;) { /* incremented in the loop */
2539:     PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2540:     PetscInt  m, nEnd;
2541:     PetscInt  groupSize;
2542:     /* for each group of dofs that have the same nodeIdx coordinate */
2543:     for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2544:       PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2545:       PetscInt  d;

2547:       /* compare the oriented permutation indices */
2548:       for (d = 0; d < nodeIdxDim; d++)
2549:         if (mind[d] != nind[d]) break;
2550:       if (d < nodeIdxDim) break;
2551:     }
2552:     /* permOrnt[[n, nEnd)] is a group of dofs that, under the symmetry are at the same location */

2554:     /* the symmetry had better map the group of dofs with the same permuted nodeIdx
2555:      * to a group of dofs with the same size, otherwise we messed up */
2556:     if (PetscDefined(USE_DEBUG)) {
2557:       PetscInt  m;
2558:       PetscInt *nind = &(intNodeIndices->nodeIdx[perm[n] * nodeIdxDim]);

2560:       for (m = n + 1; m < nEnd; m++) {
2561:         PetscInt *mind = &(intNodeIndices->nodeIdx[perm[m] * nodeIdxDim]);
2562:         PetscInt  d;

2564:         /* compare the oriented permutation indices */
2565:         for (d = 0; d < nodeIdxDim; d++)
2566:           if (mind[d] != nind[d]) break;
2567:         if (d < nodeIdxDim) break;
2568:       }
2569:       PetscCheck(m >= nEnd, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs with same index after symmetry not same block size");
2570:     }
2571:     groupSize = nEnd - n;
2572:     /* each pushforward dof vector will be expressed in a basis of the unpermuted dofs */
2573:     for (m = n; m < nEnd; m++) nnz[permOrnt[m]] = groupSize;

2575:     maxGroupSize = PetscMax(maxGroupSize, nEnd - n);
2576:     n            = nEnd;
2577:   }
2578:   PetscCheck(maxGroupSize <= nodeVecDim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs are not in blocks that can be solved");
2579:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes, nNodes, 0, nnz, &A));
2580:   PetscCall(PetscFree(nnz));
2581:   PetscCall(PetscMalloc3(maxGroupSize * nodeVecDim, &V, maxGroupSize * nodeVecDim, &W, nodeVecDim * 2, &work));
2582:   for (n = 0; n < nNodes;) { /* incremented in the loop */
2583:     PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2584:     PetscInt  nEnd;
2585:     PetscInt  m;
2586:     PetscInt  groupSize;
2587:     for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2588:       PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2589:       PetscInt  d;

2591:       /* compare the oriented permutation indices */
2592:       for (d = 0; d < nodeIdxDim; d++)
2593:         if (mind[d] != nind[d]) break;
2594:       if (d < nodeIdxDim) break;
2595:     }
2596:     groupSize = nEnd - n;
2597:     /* get all of the vectors from the original and all of the pushforward vectors */
2598:     for (m = n; m < nEnd; m++) {
2599:       PetscInt d;

2601:       for (d = 0; d < nodeVecDim; d++) {
2602:         V[(m - n) * nodeVecDim + d] = intNodeIndices->nodeVec[perm[m] * nodeVecDim + d];
2603:         W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
2604:       }
2605:     }
2606:     /* now we have to solve for W in terms of V: the systems isn't always square, but the span
2607:      * of V and W should always be the same, so the solution of the normal equations works */
2608:     {
2609:       char         transpose = 'N';
2610:       PetscBLASInt bm        = nodeVecDim;
2611:       PetscBLASInt bn        = groupSize;
2612:       PetscBLASInt bnrhs     = groupSize;
2613:       PetscBLASInt blda      = bm;
2614:       PetscBLASInt bldb      = bm;
2615:       PetscBLASInt blwork    = 2 * nodeVecDim;
2616:       PetscBLASInt info;

2618:       PetscCallBLAS("LAPACKgels", LAPACKgels_(&transpose, &bm, &bn, &bnrhs, V, &blda, W, &bldb, work, &blwork, &info));
2619:       PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELS");
2620:       /* repack */
2621:       {
2622:         PetscInt i, j;

2624:         for (i = 0; i < groupSize; i++) {
2625:           for (j = 0; j < groupSize; j++) {
2626:             /* notice the different leading dimension */
2627:             V[i * groupSize + j] = W[i * nodeVecDim + j];
2628:           }
2629:         }
2630:       }
2631:       if (PetscDefined(USE_DEBUG)) {
2632:         PetscReal res;

2634:         /* check that the normal error is 0 */
2635:         for (m = n; m < nEnd; m++) {
2636:           PetscInt d;

2638:           for (d = 0; d < nodeVecDim; d++) W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
2639:         }
2640:         res = 0.;
2641:         for (PetscInt i = 0; i < groupSize; i++) {
2642:           for (PetscInt j = 0; j < nodeVecDim; j++) {
2643:             for (PetscInt k = 0; k < groupSize; k++) W[i * nodeVecDim + j] -= V[i * groupSize + k] * intNodeIndices->nodeVec[perm[n + k] * nodeVecDim + j];
2644:             res += PetscAbsScalar(W[i * nodeVecDim + j]);
2645:           }
2646:         }
2647:         PetscCheck(res <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_LIB, "Dof block did not solve");
2648:       }
2649:     }
2650:     PetscCall(MatSetValues(A, groupSize, &permOrnt[n], groupSize, &perm[n], V, INSERT_VALUES));
2651:     n = nEnd;
2652:   }
2653:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
2654:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
2655:   *symMat = A;
2656:   PetscCall(PetscFree3(V, W, work));
2657:   PetscCall(PetscLagNodeIndicesDestroy(&ni));
2658:   PetscFunctionReturn(PETSC_SUCCESS);
2659: }

2661: #define BaryIndex(perEdge, a, b, c) (((b) * (2 * perEdge + 1 - (b))) / 2) + (c)

2663: #define CartIndex(perEdge, a, b) (perEdge * (a) + b)

2665: /* the existing interface for symmetries is insufficient for all cases:
2666:  * - it should be sufficient for form degrees that are scalar (0 and n)
2667:  * - it should be sufficient for hypercube dofs
2668:  * - it isn't sufficient for simplex cells with non-scalar form degrees if
2669:  *   there are any dofs in the interior
2670:  *
2671:  * We compute the general transformation matrices, and if they fit, we return them,
2672:  * otherwise we error (but we should probably change the interface to allow for
2673:  * these symmetries)
2674:  */
2675: static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
2676: {
2677:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2678:   PetscInt            dim, order, Nc;

2680:   PetscFunctionBegin;
2681:   PetscCall(PetscDualSpaceGetOrder(sp, &order));
2682:   PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
2683:   PetscCall(DMGetDimension(sp->dm, &dim));
2684:   if (!lag->symComputed) { /* store symmetries */
2685:     PetscInt       pStart, pEnd, p;
2686:     PetscInt       numPoints;
2687:     PetscInt       numFaces;
2688:     PetscInt       spintdim;
2689:     PetscInt    ***symperms;
2690:     PetscScalar ***symflips;

2692:     PetscCall(DMPlexGetChart(sp->dm, &pStart, &pEnd));
2693:     numPoints = pEnd - pStart;
2694:     {
2695:       DMPolytopeType ct;
2696:       /* The number of arrangements is no longer based on the number of faces */
2697:       PetscCall(DMPlexGetCellType(sp->dm, 0, &ct));
2698:       numFaces = DMPolytopeTypeGetNumArrangments(ct) / 2;
2699:     }
2700:     PetscCall(PetscCalloc1(numPoints, &symperms));
2701:     PetscCall(PetscCalloc1(numPoints, &symflips));
2702:     spintdim = sp->spintdim;
2703:     /* The nodal symmetry behavior is not present when tensorSpace != tensorCell: someone might want this for the "S"
2704:      * family of FEEC spaces.  Most used in particular are discontinuous polynomial L2 spaces in tensor cells, where
2705:      * the symmetries are not necessary for FE assembly.  So for now we assume this is the case and don't return
2706:      * symmetries if tensorSpace != tensorCell */
2707:     if (spintdim && 0 < dim && dim < 3 && (lag->tensorSpace == lag->tensorCell)) { /* compute self symmetries */
2708:       PetscInt    **cellSymperms;
2709:       PetscScalar **cellSymflips;
2710:       PetscInt      ornt;
2711:       PetscInt      nCopies = Nc / lag->intNodeIndices->nodeVecDim;
2712:       PetscInt      nNodes  = lag->intNodeIndices->nNodes;

2714:       lag->numSelfSym = 2 * numFaces;
2715:       lag->selfSymOff = numFaces;
2716:       PetscCall(PetscCalloc1(2 * numFaces, &cellSymperms));
2717:       PetscCall(PetscCalloc1(2 * numFaces, &cellSymflips));
2718:       /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
2719:       symperms[0] = &cellSymperms[numFaces];
2720:       symflips[0] = &cellSymflips[numFaces];
2721:       PetscCheck(lag->intNodeIndices->nodeVecDim * nCopies == Nc, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs");
2722:       PetscCheck(nNodes * nCopies == spintdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs");
2723:       for (ornt = -numFaces; ornt < numFaces; ornt++) { /* for every symmetry, compute the symmetry matrix, and extract rows to see if it fits in the perm + flip framework */
2724:         Mat          symMat;
2725:         PetscInt    *perm;
2726:         PetscScalar *flips;
2727:         PetscInt     i;

2729:         if (!ornt) continue;
2730:         PetscCall(PetscMalloc1(spintdim, &perm));
2731:         PetscCall(PetscCalloc1(spintdim, &flips));
2732:         for (i = 0; i < spintdim; i++) perm[i] = -1;
2733:         PetscCall(PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, ornt, &symMat));
2734:         for (i = 0; i < nNodes; i++) {
2735:           PetscInt           ncols;
2736:           PetscInt           j, k;
2737:           const PetscInt    *cols;
2738:           const PetscScalar *vals;
2739:           PetscBool          nz_seen = PETSC_FALSE;

2741:           PetscCall(MatGetRow(symMat, i, &ncols, &cols, &vals));
2742:           for (j = 0; j < ncols; j++) {
2743:             if (PetscAbsScalar(vals[j]) > PETSC_SMALL) {
2744:               PetscCheck(!nz_seen, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2745:               nz_seen = PETSC_TRUE;
2746:               PetscCheck(PetscAbsReal(PetscAbsScalar(vals[j]) - PetscRealConstant(1.)) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2747:               PetscCheck(PetscAbsReal(PetscImaginaryPart(vals[j])) <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2748:               PetscCheck(perm[cols[j] * nCopies] < 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2749:               for (k = 0; k < nCopies; k++) perm[cols[j] * nCopies + k] = i * nCopies + k;
2750:               if (PetscRealPart(vals[j]) < 0.) {
2751:                 for (k = 0; k < nCopies; k++) flips[i * nCopies + k] = -1.;
2752:               } else {
2753:                 for (k = 0; k < nCopies; k++) flips[i * nCopies + k] = 1.;
2754:               }
2755:             }
2756:           }
2757:           PetscCall(MatRestoreRow(symMat, i, &ncols, &cols, &vals));
2758:         }
2759:         PetscCall(MatDestroy(&symMat));
2760:         /* if there were no sign flips, keep NULL */
2761:         for (i = 0; i < spintdim; i++)
2762:           if (flips[i] != 1.) break;
2763:         if (i == spintdim) {
2764:           PetscCall(PetscFree(flips));
2765:           flips = NULL;
2766:         }
2767:         /* if the permutation is identity, keep NULL */
2768:         for (i = 0; i < spintdim; i++)
2769:           if (perm[i] != i) break;
2770:         if (i == spintdim) {
2771:           PetscCall(PetscFree(perm));
2772:           perm = NULL;
2773:         }
2774:         symperms[0][ornt] = perm;
2775:         symflips[0][ornt] = flips;
2776:       }
2777:       /* if no orientations produced non-identity permutations, keep NULL */
2778:       for (ornt = -numFaces; ornt < numFaces; ornt++)
2779:         if (symperms[0][ornt]) break;
2780:       if (ornt == numFaces) {
2781:         PetscCall(PetscFree(cellSymperms));
2782:         symperms[0] = NULL;
2783:       }
2784:       /* if no orientations produced sign flips, keep NULL */
2785:       for (ornt = -numFaces; ornt < numFaces; ornt++)
2786:         if (symflips[0][ornt]) break;
2787:       if (ornt == numFaces) {
2788:         PetscCall(PetscFree(cellSymflips));
2789:         symflips[0] = NULL;
2790:       }
2791:     }
2792:     { /* get the symmetries of closure points */
2793:       PetscInt  closureSize = 0;
2794:       PetscInt *closure     = NULL;
2795:       PetscInt  r;

2797:       PetscCall(DMPlexGetTransitiveClosure(sp->dm, 0, PETSC_TRUE, &closureSize, &closure));
2798:       for (r = 0; r < closureSize; r++) {
2799:         PetscDualSpace       psp;
2800:         PetscInt             point = closure[2 * r];
2801:         PetscInt             pspintdim;
2802:         const PetscInt    ***psymperms = NULL;
2803:         const PetscScalar ***psymflips = NULL;

2805:         if (!point) continue;
2806:         PetscCall(PetscDualSpaceGetPointSubspace(sp, point, &psp));
2807:         if (!psp) continue;
2808:         PetscCall(PetscDualSpaceGetInteriorDimension(psp, &pspintdim));
2809:         if (!pspintdim) continue;
2810:         PetscCall(PetscDualSpaceGetSymmetries(psp, &psymperms, &psymflips));
2811:         symperms[r] = (PetscInt **)(psymperms ? psymperms[0] : NULL);
2812:         symflips[r] = (PetscScalar **)(psymflips ? psymflips[0] : NULL);
2813:       }
2814:       PetscCall(DMPlexRestoreTransitiveClosure(sp->dm, 0, PETSC_TRUE, &closureSize, &closure));
2815:     }
2816:     for (p = 0; p < pEnd; p++)
2817:       if (symperms[p]) break;
2818:     if (p == pEnd) {
2819:       PetscCall(PetscFree(symperms));
2820:       symperms = NULL;
2821:     }
2822:     for (p = 0; p < pEnd; p++)
2823:       if (symflips[p]) break;
2824:     if (p == pEnd) {
2825:       PetscCall(PetscFree(symflips));
2826:       symflips = NULL;
2827:     }
2828:     lag->symperms    = symperms;
2829:     lag->symflips    = symflips;
2830:     lag->symComputed = PETSC_TRUE;
2831:   }
2832:   if (perms) *perms = (const PetscInt ***)lag->symperms;
2833:   if (flips) *flips = (const PetscScalar ***)lag->symflips;
2834:   PetscFunctionReturn(PETSC_SUCCESS);
2835: }

2837: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2838: {
2839:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2841:   PetscFunctionBegin;
2844:   *continuous = lag->continuous;
2845:   PetscFunctionReturn(PETSC_SUCCESS);
2846: }

2848: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2849: {
2850:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2852:   PetscFunctionBegin;
2854:   lag->continuous = continuous;
2855:   PetscFunctionReturn(PETSC_SUCCESS);
2856: }

2858: /*@
2859:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

2861:   Not Collective

2863:   Input Parameter:
2864: . sp         - the `PetscDualSpace`

2866:   Output Parameter:
2867: . continuous - flag for element continuity

2869:   Level: intermediate

2871: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDualSpaceLagrangeSetContinuity()`
2872: @*/
2873: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2874: {
2875:   PetscFunctionBegin;
2878:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace, PetscBool *), (sp, continuous));
2879:   PetscFunctionReturn(PETSC_SUCCESS);
2880: }

2882: /*@
2883:   PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous

2885:   Logically Collective

2887:   Input Parameters:
2888: + sp         - the `PetscDualSpace`
2889: - continuous - flag for element continuity

2891:   Options Database Key:
2892: . -petscdualspace_lagrange_continuity <bool> - use a continuous element

2894:   Level: intermediate

2896: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDualSpaceLagrangeGetContinuity()`
2897: @*/
2898: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2899: {
2900:   PetscFunctionBegin;
2903:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace, PetscBool), (sp, continuous));
2904:   PetscFunctionReturn(PETSC_SUCCESS);
2905: }

2907: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
2908: {
2909:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2911:   PetscFunctionBegin;
2912:   *tensor = lag->tensorSpace;
2913:   PetscFunctionReturn(PETSC_SUCCESS);
2914: }

2916: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
2917: {
2918:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2920:   PetscFunctionBegin;
2921:   lag->tensorSpace = tensor;
2922:   PetscFunctionReturn(PETSC_SUCCESS);
2923: }

2925: static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Lagrange(PetscDualSpace sp, PetscBool *trimmed)
2926: {
2927:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2929:   PetscFunctionBegin;
2930:   *trimmed = lag->trimmed;
2931:   PetscFunctionReturn(PETSC_SUCCESS);
2932: }

2934: static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp, PetscBool trimmed)
2935: {
2936:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2938:   PetscFunctionBegin;
2939:   lag->trimmed = trimmed;
2940:   PetscFunctionReturn(PETSC_SUCCESS);
2941: }

2943: static PetscErrorCode PetscDualSpaceLagrangeGetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
2944: {
2945:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2947:   PetscFunctionBegin;
2948:   if (nodeType) *nodeType = lag->nodeType;
2949:   if (boundary) *boundary = lag->endNodes;
2950:   if (exponent) *exponent = lag->nodeExponent;
2951:   PetscFunctionReturn(PETSC_SUCCESS);
2952: }

2954: static PetscErrorCode PetscDualSpaceLagrangeSetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
2955: {
2956:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2958:   PetscFunctionBegin;
2959:   PetscCheck(nodeType != PETSCDTNODES_GAUSSJACOBI || exponent > -1., PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Exponent must be > -1");
2960:   lag->nodeType     = nodeType;
2961:   lag->endNodes     = boundary;
2962:   lag->nodeExponent = exponent;
2963:   PetscFunctionReturn(PETSC_SUCCESS);
2964: }

2966: static PetscErrorCode PetscDualSpaceLagrangeGetUseMoments_Lagrange(PetscDualSpace sp, PetscBool *useMoments)
2967: {
2968:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2970:   PetscFunctionBegin;
2971:   *useMoments = lag->useMoments;
2972:   PetscFunctionReturn(PETSC_SUCCESS);
2973: }

2975: static PetscErrorCode PetscDualSpaceLagrangeSetUseMoments_Lagrange(PetscDualSpace sp, PetscBool useMoments)
2976: {
2977:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2979:   PetscFunctionBegin;
2980:   lag->useMoments = useMoments;
2981:   PetscFunctionReturn(PETSC_SUCCESS);
2982: }

2984: static PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt *momentOrder)
2985: {
2986:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2988:   PetscFunctionBegin;
2989:   *momentOrder = lag->momentOrder;
2990:   PetscFunctionReturn(PETSC_SUCCESS);
2991: }

2993: static PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt momentOrder)
2994: {
2995:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2997:   PetscFunctionBegin;
2998:   lag->momentOrder = momentOrder;
2999:   PetscFunctionReturn(PETSC_SUCCESS);
3000: }

3002: /*@
3003:   PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space

3005:   Not Collective

3007:   Input Parameter:
3008: . sp - The `PetscDualSpace`

3010:   Output Parameter:
3011: . tensor - Whether the dual space has tensor layout (vs. simplicial)

3013:   Level: intermediate

3015: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDualSpaceLagrangeSetTensor()`, `PetscDualSpaceCreate()`
3016: @*/
3017: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
3018: {
3019:   PetscFunctionBegin;
3022:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetTensor_C", (PetscDualSpace, PetscBool *), (sp, tensor));
3023:   PetscFunctionReturn(PETSC_SUCCESS);
3024: }

3026: /*@
3027:   PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space

3029:   Not Collective

3031:   Input Parameters:
3032: + sp - The `PetscDualSpace`
3033: - tensor - Whether the dual space has tensor layout (vs. simplicial)

3035:   Level: intermediate

3037: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDualSpaceLagrangeGetTensor()`, `PetscDualSpaceCreate()`
3038: @*/
3039: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
3040: {
3041:   PetscFunctionBegin;
3043:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetTensor_C", (PetscDualSpace, PetscBool), (sp, tensor));
3044:   PetscFunctionReturn(PETSC_SUCCESS);
3045: }

3047: /*@
3048:   PetscDualSpaceLagrangeGetTrimmed - Get the trimmed nature of the dual space

3050:   Not Collective

3052:   Input Parameter:
3053: . sp - The `PetscDualSpace`

3055:   Output Parameter:
3056: . trimmed - Whether the dual space represents to dual basis of a trimmed polynomial space (e.g. Raviart-Thomas and higher order / other form degree variants)

3058:   Level: intermediate

3060: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDualSpaceLagrangeSetTrimmed()`, `PetscDualSpaceCreate()`
3061: @*/
3062: PetscErrorCode PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp, PetscBool *trimmed)
3063: {
3064:   PetscFunctionBegin;
3067:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetTrimmed_C", (PetscDualSpace, PetscBool *), (sp, trimmed));
3068:   PetscFunctionReturn(PETSC_SUCCESS);
3069: }

3071: /*@
3072:   PetscDualSpaceLagrangeSetTrimmed - Set the trimmed nature of the dual space

3074:   Not Collective

3076:   Input Parameters:
3077: + sp - The `PetscDualSpace`
3078: - trimmed - Whether the dual space represents to dual basis of a trimmed polynomial space (e.g. Raviart-Thomas and higher order / other form degree variants)

3080:   Level: intermediate

3082: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDualSpaceLagrangeGetTrimmed()`, `PetscDualSpaceCreate()`
3083: @*/
3084: PetscErrorCode PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp, PetscBool trimmed)
3085: {
3086:   PetscFunctionBegin;
3088:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetTrimmed_C", (PetscDualSpace, PetscBool), (sp, trimmed));
3089:   PetscFunctionReturn(PETSC_SUCCESS);
3090: }

3092: /*@
3093:   PetscDualSpaceLagrangeGetNodeType - Get a description of how nodes are laid out for Lagrange polynomials in this
3094:   dual space

3096:   Not Collective

3098:   Input Parameter:
3099: . sp - The `PetscDualSpace`

3101:   Output Parameters:
3102: + nodeType - The type of nodes
3103: . boundary - Whether the node type is one that includes endpoints (if nodeType is `PETSCDTNODES_GAUSSJACOBI`, nodes that
3104:              include the boundary are Gauss-Lobatto-Jacobi nodes)
3105: - exponent - If nodeType is `PETSCDTNODES_GAUSJACOBI`, indicates the exponent used for both ends of the 1D Jacobi weight function
3106:              '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type

3108:   Level: advanced

3110: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDTNodeType`, `PetscDualSpaceLagrangeSetNodeType()`
3111: @*/
3112: PetscErrorCode PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
3113: {
3114:   PetscFunctionBegin;
3119:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetNodeType_C", (PetscDualSpace, PetscDTNodeType *, PetscBool *, PetscReal *), (sp, nodeType, boundary, exponent));
3120:   PetscFunctionReturn(PETSC_SUCCESS);
3121: }

3123: /*@
3124:   PetscDualSpaceLagrangeSetNodeType - Set a description of how nodes are laid out for Lagrange polynomials in this
3125:   dual space

3127:   Logically Collective

3129:   Input Parameters:
3130: + sp - The `PetscDualSpace`
3131: . nodeType - The type of nodes
3132: . boundary - Whether the node type is one that includes endpoints (if nodeType is `PETSCDTNODES_GAUSSJACOBI`, nodes that
3133:              include the boundary are Gauss-Lobatto-Jacobi nodes)
3134: - exponent - If nodeType is `PETSCDTNODES_GAUSJACOBI`, indicates the exponent used for both ends of the 1D Jacobi weight function
3135:              '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type

3137:   Level: advanced

3139: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDTNodeType`, `PetscDualSpaceLagrangeGetNodeType()`
3140: @*/
3141: PetscErrorCode PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
3142: {
3143:   PetscFunctionBegin;
3145:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetNodeType_C", (PetscDualSpace, PetscDTNodeType, PetscBool, PetscReal), (sp, nodeType, boundary, exponent));
3146:   PetscFunctionReturn(PETSC_SUCCESS);
3147: }

3149: /*@
3150:   PetscDualSpaceLagrangeGetUseMoments - Get the flag for using moment functionals

3152:   Not Collective

3154:   Input Parameter:
3155: . sp - The `PetscDualSpace`

3157:   Output Parameter:
3158: . useMoments - Moment flag

3160:   Level: advanced

3162: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDualSpaceLagrangeSetUseMoments()`
3163: @*/
3164: PetscErrorCode PetscDualSpaceLagrangeGetUseMoments(PetscDualSpace sp, PetscBool *useMoments)
3165: {
3166:   PetscFunctionBegin;
3169:   PetscUseMethod(sp, "PetscDualSpaceLagrangeGetUseMoments_C", (PetscDualSpace, PetscBool *), (sp, useMoments));
3170:   PetscFunctionReturn(PETSC_SUCCESS);
3171: }

3173: /*@
3174:   PetscDualSpaceLagrangeSetUseMoments - Set the flag for moment functionals

3176:   Logically Collective

3178:   Input Parameters:
3179: + sp - The `PetscDualSpace`
3180: - useMoments - The flag for moment functionals

3182:   Level: advanced

3184: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDualSpaceLagrangeGetUseMoments()`
3185: @*/
3186: PetscErrorCode PetscDualSpaceLagrangeSetUseMoments(PetscDualSpace sp, PetscBool useMoments)
3187: {
3188:   PetscFunctionBegin;
3190:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetUseMoments_C", (PetscDualSpace, PetscBool), (sp, useMoments));
3191:   PetscFunctionReturn(PETSC_SUCCESS);
3192: }

3194: /*@
3195:   PetscDualSpaceLagrangeGetMomentOrder - Get the order for moment integration

3197:   Not Collective

3199:   Input Parameter:
3200: . sp - The `PetscDualSpace`

3202:   Output Parameter:
3203: . order - Moment integration order

3205:   Level: advanced

3207: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDualSpaceLagrangeSetMomentOrder()`
3208: @*/
3209: PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder(PetscDualSpace sp, PetscInt *order)
3210: {
3211:   PetscFunctionBegin;
3214:   PetscUseMethod(sp, "PetscDualSpaceLagrangeGetMomentOrder_C", (PetscDualSpace, PetscInt *), (sp, order));
3215:   PetscFunctionReturn(PETSC_SUCCESS);
3216: }

3218: /*@
3219:   PetscDualSpaceLagrangeSetMomentOrder - Set the order for moment integration

3221:   Logically Collective

3223:   Input Parameters:
3224: + sp - The `PetscDualSpace`
3225: - order - The order for moment integration

3227:   Level: advanced

3229: .seealso: `PETSCDUALSPACELAGRANGE`, `PetscDualSpace`, `PetscDualSpaceLagrangeGetMomentOrder()`
3230: @*/
3231: PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder(PetscDualSpace sp, PetscInt order)
3232: {
3233:   PetscFunctionBegin;
3235:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetMomentOrder_C", (PetscDualSpace, PetscInt), (sp, order));
3236:   PetscFunctionReturn(PETSC_SUCCESS);
3237: }

3239: static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
3240: {
3241:   PetscFunctionBegin;
3242:   sp->ops->destroy              = PetscDualSpaceDestroy_Lagrange;
3243:   sp->ops->view                 = PetscDualSpaceView_Lagrange;
3244:   sp->ops->setfromoptions       = PetscDualSpaceSetFromOptions_Lagrange;
3245:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Lagrange;
3246:   sp->ops->setup                = PetscDualSpaceSetUp_Lagrange;
3247:   sp->ops->createheightsubspace = NULL;
3248:   sp->ops->createpointsubspace  = NULL;
3249:   sp->ops->getsymmetries        = PetscDualSpaceGetSymmetries_Lagrange;
3250:   sp->ops->apply                = PetscDualSpaceApplyDefault;
3251:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
3252:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
3253:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
3254:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;
3255:   PetscFunctionReturn(PETSC_SUCCESS);
3256: }

3258: /*MC
3259:   PETSCDUALSPACELAGRANGE = "lagrange" - A `PetscDualSpaceType` that encapsulates a dual space of pointwise evaluation functionals

3261:   Level: intermediate

3263:   Developer Note:
3264:   This `PetscDualSpace` seems to manage directly trimmed and untrimmed polynomials as well as tensor and non-tensor polynomials while for `PetscSpace` there seems to
3265:   be different `PetscSpaceType` for them.

3267: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`,
3268:           `PetscDualSpaceLagrangeSetMomentOrder()`, `PetscDualSpaceLagrangeGetMomentOrder()`, `PetscDualSpaceLagrangeSetUseMoments()`, `PetscDualSpaceLagrangeGetUseMoments()`,
3269:           `PetscDualSpaceLagrangeSetNodeType, PetscDualSpaceLagrangeGetNodeType, PetscDualSpaceLagrangeGetContinuity, PetscDualSpaceLagrangeSetContinuity,
3270:           `PetscDualSpaceLagrangeGetTensor()`, `PetscDualSpaceLagrangeSetTensor()`, `PetscDualSpaceLagrangeGetTrimmed()`, `PetscDualSpaceLagrangeSetTrimmed()`
3271: M*/
3272: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
3273: {
3274:   PetscDualSpace_Lag *lag;

3276:   PetscFunctionBegin;
3278:   PetscCall(PetscNew(&lag));
3279:   sp->data = lag;

3281:   lag->tensorCell  = PETSC_FALSE;
3282:   lag->tensorSpace = PETSC_FALSE;
3283:   lag->continuous  = PETSC_TRUE;
3284:   lag->numCopies   = PETSC_DEFAULT;
3285:   lag->numNodeSkip = PETSC_DEFAULT;
3286:   lag->nodeType    = PETSCDTNODES_DEFAULT;
3287:   lag->useMoments  = PETSC_FALSE;
3288:   lag->momentOrder = 0;

3290:   PetscCall(PetscDualSpaceInitialize_Lagrange(sp));
3291:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange));
3292:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange));
3293:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange));
3294:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange));
3295:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Lagrange));
3296:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Lagrange));
3297:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Lagrange));
3298:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Lagrange));
3299:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", PetscDualSpaceLagrangeGetUseMoments_Lagrange));
3300:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", PetscDualSpaceLagrangeSetUseMoments_Lagrange));
3301:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", PetscDualSpaceLagrangeGetMomentOrder_Lagrange));
3302:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", PetscDualSpaceLagrangeSetMomentOrder_Lagrange));
3303:   PetscFunctionReturn(PETSC_SUCCESS);
3304: }