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:   PetscNew(&f);
 23:   switch (family) {
 24:   case PETSCDTNODES_GAUSSJACOBI:
 25:   case PETSCDTNODES_EQUISPACED:
 26:     f->nodeFamily = family;
 27:     break;
 28:   default:
 29:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family");
 30:   }
 31:   f->endpoints      = endpoints;
 32:   f->gaussJacobiExp = 0.;
 33:   if (family == PETSCDTNODES_GAUSSJACOBI) {
 35:     f->gaussJacobiExp = gaussJacobiExp;
 36:   }
 37:   f->refct = 1;
 38:   *nf      = f;
 39:   return 0;
 40: }

 42: static PetscErrorCode Petsc1DNodeFamilyReference(Petsc1DNodeFamily nf)
 43: {
 44:   if (nf) nf->refct++;
 45:   return 0;
 46: }

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

 52:   if (!(*nf)) return 0;
 53:   if (--(*nf)->refct > 0) {
 54:     *nf = NULL;
 55:     return 0;
 56:   }
 57:   nc = (*nf)->nComputed;
 58:   for (i = 0; i < nc; i++) PetscFree((*nf)->nodesets[i]);
 59:   PetscFree((*nf)->nodesets);
 60:   PetscFree(*nf);
 61:   *nf = NULL;
 62:   return 0;
 63: }

 65: static PetscErrorCode Petsc1DNodeFamilyGetNodeSets(Petsc1DNodeFamily f, PetscInt degree, PetscReal ***nodesets)
 66: {
 67:   PetscInt nc;

 69:   nc = f->nComputed;
 70:   if (degree >= nc) {
 71:     PetscInt    i, j;
 72:     PetscReal **new_nodesets;
 73:     PetscReal  *w;

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

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

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

130:       for (j = 0; j < dim + 1; j++) tup[dim + 1 + j] = tup[j + (j >= i)];
131:       PetscNodeRecursive_Internal(dim - 1, degree - tup[i], nodesets, &tup[dim + 1], &node[dim + 1]);
132:       for (j = 0; j < dim + 1; j++) node[j + (j >= i)] += wi * node[dim + 1 + j];
133:       w += wi;
134:     }
135:     for (i = 0; i < dim + 1; i++) node[i] /= w;
136:   }
137:   return 0;
138: }

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

153:   if (!dim) return 0;
154:   PetscCalloc1(dim + 2, &tup);
155:   k = 0;
156:   PetscDTBinomialInt(degree + dim, dim, &npoints);
157:   Petsc1DNodeFamilyGetNodeSets(f, degree, &nodesets);
158:   worksize = ((dim + 2) * (dim + 3)) / 2;
159:   PetscMalloc2(worksize, &nodework, worksize, &tupwork);
160:   /* loop over the tuples of length dim with sum at most degree */
161:   for (k = 0; k < npoints; k++) {
162:     PetscInt i;

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

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

315: /* this is just here so I can access the values in tests/ex1.c outside the library */
316: PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices ni, PetscInt *nodeIdxDim, PetscInt *nodeVecDim, PetscInt *nNodes, const PetscInt *nodeIdx[], const PetscReal *nodeVec[])
317: {
318:   *nodeIdxDim = ni->nodeIdxDim;
319:   *nodeVecDim = ni->nodeVecDim;
320:   *nNodes     = ni->nNodes;
321:   *nodeIdx    = ni->nodeIdx;
322:   *nodeVec    = ni->nodeVec;
323:   return 0;
324: }

326: static PetscErrorCode PetscLagNodeIndicesReference(PetscLagNodeIndices ni)
327: {
328:   if (ni) ni->refct++;
329:   return 0;
330: }

332: static PetscErrorCode PetscLagNodeIndicesDuplicate(PetscLagNodeIndices ni, PetscLagNodeIndices *niNew)
333: {
334:   PetscNew(niNew);
335:   (*niNew)->refct      = 1;
336:   (*niNew)->nodeIdxDim = ni->nodeIdxDim;
337:   (*niNew)->nodeVecDim = ni->nodeVecDim;
338:   (*niNew)->nNodes     = ni->nNodes;
339:   PetscMalloc1(ni->nNodes * ni->nodeIdxDim, &((*niNew)->nodeIdx));
340:   PetscArraycpy((*niNew)->nodeIdx, ni->nodeIdx, ni->nNodes * ni->nodeIdxDim);
341:   PetscMalloc1(ni->nNodes * ni->nodeVecDim, &((*niNew)->nodeVec));
342:   PetscArraycpy((*niNew)->nodeVec, ni->nodeVec, ni->nNodes * ni->nodeVecDim);
343:   (*niNew)->perm = NULL;
344:   return 0;
345: }

347: static PetscErrorCode PetscLagNodeIndicesDestroy(PetscLagNodeIndices *ni)
348: {
349:   if (!(*ni)) return 0;
350:   if (--(*ni)->refct > 0) {
351:     *ni = NULL;
352:     return 0;
353:   }
354:   PetscFree((*ni)->nodeIdx);
355:   PetscFree((*ni)->nodeVec);
356:   PetscFree((*ni)->perm);
357:   PetscFree(*ni);
358:   *ni = NULL;
359:   return 0;
360: }

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

386:   DMGetDimension(dm, &dim);
387:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
388:   nVerts = vEnd - vStart;
389:   PetscMalloc1(nVerts, &closureOrder);
390:   PetscMalloc1(nVerts, &invClosureOrder);
391:   PetscMalloc1(nVerts, &revlexOrder);
392:   if (sortIdx) { /* bubble sort nodeIdx into revlex order */
393:     PetscInt  nodeIdxDim = ni->nodeIdxDim;
394:     PetscInt *idxOrder;

396:     PetscMalloc1(nVerts * nodeIdxDim, &newNodeIdx);
397:     PetscMalloc1(nVerts, &idxOrder);
398:     for (v = 0; v < nVerts; v++) idxOrder[v] = v;
399:     for (v = 0; v < nVerts; v++) {
400:       for (w = v + 1; w < nVerts; w++) {
401:         const PetscInt *iv   = &(ni->nodeIdx[idxOrder[v] * nodeIdxDim]);
402:         const PetscInt *iw   = &(ni->nodeIdx[idxOrder[w] * nodeIdxDim]);
403:         PetscInt        diff = 0;

405:         for (d = nodeIdxDim - 1; d >= 0; d--)
406:           if ((diff = (iv[d] - iw[d]))) break;
407:         if (diff > 0) {
408:           PetscInt swap = idxOrder[v];

410:           idxOrder[v] = idxOrder[w];
411:           idxOrder[w] = swap;
412:         }
413:       }
414:     }
415:     for (v = 0; v < nVerts; v++) {
416:       for (d = 0; d < nodeIdxDim; d++) newNodeIdx[v * ni->nodeIdxDim + d] = ni->nodeIdx[idxOrder[v] * nodeIdxDim + d];
417:     }
418:     PetscFree(ni->nodeIdx);
419:     ni->nodeIdx = newNodeIdx;
420:     newNodeIdx  = NULL;
421:     PetscFree(idxOrder);
422:   }
423:   DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);
424:   c = closureSize - nVerts;
425:   for (v = 0; v < nVerts; v++) closureOrder[v] = closure[2 * (c + v)] - vStart;
426:   for (v = 0; v < nVerts; v++) invClosureOrder[closureOrder[v]] = v;
427:   DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);
428:   DMGetCoordinatesLocal(dm, &coordVec);
429:   VecGetArrayRead(coordVec, &coords);
430:   /* bubble sort closure vertices by coordinates in revlex order */
431:   for (v = 0; v < nVerts; v++) revlexOrder[v] = v;
432:   for (v = 0; v < nVerts; v++) {
433:     for (w = v + 1; w < nVerts; w++) {
434:       const PetscScalar *cv   = &coords[closureOrder[revlexOrder[v]] * dim];
435:       const PetscScalar *cw   = &coords[closureOrder[revlexOrder[w]] * dim];
436:       PetscReal          diff = 0;

438:       for (d = dim - 1; d >= 0; d--)
439:         if ((diff = PetscRealPart(cv[d] - cw[d])) != 0.) break;
440:       if (diff > 0.) {
441:         PetscInt swap = revlexOrder[v];

443:         revlexOrder[v] = revlexOrder[w];
444:         revlexOrder[w] = swap;
445:       }
446:     }
447:   }
448:   VecRestoreArrayRead(coordVec, &coords);
449:   PetscMalloc1(ni->nodeIdxDim * nVerts, &newNodeIdx);
450:   /* reorder nodeIdx to be in closure order */
451:   for (v = 0; v < nVerts; v++) {
452:     for (d = 0; d < ni->nodeIdxDim; d++) newNodeIdx[revlexOrder[v] * ni->nodeIdxDim + d] = ni->nodeIdx[v * ni->nodeIdxDim + d];
453:   }
454:   PetscFree(ni->nodeIdx);
455:   ni->nodeIdx = newNodeIdx;
456:   ni->perm    = invClosureOrder;
457:   PetscFree(revlexOrder);
458:   PetscFree(closureOrder);
459:   return 0;
460: }

462: /* the coordinates of the simplex vertices are the corners of the barycentric simplex.
463:  * When we stack them on top of each other in revlex order, they look like the identity matrix */
464: static PetscErrorCode PetscLagNodeIndicesCreateSimplexVertices(DM dm, PetscLagNodeIndices *nodeIndices)
465: {
466:   PetscLagNodeIndices ni;
467:   PetscInt            dim, d;

469:   PetscNew(&ni);
470:   DMGetDimension(dm, &dim);
471:   ni->nodeIdxDim = dim + 1;
472:   ni->nodeVecDim = 0;
473:   ni->nNodes     = dim + 1;
474:   ni->refct      = 1;
475:   PetscCalloc1((dim + 1) * (dim + 1), &(ni->nodeIdx));
476:   for (d = 0; d < dim + 1; d++) ni->nodeIdx[d * (dim + 2)] = 1;
477:   PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_FALSE);
478:   *nodeIndices = ni;
479:   return 0;
480: }

482: /* A polytope that is a tensor product of a facet and a segment.
483:  * We take whatever coordinate system was being used for the facet
484:  * and we concatenate the barycentric coordinates for the vertices
485:  * at the end of the segment, (1,0) and (0,1), to get a coordinate
486:  * system for the tensor product element */
487: static PetscErrorCode PetscLagNodeIndicesCreateTensorVertices(DM dm, PetscLagNodeIndices facetni, PetscLagNodeIndices *nodeIndices)
488: {
489:   PetscLagNodeIndices ni;
490:   PetscInt            nodeIdxDim, subNodeIdxDim = facetni->nodeIdxDim;
491:   PetscInt            nVerts, nSubVerts         = facetni->nNodes;
492:   PetscInt            dim, d, e, f, g;

494:   PetscNew(&ni);
495:   DMGetDimension(dm, &dim);
496:   ni->nodeIdxDim = nodeIdxDim = subNodeIdxDim + 2;
497:   ni->nodeVecDim              = 0;
498:   ni->nNodes = nVerts = 2 * nSubVerts;
499:   ni->refct           = 1;
500:   PetscCalloc1(nodeIdxDim * nVerts, &(ni->nodeIdx));
501:   for (f = 0, d = 0; d < 2; d++) {
502:     for (e = 0; e < nSubVerts; e++, f++) {
503:       for (g = 0; g < subNodeIdxDim; g++) ni->nodeIdx[f * nodeIdxDim + g] = facetni->nodeIdx[e * subNodeIdxDim + g];
504:       ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim]     = (1 - d);
505:       ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim + 1] = d;
506:     }
507:   }
508:   PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_TRUE);
509:   *nodeIndices = ni;
510:   return 0;
511: }

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

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

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

636: /* given to sets of nodes, take the tensor product, where the product of the dof indices is concatenation and the
637:  * product of the dof vectors is the wedge product */
638: static PetscErrorCode PetscLagNodeIndicesTensor(PetscLagNodeIndices tracei, PetscInt dimT, PetscInt kT, PetscLagNodeIndices fiberi, PetscInt dimF, PetscInt kF, PetscLagNodeIndices *nodeIndices)
639: {
640:   PetscInt            dim = dimT + dimF;
641:   PetscInt            nodeIdxDim, nNodes;
642:   PetscInt            formDegree = kT + kF;
643:   PetscInt            Nk, NkT, NkF;
644:   PetscInt            MkT, MkF;
645:   PetscLagNodeIndices ni;
646:   PetscInt            i, j, l;
647:   PetscReal          *projF, *projT;
648:   PetscReal          *projFstar, *projTstar;
649:   PetscReal          *workF, *workF2, *workT, *workT2, *work, *work2;
650:   PetscReal          *wedgeMat;
651:   PetscReal           sign;

653:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
654:   PetscDTBinomialInt(dimT, PetscAbsInt(kT), &NkT);
655:   PetscDTBinomialInt(dimF, PetscAbsInt(kF), &NkF);
656:   PetscDTBinomialInt(dim, PetscAbsInt(kT), &MkT);
657:   PetscDTBinomialInt(dim, PetscAbsInt(kF), &MkF);
658:   PetscNew(&ni);
659:   ni->nodeIdxDim = nodeIdxDim = tracei->nodeIdxDim + fiberi->nodeIdxDim;
660:   ni->nodeVecDim              = Nk;
661:   ni->nNodes = nNodes = tracei->nNodes * fiberi->nNodes;
662:   ni->refct           = 1;
663:   PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
664:   /* first concatenate the indices */
665:   for (l = 0, j = 0; j < fiberi->nNodes; j++) {
666:     for (i = 0; i < tracei->nNodes; i++, l++) {
667:       PetscInt m, n = 0;

669:       for (m = 0; m < tracei->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = tracei->nodeIdx[i * tracei->nodeIdxDim + m];
670:       for (m = 0; m < fiberi->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = fiberi->nodeIdx[j * fiberi->nodeIdxDim + m];
671:     }
672:   }

674:   /* now wedge together the push-forward vectors */
675:   PetscMalloc1(nNodes * Nk, &(ni->nodeVec));
676:   PetscCalloc2(dimT * dim, &projT, dimF * dim, &projF);
677:   for (i = 0; i < dimT; i++) projT[i * (dim + 1)] = 1.;
678:   for (i = 0; i < dimF; i++) projF[i * (dim + dimT + 1) + dimT] = 1.;
679:   PetscMalloc2(MkT * NkT, &projTstar, MkF * NkF, &projFstar);
680:   PetscDTAltVPullbackMatrix(dim, dimT, projT, kT, projTstar);
681:   PetscDTAltVPullbackMatrix(dim, dimF, projF, kF, projFstar);
682:   PetscMalloc6(MkT, &workT, MkT, &workT2, MkF, &workF, MkF, &workF2, Nk, &work, Nk, &work2);
683:   PetscMalloc1(Nk * MkT, &wedgeMat);
684:   sign = (PetscAbsInt(kT * kF) & 1) ? -1. : 1.;
685:   for (l = 0, j = 0; j < fiberi->nNodes; j++) {
686:     PetscInt d, e;

688:     /* push forward fiber k-form */
689:     for (d = 0; d < MkF; d++) {
690:       PetscReal val = 0.;
691:       for (e = 0; e < NkF; e++) val += projFstar[d * NkF + e] * fiberi->nodeVec[j * NkF + e];
692:       workF[d] = val;
693:     }
694:     /* Hodge star to proper form if necessary */
695:     if (kF < 0) {
696:       for (d = 0; d < MkF; d++) workF2[d] = workF[d];
697:       PetscDTAltVStar(dim, PetscAbsInt(kF), 1, workF2, workF);
698:     }
699:     /* Compute the matrix that wedges this form with one of the trace k-form */
700:     PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kF), PetscAbsInt(kT), workF, wedgeMat);
701:     for (i = 0; i < tracei->nNodes; i++, l++) {
702:       /* push forward trace k-form */
703:       for (d = 0; d < MkT; d++) {
704:         PetscReal val = 0.;
705:         for (e = 0; e < NkT; e++) val += projTstar[d * NkT + e] * tracei->nodeVec[i * NkT + e];
706:         workT[d] = val;
707:       }
708:       /* Hodge star to proper form if necessary */
709:       if (kT < 0) {
710:         for (d = 0; d < MkT; d++) workT2[d] = workT[d];
711:         PetscDTAltVStar(dim, PetscAbsInt(kT), 1, workT2, workT);
712:       }
713:       /* compute the wedge product of the push-forward trace form and firer forms */
714:       for (d = 0; d < Nk; d++) {
715:         PetscReal val = 0.;
716:         for (e = 0; e < MkT; e++) val += wedgeMat[d * MkT + e] * workT[e];
717:         work[d] = val;
718:       }
719:       /* inverse Hodge star from proper form if necessary */
720:       if (formDegree < 0) {
721:         for (d = 0; d < Nk; d++) work2[d] = work[d];
722:         PetscDTAltVStar(dim, PetscAbsInt(formDegree), -1, work2, work);
723:       }
724:       /* insert into the array (adjusting for sign) */
725:       for (d = 0; d < Nk; d++) ni->nodeVec[l * Nk + d] = sign * work[d];
726:     }
727:   }
728:   PetscFree(wedgeMat);
729:   PetscFree6(workT, workT2, workF, workF2, work, work2);
730:   PetscFree2(projTstar, projFstar);
731:   PetscFree2(projT, projF);
732:   *nodeIndices = ni;
733:   return 0;
734: }

736: /* simple union of two sets of nodes */
737: static PetscErrorCode PetscLagNodeIndicesMerge(PetscLagNodeIndices niA, PetscLagNodeIndices niB, PetscLagNodeIndices *nodeIndices)
738: {
739:   PetscLagNodeIndices ni;
740:   PetscInt            nodeIdxDim, nodeVecDim, nNodes;

742:   PetscNew(&ni);
743:   ni->nodeIdxDim = nodeIdxDim = niA->nodeIdxDim;
745:   ni->nodeVecDim = nodeVecDim = niA->nodeVecDim;
747:   ni->nNodes = nNodes = niA->nNodes + niB->nNodes;
748:   ni->refct           = 1;
749:   PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
750:   PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
751:   PetscArraycpy(ni->nodeIdx, niA->nodeIdx, niA->nNodes * nodeIdxDim);
752:   PetscArraycpy(ni->nodeVec, niA->nodeVec, niA->nNodes * nodeVecDim);
753:   PetscArraycpy(&(ni->nodeIdx[niA->nNodes * nodeIdxDim]), niB->nodeIdx, niB->nNodes * nodeIdxDim);
754:   PetscArraycpy(&(ni->nodeVec[niA->nNodes * nodeVecDim]), niB->nodeVec, niB->nNodes * nodeVecDim);
755:   *nodeIndices = ni;
756:   return 0;
757: }

759: #define PETSCTUPINTCOMPREVLEX(N) \
760:   static int PetscConcat_(PetscTupIntCompRevlex_, N)(const void *a, const void *b) \
761:   { \
762:     const PetscInt *A = (const PetscInt *)a; \
763:     const PetscInt *B = (const PetscInt *)b; \
764:     int             i; \
765:     PetscInt        diff = 0; \
766:     for (i = 0; i < N; i++) { \
767:       diff = A[N - i] - B[N - i]; \
768:       if (diff) break; \
769:     } \
770:     return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1; \
771:   }

773: PETSCTUPINTCOMPREVLEX(3)
774: PETSCTUPINTCOMPREVLEX(4)
775: PETSCTUPINTCOMPREVLEX(5)
776: PETSCTUPINTCOMPREVLEX(6)
777: PETSCTUPINTCOMPREVLEX(7)

779: static int PetscTupIntCompRevlex_N(const void *a, const void *b)
780: {
781:   const PetscInt *A = (const PetscInt *)a;
782:   const PetscInt *B = (const PetscInt *)b;
783:   int             i;
784:   int             N    = A[0];
785:   PetscInt        diff = 0;
786:   for (i = 0; i < N; i++) {
787:     diff = A[N - i] - B[N - i];
788:     if (diff) break;
789:   }
790:   return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1;
791: }

793: /* The nodes are not necessarily in revlex order wrt nodeIdx: get the permutation
794:  * that puts them in that order */
795: static PetscErrorCode PetscLagNodeIndicesGetPermutation(PetscLagNodeIndices ni, PetscInt *perm[])
796: {
797:   if (!(ni->perm)) {
798:     PetscInt *sorter;
799:     PetscInt  m          = ni->nNodes;
800:     PetscInt  nodeIdxDim = ni->nodeIdxDim;
801:     PetscInt  i, j, k, l;
802:     PetscInt *prm;
803:     int (*comp)(const void *, const void *);

805:     PetscMalloc1((nodeIdxDim + 2) * m, &sorter);
806:     for (k = 0, l = 0, i = 0; i < m; i++) {
807:       sorter[k++] = nodeIdxDim + 1;
808:       sorter[k++] = i;
809:       for (j = 0; j < nodeIdxDim; j++) sorter[k++] = ni->nodeIdx[l++];
810:     }
811:     switch (nodeIdxDim) {
812:     case 2:
813:       comp = PetscTupIntCompRevlex_3;
814:       break;
815:     case 3:
816:       comp = PetscTupIntCompRevlex_4;
817:       break;
818:     case 4:
819:       comp = PetscTupIntCompRevlex_5;
820:       break;
821:     case 5:
822:       comp = PetscTupIntCompRevlex_6;
823:       break;
824:     case 6:
825:       comp = PetscTupIntCompRevlex_7;
826:       break;
827:     default:
828:       comp = PetscTupIntCompRevlex_N;
829:       break;
830:     }
831:     qsort(sorter, m, (nodeIdxDim + 2) * sizeof(PetscInt), comp);
832:     PetscMalloc1(m, &prm);
833:     for (i = 0; i < m; i++) prm[i] = sorter[(nodeIdxDim + 2) * i + 1];
834:     ni->perm = prm;
835:     PetscFree(sorter);
836:   }
837:   *perm = ni->perm;
838:   return 0;
839: }

841: static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
842: {
843:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

845:   if (lag->symperms) {
846:     PetscInt **selfSyms = lag->symperms[0];

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

851:       for (i = 0; i < lag->numSelfSym; i++) PetscFree(allocated[i]);
852:       PetscFree(allocated);
853:     }
854:     PetscFree(lag->symperms);
855:   }
856:   if (lag->symflips) {
857:     PetscScalar **selfSyms = lag->symflips[0];

859:     if (selfSyms) {
860:       PetscInt      i;
861:       PetscScalar **allocated = &selfSyms[-lag->selfSymOff];

863:       for (i = 0; i < lag->numSelfSym; i++) PetscFree(allocated[i]);
864:       PetscFree(allocated);
865:     }
866:     PetscFree(lag->symflips);
867:   }
868:   Petsc1DNodeFamilyDestroy(&(lag->nodeFamily));
869:   PetscLagNodeIndicesDestroy(&(lag->vertIndices));
870:   PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));
871:   PetscLagNodeIndicesDestroy(&(lag->allNodeIndices));
872:   PetscFree(lag);
873:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
874:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
875:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);
876:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);
877:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL);
878:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL);
879:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL);
880:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL);
881:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", NULL);
882:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", NULL);
883:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", NULL);
884:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", NULL);
885:   return 0;
886: }

888: static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
889: {
890:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

892:   PetscViewerASCIIPrintf(viewer, "%s %s%sLagrange dual space\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "", lag->trimmed ? "trimmed " : "");
893:   return 0;
894: }

896: static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
897: {
898:   PetscBool iascii;

902:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
903:   if (iascii) PetscDualSpaceLagrangeView_Ascii(sp, viewer);
904:   return 0;
905: }

907: static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscDualSpace sp, PetscOptionItems *PetscOptionsObject)
908: {
909:   PetscBool       continuous, tensor, trimmed, flg, flg2, flg3;
910:   PetscDTNodeType nodeType;
911:   PetscReal       nodeExponent;
912:   PetscInt        momentOrder;
913:   PetscBool       nodeEndpoints, useMoments;

915:   PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
916:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
917:   PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
918:   PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &nodeEndpoints, &nodeExponent);
919:   if (nodeType == PETSCDTNODES_DEFAULT) nodeType = PETSCDTNODES_GAUSSJACOBI;
920:   PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments);
921:   PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder);
922:   PetscOptionsHeadBegin(PetscOptionsObject, "PetscDualSpace Lagrange Options");
923:   PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
924:   if (flg) PetscDualSpaceLagrangeSetContinuity(sp, continuous);
925:   PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetTensor", tensor, &tensor, &flg);
926:   if (flg) PetscDualSpaceLagrangeSetTensor(sp, tensor);
927:   PetscOptionsBool("-petscdualspace_lagrange_trimmed", "Flag for trimmed dual space", "PetscDualSpaceLagrangeSetTrimmed", trimmed, &trimmed, &flg);
928:   if (flg) PetscDualSpaceLagrangeSetTrimmed(sp, trimmed);
929:   PetscOptionsEnum("-petscdualspace_lagrange_node_type", "Lagrange node location type", "PetscDualSpaceLagrangeSetNodeType", PetscDTNodeTypes, (PetscEnum)nodeType, (PetscEnum *)&nodeType, &flg);
930:   PetscOptionsBool("-petscdualspace_lagrange_node_endpoints", "Flag for nodes that include endpoints", "PetscDualSpaceLagrangeSetNodeType", nodeEndpoints, &nodeEndpoints, &flg2);
931:   flg3 = PETSC_FALSE;
932:   if (nodeType == PETSCDTNODES_GAUSSJACOBI) PetscOptionsReal("-petscdualspace_lagrange_node_exponent", "Gauss-Jacobi weight function exponent", "PetscDualSpaceLagrangeSetNodeType", nodeExponent, &nodeExponent, &flg3);
933:   if (flg || flg2 || flg3) PetscDualSpaceLagrangeSetNodeType(sp, nodeType, nodeEndpoints, nodeExponent);
934:   PetscOptionsBool("-petscdualspace_lagrange_use_moments", "Use moments (where appropriate) for functionals", "PetscDualSpaceLagrangeSetUseMoments", useMoments, &useMoments, &flg);
935:   if (flg) PetscDualSpaceLagrangeSetUseMoments(sp, useMoments);
936:   PetscOptionsInt("-petscdualspace_lagrange_moment_order", "Quadrature order for moment functionals", "PetscDualSpaceLagrangeSetMomentOrder", momentOrder, &momentOrder, &flg);
937:   if (flg) PetscDualSpaceLagrangeSetMomentOrder(sp, momentOrder);
938:   PetscOptionsHeadEnd();
939:   return 0;
940: }

942: static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace spNew)
943: {
944:   PetscBool           cont, tensor, trimmed, boundary;
945:   PetscDTNodeType     nodeType;
946:   PetscReal           exponent;
947:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

949:   PetscDualSpaceLagrangeGetContinuity(sp, &cont);
950:   PetscDualSpaceLagrangeSetContinuity(spNew, cont);
951:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
952:   PetscDualSpaceLagrangeSetTensor(spNew, tensor);
953:   PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
954:   PetscDualSpaceLagrangeSetTrimmed(spNew, trimmed);
955:   PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &boundary, &exponent);
956:   PetscDualSpaceLagrangeSetNodeType(spNew, nodeType, boundary, exponent);
957:   if (lag->nodeFamily) {
958:     PetscDualSpace_Lag *lagnew = (PetscDualSpace_Lag *)spNew->data;

960:     Petsc1DNodeFamilyReference(lag->nodeFamily);
961:     lagnew->nodeFamily = lag->nodeFamily;
962:   }
963:   return 0;
964: }

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

973:   PetscDualSpaceDuplicate(sp, bdsp);
974:   PetscDualSpaceSetFormDegree(*bdsp, k);
975:   DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(1, PETSC_TRUE), &K);
976:   PetscDualSpaceSetDM(*bdsp, K);
977:   DMDestroy(&K);
978:   PetscDualSpaceSetOrder(*bdsp, order);
979:   PetscDualSpaceSetNumComponents(*bdsp, Nc);
980:   newlag               = (PetscDualSpace_Lag *)(*bdsp)->data;
981:   newlag->interiorOnly = interiorOnly;
982:   PetscDualSpaceSetUp(*bdsp);
983:   return 0;
984: }

986: /* just the points, weights aren't handled */
987: static PetscErrorCode PetscQuadratureCreateTensor(PetscQuadrature trace, PetscQuadrature fiber, PetscQuadrature *product)
988: {
989:   PetscInt         dimTrace, dimFiber;
990:   PetscInt         numPointsTrace, numPointsFiber;
991:   PetscInt         dim, numPoints;
992:   const PetscReal *pointsTrace;
993:   const PetscReal *pointsFiber;
994:   PetscReal       *points;
995:   PetscInt         i, j, k, p;

997:   PetscQuadratureGetData(trace, &dimTrace, NULL, &numPointsTrace, &pointsTrace, NULL);
998:   PetscQuadratureGetData(fiber, &dimFiber, NULL, &numPointsFiber, &pointsFiber, NULL);
999:   dim       = dimTrace + dimFiber;
1000:   numPoints = numPointsFiber * numPointsTrace;
1001:   PetscMalloc1(numPoints * dim, &points);
1002:   for (p = 0, j = 0; j < numPointsFiber; j++) {
1003:     for (i = 0; i < numPointsTrace; i++, p++) {
1004:       for (k = 0; k < dimTrace; k++) points[p * dim + k] = pointsTrace[i * dimTrace + k];
1005:       for (k = 0; k < dimFiber; k++) points[p * dim + dimTrace + k] = pointsFiber[j * dimFiber + k];
1006:     }
1007:   }
1008:   PetscQuadratureCreate(PETSC_COMM_SELF, product);
1009:   PetscQuadratureSetData(*product, dim, 0, numPoints, points, NULL);
1010:   return 0;
1011: }

1013: /* Kronecker tensor product where matrix is considered a matrix of k-forms, so that
1014:  * the entries in the product matrix are wedge products of the entries in the original matrices */
1015: static PetscErrorCode MatTensorAltV(Mat trace, Mat fiber, PetscInt dimTrace, PetscInt kTrace, PetscInt dimFiber, PetscInt kFiber, Mat *product)
1016: {
1017:   PetscInt     mTrace, nTrace, mFiber, nFiber, m, n, k, i, j, l;
1018:   PetscInt     dim, NkTrace, NkFiber, Nk;
1019:   PetscInt     dT, dF;
1020:   PetscInt    *nnzTrace, *nnzFiber, *nnz;
1021:   PetscInt     iT, iF, jT, jF, il, jl;
1022:   PetscReal   *workT, *workT2, *workF, *workF2, *work, *workstar;
1023:   PetscReal   *projT, *projF;
1024:   PetscReal   *projTstar, *projFstar;
1025:   PetscReal   *wedgeMat;
1026:   PetscReal    sign;
1027:   PetscScalar *workS;
1028:   Mat          prod;
1029:   /* this produces dof groups that look like the identity */

1031:   MatGetSize(trace, &mTrace, &nTrace);
1032:   PetscDTBinomialInt(dimTrace, PetscAbsInt(kTrace), &NkTrace);
1034:   MatGetSize(fiber, &mFiber, &nFiber);
1035:   PetscDTBinomialInt(dimFiber, PetscAbsInt(kFiber), &NkFiber);
1037:   PetscMalloc2(mTrace, &nnzTrace, mFiber, &nnzFiber);
1038:   for (i = 0; i < mTrace; i++) {
1039:     MatGetRow(trace, i, &(nnzTrace[i]), NULL, NULL);
1041:   }
1042:   for (i = 0; i < mFiber; i++) {
1043:     MatGetRow(fiber, i, &(nnzFiber[i]), NULL, NULL);
1045:   }
1046:   dim = dimTrace + dimFiber;
1047:   k   = kFiber + kTrace;
1048:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1049:   m = mTrace * mFiber;
1050:   PetscMalloc1(m, &nnz);
1051:   for (l = 0, j = 0; j < mFiber; j++)
1052:     for (i = 0; i < mTrace; i++, l++) nnz[l] = (nnzTrace[i] / NkTrace) * (nnzFiber[j] / NkFiber) * Nk;
1053:   n = (nTrace / NkTrace) * (nFiber / NkFiber) * Nk;
1054:   MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &prod);
1055:   PetscFree(nnz);
1056:   PetscFree2(nnzTrace, nnzFiber);
1057:   /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1058:   MatSetOption(prod, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1059:   /* compute pullbacks */
1060:   PetscDTBinomialInt(dim, PetscAbsInt(kTrace), &dT);
1061:   PetscDTBinomialInt(dim, PetscAbsInt(kFiber), &dF);
1062:   PetscMalloc4(dimTrace * dim, &projT, dimFiber * dim, &projF, dT * NkTrace, &projTstar, dF * NkFiber, &projFstar);
1063:   PetscArrayzero(projT, dimTrace * dim);
1064:   for (i = 0; i < dimTrace; i++) projT[i * (dim + 1)] = 1.;
1065:   PetscArrayzero(projF, dimFiber * dim);
1066:   for (i = 0; i < dimFiber; i++) projF[i * (dim + 1) + dimTrace] = 1.;
1067:   PetscDTAltVPullbackMatrix(dim, dimTrace, projT, kTrace, projTstar);
1068:   PetscDTAltVPullbackMatrix(dim, dimFiber, projF, kFiber, projFstar);
1069:   PetscMalloc5(dT, &workT, dF, &workF, Nk, &work, Nk, &workstar, Nk, &workS);
1070:   PetscMalloc2(dT, &workT2, dF, &workF2);
1071:   PetscMalloc1(Nk * dT, &wedgeMat);
1072:   sign = (PetscAbsInt(kTrace * kFiber) & 1) ? -1. : 1.;
1073:   for (i = 0, iF = 0; iF < mFiber; iF++) {
1074:     PetscInt           ncolsF, nformsF;
1075:     const PetscInt    *colsF;
1076:     const PetscScalar *valsF;

1078:     MatGetRow(fiber, iF, &ncolsF, &colsF, &valsF);
1079:     nformsF = ncolsF / NkFiber;
1080:     for (iT = 0; iT < mTrace; iT++, i++) {
1081:       PetscInt           ncolsT, nformsT;
1082:       const PetscInt    *colsT;
1083:       const PetscScalar *valsT;

1085:       MatGetRow(trace, iT, &ncolsT, &colsT, &valsT);
1086:       nformsT = ncolsT / NkTrace;
1087:       for (j = 0, jF = 0; jF < nformsF; jF++) {
1088:         PetscInt colF = colsF[jF * NkFiber] / NkFiber;

1090:         for (il = 0; il < dF; il++) {
1091:           PetscReal val = 0.;
1092:           for (jl = 0; jl < NkFiber; jl++) val += projFstar[il * NkFiber + jl] * PetscRealPart(valsF[jF * NkFiber + jl]);
1093:           workF[il] = val;
1094:         }
1095:         if (kFiber < 0) {
1096:           for (il = 0; il < dF; il++) workF2[il] = workF[il];
1097:           PetscDTAltVStar(dim, PetscAbsInt(kFiber), 1, workF2, workF);
1098:         }
1099:         PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kFiber), PetscAbsInt(kTrace), workF, wedgeMat);
1100:         for (jT = 0; jT < nformsT; jT++, j++) {
1101:           PetscInt           colT = colsT[jT * NkTrace] / NkTrace;
1102:           PetscInt           col  = colF * (nTrace / NkTrace) + colT;
1103:           const PetscScalar *vals;

1105:           for (il = 0; il < dT; il++) {
1106:             PetscReal val = 0.;
1107:             for (jl = 0; jl < NkTrace; jl++) val += projTstar[il * NkTrace + jl] * PetscRealPart(valsT[jT * NkTrace + jl]);
1108:             workT[il] = val;
1109:           }
1110:           if (kTrace < 0) {
1111:             for (il = 0; il < dT; il++) workT2[il] = workT[il];
1112:             PetscDTAltVStar(dim, PetscAbsInt(kTrace), 1, workT2, workT);
1113:           }

1115:           for (il = 0; il < Nk; il++) {
1116:             PetscReal val = 0.;
1117:             for (jl = 0; jl < dT; jl++) val += sign * wedgeMat[il * dT + jl] * workT[jl];
1118:             work[il] = val;
1119:           }
1120:           if (k < 0) {
1121:             PetscDTAltVStar(dim, PetscAbsInt(k), -1, work, workstar);
1122: #if defined(PETSC_USE_COMPLEX)
1123:             for (l = 0; l < Nk; l++) workS[l] = workstar[l];
1124:             vals = &workS[0];
1125: #else
1126:             vals = &workstar[0];
1127: #endif
1128:           } else {
1129: #if defined(PETSC_USE_COMPLEX)
1130:             for (l = 0; l < Nk; l++) workS[l] = work[l];
1131:             vals = &workS[0];
1132: #else
1133:             vals = &work[0];
1134: #endif
1135:           }
1136:           for (l = 0; l < Nk; l++) MatSetValue(prod, i, col * Nk + l, vals[l], INSERT_VALUES); /* Nk */
1137:         }                                                                                                 /* jT */
1138:       }                                                                                                   /* jF */
1139:       MatRestoreRow(trace, iT, &ncolsT, &colsT, &valsT);
1140:     } /* iT */
1141:     MatRestoreRow(fiber, iF, &ncolsF, &colsF, &valsF);
1142:   } /* iF */
1143:   PetscFree(wedgeMat);
1144:   PetscFree4(projT, projF, projTstar, projFstar);
1145:   PetscFree2(workT2, workF2);
1146:   PetscFree5(workT, workF, work, workstar, workS);
1147:   MatAssemblyBegin(prod, MAT_FINAL_ASSEMBLY);
1148:   MatAssemblyEnd(prod, MAT_FINAL_ASSEMBLY);
1149:   *product = prod;
1150:   return 0;
1151: }

1153: /* Union of quadrature points, with an attempt to identify commont points in the two sets */
1154: static PetscErrorCode PetscQuadraturePointsMerge(PetscQuadrature quadA, PetscQuadrature quadB, PetscQuadrature *quadJoint, PetscInt *aToJoint[], PetscInt *bToJoint[])
1155: {
1156:   PetscInt         dimA, dimB;
1157:   PetscInt         nA, nB, nJoint, i, j, d;
1158:   const PetscReal *pointsA;
1159:   const PetscReal *pointsB;
1160:   PetscReal       *pointsJoint;
1161:   PetscInt        *aToJ, *bToJ;
1162:   PetscQuadrature  qJ;

1164:   PetscQuadratureGetData(quadA, &dimA, NULL, &nA, &pointsA, NULL);
1165:   PetscQuadratureGetData(quadB, &dimB, NULL, &nB, &pointsB, NULL);
1167:   nJoint = nA;
1168:   PetscMalloc1(nA, &aToJ);
1169:   for (i = 0; i < nA; i++) aToJ[i] = i;
1170:   PetscMalloc1(nB, &bToJ);
1171:   for (i = 0; i < nB; i++) {
1172:     for (j = 0; j < nA; j++) {
1173:       bToJ[i] = -1;
1174:       for (d = 0; d < dimA; d++)
1175:         if (PetscAbsReal(pointsB[i * dimA + d] - pointsA[j * dimA + d]) > PETSC_SMALL) break;
1176:       if (d == dimA) {
1177:         bToJ[i] = j;
1178:         break;
1179:       }
1180:     }
1181:     if (bToJ[i] == -1) bToJ[i] = nJoint++;
1182:   }
1183:   *aToJoint = aToJ;
1184:   *bToJoint = bToJ;
1185:   PetscMalloc1(nJoint * dimA, &pointsJoint);
1186:   PetscArraycpy(pointsJoint, pointsA, nA * dimA);
1187:   for (i = 0; i < nB; i++) {
1188:     if (bToJ[i] >= nA) {
1189:       for (d = 0; d < dimA; d++) pointsJoint[bToJ[i] * dimA + d] = pointsB[i * dimA + d];
1190:     }
1191:   }
1192:   PetscQuadratureCreate(PETSC_COMM_SELF, &qJ);
1193:   PetscQuadratureSetData(qJ, dimA, 0, nJoint, pointsJoint, NULL);
1194:   *quadJoint = qJ;
1195:   return 0;
1196: }

1198: /* Matrices matA and matB are both quadrature -> dof matrices: produce a matrix that is joint quadrature -> union of
1199:  * dofs, where the joint quadrature was produced by PetscQuadraturePointsMerge */
1200: static PetscErrorCode MatricesMerge(Mat matA, Mat matB, PetscInt dim, PetscInt k, PetscInt numMerged, const PetscInt aToMerged[], const PetscInt bToMerged[], Mat *matMerged)
1201: {
1202:   PetscInt  m, n, mA, nA, mB, nB, Nk, i, j, l;
1203:   Mat       M;
1204:   PetscInt *nnz;
1205:   PetscInt  maxnnz;
1206:   PetscInt *work;

1208:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1209:   MatGetSize(matA, &mA, &nA);
1211:   MatGetSize(matB, &mB, &nB);
1213:   m = mA + mB;
1214:   n = numMerged * Nk;
1215:   PetscMalloc1(m, &nnz);
1216:   maxnnz = 0;
1217:   for (i = 0; i < mA; i++) {
1218:     MatGetRow(matA, i, &(nnz[i]), NULL, NULL);
1220:     maxnnz = PetscMax(maxnnz, nnz[i]);
1221:   }
1222:   for (i = 0; i < mB; i++) {
1223:     MatGetRow(matB, i, &(nnz[i + mA]), NULL, NULL);
1225:     maxnnz = PetscMax(maxnnz, nnz[i + mA]);
1226:   }
1227:   MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &M);
1228:   PetscFree(nnz);
1229:   /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1230:   MatSetOption(M, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1231:   PetscMalloc1(maxnnz, &work);
1232:   for (i = 0; i < mA; i++) {
1233:     const PetscInt    *cols;
1234:     const PetscScalar *vals;
1235:     PetscInt           nCols;
1236:     MatGetRow(matA, i, &nCols, &cols, &vals);
1237:     for (j = 0; j < nCols / Nk; j++) {
1238:       PetscInt newCol = aToMerged[cols[j * Nk] / Nk];
1239:       for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1240:     }
1241:     MatSetValuesBlocked(M, 1, &i, nCols, work, vals, INSERT_VALUES);
1242:     MatRestoreRow(matA, i, &nCols, &cols, &vals);
1243:   }
1244:   for (i = 0; i < mB; i++) {
1245:     const PetscInt    *cols;
1246:     const PetscScalar *vals;

1248:     PetscInt row = i + mA;
1249:     PetscInt nCols;
1250:     MatGetRow(matB, i, &nCols, &cols, &vals);
1251:     for (j = 0; j < nCols / Nk; j++) {
1252:       PetscInt newCol = bToMerged[cols[j * Nk] / Nk];
1253:       for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1254:     }
1255:     MatSetValuesBlocked(M, 1, &row, nCols, work, vals, INSERT_VALUES);
1256:     MatRestoreRow(matB, i, &nCols, &cols, &vals);
1257:   }
1258:   PetscFree(work);
1259:   MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
1260:   MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
1261:   *matMerged = M;
1262:   return 0;
1263: }

1265: /* Take a dual space and product a segment space that has all the same specifications (trimmed, continuous, order,
1266:  * node set), except for the form degree.  For computing boundary dofs and for making tensor product spaces */
1267: static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Lagrange(PetscDualSpace sp, DM K, PetscInt f, PetscInt k, PetscInt Ncopies, PetscBool interiorOnly, PetscDualSpace *bdsp)
1268: {
1269:   PetscInt            Nknew, Ncnew;
1270:   PetscInt            dim, pointDim = -1;
1271:   PetscInt            depth;
1272:   DM                  dm;
1273:   PetscDualSpace_Lag *newlag;

1275:   PetscDualSpaceGetDM(sp, &dm);
1276:   DMGetDimension(dm, &dim);
1277:   DMPlexGetDepth(dm, &depth);
1278:   PetscDualSpaceDuplicate(sp, bdsp);
1279:   PetscDualSpaceSetFormDegree(*bdsp, k);
1280:   if (!K) {
1281:     if (depth == dim) {
1282:       DMPolytopeType ct;

1284:       pointDim = dim - 1;
1285:       DMPlexGetCellType(dm, f, &ct);
1286:       DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K);
1287:     } else if (depth == 1) {
1288:       pointDim = 0;
1289:       DMPlexCreateReferenceCell(PETSC_COMM_SELF, DM_POLYTOPE_POINT, &K);
1290:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
1291:   } else {
1292:     PetscObjectReference((PetscObject)K);
1293:     DMGetDimension(K, &pointDim);
1294:   }
1295:   PetscDualSpaceSetDM(*bdsp, K);
1296:   DMDestroy(&K);
1297:   PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew);
1298:   Ncnew = Nknew * Ncopies;
1299:   PetscDualSpaceSetNumComponents(*bdsp, Ncnew);
1300:   newlag               = (PetscDualSpace_Lag *)(*bdsp)->data;
1301:   newlag->interiorOnly = interiorOnly;
1302:   PetscDualSpaceSetUp(*bdsp);
1303:   return 0;
1304: }

1306: /* Construct simplex nodes from a nodefamily, add Nk dof vectors of length Nk at each node.
1307:  * Return the (quadrature, matrix) form of the dofs and the nodeIndices form as well.
1308:  *
1309:  * Sometimes we want a set of nodes to be contained in the interior of the element,
1310:  * even when the node scheme puts nodes on the boundaries.  numNodeSkip tells
1311:  * the routine how many "layers" of nodes need to be skipped.
1312:  * */
1313: static PetscErrorCode PetscDualSpaceLagrangeCreateSimplexNodeMat(Petsc1DNodeFamily nodeFamily, PetscInt dim, PetscInt sum, PetscInt Nk, PetscInt numNodeSkip, PetscQuadrature *iNodes, Mat *iMat, PetscLagNodeIndices *nodeIndices)
1314: {
1315:   PetscReal          *extraNodeCoords, *nodeCoords;
1316:   PetscInt            nNodes, nExtraNodes;
1317:   PetscInt            i, j, k, extraSum = sum + numNodeSkip * (1 + dim);
1318:   PetscQuadrature     intNodes;
1319:   Mat                 intMat;
1320:   PetscLagNodeIndices ni;

1322:   PetscDTBinomialInt(dim + sum, dim, &nNodes);
1323:   PetscDTBinomialInt(dim + extraSum, dim, &nExtraNodes);

1325:   PetscMalloc1(dim * nExtraNodes, &extraNodeCoords);
1326:   PetscNew(&ni);
1327:   ni->nodeIdxDim = dim + 1;
1328:   ni->nodeVecDim = Nk;
1329:   ni->nNodes     = nNodes * Nk;
1330:   ni->refct      = 1;
1331:   PetscMalloc1(nNodes * Nk * (dim + 1), &(ni->nodeIdx));
1332:   PetscMalloc1(nNodes * Nk * Nk, &(ni->nodeVec));
1333:   for (i = 0; i < nNodes; i++)
1334:     for (j = 0; j < Nk; j++)
1335:       for (k = 0; k < Nk; k++) ni->nodeVec[(i * Nk + j) * Nk + k] = (j == k) ? 1. : 0.;
1336:   Petsc1DNodeFamilyComputeSimplexNodes(nodeFamily, dim, extraSum, extraNodeCoords);
1337:   if (numNodeSkip) {
1338:     PetscInt  k;
1339:     PetscInt *tup;

1341:     PetscMalloc1(dim * nNodes, &nodeCoords);
1342:     PetscMalloc1(dim + 1, &tup);
1343:     for (k = 0; k < nNodes; k++) {
1344:       PetscInt j, c;
1345:       PetscInt index;

1347:       PetscDTIndexToBary(dim + 1, sum, k, tup);
1348:       for (j = 0; j < dim + 1; j++) tup[j] += numNodeSkip;
1349:       for (c = 0; c < Nk; c++) {
1350:         for (j = 0; j < dim + 1; j++) ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1351:       }
1352:       PetscDTBaryToIndex(dim + 1, extraSum, tup, &index);
1353:       for (j = 0; j < dim; j++) nodeCoords[k * dim + j] = extraNodeCoords[index * dim + j];
1354:     }
1355:     PetscFree(tup);
1356:     PetscFree(extraNodeCoords);
1357:   } else {
1358:     PetscInt  k;
1359:     PetscInt *tup;

1361:     nodeCoords = extraNodeCoords;
1362:     PetscMalloc1(dim + 1, &tup);
1363:     for (k = 0; k < nNodes; k++) {
1364:       PetscInt j, c;

1366:       PetscDTIndexToBary(dim + 1, sum, k, tup);
1367:       for (c = 0; c < Nk; c++) {
1368:         for (j = 0; j < dim + 1; j++) {
1369:           /* barycentric indices can have zeros, but we don't want to push forward zeros because it makes it harder to
1370:            * determine which nodes correspond to which under symmetries, so we increase by 1.  This is fine
1371:            * because the nodeIdx coordinates don't have any meaning other than helping to identify symmetries */
1372:           ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1373:         }
1374:       }
1375:     }
1376:     PetscFree(tup);
1377:   }
1378:   PetscQuadratureCreate(PETSC_COMM_SELF, &intNodes);
1379:   PetscQuadratureSetData(intNodes, dim, 0, nNodes, nodeCoords, NULL);
1380:   MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes * Nk, nNodes * Nk, Nk, NULL, &intMat);
1381:   MatSetOption(intMat, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1382:   for (j = 0; j < nNodes * Nk; j++) {
1383:     PetscInt rem = j % Nk;
1384:     PetscInt a, aprev = j - rem;
1385:     PetscInt anext = aprev + Nk;

1387:     for (a = aprev; a < anext; a++) MatSetValue(intMat, j, a, (a == j) ? 1. : 0., INSERT_VALUES);
1388:   }
1389:   MatAssemblyBegin(intMat, MAT_FINAL_ASSEMBLY);
1390:   MatAssemblyEnd(intMat, MAT_FINAL_ASSEMBLY);
1391:   *iNodes      = intNodes;
1392:   *iMat        = intMat;
1393:   *nodeIndices = ni;
1394:   return 0;
1395: }

1397: /* once the nodeIndices have been created for the interior of the reference cell, and for all of the boundary cells,
1398:  * push forward the boundary dofs and concatenate them into the full node indices for the dual space */
1399: static PetscErrorCode PetscDualSpaceLagrangeCreateAllNodeIdx(PetscDualSpace sp)
1400: {
1401:   DM                  dm;
1402:   PetscInt            dim, nDofs;
1403:   PetscSection        section;
1404:   PetscInt            pStart, pEnd, p;
1405:   PetscInt            formDegree, Nk;
1406:   PetscInt            nodeIdxDim, spintdim;
1407:   PetscDualSpace_Lag *lag;
1408:   PetscLagNodeIndices ni, verti;

1410:   lag   = (PetscDualSpace_Lag *)sp->data;
1411:   verti = lag->vertIndices;
1412:   PetscDualSpaceGetDM(sp, &dm);
1413:   DMGetDimension(dm, &dim);
1414:   PetscDualSpaceGetFormDegree(sp, &formDegree);
1415:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
1416:   PetscDualSpaceGetSection(sp, &section);
1417:   PetscSectionGetStorageSize(section, &nDofs);
1418:   PetscNew(&ni);
1419:   ni->nodeIdxDim = nodeIdxDim = verti->nodeIdxDim;
1420:   ni->nodeVecDim              = Nk;
1421:   ni->nNodes                  = nDofs;
1422:   ni->refct                   = 1;
1423:   PetscMalloc1(nodeIdxDim * nDofs, &(ni->nodeIdx));
1424:   PetscMalloc1(Nk * nDofs, &(ni->nodeVec));
1425:   DMPlexGetChart(dm, &pStart, &pEnd);
1426:   PetscSectionGetDof(section, 0, &spintdim);
1427:   if (spintdim) {
1428:     PetscArraycpy(ni->nodeIdx, lag->intNodeIndices->nodeIdx, spintdim * nodeIdxDim);
1429:     PetscArraycpy(ni->nodeVec, lag->intNodeIndices->nodeVec, spintdim * Nk);
1430:   }
1431:   for (p = pStart + 1; p < pEnd; p++) {
1432:     PetscDualSpace      psp = sp->pointSpaces[p];
1433:     PetscDualSpace_Lag *plag;
1434:     PetscInt            dof, off;

1436:     PetscSectionGetDof(section, p, &dof);
1437:     if (!dof) continue;
1438:     plag = (PetscDualSpace_Lag *)psp->data;
1439:     PetscSectionGetOffset(section, p, &off);
1440:     PetscLagNodeIndicesPushForward(dm, verti, p, plag->vertIndices, plag->intNodeIndices, 0, formDegree, &(ni->nodeIdx[off * nodeIdxDim]), &(ni->nodeVec[off * Nk]));
1441:   }
1442:   lag->allNodeIndices = ni;
1443:   return 0;
1444: }

1446: /* once the (quadrature, Matrix) forms of the dofs have been created for the interior of the
1447:  * reference cell and for the boundary cells, jk
1448:  * push forward the boundary data and concatenate them into the full (quadrature, matrix) data
1449:  * for the dual space */
1450: static PetscErrorCode PetscDualSpaceCreateAllDataFromInteriorData(PetscDualSpace sp)
1451: {
1452:   DM              dm;
1453:   PetscSection    section;
1454:   PetscInt        pStart, pEnd, p, k, Nk, dim, Nc;
1455:   PetscInt        nNodes;
1456:   PetscInt        countNodes;
1457:   Mat             allMat;
1458:   PetscQuadrature allNodes;
1459:   PetscInt        nDofs;
1460:   PetscInt        maxNzforms, j;
1461:   PetscScalar    *work;
1462:   PetscReal      *L, *J, *Jinv, *v0, *pv0;
1463:   PetscInt       *iwork;
1464:   PetscReal      *nodes;

1466:   PetscDualSpaceGetDM(sp, &dm);
1467:   DMGetDimension(dm, &dim);
1468:   PetscDualSpaceGetSection(sp, &section);
1469:   PetscSectionGetStorageSize(section, &nDofs);
1470:   DMPlexGetChart(dm, &pStart, &pEnd);
1471:   PetscDualSpaceGetFormDegree(sp, &k);
1472:   PetscDualSpaceGetNumComponents(sp, &Nc);
1473:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1474:   for (p = pStart, nNodes = 0, maxNzforms = 0; p < pEnd; p++) {
1475:     PetscDualSpace  psp;
1476:     DM              pdm;
1477:     PetscInt        pdim, pNk;
1478:     PetscQuadrature intNodes;
1479:     Mat             intMat;

1481:     PetscDualSpaceGetPointSubspace(sp, p, &psp);
1482:     if (!psp) continue;
1483:     PetscDualSpaceGetDM(psp, &pdm);
1484:     DMGetDimension(pdm, &pdim);
1485:     if (pdim < PetscAbsInt(k)) continue;
1486:     PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
1487:     PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
1488:     if (intNodes) {
1489:       PetscInt nNodesp;

1491:       PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, NULL, NULL);
1492:       nNodes += nNodesp;
1493:     }
1494:     if (intMat) {
1495:       PetscInt maxNzsp;
1496:       PetscInt maxNzformsp;

1498:       MatSeqAIJGetMaxRowNonzeros(intMat, &maxNzsp);
1500:       maxNzformsp = maxNzsp / pNk;
1501:       maxNzforms  = PetscMax(maxNzforms, maxNzformsp);
1502:     }
1503:   }
1504:   MatCreateSeqAIJ(PETSC_COMM_SELF, nDofs, nNodes * Nc, maxNzforms * Nk, NULL, &allMat);
1505:   MatSetOption(allMat, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1506:   PetscMalloc7(dim, &v0, dim, &pv0, dim * dim, &J, dim * dim, &Jinv, Nk * Nk, &L, maxNzforms * Nk, &work, maxNzforms * Nk, &iwork);
1507:   for (j = 0; j < dim; j++) pv0[j] = -1.;
1508:   PetscMalloc1(dim * nNodes, &nodes);
1509:   for (p = pStart, countNodes = 0; p < pEnd; p++) {
1510:     PetscDualSpace  psp;
1511:     PetscQuadrature intNodes;
1512:     DM              pdm;
1513:     PetscInt        pdim, pNk;
1514:     PetscInt        countNodesIn = countNodes;
1515:     PetscReal       detJ;
1516:     Mat             intMat;

1518:     PetscDualSpaceGetPointSubspace(sp, p, &psp);
1519:     if (!psp) continue;
1520:     PetscDualSpaceGetDM(psp, &pdm);
1521:     DMGetDimension(pdm, &pdim);
1522:     if (pdim < PetscAbsInt(k)) continue;
1523:     PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
1524:     if (intNodes == NULL && intMat == NULL) continue;
1525:     PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
1526:     if (p) {
1527:       DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, Jinv, &detJ);
1528:     } else { /* identity */
1529:       PetscInt i, j;

1531:       for (i = 0; i < dim; i++)
1532:         for (j = 0; j < dim; j++) J[i * dim + j] = Jinv[i * dim + j] = 0.;
1533:       for (i = 0; i < dim; i++) J[i * dim + i] = Jinv[i * dim + i] = 1.;
1534:       for (i = 0; i < dim; i++) v0[i] = -1.;
1535:     }
1536:     if (pdim != dim) { /* compactify Jacobian */
1537:       PetscInt i, j;

1539:       for (i = 0; i < dim; i++)
1540:         for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
1541:     }
1542:     PetscDTAltVPullbackMatrix(pdim, dim, J, k, L);
1543:     if (intNodes) { /* push forward quadrature locations by the affine transformation */
1544:       PetscInt         nNodesp;
1545:       const PetscReal *nodesp;
1546:       PetscInt         j;

1548:       PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, &nodesp, NULL);
1549:       for (j = 0; j < nNodesp; j++, countNodes++) {
1550:         PetscInt d, e;

1552:         for (d = 0; d < dim; d++) {
1553:           nodes[countNodes * dim + d] = v0[d];
1554:           for (e = 0; e < pdim; e++) nodes[countNodes * dim + d] += J[d * pdim + e] * (nodesp[j * pdim + e] - pv0[e]);
1555:         }
1556:       }
1557:     }
1558:     if (intMat) {
1559:       PetscInt nrows;
1560:       PetscInt off;

1562:       PetscSectionGetDof(section, p, &nrows);
1563:       PetscSectionGetOffset(section, p, &off);
1564:       for (j = 0; j < nrows; j++) {
1565:         PetscInt           ncols;
1566:         const PetscInt    *cols;
1567:         const PetscScalar *vals;
1568:         PetscInt           l, d, e;
1569:         PetscInt           row = j + off;

1571:         MatGetRow(intMat, j, &ncols, &cols, &vals);
1573:         for (l = 0; l < ncols / pNk; l++) {
1574:           PetscInt blockcol;

1577:           blockcol = cols[l * pNk] / pNk;
1578:           for (d = 0; d < Nk; d++) iwork[l * Nk + d] = (blockcol + countNodesIn) * Nk + d;
1579:           for (d = 0; d < Nk; d++) work[l * Nk + d] = 0.;
1580:           for (d = 0; d < Nk; d++) {
1581:             for (e = 0; e < pNk; e++) {
1582:               /* "push forward" dof by pulling back a k-form to be evaluated on the point: multiply on the right by L */
1583:               work[l * Nk + d] += vals[l * pNk + e] * L[e * Nk + d];
1584:             }
1585:           }
1586:         }
1587:         MatSetValues(allMat, 1, &row, (ncols / pNk) * Nk, iwork, work, INSERT_VALUES);
1588:         MatRestoreRow(intMat, j, &ncols, &cols, &vals);
1589:       }
1590:     }
1591:   }
1592:   MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);
1593:   MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);
1594:   PetscQuadratureCreate(PETSC_COMM_SELF, &allNodes);
1595:   PetscQuadratureSetData(allNodes, dim, 0, nNodes, nodes, NULL);
1596:   PetscFree7(v0, pv0, J, Jinv, L, work, iwork);
1597:   MatDestroy(&(sp->allMat));
1598:   sp->allMat = allMat;
1599:   PetscQuadratureDestroy(&(sp->allNodes));
1600:   sp->allNodes = allNodes;
1601:   return 0;
1602: }

1604: /* rather than trying to get all data from the functionals, we create
1605:  * the functionals from rows of the quadrature -> dof matrix.
1606:  *
1607:  * Ideally most of the uses of PetscDualSpace in PetscFE will switch
1608:  * to using intMat and allMat, so that the individual functionals
1609:  * don't need to be constructed at all */
1610: static PetscErrorCode PetscDualSpaceComputeFunctionalsFromAllData(PetscDualSpace sp)
1611: {
1612:   PetscQuadrature  allNodes;
1613:   Mat              allMat;
1614:   PetscInt         nDofs;
1615:   PetscInt         dim, k, Nk, Nc, f;
1616:   DM               dm;
1617:   PetscInt         nNodes, spdim;
1618:   const PetscReal *nodes = NULL;
1619:   PetscSection     section;
1620:   PetscBool        useMoments;

1622:   PetscDualSpaceGetDM(sp, &dm);
1623:   DMGetDimension(dm, &dim);
1624:   PetscDualSpaceGetNumComponents(sp, &Nc);
1625:   PetscDualSpaceGetFormDegree(sp, &k);
1626:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1627:   PetscDualSpaceGetAllData(sp, &allNodes, &allMat);
1628:   nNodes = 0;
1629:   if (allNodes) PetscQuadratureGetData(allNodes, NULL, NULL, &nNodes, &nodes, NULL);
1630:   MatGetSize(allMat, &nDofs, NULL);
1631:   PetscDualSpaceGetSection(sp, &section);
1632:   PetscSectionGetStorageSize(section, &spdim);
1634:   PetscMalloc1(nDofs, &(sp->functional));
1635:   PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments);
1636:   if (useMoments) {
1637:     Mat              allMat;
1638:     PetscInt         momentOrder, i;
1639:     PetscBool        tensor;
1640:     const PetscReal *weights;
1641:     PetscScalar     *array;

1644:     PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder);
1645:     PetscDualSpaceLagrangeGetTensor(sp, &tensor);
1646:     if (!tensor) PetscDTStroudConicalQuadrature(dim, Nc, PetscMax(momentOrder + 1, 1), -1.0, 1.0, &(sp->functional[0]));
1647:     else PetscDTGaussTensorQuadrature(dim, Nc, PetscMax(momentOrder + 1, 1), -1.0, 1.0, &(sp->functional[0]));
1648:     /* Need to replace allNodes and allMat */
1649:     PetscObjectReference((PetscObject)sp->functional[0]);
1650:     PetscQuadratureDestroy(&(sp->allNodes));
1651:     sp->allNodes = sp->functional[0];
1652:     PetscQuadratureGetData(sp->allNodes, NULL, NULL, &nNodes, NULL, &weights);
1653:     MatCreateSeqDense(PETSC_COMM_SELF, nDofs, nNodes * Nc, NULL, &allMat);
1654:     MatDenseGetArrayWrite(allMat, &array);
1655:     for (i = 0; i < nNodes * Nc; ++i) array[i] = weights[i];
1656:     MatDenseRestoreArrayWrite(allMat, &array);
1657:     MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);
1658:     MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);
1659:     MatDestroy(&(sp->allMat));
1660:     sp->allMat = allMat;
1661:     return 0;
1662:   }
1663:   for (f = 0; f < nDofs; f++) {
1664:     PetscInt           ncols, c;
1665:     const PetscInt    *cols;
1666:     const PetscScalar *vals;
1667:     PetscReal         *nodesf;
1668:     PetscReal         *weightsf;
1669:     PetscInt           nNodesf;
1670:     PetscInt           countNodes;

1672:     MatGetRow(allMat, f, &ncols, &cols, &vals);
1674:     for (c = 1, nNodesf = 1; c < ncols; c++) {
1675:       if ((cols[c] / Nc) != (cols[c - 1] / Nc)) nNodesf++;
1676:     }
1677:     PetscMalloc1(dim * nNodesf, &nodesf);
1678:     PetscMalloc1(Nc * nNodesf, &weightsf);
1679:     for (c = 0, countNodes = 0; c < ncols; c++) {
1680:       if (!c || ((cols[c] / Nc) != (cols[c - 1] / Nc))) {
1681:         PetscInt d;

1683:         for (d = 0; d < Nc; d++) weightsf[countNodes * Nc + d] = 0.;
1684:         for (d = 0; d < dim; d++) nodesf[countNodes * dim + d] = nodes[(cols[c] / Nc) * dim + d];
1685:         countNodes++;
1686:       }
1687:       weightsf[(countNodes - 1) * Nc + (cols[c] % Nc)] = PetscRealPart(vals[c]);
1688:     }
1689:     PetscQuadratureCreate(PETSC_COMM_SELF, &(sp->functional[f]));
1690:     PetscQuadratureSetData(sp->functional[f], dim, Nc, nNodesf, nodesf, weightsf);
1691:     MatRestoreRow(allMat, f, &ncols, &cols, &vals);
1692:   }
1693:   return 0;
1694: }

1696: /* take a matrix meant for k-forms and expand it to one for Ncopies */
1697: static PetscErrorCode PetscDualSpaceLagrangeMatrixCreateCopies(Mat A, PetscInt Nk, PetscInt Ncopies, Mat *Abs)
1698: {
1699:   PetscInt m, n, i, j, k;
1700:   PetscInt maxnnz, *nnz, *iwork;
1701:   Mat      Ac;

1703:   MatGetSize(A, &m, &n);
1705:   PetscMalloc1(m * Ncopies, &nnz);
1706:   for (i = 0, maxnnz = 0; i < m; i++) {
1707:     PetscInt innz;
1708:     MatGetRow(A, i, &innz, NULL, NULL);
1710:     for (j = 0; j < Ncopies; j++) nnz[i * Ncopies + j] = innz;
1711:     maxnnz = PetscMax(maxnnz, innz);
1712:   }
1713:   MatCreateSeqAIJ(PETSC_COMM_SELF, m * Ncopies, n * Ncopies, 0, nnz, &Ac);
1714:   MatSetOption(Ac, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1715:   PetscFree(nnz);
1716:   PetscMalloc1(maxnnz, &iwork);
1717:   for (i = 0; i < m; i++) {
1718:     PetscInt           innz;
1719:     const PetscInt    *cols;
1720:     const PetscScalar *vals;

1722:     MatGetRow(A, i, &innz, &cols, &vals);
1723:     for (j = 0; j < innz; j++) iwork[j] = (cols[j] / Nk) * (Nk * Ncopies) + (cols[j] % Nk);
1724:     for (j = 0; j < Ncopies; j++) {
1725:       PetscInt row = i * Ncopies + j;

1727:       MatSetValues(Ac, 1, &row, innz, iwork, vals, INSERT_VALUES);
1728:       for (k = 0; k < innz; k++) iwork[k] += Nk;
1729:     }
1730:     MatRestoreRow(A, i, &innz, &cols, &vals);
1731:   }
1732:   PetscFree(iwork);
1733:   MatAssemblyBegin(Ac, MAT_FINAL_ASSEMBLY);
1734:   MatAssemblyEnd(Ac, MAT_FINAL_ASSEMBLY);
1735:   *Abs = Ac;
1736:   return 0;
1737: }

1739: /* check if a cell is a tensor product of the segment with a facet,
1740:  * specifically checking if f and f2 can be the "endpoints" (like the triangles
1741:  * at either end of a wedge) */
1742: static PetscErrorCode DMPlexPointIsTensor_Internal_Given(DM dm, PetscInt p, PetscInt f, PetscInt f2, PetscBool *isTensor)
1743: {
1744:   PetscInt        coneSize, c;
1745:   const PetscInt *cone;
1746:   const PetscInt *fCone;
1747:   const PetscInt *f2Cone;
1748:   PetscInt        fs[2];
1749:   PetscInt        meetSize, nmeet;
1750:   const PetscInt *meet;

1752:   fs[0] = f;
1753:   fs[1] = f2;
1754:   DMPlexGetMeet(dm, 2, fs, &meetSize, &meet);
1755:   nmeet = meetSize;
1756:   DMPlexRestoreMeet(dm, 2, fs, &meetSize, &meet);
1757:   /* two points that have a non-empty meet cannot be at opposite ends of a cell */
1758:   if (nmeet) {
1759:     *isTensor = PETSC_FALSE;
1760:     return 0;
1761:   }
1762:   DMPlexGetConeSize(dm, p, &coneSize);
1763:   DMPlexGetCone(dm, p, &cone);
1764:   DMPlexGetCone(dm, f, &fCone);
1765:   DMPlexGetCone(dm, f2, &f2Cone);
1766:   for (c = 0; c < coneSize; c++) {
1767:     PetscInt        e, ef;
1768:     PetscInt        d = -1, d2 = -1;
1769:     PetscInt        dcount, d2count;
1770:     PetscInt        t = cone[c];
1771:     PetscInt        tConeSize;
1772:     PetscBool       tIsTensor;
1773:     const PetscInt *tCone;

1775:     if (t == f || t == f2) continue;
1776:     /* for every other facet in the cone, check that is has
1777:      * one ridge in common with each end */
1778:     DMPlexGetConeSize(dm, t, &tConeSize);
1779:     DMPlexGetCone(dm, t, &tCone);

1781:     dcount  = 0;
1782:     d2count = 0;
1783:     for (e = 0; e < tConeSize; e++) {
1784:       PetscInt q = tCone[e];
1785:       for (ef = 0; ef < coneSize - 2; ef++) {
1786:         if (fCone[ef] == q) {
1787:           if (dcount) {
1788:             *isTensor = PETSC_FALSE;
1789:             return 0;
1790:           }
1791:           d = q;
1792:           dcount++;
1793:         } else if (f2Cone[ef] == q) {
1794:           if (d2count) {
1795:             *isTensor = PETSC_FALSE;
1796:             return 0;
1797:           }
1798:           d2 = q;
1799:           d2count++;
1800:         }
1801:       }
1802:     }
1803:     /* if the whole cell is a tensor with the segment, then this
1804:      * facet should be a tensor with the segment */
1805:     DMPlexPointIsTensor_Internal_Given(dm, t, d, d2, &tIsTensor);
1806:     if (!tIsTensor) {
1807:       *isTensor = PETSC_FALSE;
1808:       return 0;
1809:     }
1810:   }
1811:   *isTensor = PETSC_TRUE;
1812:   return 0;
1813: }

1815: /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair
1816:  * that could be the opposite ends */
1817: static PetscErrorCode DMPlexPointIsTensor_Internal(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB)
1818: {
1819:   PetscInt        coneSize, c, c2;
1820:   const PetscInt *cone;

1822:   DMPlexGetConeSize(dm, p, &coneSize);
1823:   if (!coneSize) {
1824:     if (isTensor) *isTensor = PETSC_FALSE;
1825:     if (endA) *endA = -1;
1826:     if (endB) *endB = -1;
1827:   }
1828:   DMPlexGetCone(dm, p, &cone);
1829:   for (c = 0; c < coneSize; c++) {
1830:     PetscInt f = cone[c];
1831:     PetscInt fConeSize;

1833:     DMPlexGetConeSize(dm, f, &fConeSize);
1834:     if (fConeSize != coneSize - 2) continue;

1836:     for (c2 = c + 1; c2 < coneSize; c2++) {
1837:       PetscInt  f2 = cone[c2];
1838:       PetscBool isTensorff2;
1839:       PetscInt  f2ConeSize;

1841:       DMPlexGetConeSize(dm, f2, &f2ConeSize);
1842:       if (f2ConeSize != coneSize - 2) continue;

1844:       DMPlexPointIsTensor_Internal_Given(dm, p, f, f2, &isTensorff2);
1845:       if (isTensorff2) {
1846:         if (isTensor) *isTensor = PETSC_TRUE;
1847:         if (endA) *endA = f;
1848:         if (endB) *endB = f2;
1849:         return 0;
1850:       }
1851:     }
1852:   }
1853:   if (isTensor) *isTensor = PETSC_FALSE;
1854:   if (endA) *endA = -1;
1855:   if (endB) *endB = -1;
1856:   return 0;
1857: }

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

1865:   DMPlexIsInterpolated(dm, &interpolated);
1867:   DMPlexPointIsTensor_Internal(dm, p, isTensor, endA, endB);
1868:   return 0;
1869: }

1871: /* Let k = formDegree and k' = -sign(k) * dim + k.  Transform a symmetric frame for k-forms on the biunit simplex into
1872:  * a symmetric frame for k'-forms on the biunit simplex.
1873:  *
1874:  * A frame is "symmetric" if the pullback of every symmetry of the biunit simplex is a permutation of the frame.
1875:  *
1876:  * forms in the symmetric frame are used as dofs in the untrimmed simplex spaces.  This way, symmetries of the
1877:  * reference cell result in permutations of dofs grouped by node.
1878:  *
1879:  * Use T to transform dof matrices for k'-forms into dof matrices for k-forms as a block diagonal transformation on
1880:  * the right.
1881:  */
1882: static PetscErrorCode BiunitSimplexSymmetricFormTransformation(PetscInt dim, PetscInt formDegree, PetscReal T[])
1883: {
1884:   PetscInt   k  = formDegree;
1885:   PetscInt   kd = k < 0 ? dim + k : k - dim;
1886:   PetscInt   Nk;
1887:   PetscReal *biToEq, *eqToBi, *biToEqStar, *eqToBiStar;
1888:   PetscInt   fact;

1890:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1891:   PetscCalloc4(dim * dim, &biToEq, dim * dim, &eqToBi, Nk * Nk, &biToEqStar, Nk * Nk, &eqToBiStar);
1892:   /* fill in biToEq: Jacobian of the transformation from the biunit simplex to the equilateral simplex */
1893:   fact = 0;
1894:   for (PetscInt i = 0; i < dim; i++) {
1895:     biToEq[i * dim + i] = PetscSqrtReal(((PetscReal)i + 2.) / (2. * ((PetscReal)i + 1.)));
1896:     fact += 4 * (i + 1);
1897:     for (PetscInt j = i + 1; j < dim; j++) biToEq[i * dim + j] = PetscSqrtReal(1. / (PetscReal)fact);
1898:   }
1899:   /* fill in eqToBi: Jacobian of the transformation from the equilateral simplex to the biunit simplex */
1900:   fact = 0;
1901:   for (PetscInt j = 0; j < dim; j++) {
1902:     eqToBi[j * dim + j] = PetscSqrtReal(2. * ((PetscReal)j + 1.) / ((PetscReal)j + 2));
1903:     fact += j + 1;
1904:     for (PetscInt i = 0; i < j; i++) eqToBi[i * dim + j] = -PetscSqrtReal(1. / (PetscReal)fact);
1905:   }
1906:   PetscDTAltVPullbackMatrix(dim, dim, biToEq, kd, biToEqStar);
1907:   PetscDTAltVPullbackMatrix(dim, dim, eqToBi, k, eqToBiStar);
1908:   /* product of pullbacks simulates the following steps
1909:    *
1910:    * 1. start with frame W = [w_1, w_2, ..., w_m] of k forms that is symmetric on the biunit simplex:
1911:           if J is the Jacobian of a symmetry of the biunit simplex, then J_k* W = [J_k*w_1, ..., J_k*w_m]
1912:           is a permutation of W.
1913:           Even though a k' form --- a (dim - k) form represented by its Hodge star --- has the same geometric
1914:           content as a k form, W is not a symmetric frame of k' forms on the biunit simplex.  That's because,
1915:           for general Jacobian J, J_k* != J_k'*.
1916:    * 2. pullback W to the equilateral triangle using the k pullback, W_eq = eqToBi_k* W.  All symmetries of the
1917:           equilateral simplex have orthonormal Jacobians.  For an orthonormal Jacobian O, J_k* = J_k'*, so W_eq is
1918:           also a symmetric frame for k' forms on the equilateral simplex.
1919:      3. pullback W_eq back to the biunit simplex using the k' pulback, V = biToEq_k'* W_eq = biToEq_k'* eqToBi_k* W.
1920:           V is a symmetric frame for k' forms on the biunit simplex.
1921:    */
1922:   for (PetscInt i = 0; i < Nk; i++) {
1923:     for (PetscInt j = 0; j < Nk; j++) {
1924:       PetscReal val = 0.;
1925:       for (PetscInt k = 0; k < Nk; k++) val += biToEqStar[i * Nk + k] * eqToBiStar[k * Nk + j];
1926:       T[i * Nk + j] = val;
1927:     }
1928:   }
1929:   PetscFree4(biToEq, eqToBi, biToEqStar, eqToBiStar);
1930:   return 0;
1931: }

1933: /* permute a quadrature -> dof matrix so that its rows are in revlex order by nodeIdx */
1934: static PetscErrorCode MatPermuteByNodeIdx(Mat A, PetscLagNodeIndices ni, Mat *Aperm)
1935: {
1936:   PetscInt   m, n, i, j;
1937:   PetscInt   nodeIdxDim = ni->nodeIdxDim;
1938:   PetscInt   nodeVecDim = ni->nodeVecDim;
1939:   PetscInt  *perm;
1940:   IS         permIS;
1941:   IS         id;
1942:   PetscInt  *nIdxPerm;
1943:   PetscReal *nVecPerm;

1945:   PetscLagNodeIndicesGetPermutation(ni, &perm);
1946:   MatGetSize(A, &m, &n);
1947:   PetscMalloc1(nodeIdxDim * m, &nIdxPerm);
1948:   PetscMalloc1(nodeVecDim * m, &nVecPerm);
1949:   for (i = 0; i < m; i++)
1950:     for (j = 0; j < nodeIdxDim; j++) nIdxPerm[i * nodeIdxDim + j] = ni->nodeIdx[perm[i] * nodeIdxDim + j];
1951:   for (i = 0; i < m; i++)
1952:     for (j = 0; j < nodeVecDim; j++) nVecPerm[i * nodeVecDim + j] = ni->nodeVec[perm[i] * nodeVecDim + j];
1953:   ISCreateGeneral(PETSC_COMM_SELF, m, perm, PETSC_USE_POINTER, &permIS);
1954:   ISSetPermutation(permIS);
1955:   ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &id);
1956:   ISSetPermutation(id);
1957:   MatPermute(A, permIS, id, Aperm);
1958:   ISDestroy(&permIS);
1959:   ISDestroy(&id);
1960:   for (i = 0; i < m; i++) perm[i] = i;
1961:   PetscFree(ni->nodeIdx);
1962:   PetscFree(ni->nodeVec);
1963:   ni->nodeIdx = nIdxPerm;
1964:   ni->nodeVec = nVecPerm;
1965:   return 0;
1966: }

1968: static PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1969: {
1970:   PetscDualSpace_Lag    *lag   = (PetscDualSpace_Lag *)sp->data;
1971:   DM                     dm    = sp->dm;
1972:   DM                     dmint = NULL;
1973:   PetscInt               order;
1974:   PetscInt               Nc = sp->Nc;
1975:   MPI_Comm               comm;
1976:   PetscBool              continuous;
1977:   PetscSection           section;
1978:   PetscInt               depth, dim, pStart, pEnd, cStart, cEnd, p, *pStratStart, *pStratEnd, d;
1979:   PetscInt               formDegree, Nk, Ncopies;
1980:   PetscInt               tensorf = -1, tensorf2 = -1;
1981:   PetscBool              tensorCell, tensorSpace;
1982:   PetscBool              uniform, trimmed;
1983:   Petsc1DNodeFamily      nodeFamily;
1984:   PetscInt               numNodeSkip;
1985:   DMPlexInterpolatedFlag interpolated;
1986:   PetscBool              isbdm;

1988:   /* step 1: sanitize input */
1989:   PetscObjectGetComm((PetscObject)sp, &comm);
1990:   DMGetDimension(dm, &dim);
1991:   PetscObjectTypeCompare((PetscObject)sp, PETSCDUALSPACEBDM, &isbdm);
1992:   if (isbdm) {
1993:     sp->k = -(dim - 1); /* form degree of H-div */
1994:     PetscObjectChangeTypeName((PetscObject)sp, PETSCDUALSPACELAGRANGE);
1995:   }
1996:   PetscDualSpaceGetFormDegree(sp, &formDegree);
1998:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
1999:   if (sp->Nc <= 0 && lag->numCopies > 0) sp->Nc = Nk * lag->numCopies;
2000:   Nc = sp->Nc;
2002:   if (lag->numCopies <= 0) lag->numCopies = Nc / Nk;
2003:   Ncopies = lag->numCopies;
2005:   if (!dim) sp->order = 0;
2006:   order   = sp->order;
2007:   uniform = sp->uniform;
2009:   if (lag->trimmed && !formDegree) lag->trimmed = PETSC_FALSE; /* trimmed spaces are the same as full spaces for 0-forms */
2010:   if (lag->nodeType == PETSCDTNODES_DEFAULT) {
2011:     lag->nodeType     = PETSCDTNODES_GAUSSJACOBI;
2012:     lag->nodeExponent = 0.;
2013:     /* trimmed spaces don't include corner vertices, so don't use end nodes by default */
2014:     lag->endNodes = lag->trimmed ? PETSC_FALSE : PETSC_TRUE;
2015:   }
2016:   /* If a trimmed space and the user did choose nodes with endpoints, skip them by default */
2017:   if (lag->numNodeSkip < 0) lag->numNodeSkip = (lag->trimmed && lag->endNodes) ? 1 : 0;
2018:   numNodeSkip = lag->numNodeSkip;
2020:   if (lag->trimmed && PetscAbsInt(formDegree) == dim) { /* convert trimmed n-forms to untrimmed of one polynomial order less */
2021:     sp->order--;
2022:     order--;
2023:     lag->trimmed = PETSC_FALSE;
2024:   }
2025:   trimmed = lag->trimmed;
2026:   if (!order || PetscAbsInt(formDegree) == dim) lag->continuous = PETSC_FALSE;
2027:   continuous = lag->continuous;
2028:   DMPlexGetDepth(dm, &depth);
2029:   DMPlexGetChart(dm, &pStart, &pEnd);
2030:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2033:   DMPlexIsInterpolated(dm, &interpolated);
2034:   if (interpolated != DMPLEX_INTERPOLATED_FULL) {
2035:     DMPlexInterpolate(dm, &dmint);
2036:   } else {
2037:     PetscObjectReference((PetscObject)dm);
2038:     dmint = dm;
2039:   }
2040:   tensorCell = PETSC_FALSE;
2041:   if (dim > 1) DMPlexPointIsTensor(dmint, 0, &tensorCell, &tensorf, &tensorf2);
2042:   lag->tensorCell = tensorCell;
2043:   if (dim < 2 || !lag->tensorCell) lag->tensorSpace = PETSC_FALSE;
2044:   tensorSpace = lag->tensorSpace;
2045:   if (!lag->nodeFamily) Petsc1DNodeFamilyCreate(lag->nodeType, lag->nodeExponent, lag->endNodes, &lag->nodeFamily);
2046:   nodeFamily = lag->nodeFamily;

2049:   /* step 2: construct the boundary spaces */
2050:   PetscMalloc2(depth + 1, &pStratStart, depth + 1, &pStratEnd);
2051:   PetscCalloc1(pEnd, &(sp->pointSpaces));
2052:   for (d = 0; d <= depth; ++d) DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);
2053:   PetscDualSpaceSectionCreate_Internal(sp, &section);
2054:   sp->pointSection = section;
2055:   if (continuous && !(lag->interiorOnly)) {
2056:     PetscInt h;

2058:     for (p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */
2059:       PetscReal      v0[3];
2060:       DMPolytopeType ptype;
2061:       PetscReal      J[9], detJ;
2062:       PetscInt       q;

2064:       DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, NULL, &detJ);
2065:       DMPlexGetCellType(dm, p, &ptype);

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

2071:         DMPlexGetCellType(dm, q, &qtype);
2072:         if (qtype == ptype) break;
2073:       }
2074:       if (q < p) { /* this facet has the same dual space as that one */
2075:         PetscObjectReference((PetscObject)sp->pointSpaces[q]);
2076:         sp->pointSpaces[p] = sp->pointSpaces[q];
2077:         continue;
2078:       }
2079:       /* if not, recursively compute this dual space */
2080:       PetscDualSpaceCreateFacetSubspace_Lagrange(sp, NULL, p, formDegree, Ncopies, PETSC_FALSE, &sp->pointSpaces[p]);
2081:     }
2082:     for (h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */
2083:       PetscInt hd   = depth - h;
2084:       PetscInt hdim = dim - h;

2086:       if (hdim < PetscAbsInt(formDegree)) break;
2087:       for (p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
2088:         PetscInt        suppSize, s;
2089:         const PetscInt *supp;

2091:         DMPlexGetSupportSize(dm, p, &suppSize);
2092:         DMPlexGetSupport(dm, p, &supp);
2093:         for (s = 0; s < suppSize; s++) {
2094:           DM              qdm;
2095:           PetscDualSpace  qsp, psp;
2096:           PetscInt        c, coneSize, q;
2097:           const PetscInt *cone;
2098:           const PetscInt *refCone;

2100:           q   = supp[0];
2101:           qsp = sp->pointSpaces[q];
2102:           DMPlexGetConeSize(dm, q, &coneSize);
2103:           DMPlexGetCone(dm, q, &cone);
2104:           for (c = 0; c < coneSize; c++)
2105:             if (cone[c] == p) break;
2107:           PetscDualSpaceGetDM(qsp, &qdm);
2108:           DMPlexGetCone(qdm, 0, &refCone);
2109:           /* get the equivalent dual space from the support dual space */
2110:           PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp);
2111:           if (!s) {
2112:             PetscObjectReference((PetscObject)psp);
2113:             sp->pointSpaces[p] = psp;
2114:           }
2115:         }
2116:       }
2117:     }
2118:     for (p = 1; p < pEnd; p++) {
2119:       PetscInt pspdim;
2120:       if (!sp->pointSpaces[p]) continue;
2121:       PetscDualSpaceGetInteriorDimension(sp->pointSpaces[p], &pspdim);
2122:       PetscSectionSetDof(section, p, pspdim);
2123:     }
2124:   }

2126:   if (Ncopies > 1) {
2127:     Mat                 intMatScalar, allMatScalar;
2128:     PetscDualSpace      scalarsp;
2129:     PetscDualSpace_Lag *scalarlag;

2131:     PetscDualSpaceDuplicate(sp, &scalarsp);
2132:     /* Setting the number of components to Nk is a space with 1 copy of each k-form */
2133:     PetscDualSpaceSetNumComponents(scalarsp, Nk);
2134:     PetscDualSpaceSetUp(scalarsp);
2135:     PetscDualSpaceGetInteriorData(scalarsp, &(sp->intNodes), &intMatScalar);
2136:     PetscObjectReference((PetscObject)(sp->intNodes));
2137:     if (intMatScalar) PetscDualSpaceLagrangeMatrixCreateCopies(intMatScalar, Nk, Ncopies, &(sp->intMat));
2138:     PetscDualSpaceGetAllData(scalarsp, &(sp->allNodes), &allMatScalar);
2139:     PetscObjectReference((PetscObject)(sp->allNodes));
2140:     PetscDualSpaceLagrangeMatrixCreateCopies(allMatScalar, Nk, Ncopies, &(sp->allMat));
2141:     sp->spdim    = scalarsp->spdim * Ncopies;
2142:     sp->spintdim = scalarsp->spintdim * Ncopies;
2143:     scalarlag    = (PetscDualSpace_Lag *)scalarsp->data;
2144:     PetscLagNodeIndicesReference(scalarlag->vertIndices);
2145:     lag->vertIndices = scalarlag->vertIndices;
2146:     PetscLagNodeIndicesReference(scalarlag->intNodeIndices);
2147:     lag->intNodeIndices = scalarlag->intNodeIndices;
2148:     PetscLagNodeIndicesReference(scalarlag->allNodeIndices);
2149:     lag->allNodeIndices = scalarlag->allNodeIndices;
2150:     PetscDualSpaceDestroy(&scalarsp);
2151:     PetscSectionSetDof(section, 0, sp->spintdim);
2152:     PetscDualSpaceSectionSetUp_Internal(sp, section);
2153:     PetscDualSpaceComputeFunctionalsFromAllData(sp);
2154:     PetscFree2(pStratStart, pStratEnd);
2155:     DMDestroy(&dmint);
2156:     return 0;
2157:   }

2159:   if (trimmed && !continuous) {
2160:     /* the dofs of a trimmed space don't have a nice tensor/lattice structure:
2161:      * just construct the continuous dual space and copy all of the data over,
2162:      * allocating it all to the cell instead of splitting it up between the boundaries */
2163:     PetscDualSpace      spcont;
2164:     PetscInt            spdim, f;
2165:     PetscQuadrature     allNodes;
2166:     PetscDualSpace_Lag *lagc;
2167:     Mat                 allMat;

2169:     PetscDualSpaceDuplicate(sp, &spcont);
2170:     PetscDualSpaceLagrangeSetContinuity(spcont, PETSC_TRUE);
2171:     PetscDualSpaceSetUp(spcont);
2172:     PetscDualSpaceGetDimension(spcont, &spdim);
2173:     sp->spdim = sp->spintdim = spdim;
2174:     PetscSectionSetDof(section, 0, spdim);
2175:     PetscDualSpaceSectionSetUp_Internal(sp, section);
2176:     PetscMalloc1(spdim, &(sp->functional));
2177:     for (f = 0; f < spdim; f++) {
2178:       PetscQuadrature fn;

2180:       PetscDualSpaceGetFunctional(spcont, f, &fn);
2181:       PetscObjectReference((PetscObject)fn);
2182:       sp->functional[f] = fn;
2183:     }
2184:     PetscDualSpaceGetAllData(spcont, &allNodes, &allMat);
2185:     PetscObjectReference((PetscObject)allNodes);
2186:     PetscObjectReference((PetscObject)allNodes);
2187:     sp->allNodes = sp->intNodes = allNodes;
2188:     PetscObjectReference((PetscObject)allMat);
2189:     PetscObjectReference((PetscObject)allMat);
2190:     sp->allMat = sp->intMat = allMat;
2191:     lagc                    = (PetscDualSpace_Lag *)spcont->data;
2192:     PetscLagNodeIndicesReference(lagc->vertIndices);
2193:     lag->vertIndices = lagc->vertIndices;
2194:     PetscLagNodeIndicesReference(lagc->allNodeIndices);
2195:     PetscLagNodeIndicesReference(lagc->allNodeIndices);
2196:     lag->intNodeIndices = lagc->allNodeIndices;
2197:     lag->allNodeIndices = lagc->allNodeIndices;
2198:     PetscDualSpaceDestroy(&spcont);
2199:     PetscFree2(pStratStart, pStratEnd);
2200:     DMDestroy(&dmint);
2201:     return 0;
2202:   }

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

2208:     if (trimmed) {
2209:       /* there is one dof in the interior of the a trimmed element for each full polynomial of with degree at most
2210:        * order + k - dim - 1 */
2211:       if (order + PetscAbsInt(formDegree) > dim) {
2212:         PetscInt sum = order + PetscAbsInt(formDegree) - dim - 1;
2213:         PetscInt nDofs;

2215:         PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));
2216:         MatGetSize(sp->intMat, &nDofs, NULL);
2217:         PetscSectionSetDof(section, 0, nDofs);
2218:       }
2219:       PetscDualSpaceSectionSetUp_Internal(sp, section);
2220:       PetscDualSpaceCreateAllDataFromInteriorData(sp);
2221:       PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2222:     } else {
2223:       if (!continuous) {
2224:         /* if discontinuous just construct one node for each set of dofs (a set of dofs is a basis for the k-form
2225:          * space) */
2226:         PetscInt sum = order;
2227:         PetscInt nDofs;

2229:         PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));
2230:         MatGetSize(sp->intMat, &nDofs, NULL);
2231:         PetscSectionSetDof(section, 0, nDofs);
2232:         PetscDualSpaceSectionSetUp_Internal(sp, section);
2233:         PetscObjectReference((PetscObject)(sp->intNodes));
2234:         sp->allNodes = sp->intNodes;
2235:         PetscObjectReference((PetscObject)(sp->intMat));
2236:         sp->allMat = sp->intMat;
2237:         PetscLagNodeIndicesReference(lag->intNodeIndices);
2238:         lag->allNodeIndices = lag->intNodeIndices;
2239:       } else {
2240:         /* there is one dof in the interior of the a full element for each trimmed polynomial of with degree at most
2241:          * order + k - dim, but with complementary form degree */
2242:         if (order + PetscAbsInt(formDegree) > dim) {
2243:           PetscDualSpace      trimmedsp;
2244:           PetscDualSpace_Lag *trimmedlag;
2245:           PetscQuadrature     intNodes;
2246:           PetscInt            trFormDegree = formDegree >= 0 ? formDegree - dim : dim - PetscAbsInt(formDegree);
2247:           PetscInt            nDofs;
2248:           Mat                 intMat;

2250:           PetscDualSpaceDuplicate(sp, &trimmedsp);
2251:           PetscDualSpaceLagrangeSetTrimmed(trimmedsp, PETSC_TRUE);
2252:           PetscDualSpaceSetOrder(trimmedsp, order + PetscAbsInt(formDegree) - dim);
2253:           PetscDualSpaceSetFormDegree(trimmedsp, trFormDegree);
2254:           trimmedlag              = (PetscDualSpace_Lag *)trimmedsp->data;
2255:           trimmedlag->numNodeSkip = numNodeSkip + 1;
2256:           PetscDualSpaceSetUp(trimmedsp);
2257:           PetscDualSpaceGetAllData(trimmedsp, &intNodes, &intMat);
2258:           PetscObjectReference((PetscObject)intNodes);
2259:           sp->intNodes = intNodes;
2260:           PetscLagNodeIndicesReference(trimmedlag->allNodeIndices);
2261:           lag->intNodeIndices = trimmedlag->allNodeIndices;
2262:           PetscObjectReference((PetscObject)intMat);
2263:           if (PetscAbsInt(formDegree) > 0 && PetscAbsInt(formDegree) < dim) {
2264:             PetscReal   *T;
2265:             PetscScalar *work;
2266:             PetscInt     nCols, nRows;
2267:             Mat          intMatT;

2269:             MatDuplicate(intMat, MAT_COPY_VALUES, &intMatT);
2270:             MatGetSize(intMat, &nRows, &nCols);
2271:             PetscMalloc2(Nk * Nk, &T, nCols, &work);
2272:             BiunitSimplexSymmetricFormTransformation(dim, formDegree, T);
2273:             for (PetscInt row = 0; row < nRows; row++) {
2274:               PetscInt           nrCols;
2275:               const PetscInt    *rCols;
2276:               const PetscScalar *rVals;

2278:               MatGetRow(intMat, row, &nrCols, &rCols, &rVals);
2280:               for (PetscInt b = 0; b < nrCols; b += Nk) {
2281:                 const PetscScalar *v = &rVals[b];
2282:                 PetscScalar       *w = &work[b];
2283:                 for (PetscInt j = 0; j < Nk; j++) {
2284:                   w[j] = 0.;
2285:                   for (PetscInt i = 0; i < Nk; i++) w[j] += v[i] * T[i * Nk + j];
2286:                 }
2287:               }
2288:               MatSetValuesBlocked(intMatT, 1, &row, nrCols, rCols, work, INSERT_VALUES);
2289:               MatRestoreRow(intMat, row, &nrCols, &rCols, &rVals);
2290:             }
2291:             MatAssemblyBegin(intMatT, MAT_FINAL_ASSEMBLY);
2292:             MatAssemblyEnd(intMatT, MAT_FINAL_ASSEMBLY);
2293:             MatDestroy(&intMat);
2294:             intMat = intMatT;
2295:             PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));
2296:             PetscLagNodeIndicesDuplicate(trimmedlag->allNodeIndices, &(lag->intNodeIndices));
2297:             {
2298:               PetscInt         nNodes     = lag->intNodeIndices->nNodes;
2299:               PetscReal       *newNodeVec = lag->intNodeIndices->nodeVec;
2300:               const PetscReal *oldNodeVec = trimmedlag->allNodeIndices->nodeVec;

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

2306:                 for (PetscInt j = 0; j < Nk; j++) {
2307:                   w[j] = 0.;
2308:                   for (PetscInt i = 0; i < Nk; i++) w[j] += v[i] * T[i * Nk + j];
2309:                 }
2310:               }
2311:             }
2312:             PetscFree2(T, work);
2313:           }
2314:           sp->intMat = intMat;
2315:           MatGetSize(sp->intMat, &nDofs, NULL);
2316:           PetscDualSpaceDestroy(&trimmedsp);
2317:           PetscSectionSetDof(section, 0, nDofs);
2318:         }
2319:         PetscDualSpaceSectionSetUp_Internal(sp, section);
2320:         PetscDualSpaceCreateAllDataFromInteriorData(sp);
2321:         PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2322:       }
2323:     }
2324:   } else {
2325:     PetscQuadrature     intNodesTrace  = NULL;
2326:     PetscQuadrature     intNodesFiber  = NULL;
2327:     PetscQuadrature     intNodes       = NULL;
2328:     PetscLagNodeIndices intNodeIndices = NULL;
2329:     Mat                 intMat         = NULL;

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

2337:       if (sp->pointSpaces[tensorf]) {
2338:         PetscObjectReference((PetscObject)(sp->pointSpaces[tensorf]));
2339:         trace = sp->pointSpaces[tensorf];
2340:       } else {
2341:         PetscDualSpaceCreateFacetSubspace_Lagrange(sp, NULL, tensorf, formDegree, Ncopies, PETSC_TRUE, &trace);
2342:       }
2343:       PetscDualSpaceCreateEdgeSubspace_Lagrange(sp, order, 0, 1, PETSC_TRUE, &fiber);
2344:       tracel = (PetscDualSpace_Lag *)trace->data;
2345:       fiberl = (PetscDualSpace_Lag *)fiber->data;
2346:       PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));
2347:       PetscDualSpaceGetInteriorData(trace, &intNodesTrace, &intMatTrace);
2348:       PetscDualSpaceGetInteriorData(fiber, &intNodesFiber, &intMatFiber);
2349:       if (intNodesTrace && intNodesFiber) {
2350:         PetscQuadratureCreateTensor(intNodesTrace, intNodesFiber, &intNodes);
2351:         MatTensorAltV(intMatTrace, intMatFiber, dim - 1, formDegree, 1, 0, &intMat);
2352:         PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, formDegree, fiberl->intNodeIndices, 1, 0, &intNodeIndices);
2353:       }
2354:       PetscObjectReference((PetscObject)intNodesTrace);
2355:       PetscObjectReference((PetscObject)intNodesFiber);
2356:       PetscDualSpaceDestroy(&fiber);
2357:       PetscDualSpaceDestroy(&trace);
2358:     }
2359:     if (PetscAbsInt(formDegree) > 0) { /* get the trace (k-1)-forms on the first facet, and the 1-forms on the edge,
2360:                                           and wedge them together to create the remaining k-form dofs */
2361:       PetscDualSpace      trace, fiber;
2362:       PetscDualSpace_Lag *tracel, *fiberl;
2363:       PetscQuadrature     intNodesTrace2, intNodesFiber2, intNodes2;
2364:       PetscLagNodeIndices intNodeIndices2;
2365:       Mat                 intMatTrace, intMatFiber, intMat2;
2366:       PetscInt            traceDegree = formDegree > 0 ? formDegree - 1 : formDegree + 1;
2367:       PetscInt            fiberDegree = formDegree > 0 ? 1 : -1;

2369:       PetscDualSpaceCreateFacetSubspace_Lagrange(sp, NULL, tensorf, traceDegree, Ncopies, PETSC_TRUE, &trace);
2370:       PetscDualSpaceCreateEdgeSubspace_Lagrange(sp, order, fiberDegree, 1, PETSC_TRUE, &fiber);
2371:       tracel = (PetscDualSpace_Lag *)trace->data;
2372:       fiberl = (PetscDualSpace_Lag *)fiber->data;
2373:       if (!lag->vertIndices) PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));
2374:       PetscDualSpaceGetInteriorData(trace, &intNodesTrace2, &intMatTrace);
2375:       PetscDualSpaceGetInteriorData(fiber, &intNodesFiber2, &intMatFiber);
2376:       if (intNodesTrace2 && intNodesFiber2) {
2377:         PetscQuadratureCreateTensor(intNodesTrace2, intNodesFiber2, &intNodes2);
2378:         MatTensorAltV(intMatTrace, intMatFiber, dim - 1, traceDegree, 1, fiberDegree, &intMat2);
2379:         PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, traceDegree, fiberl->intNodeIndices, 1, fiberDegree, &intNodeIndices2);
2380:         if (!intMat) {
2381:           intMat         = intMat2;
2382:           intNodes       = intNodes2;
2383:           intNodeIndices = intNodeIndices2;
2384:         } else {
2385:           /* merge the matrices, quadrature points, and nodes */
2386:           PetscInt            nM;
2387:           PetscInt            nDof, nDof2;
2388:           PetscInt           *toMerged = NULL, *toMerged2 = NULL;
2389:           PetscQuadrature     merged               = NULL;
2390:           PetscLagNodeIndices intNodeIndicesMerged = NULL;
2391:           Mat                 matMerged            = NULL;

2393:           MatGetSize(intMat, &nDof, NULL);
2394:           MatGetSize(intMat2, &nDof2, NULL);
2395:           PetscQuadraturePointsMerge(intNodes, intNodes2, &merged, &toMerged, &toMerged2);
2396:           PetscQuadratureGetData(merged, NULL, NULL, &nM, NULL, NULL);
2397:           MatricesMerge(intMat, intMat2, dim, formDegree, nM, toMerged, toMerged2, &matMerged);
2398:           PetscLagNodeIndicesMerge(intNodeIndices, intNodeIndices2, &intNodeIndicesMerged);
2399:           PetscFree(toMerged);
2400:           PetscFree(toMerged2);
2401:           MatDestroy(&intMat);
2402:           MatDestroy(&intMat2);
2403:           PetscQuadratureDestroy(&intNodes);
2404:           PetscQuadratureDestroy(&intNodes2);
2405:           PetscLagNodeIndicesDestroy(&intNodeIndices);
2406:           PetscLagNodeIndicesDestroy(&intNodeIndices2);
2407:           intNodes       = merged;
2408:           intMat         = matMerged;
2409:           intNodeIndices = intNodeIndicesMerged;
2410:           if (!trimmed) {
2411:             /* I think users expect that, when a node has a full basis for the k-forms,
2412:              * they should be consecutive dofs.  That isn't the case for trimmed spaces,
2413:              * but is for some of the nodes in untrimmed spaces, so in that case we
2414:              * sort them to group them by node */
2415:             Mat intMatPerm;

2417:             MatPermuteByNodeIdx(intMat, intNodeIndices, &intMatPerm);
2418:             MatDestroy(&intMat);
2419:             intMat = intMatPerm;
2420:           }
2421:         }
2422:       }
2423:       PetscDualSpaceDestroy(&fiber);
2424:       PetscDualSpaceDestroy(&trace);
2425:     }
2426:     PetscQuadratureDestroy(&intNodesTrace);
2427:     PetscQuadratureDestroy(&intNodesFiber);
2428:     sp->intNodes        = intNodes;
2429:     sp->intMat          = intMat;
2430:     lag->intNodeIndices = intNodeIndices;
2431:     {
2432:       PetscInt nDofs = 0;

2434:       if (intMat) MatGetSize(intMat, &nDofs, NULL);
2435:       PetscSectionSetDof(section, 0, nDofs);
2436:     }
2437:     PetscDualSpaceSectionSetUp_Internal(sp, section);
2438:     if (continuous) {
2439:       PetscDualSpaceCreateAllDataFromInteriorData(sp);
2440:       PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2441:     } else {
2442:       PetscObjectReference((PetscObject)intNodes);
2443:       sp->allNodes = intNodes;
2444:       PetscObjectReference((PetscObject)intMat);
2445:       sp->allMat = intMat;
2446:       PetscLagNodeIndicesReference(intNodeIndices);
2447:       lag->allNodeIndices = intNodeIndices;
2448:     }
2449:   }
2450:   PetscSectionGetStorageSize(section, &sp->spdim);
2451:   PetscSectionGetConstrainedStorageSize(section, &sp->spintdim);
2452:   PetscDualSpaceComputeFunctionalsFromAllData(sp);
2453:   PetscFree2(pStratStart, pStratEnd);
2454:   DMDestroy(&dmint);
2455:   return 0;
2456: }

2458: /* Create a matrix that represents the transformation that DMPlexVecGetClosure() would need
2459:  * to get the representation of the dofs for a mesh point if the mesh point had this orientation
2460:  * relative to the cell */
2461: PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat)
2462: {
2463:   PetscDualSpace_Lag *lag;
2464:   DM                  dm;
2465:   PetscLagNodeIndices vertIndices, intNodeIndices;
2466:   PetscLagNodeIndices ni;
2467:   PetscInt            nodeIdxDim, nodeVecDim, nNodes;
2468:   PetscInt            formDegree;
2469:   PetscInt           *perm, *permOrnt;
2470:   PetscInt           *nnz;
2471:   PetscInt            n;
2472:   PetscInt            maxGroupSize;
2473:   PetscScalar        *V, *W, *work;
2474:   Mat                 A;

2476:   if (!sp->spintdim) {
2477:     *symMat = NULL;
2478:     return 0;
2479:   }
2480:   lag            = (PetscDualSpace_Lag *)sp->data;
2481:   vertIndices    = lag->vertIndices;
2482:   intNodeIndices = lag->intNodeIndices;
2483:   PetscDualSpaceGetDM(sp, &dm);
2484:   PetscDualSpaceGetFormDegree(sp, &formDegree);
2485:   PetscNew(&ni);
2486:   ni->refct      = 1;
2487:   ni->nodeIdxDim = nodeIdxDim = intNodeIndices->nodeIdxDim;
2488:   ni->nodeVecDim = nodeVecDim = intNodeIndices->nodeVecDim;
2489:   ni->nNodes = nNodes = intNodeIndices->nNodes;
2490:   PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
2491:   PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
2492:   /* push forward the dofs by the symmetry of the reference element induced by ornt */
2493:   PetscLagNodeIndicesPushForward(dm, vertIndices, 0, vertIndices, intNodeIndices, ornt, formDegree, ni->nodeIdx, ni->nodeVec);
2494:   /* get the revlex order for both the original and transformed dofs */
2495:   PetscLagNodeIndicesGetPermutation(intNodeIndices, &perm);
2496:   PetscLagNodeIndicesGetPermutation(ni, &permOrnt);
2497:   PetscMalloc1(nNodes, &nnz);
2498:   for (n = 0, maxGroupSize = 0; n < nNodes;) { /* incremented in the loop */
2499:     PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2500:     PetscInt  m, nEnd;
2501:     PetscInt  groupSize;
2502:     /* for each group of dofs that have the same nodeIdx coordinate */
2503:     for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2504:       PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2505:       PetscInt  d;

2507:       /* compare the oriented permutation indices */
2508:       for (d = 0; d < nodeIdxDim; d++)
2509:         if (mind[d] != nind[d]) break;
2510:       if (d < nodeIdxDim) break;
2511:     }
2512:     /* permOrnt[[n, nEnd)] is a group of dofs that, under the symmetry are at the same location */

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

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

2524:         /* compare the oriented permutation indices */
2525:         for (d = 0; d < nodeIdxDim; d++)
2526:           if (mind[d] != nind[d]) break;
2527:         if (d < nodeIdxDim) break;
2528:       }
2530:     }
2531:     groupSize = nEnd - n;
2532:     /* each pushforward dof vector will be expressed in a basis of the unpermuted dofs */
2533:     for (m = n; m < nEnd; m++) nnz[permOrnt[m]] = groupSize;

2535:     maxGroupSize = PetscMax(maxGroupSize, nEnd - n);
2536:     n            = nEnd;
2537:   }
2539:   MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes, nNodes, 0, nnz, &A);
2540:   PetscFree(nnz);
2541:   PetscMalloc3(maxGroupSize * nodeVecDim, &V, maxGroupSize * nodeVecDim, &W, nodeVecDim * 2, &work);
2542:   for (n = 0; n < nNodes;) { /* incremented in the loop */
2543:     PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2544:     PetscInt  nEnd;
2545:     PetscInt  m;
2546:     PetscInt  groupSize;
2547:     for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2548:       PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2549:       PetscInt  d;

2551:       /* compare the oriented permutation indices */
2552:       for (d = 0; d < nodeIdxDim; d++)
2553:         if (mind[d] != nind[d]) break;
2554:       if (d < nodeIdxDim) break;
2555:     }
2556:     groupSize = nEnd - n;
2557:     /* get all of the vectors from the original and all of the pushforward vectors */
2558:     for (m = n; m < nEnd; m++) {
2559:       PetscInt d;

2561:       for (d = 0; d < nodeVecDim; d++) {
2562:         V[(m - n) * nodeVecDim + d] = intNodeIndices->nodeVec[perm[m] * nodeVecDim + d];
2563:         W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
2564:       }
2565:     }
2566:     /* now we have to solve for W in terms of V: the systems isn't always square, but the span
2567:      * of V and W should always be the same, so the solution of the normal equations works */
2568:     {
2569:       char         transpose = 'N';
2570:       PetscBLASInt bm        = nodeVecDim;
2571:       PetscBLASInt bn        = groupSize;
2572:       PetscBLASInt bnrhs     = groupSize;
2573:       PetscBLASInt blda      = bm;
2574:       PetscBLASInt bldb      = bm;
2575:       PetscBLASInt blwork    = 2 * nodeVecDim;
2576:       PetscBLASInt info;

2578:       PetscCallBLAS("LAPACKgels", LAPACKgels_(&transpose, &bm, &bn, &bnrhs, V, &blda, W, &bldb, work, &blwork, &info));
2580:       /* repack */
2581:       {
2582:         PetscInt i, j;

2584:         for (i = 0; i < groupSize; i++) {
2585:           for (j = 0; j < groupSize; j++) {
2586:             /* notice the different leading dimension */
2587:             V[i * groupSize + j] = W[i * nodeVecDim + j];
2588:           }
2589:         }
2590:       }
2591:       if (PetscDefined(USE_DEBUG)) {
2592:         PetscReal res;

2594:         /* check that the normal error is 0 */
2595:         for (m = n; m < nEnd; m++) {
2596:           PetscInt d;

2598:           for (d = 0; d < nodeVecDim; d++) W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
2599:         }
2600:         res = 0.;
2601:         for (PetscInt i = 0; i < groupSize; i++) {
2602:           for (PetscInt j = 0; j < nodeVecDim; j++) {
2603:             for (PetscInt k = 0; k < groupSize; k++) W[i * nodeVecDim + j] -= V[i * groupSize + k] * intNodeIndices->nodeVec[perm[n + k] * nodeVecDim + j];
2604:             res += PetscAbsScalar(W[i * nodeVecDim + j]);
2605:           }
2606:         }
2608:       }
2609:     }
2610:     MatSetValues(A, groupSize, &permOrnt[n], groupSize, &perm[n], V, INSERT_VALUES);
2611:     n = nEnd;
2612:   }
2613:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
2614:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
2615:   *symMat = A;
2616:   PetscFree3(V, W, work);
2617:   PetscLagNodeIndicesDestroy(&ni);
2618:   return 0;
2619: }

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

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

2625: /* the existing interface for symmetries is insufficient for all cases:
2626:  * - it should be sufficient for form degrees that are scalar (0 and n)
2627:  * - it should be sufficient for hypercube dofs
2628:  * - it isn't sufficient for simplex cells with non-scalar form degrees if
2629:  *   there are any dofs in the interior
2630:  *
2631:  * We compute the general transformation matrices, and if they fit, we return them,
2632:  * otherwise we error (but we should probably change the interface to allow for
2633:  * these symmetries)
2634:  */
2635: static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
2636: {
2637:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2638:   PetscInt            dim, order, Nc;

2640:   PetscDualSpaceGetOrder(sp, &order);
2641:   PetscDualSpaceGetNumComponents(sp, &Nc);
2642:   DMGetDimension(sp->dm, &dim);
2643:   if (!lag->symComputed) { /* store symmetries */
2644:     PetscInt       pStart, pEnd, p;
2645:     PetscInt       numPoints;
2646:     PetscInt       numFaces;
2647:     PetscInt       spintdim;
2648:     PetscInt    ***symperms;
2649:     PetscScalar ***symflips;

2651:     DMPlexGetChart(sp->dm, &pStart, &pEnd);
2652:     numPoints = pEnd - pStart;
2653:     {
2654:       DMPolytopeType ct;
2655:       /* The number of arrangements is no longer based on the number of faces */
2656:       DMPlexGetCellType(sp->dm, 0, &ct);
2657:       numFaces = DMPolytopeTypeGetNumArrangments(ct) / 2;
2658:     }
2659:     PetscCalloc1(numPoints, &symperms);
2660:     PetscCalloc1(numPoints, &symflips);
2661:     spintdim = sp->spintdim;
2662:     /* The nodal symmetry behavior is not present when tensorSpace != tensorCell: someone might want this for the "S"
2663:      * family of FEEC spaces.  Most used in particular are discontinuous polynomial L2 spaces in tensor cells, where
2664:      * the symmetries are not necessary for FE assembly.  So for now we assume this is the case and don't return
2665:      * symmetries if tensorSpace != tensorCell */
2666:     if (spintdim && 0 < dim && dim < 3 && (lag->tensorSpace == lag->tensorCell)) { /* compute self symmetries */
2667:       PetscInt    **cellSymperms;
2668:       PetscScalar **cellSymflips;
2669:       PetscInt      ornt;
2670:       PetscInt      nCopies = Nc / lag->intNodeIndices->nodeVecDim;
2671:       PetscInt      nNodes  = lag->intNodeIndices->nNodes;

2673:       lag->numSelfSym = 2 * numFaces;
2674:       lag->selfSymOff = numFaces;
2675:       PetscCalloc1(2 * numFaces, &cellSymperms);
2676:       PetscCalloc1(2 * numFaces, &cellSymflips);
2677:       /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
2678:       symperms[0] = &cellSymperms[numFaces];
2679:       symflips[0] = &cellSymflips[numFaces];
2682:       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 */
2683:         Mat          symMat;
2684:         PetscInt    *perm;
2685:         PetscScalar *flips;
2686:         PetscInt     i;

2688:         if (!ornt) continue;
2689:         PetscMalloc1(spintdim, &perm);
2690:         PetscCalloc1(spintdim, &flips);
2691:         for (i = 0; i < spintdim; i++) perm[i] = -1;
2692:         PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, ornt, &symMat);
2693:         for (i = 0; i < nNodes; i++) {
2694:           PetscInt           ncols;
2695:           PetscInt           j, k;
2696:           const PetscInt    *cols;
2697:           const PetscScalar *vals;
2698:           PetscBool          nz_seen = PETSC_FALSE;

2700:           MatGetRow(symMat, i, &ncols, &cols, &vals);
2701:           for (j = 0; j < ncols; j++) {
2702:             if (PetscAbsScalar(vals[j]) > PETSC_SMALL) {
2704:               nz_seen = PETSC_TRUE;
2708:               for (k = 0; k < nCopies; k++) perm[cols[j] * nCopies + k] = i * nCopies + k;
2709:               if (PetscRealPart(vals[j]) < 0.) {
2710:                 for (k = 0; k < nCopies; k++) flips[i * nCopies + k] = -1.;
2711:               } else {
2712:                 for (k = 0; k < nCopies; k++) flips[i * nCopies + k] = 1.;
2713:               }
2714:             }
2715:           }
2716:           MatRestoreRow(symMat, i, &ncols, &cols, &vals);
2717:         }
2718:         MatDestroy(&symMat);
2719:         /* if there were no sign flips, keep NULL */
2720:         for (i = 0; i < spintdim; i++)
2721:           if (flips[i] != 1.) break;
2722:         if (i == spintdim) {
2723:           PetscFree(flips);
2724:           flips = NULL;
2725:         }
2726:         /* if the permutation is identity, keep NULL */
2727:         for (i = 0; i < spintdim; i++)
2728:           if (perm[i] != i) break;
2729:         if (i == spintdim) {
2730:           PetscFree(perm);
2731:           perm = NULL;
2732:         }
2733:         symperms[0][ornt] = perm;
2734:         symflips[0][ornt] = flips;
2735:       }
2736:       /* if no orientations produced non-identity permutations, keep NULL */
2737:       for (ornt = -numFaces; ornt < numFaces; ornt++)
2738:         if (symperms[0][ornt]) break;
2739:       if (ornt == numFaces) {
2740:         PetscFree(cellSymperms);
2741:         symperms[0] = NULL;
2742:       }
2743:       /* if no orientations produced sign flips, keep NULL */
2744:       for (ornt = -numFaces; ornt < numFaces; ornt++)
2745:         if (symflips[0][ornt]) break;
2746:       if (ornt == numFaces) {
2747:         PetscFree(cellSymflips);
2748:         symflips[0] = NULL;
2749:       }
2750:     }
2751:     { /* get the symmetries of closure points */
2752:       PetscInt  closureSize = 0;
2753:       PetscInt *closure     = NULL;
2754:       PetscInt  r;

2756:       DMPlexGetTransitiveClosure(sp->dm, 0, PETSC_TRUE, &closureSize, &closure);
2757:       for (r = 0; r < closureSize; r++) {
2758:         PetscDualSpace       psp;
2759:         PetscInt             point = closure[2 * r];
2760:         PetscInt             pspintdim;
2761:         const PetscInt    ***psymperms = NULL;
2762:         const PetscScalar ***psymflips = NULL;

2764:         if (!point) continue;
2765:         PetscDualSpaceGetPointSubspace(sp, point, &psp);
2766:         if (!psp) continue;
2767:         PetscDualSpaceGetInteriorDimension(psp, &pspintdim);
2768:         if (!pspintdim) continue;
2769:         PetscDualSpaceGetSymmetries(psp, &psymperms, &psymflips);
2770:         symperms[r] = (PetscInt **)(psymperms ? psymperms[0] : NULL);
2771:         symflips[r] = (PetscScalar **)(psymflips ? psymflips[0] : NULL);
2772:       }
2773:       DMPlexRestoreTransitiveClosure(sp->dm, 0, PETSC_TRUE, &closureSize, &closure);
2774:     }
2775:     for (p = 0; p < pEnd; p++)
2776:       if (symperms[p]) break;
2777:     if (p == pEnd) {
2778:       PetscFree(symperms);
2779:       symperms = NULL;
2780:     }
2781:     for (p = 0; p < pEnd; p++)
2782:       if (symflips[p]) break;
2783:     if (p == pEnd) {
2784:       PetscFree(symflips);
2785:       symflips = NULL;
2786:     }
2787:     lag->symperms    = symperms;
2788:     lag->symflips    = symflips;
2789:     lag->symComputed = PETSC_TRUE;
2790:   }
2791:   if (perms) *perms = (const PetscInt ***)lag->symperms;
2792:   if (flips) *flips = (const PetscScalar ***)lag->symflips;
2793:   return 0;
2794: }

2796: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2797: {
2798:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2802:   *continuous = lag->continuous;
2803:   return 0;
2804: }

2806: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2807: {
2808:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2811:   lag->continuous = continuous;
2812:   return 0;
2813: }

2815: /*@
2816:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

2818:   Not Collective

2820:   Input Parameter:
2821: . sp         - the `PetscDualSpace`

2823:   Output Parameter:
2824: . continuous - flag for element continuity

2826:   Level: intermediate

2828: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeSetContinuity()`
2829: @*/
2830: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2831: {
2834:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace, PetscBool *), (sp, continuous));
2835:   return 0;
2836: }

2838: /*@
2839:   PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous

2841:   Logically Collective on sp

2843:   Input Parameters:
2844: + sp         - the `PetscDualSpace`
2845: - continuous - flag for element continuity

2847:   Options Database:
2848: . -petscdualspace_lagrange_continuity <bool> - use a continuous element

2850:   Level: intermediate

2852: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeGetContinuity()`
2853: @*/
2854: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2855: {
2858:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace, PetscBool), (sp, continuous));
2859:   return 0;
2860: }

2862: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
2863: {
2864:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2866:   *tensor = lag->tensorSpace;
2867:   return 0;
2868: }

2870: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
2871: {
2872:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2874:   lag->tensorSpace = tensor;
2875:   return 0;
2876: }

2878: static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Lagrange(PetscDualSpace sp, PetscBool *trimmed)
2879: {
2880:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2882:   *trimmed = lag->trimmed;
2883:   return 0;
2884: }

2886: static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp, PetscBool trimmed)
2887: {
2888:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2890:   lag->trimmed = trimmed;
2891:   return 0;
2892: }

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

2898:   if (nodeType) *nodeType = lag->nodeType;
2899:   if (boundary) *boundary = lag->endNodes;
2900:   if (exponent) *exponent = lag->nodeExponent;
2901:   return 0;
2902: }

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

2909:   lag->nodeType     = nodeType;
2910:   lag->endNodes     = boundary;
2911:   lag->nodeExponent = exponent;
2912:   return 0;
2913: }

2915: static PetscErrorCode PetscDualSpaceLagrangeGetUseMoments_Lagrange(PetscDualSpace sp, PetscBool *useMoments)
2916: {
2917:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2919:   *useMoments = lag->useMoments;
2920:   return 0;
2921: }

2923: static PetscErrorCode PetscDualSpaceLagrangeSetUseMoments_Lagrange(PetscDualSpace sp, PetscBool useMoments)
2924: {
2925:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2927:   lag->useMoments = useMoments;
2928:   return 0;
2929: }

2931: static PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt *momentOrder)
2932: {
2933:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2935:   *momentOrder = lag->momentOrder;
2936:   return 0;
2937: }

2939: static PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt momentOrder)
2940: {
2941:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2943:   lag->momentOrder = momentOrder;
2944:   return 0;
2945: }

2947: /*@
2948:   PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space

2950:   Not collective

2952:   Input Parameter:
2953: . sp - The `PetscDualSpace`

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

2958:   Level: intermediate

2960: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeSetTensor()`, `PetscDualSpaceCreate()`
2961: @*/
2962: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
2963: {
2966:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetTensor_C", (PetscDualSpace, PetscBool *), (sp, tensor));
2967:   return 0;
2968: }

2970: /*@
2971:   PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space

2973:   Not collective

2975:   Input Parameters:
2976: + sp - The `PetscDualSpace`
2977: - tensor - Whether the dual space has tensor layout (vs. simplicial)

2979:   Level: intermediate

2981: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeGetTensor()`, `PetscDualSpaceCreate()`
2982: @*/
2983: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
2984: {
2986:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetTensor_C", (PetscDualSpace, PetscBool), (sp, tensor));
2987:   return 0;
2988: }

2990: /*@
2991:   PetscDualSpaceLagrangeGetTrimmed - Get the trimmed nature of the dual space

2993:   Not collective

2995:   Input Parameter:
2996: . sp - The `PetscDualSpace`

2998:   Output Parameter:
2999: . 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)

3001:   Level: intermediate

3003: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeSetTrimmed()`, `PetscDualSpaceCreate()`
3004: @*/
3005: PetscErrorCode PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp, PetscBool *trimmed)
3006: {
3009:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetTrimmed_C", (PetscDualSpace, PetscBool *), (sp, trimmed));
3010:   return 0;
3011: }

3013: /*@
3014:   PetscDualSpaceLagrangeSetTrimmed - Set the trimmed nature of the dual space

3016:   Not collective

3018:   Input Parameters:
3019: + sp - The `PetscDualSpace`
3020: - 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)

3022:   Level: intermediate

3024: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeGetTrimmed()`, `PetscDualSpaceCreate()`
3025: @*/
3026: PetscErrorCode PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp, PetscBool trimmed)
3027: {
3029:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetTrimmed_C", (PetscDualSpace, PetscBool), (sp, trimmed));
3030:   return 0;
3031: }

3033: /*@
3034:   PetscDualSpaceLagrangeGetNodeType - Get a description of how nodes are laid out for Lagrange polynomials in this
3035:   dual space

3037:   Not collective

3039:   Input Parameter:
3040: . sp - The `PetscDualSpace`

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

3049:   Level: advanced

3051: .seealso: `PetscDualSpace`, `PetscDTNodeType`, `PetscDualSpaceLagrangeSetNodeType()`
3052: @*/
3053: PetscErrorCode PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
3054: {
3059:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetNodeType_C", (PetscDualSpace, PetscDTNodeType *, PetscBool *, PetscReal *), (sp, nodeType, boundary, exponent));
3060:   return 0;
3061: }

3063: /*@
3064:   PetscDualSpaceLagrangeSetNodeType - Set a description of how nodes are laid out for Lagrange polynomials in this
3065:   dual space

3067:   Logically collective

3069:   Input Parameters:
3070: + sp - The `PetscDualSpace`
3071: . nodeType - The type of nodes
3072: . boundary - Whether the node type is one that includes endpoints (if nodeType is `PETSCDTNODES_GAUSSJACOBI`, nodes that
3073:              include the boundary are Gauss-Lobatto-Jacobi nodes)
3074: - exponent - If nodeType is `PETSCDTNODES_GAUSJACOBI`, indicates the exponent used for both ends of the 1D Jacobi weight function
3075:              '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type

3077:   Level: advanced

3079: .seealso: `PetscDualSpace`, `PetscDTNodeType`, `PetscDualSpaceLagrangeGetNodeType()`
3080: @*/
3081: PetscErrorCode PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
3082: {
3084:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetNodeType_C", (PetscDualSpace, PetscDTNodeType, PetscBool, PetscReal), (sp, nodeType, boundary, exponent));
3085:   return 0;
3086: }

3088: /*@
3089:   PetscDualSpaceLagrangeGetUseMoments - Get the flag for using moment functionals

3091:   Not collective

3093:   Input Parameter:
3094: . sp - The `PetscDualSpace`

3096:   Output Parameter:
3097: . useMoments - Moment flag

3099:   Level: advanced

3101: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeSetUseMoments()`
3102: @*/
3103: PetscErrorCode PetscDualSpaceLagrangeGetUseMoments(PetscDualSpace sp, PetscBool *useMoments)
3104: {
3107:   PetscUseMethod(sp, "PetscDualSpaceLagrangeGetUseMoments_C", (PetscDualSpace, PetscBool *), (sp, useMoments));
3108:   return 0;
3109: }

3111: /*@
3112:   PetscDualSpaceLagrangeSetUseMoments - Set the flag for moment functionals

3114:   Logically collective

3116:   Input Parameters:
3117: + sp - The `PetscDualSpace`
3118: - useMoments - The flag for moment functionals

3120:   Level: advanced

3122: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeGetUseMoments()`
3123: @*/
3124: PetscErrorCode PetscDualSpaceLagrangeSetUseMoments(PetscDualSpace sp, PetscBool useMoments)
3125: {
3127:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetUseMoments_C", (PetscDualSpace, PetscBool), (sp, useMoments));
3128:   return 0;
3129: }

3131: /*@
3132:   PetscDualSpaceLagrangeGetMomentOrder - Get the order for moment integration

3134:   Not collective

3136:   Input Parameter:
3137: . sp - The `PetscDualSpace`

3139:   Output Parameter:
3140: . order - Moment integration order

3142:   Level: advanced

3144: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeSetMomentOrder()`
3145: @*/
3146: PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder(PetscDualSpace sp, PetscInt *order)
3147: {
3150:   PetscUseMethod(sp, "PetscDualSpaceLagrangeGetMomentOrder_C", (PetscDualSpace, PetscInt *), (sp, order));
3151:   return 0;
3152: }

3154: /*@
3155:   PetscDualSpaceLagrangeSetMomentOrder - Set the order for moment integration

3157:   Logically collective

3159:   Input Parameters:
3160: + sp - The `PetscDualSpace`
3161: - order - The order for moment integration

3163:   Level: advanced

3165: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeGetMomentOrder()`
3166: @*/
3167: PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder(PetscDualSpace sp, PetscInt order)
3168: {
3170:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetMomentOrder_C", (PetscDualSpace, PetscInt), (sp, order));
3171:   return 0;
3172: }

3174: static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
3175: {
3176:   sp->ops->destroy              = PetscDualSpaceDestroy_Lagrange;
3177:   sp->ops->view                 = PetscDualSpaceView_Lagrange;
3178:   sp->ops->setfromoptions       = PetscDualSpaceSetFromOptions_Lagrange;
3179:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Lagrange;
3180:   sp->ops->setup                = PetscDualSpaceSetUp_Lagrange;
3181:   sp->ops->createheightsubspace = NULL;
3182:   sp->ops->createpointsubspace  = NULL;
3183:   sp->ops->getsymmetries        = PetscDualSpaceGetSymmetries_Lagrange;
3184:   sp->ops->apply                = PetscDualSpaceApplyDefault;
3185:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
3186:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
3187:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
3188:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;
3189:   return 0;
3190: }

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

3195:   Level: intermediate

3197: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
3198: M*/
3199: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
3200: {
3201:   PetscDualSpace_Lag *lag;

3204:   PetscNew(&lag);
3205:   sp->data = lag;

3207:   lag->tensorCell  = PETSC_FALSE;
3208:   lag->tensorSpace = PETSC_FALSE;
3209:   lag->continuous  = PETSC_TRUE;
3210:   lag->numCopies   = PETSC_DEFAULT;
3211:   lag->numNodeSkip = PETSC_DEFAULT;
3212:   lag->nodeType    = PETSCDTNODES_DEFAULT;
3213:   lag->useMoments  = PETSC_FALSE;
3214:   lag->momentOrder = 0;

3216:   PetscDualSpaceInitialize_Lagrange(sp);
3217:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
3218:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
3219:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
3220:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
3221:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Lagrange);
3222:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Lagrange);
3223:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Lagrange);
3224:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Lagrange);
3225:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", PetscDualSpaceLagrangeGetUseMoments_Lagrange);
3226:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", PetscDualSpaceLagrangeSetUseMoments_Lagrange);
3227:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", PetscDualSpaceLagrangeGetMomentOrder_Lagrange);
3228:   PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", PetscDualSpaceLagrangeSetMomentOrder_Lagrange);
3229:   return 0;
3230: }