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: {
  9:   PetscInt         refct;
 10:   PetscDTNodeType  nodeFamily;
 11:   PetscReal        gaussJacobiExp;
 12:   PetscInt         nComputed;
 13:   PetscReal      **nodesets;
 14:   PetscBool        endpoints;
 15: };

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

 23:   PetscNew(&f);
 24:   switch (family) {
 25:   case PETSCDTNODES_GAUSSJACOBI:
 26:   case PETSCDTNODES_EQUISPACED:
 27:     f->nodeFamily = family;
 28:     break;
 29:   default: 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++) {
 59:     PetscFree((*nf)->nodesets[i]);
 60:   }
 61:   PetscFree((*nf)->nodesets);
 62:   PetscFree(*nf);
 63:   *nf = NULL;
 64:   return 0;
 65: }

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

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

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

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

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

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

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

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

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

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

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

332: static PetscErrorCode PetscLagNodeIndicesReference(PetscLagNodeIndices ni)
333: {
334:   if (ni) ni->refct++;
335:   return 0;
336: }

338: static PetscErrorCode PetscLagNodeIndicesDuplicate(PetscLagNodeIndices ni, PetscLagNodeIndices *niNew)
339: {
340:   PetscNew(niNew);
341:   (*niNew)->refct = 1;
342:   (*niNew)->nodeIdxDim = ni->nodeIdxDim;
343:   (*niNew)->nodeVecDim = ni->nodeVecDim;
344:   (*niNew)->nNodes = ni->nNodes;
345:   PetscMalloc1(ni->nNodes * ni->nodeIdxDim, &((*niNew)->nodeIdx));
346:   PetscArraycpy((*niNew)->nodeIdx, ni->nodeIdx, ni->nNodes * ni->nodeIdxDim);
347:   PetscMalloc1(ni->nNodes * ni->nodeVecDim, &((*niNew)->nodeVec));
348:   PetscArraycpy((*niNew)->nodeVec, ni->nodeVec, ni->nNodes * ni->nodeVecDim);
349:   (*niNew)->perm = NULL;
350:   return 0;
351: }

353: static PetscErrorCode PetscLagNodeIndicesDestroy(PetscLagNodeIndices *ni)
354: {
355:   if (!(*ni)) return 0;
356:   if (--(*ni)->refct > 0) {
357:     *ni = NULL;
358:     return 0;
359:   }
360:   PetscFree((*ni)->nodeIdx);
361:   PetscFree((*ni)->nodeVec);
362:   PetscFree((*ni)->perm);
363:   PetscFree(*ni);
364:   *ni = NULL;
365:   return 0;
366: }

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

392:   DMGetDimension(dm, &dim);
393:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
394:   nVerts = vEnd - vStart;
395:   PetscMalloc1(nVerts, &closureOrder);
396:   PetscMalloc1(nVerts, &invClosureOrder);
397:   PetscMalloc1(nVerts, &revlexOrder);
398:   if (sortIdx) { /* bubble sort nodeIdx into revlex order */
399:     PetscInt nodeIdxDim = ni->nodeIdxDim;
400:     PetscInt *idxOrder;

402:     PetscMalloc1(nVerts * nodeIdxDim, &newNodeIdx);
403:     PetscMalloc1(nVerts, &idxOrder);
404:     for (v = 0; v < nVerts; v++) idxOrder[v] = v;
405:     for (v = 0; v < nVerts; v++) {
406:       for (w = v + 1; w < nVerts; w++) {
407:         const PetscInt *iv = &(ni->nodeIdx[idxOrder[v] * nodeIdxDim]);
408:         const PetscInt *iw = &(ni->nodeIdx[idxOrder[w] * nodeIdxDim]);
409:         PetscInt diff = 0;

411:         for (d = nodeIdxDim - 1; d >= 0; d--) if ((diff = (iv[d] - iw[d]))) break;
412:         if (diff > 0) {
413:           PetscInt swap = idxOrder[v];

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

445:       for (d = dim - 1; d >= 0; d--) if ((diff = PetscRealPart(cv[d] - cw[d])) != 0.) break;
446:       if (diff > 0.) {
447:         PetscInt swap = revlexOrder[v];

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

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

477:   PetscNew(&ni);
478:   DMGetDimension(dm, &dim);
479:   ni->nodeIdxDim = dim + 1;
480:   ni->nodeVecDim = 0;
481:   ni->nNodes = dim + 1;
482:   ni->refct = 1;
483:   PetscCalloc1((dim + 1)*(dim + 1), &(ni->nodeIdx));
484:   for (d = 0; d < dim + 1; d++) ni->nodeIdx[d*(dim + 2)] = 1;
485:   PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_FALSE);
486:   *nodeIndices = ni;
487:   return 0;
488: }

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

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

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

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

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

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

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

680:       for (m = 0; m < tracei->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = tracei->nodeIdx[i * tracei->nodeIdxDim + m];
681:       for (m = 0; m < fiberi->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = fiberi->nodeIdx[j * fiberi->nodeIdxDim + m];
682:     }
683:   }

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

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

747: /* simple union of two sets of nodes */
748: static PetscErrorCode PetscLagNodeIndicesMerge(PetscLagNodeIndices niA, PetscLagNodeIndices niB, PetscLagNodeIndices *nodeIndices)
749: {
750:   PetscLagNodeIndices ni;
751:   PetscInt            nodeIdxDim, nodeVecDim, nNodes;

753:   PetscNew(&ni);
754:   ni->nodeIdxDim = nodeIdxDim = niA->nodeIdxDim;
756:   ni->nodeVecDim = nodeVecDim = niA->nodeVecDim;
758:   ni->nNodes = nNodes = niA->nNodes + niB->nNodes;
759:   ni->refct = 1;
760:   PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
761:   PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
762:   PetscArraycpy(ni->nodeIdx, niA->nodeIdx, niA->nNodes * nodeIdxDim);
763:   PetscArraycpy(ni->nodeVec, niA->nodeVec, niA->nNodes * nodeVecDim);
764:   PetscArraycpy(&(ni->nodeIdx[niA->nNodes * nodeIdxDim]), niB->nodeIdx, niB->nNodes * nodeIdxDim);
765:   PetscArraycpy(&(ni->nodeVec[niA->nNodes * nodeVecDim]), niB->nodeVec, niB->nNodes * nodeVecDim);
766:   *nodeIndices = ni;
767:   return 0;
768: }

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

784: PETSCTUPINTCOMPREVLEX(3)
785: PETSCTUPINTCOMPREVLEX(4)
786: PETSCTUPINTCOMPREVLEX(5)
787: PETSCTUPINTCOMPREVLEX(6)
788: PETSCTUPINTCOMPREVLEX(7)

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

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

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

852: static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
853: {
854:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

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

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

862:       for (i = 0; i < lag->numSelfSym; i++) {
863:         PetscFree(allocated[i]);
864:       }
865:       PetscFree(allocated);
866:     }
867:     PetscFree(lag->symperms);
868:   }
869:   if (lag->symflips) {
870:     PetscScalar **selfSyms = lag->symflips[0];

872:     if (selfSyms) {
873:       PetscInt i;
874:       PetscScalar **allocated = &selfSyms[-lag->selfSymOff];

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

903: static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
904: {
905:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

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

911: static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
912: {
913:   PetscBool      iascii;

917:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
918:   if (iascii) PetscDualSpaceLagrangeView_Ascii(sp, viewer);
919:   return 0;
920: }

922: static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
923: {
924:   PetscBool      continuous, tensor, trimmed, flg, flg2, flg3;
925:   PetscDTNodeType nodeType;
926:   PetscReal      nodeExponent;
927:   PetscInt       momentOrder;
928:   PetscBool      nodeEndpoints, useMoments;

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

959: static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace spNew)
960: {
961:   PetscBool           cont, tensor, trimmed, boundary;
962:   PetscDTNodeType     nodeType;
963:   PetscReal           exponent;
964:   PetscDualSpace_Lag *lag    = (PetscDualSpace_Lag *) sp->data;

966:   PetscDualSpaceLagrangeGetContinuity(sp, &cont);
967:   PetscDualSpaceLagrangeSetContinuity(spNew, cont);
968:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
969:   PetscDualSpaceLagrangeSetTensor(spNew, tensor);
970:   PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
971:   PetscDualSpaceLagrangeSetTrimmed(spNew, trimmed);
972:   PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &boundary, &exponent);
973:   PetscDualSpaceLagrangeSetNodeType(spNew, nodeType, boundary, exponent);
974:   if (lag->nodeFamily) {
975:     PetscDualSpace_Lag *lagnew = (PetscDualSpace_Lag *) spNew->data;

977:     Petsc1DNodeFamilyReference(lag->nodeFamily);
978:     lagnew->nodeFamily = lag->nodeFamily;
979:   }
980:   return 0;
981: }

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

990:   PetscDualSpaceDuplicate(sp,bdsp);
991:   PetscDualSpaceSetFormDegree(*bdsp, k);
992:   DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(1, PETSC_TRUE), &K);
993:   PetscDualSpaceSetDM(*bdsp, K);
994:   DMDestroy(&K);
995:   PetscDualSpaceSetOrder(*bdsp, order);
996:   PetscDualSpaceSetNumComponents(*bdsp, Nc);
997:   newlag = (PetscDualSpace_Lag *) (*bdsp)->data;
998:   newlag->interiorOnly = interiorOnly;
999:   PetscDualSpaceSetUp(*bdsp);
1000:   return 0;
1001: }

1003: /* just the points, weights aren't handled */
1004: static PetscErrorCode PetscQuadratureCreateTensor(PetscQuadrature trace, PetscQuadrature fiber, PetscQuadrature *product)
1005: {
1006:   PetscInt         dimTrace, dimFiber;
1007:   PetscInt         numPointsTrace, numPointsFiber;
1008:   PetscInt         dim, numPoints;
1009:   const PetscReal *pointsTrace;
1010:   const PetscReal *pointsFiber;
1011:   PetscReal       *points;
1012:   PetscInt         i, j, k, p;

1014:   PetscQuadratureGetData(trace, &dimTrace, NULL, &numPointsTrace, &pointsTrace, NULL);
1015:   PetscQuadratureGetData(fiber, &dimFiber, NULL, &numPointsFiber, &pointsFiber, NULL);
1016:   dim = dimTrace + dimFiber;
1017:   numPoints = numPointsFiber * numPointsTrace;
1018:   PetscMalloc1(numPoints * dim, &points);
1019:   for (p = 0, j = 0; j < numPointsFiber; j++) {
1020:     for (i = 0; i < numPointsTrace; i++, p++) {
1021:       for (k = 0; k < dimTrace; k++) points[p * dim +            k] = pointsTrace[i * dimTrace + k];
1022:       for (k = 0; k < dimFiber; k++) points[p * dim + dimTrace + k] = pointsFiber[j * dimFiber + k];
1023:     }
1024:   }
1025:   PetscQuadratureCreate(PETSC_COMM_SELF, product);
1026:   PetscQuadratureSetData(*product, dim, 0, numPoints, points, NULL);
1027:   return 0;
1028: }

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

1048:   MatGetSize(trace, &mTrace, &nTrace);
1049:   PetscDTBinomialInt(dimTrace, PetscAbsInt(kTrace), &NkTrace);
1051:   MatGetSize(fiber, &mFiber, &nFiber);
1052:   PetscDTBinomialInt(dimFiber, PetscAbsInt(kFiber), &NkFiber);
1054:   PetscMalloc2(mTrace, &nnzTrace, mFiber, &nnzFiber);
1055:   for (i = 0; i < mTrace; i++) {
1056:     MatGetRow(trace, i, &(nnzTrace[i]), NULL, NULL);
1058:   }
1059:   for (i = 0; i < mFiber; i++) {
1060:     MatGetRow(fiber, i, &(nnzFiber[i]), NULL, NULL);
1062:   }
1063:   dim = dimTrace + dimFiber;
1064:   k = kFiber + kTrace;
1065:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1066:   m = mTrace * mFiber;
1067:   PetscMalloc1(m, &nnz);
1068:   for (l = 0, j = 0; j < mFiber; j++) for (i = 0; i < mTrace; i++, l++) nnz[l] = (nnzTrace[i] / NkTrace) * (nnzFiber[j] / NkFiber) * Nk;
1069:   n = (nTrace / NkTrace) * (nFiber / NkFiber) * Nk;
1070:   MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &prod);
1071:   PetscFree(nnz);
1072:   PetscFree2(nnzTrace,nnzFiber);
1073:   /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1074:   MatSetOption(prod, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1075:   /* compute pullbacks */
1076:   PetscDTBinomialInt(dim, PetscAbsInt(kTrace), &dT);
1077:   PetscDTBinomialInt(dim, PetscAbsInt(kFiber), &dF);
1078:   PetscMalloc4(dimTrace * dim, &projT, dimFiber * dim, &projF, dT * NkTrace, &projTstar, dF * NkFiber, &projFstar);
1079:   PetscArrayzero(projT, dimTrace * dim);
1080:   for (i = 0; i < dimTrace; i++) projT[i * (dim + 1)] = 1.;
1081:   PetscArrayzero(projF, dimFiber * dim);
1082:   for (i = 0; i < dimFiber; i++) projF[i * (dim + 1) + dimTrace] = 1.;
1083:   PetscDTAltVPullbackMatrix(dim, dimTrace, projT, kTrace, projTstar);
1084:   PetscDTAltVPullbackMatrix(dim, dimFiber, projF, kFiber, projFstar);
1085:   PetscMalloc5(dT, &workT, dF, &workF, Nk, &work, Nk, &workstar, Nk, &workS);
1086:   PetscMalloc2(dT, &workT2, dF, &workF2);
1087:   PetscMalloc1(Nk * dT, &wedgeMat);
1088:   sign = (PetscAbsInt(kTrace * kFiber) & 1) ? -1. : 1.;
1089:   for (i = 0, iF = 0; iF < mFiber; iF++) {
1090:     PetscInt           ncolsF, nformsF;
1091:     const PetscInt    *colsF;
1092:     const PetscScalar *valsF;

1094:     MatGetRow(fiber, iF, &ncolsF, &colsF, &valsF);
1095:     nformsF = ncolsF / NkFiber;
1096:     for (iT = 0; iT < mTrace; iT++, i++) {
1097:       PetscInt           ncolsT, nformsT;
1098:       const PetscInt    *colsT;
1099:       const PetscScalar *valsT;

1101:       MatGetRow(trace, iT, &ncolsT, &colsT, &valsT);
1102:       nformsT = ncolsT / NkTrace;
1103:       for (j = 0, jF = 0; jF < nformsF; jF++) {
1104:         PetscInt colF = colsF[jF * NkFiber] / NkFiber;

1106:         for (il = 0; il < dF; il++) {
1107:           PetscReal val = 0.;
1108:           for (jl = 0; jl < NkFiber; jl++) val += projFstar[il * NkFiber + jl] * PetscRealPart(valsF[jF * NkFiber + jl]);
1109:           workF[il] = val;
1110:         }
1111:         if (kFiber < 0) {
1112:           for (il = 0; il < dF; il++) workF2[il] = workF[il];
1113:           PetscDTAltVStar(dim, PetscAbsInt(kFiber), 1, workF2, workF);
1114:         }
1115:         PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kFiber), PetscAbsInt(kTrace), workF, wedgeMat);
1116:         for (jT = 0; jT < nformsT; jT++, j++) {
1117:           PetscInt colT = colsT[jT * NkTrace] / NkTrace;
1118:           PetscInt col = colF * (nTrace / NkTrace) + colT;
1119:           const PetscScalar *vals;

1121:           for (il = 0; il < dT; il++) {
1122:             PetscReal val = 0.;
1123:             for (jl = 0; jl < NkTrace; jl++) val += projTstar[il * NkTrace + jl] * PetscRealPart(valsT[jT * NkTrace + jl]);
1124:             workT[il] = val;
1125:           }
1126:           if (kTrace < 0) {
1127:             for (il = 0; il < dT; il++) workT2[il] = workT[il];
1128:             PetscDTAltVStar(dim, PetscAbsInt(kTrace), 1, workT2, workT);
1129:           }

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

1171: /* Union of quadrature points, with an attempt to identify commont points in the two sets */
1172: static PetscErrorCode PetscQuadraturePointsMerge(PetscQuadrature quadA, PetscQuadrature quadB, PetscQuadrature *quadJoint, PetscInt *aToJoint[], PetscInt *bToJoint[])
1173: {
1174:   PetscInt         dimA, dimB;
1175:   PetscInt         nA, nB, nJoint, i, j, d;
1176:   const PetscReal *pointsA;
1177:   const PetscReal *pointsB;
1178:   PetscReal       *pointsJoint;
1179:   PetscInt        *aToJ, *bToJ;
1180:   PetscQuadrature  qJ;

1182:   PetscQuadratureGetData(quadA, &dimA, NULL, &nA, &pointsA, NULL);
1183:   PetscQuadratureGetData(quadB, &dimB, NULL, &nB, &pointsB, NULL);
1185:   nJoint = nA;
1186:   PetscMalloc1(nA, &aToJ);
1187:   for (i = 0; i < nA; i++) aToJ[i] = i;
1188:   PetscMalloc1(nB, &bToJ);
1189:   for (i = 0; i < nB; i++) {
1190:     for (j = 0; j < nA; j++) {
1191:       bToJ[i] = -1;
1192:       for (d = 0; d < dimA; d++) if (PetscAbsReal(pointsB[i * dimA + d] - pointsA[j * dimA + d]) > PETSC_SMALL) break;
1193:       if (d == dimA) {
1194:         bToJ[i] = j;
1195:         break;
1196:       }
1197:     }
1198:     if (bToJ[i] == -1) {
1199:       bToJ[i] = nJoint++;
1200:     }
1201:   }
1202:   *aToJoint = aToJ;
1203:   *bToJoint = bToJ;
1204:   PetscMalloc1(nJoint * dimA, &pointsJoint);
1205:   PetscArraycpy(pointsJoint, pointsA, nA * dimA);
1206:   for (i = 0; i < nB; i++) {
1207:     if (bToJ[i] >= nA) {
1208:       for (d = 0; d < dimA; d++) pointsJoint[bToJ[i] * dimA + d] = pointsB[i * dimA + d];
1209:     }
1210:   }
1211:   PetscQuadratureCreate(PETSC_COMM_SELF, &qJ);
1212:   PetscQuadratureSetData(qJ, dimA, 0, nJoint, pointsJoint, NULL);
1213:   *quadJoint = qJ;
1214:   return 0;
1215: }

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

1227:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1228:   MatGetSize(matA, &mA, &nA);
1230:   MatGetSize(matB, &mB, &nB);
1232:   m = mA + mB;
1233:   n = numMerged * Nk;
1234:   PetscMalloc1(m, &nnz);
1235:   maxnnz = 0;
1236:   for (i = 0; i < mA; i++) {
1237:     MatGetRow(matA, i, &(nnz[i]), NULL, NULL);
1239:     maxnnz = PetscMax(maxnnz, nnz[i]);
1240:   }
1241:   for (i = 0; i < mB; i++) {
1242:     MatGetRow(matB, i, &(nnz[i+mA]), NULL, NULL);
1244:     maxnnz = PetscMax(maxnnz, nnz[i+mA]);
1245:   }
1246:   MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &M);
1247:   PetscFree(nnz);
1248:   /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1249:   MatSetOption(M, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1250:   PetscMalloc1(maxnnz, &work);
1251:   for (i = 0; i < mA; i++) {
1252:     const PetscInt *cols;
1253:     const PetscScalar *vals;
1254:     PetscInt nCols;
1255:     MatGetRow(matA, i, &nCols, &cols, &vals);
1256:     for (j = 0; j < nCols / Nk; j++) {
1257:       PetscInt newCol = aToMerged[cols[j * Nk] / Nk];
1258:       for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1259:     }
1260:     MatSetValuesBlocked(M, 1, &i, nCols, work, vals, INSERT_VALUES);
1261:     MatRestoreRow(matA, i, &nCols, &cols, &vals);
1262:   }
1263:   for (i = 0; i < mB; i++) {
1264:     const PetscInt *cols;
1265:     const PetscScalar *vals;

1267:     PetscInt row = i + mA;
1268:     PetscInt nCols;
1269:     MatGetRow(matB, i, &nCols, &cols, &vals);
1270:     for (j = 0; j < nCols / Nk; j++) {
1271:       PetscInt newCol = bToMerged[cols[j * Nk] / Nk];
1272:       for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1273:     }
1274:     MatSetValuesBlocked(M, 1, &row, nCols, work, vals, INSERT_VALUES);
1275:     MatRestoreRow(matB, i, &nCols, &cols, &vals);
1276:   }
1277:   PetscFree(work);
1278:   MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
1279:   MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
1280:   *matMerged = M;
1281:   return 0;
1282: }

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

1294:   PetscDualSpaceGetDM(sp,&dm);
1295:   DMGetDimension(dm,&dim);
1296:   DMPlexGetDepth(dm,&depth);
1297:   PetscDualSpaceDuplicate(sp,bdsp);
1298:   PetscDualSpaceSetFormDegree(*bdsp,k);
1299:   if (!K) {
1300:     if (depth == dim) {
1301:       DMPolytopeType ct;

1303:       pointDim = dim - 1;
1304:       DMPlexGetCellType(dm, f, &ct);
1305:       DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K);
1306:     } else if (depth == 1) {
1307:       pointDim = 0;
1308:       DMPlexCreateReferenceCell(PETSC_COMM_SELF, DM_POLYTOPE_POINT, &K);
1309:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
1310:   } else {
1311:     PetscObjectReference((PetscObject)K);
1312:     DMGetDimension(K, &pointDim);
1313:   }
1314:   PetscDualSpaceSetDM(*bdsp, K);
1315:   DMDestroy(&K);
1316:   PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew);
1317:   Ncnew = Nknew * Ncopies;
1318:   PetscDualSpaceSetNumComponents(*bdsp, Ncnew);
1319:   newlag = (PetscDualSpace_Lag *) (*bdsp)->data;
1320:   newlag->interiorOnly = interiorOnly;
1321:   PetscDualSpaceSetUp(*bdsp);
1322:   return 0;
1323: }

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

1341:   PetscDTBinomialInt(dim + sum, dim, &nNodes);
1342:   PetscDTBinomialInt(dim + extraSum, dim, &nExtraNodes);

1344:   PetscMalloc1(dim * nExtraNodes, &extraNodeCoords);
1345:   PetscNew(&ni);
1346:   ni->nodeIdxDim = dim + 1;
1347:   ni->nodeVecDim = Nk;
1348:   ni->nNodes = nNodes * Nk;
1349:   ni->refct = 1;
1350:   PetscMalloc1(nNodes * Nk * (dim + 1), &(ni->nodeIdx));
1351:   PetscMalloc1(nNodes * Nk * Nk, &(ni->nodeVec));
1352:   for (i = 0; i < nNodes; i++) for (j = 0; j < Nk; j++) for (k = 0; k < Nk; k++) ni->nodeVec[(i * Nk + j) * Nk + k] = (j == k) ? 1. : 0.;
1353:   Petsc1DNodeFamilyComputeSimplexNodes(nodeFamily, dim, extraSum, extraNodeCoords);
1354:   if (numNodeSkip) {
1355:     PetscInt k;
1356:     PetscInt *tup;

1358:     PetscMalloc1(dim * nNodes, &nodeCoords);
1359:     PetscMalloc1(dim + 1, &tup);
1360:     for (k = 0; k < nNodes; k++) {
1361:       PetscInt j, c;
1362:       PetscInt index;

1364:       PetscDTIndexToBary(dim + 1, sum, k, tup);
1365:       for (j = 0; j < dim + 1; j++) tup[j] += numNodeSkip;
1366:       for (c = 0; c < Nk; c++) {
1367:         for (j = 0; j < dim + 1; j++) {
1368:           ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1369:         }
1370:       }
1371:       PetscDTBaryToIndex(dim + 1, extraSum, tup, &index);
1372:       for (j = 0; j < dim; j++) nodeCoords[k * dim + j] = extraNodeCoords[index * dim + j];
1373:     }
1374:     PetscFree(tup);
1375:     PetscFree(extraNodeCoords);
1376:   } else {
1377:     PetscInt k;
1378:     PetscInt *tup;

1380:     nodeCoords = extraNodeCoords;
1381:     PetscMalloc1(dim + 1, &tup);
1382:     for (k = 0; k < nNodes; k++) {
1383:       PetscInt j, c;

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

1406:     for (a = aprev; a < anext; a++) {
1407:       MatSetValue(intMat, j, a, (a == j) ? 1. : 0., INSERT_VALUES);
1408:     }
1409:   }
1410:   MatAssemblyBegin(intMat, MAT_FINAL_ASSEMBLY);
1411:   MatAssemblyEnd(intMat, MAT_FINAL_ASSEMBLY);
1412:   *iNodes = intNodes;
1413:   *iMat = intMat;
1414:   *nodeIndices = ni;
1415:   return 0;
1416: }

1418: /* once the nodeIndices have been created for the interior of the reference cell, and for all of the boundary cells,
1419:  * push forward the boundary dofs and concatenate them into the full node indices for the dual space */
1420: static PetscErrorCode PetscDualSpaceLagrangeCreateAllNodeIdx(PetscDualSpace sp)
1421: {
1422:   DM             dm;
1423:   PetscInt       dim, nDofs;
1424:   PetscSection   section;
1425:   PetscInt       pStart, pEnd, p;
1426:   PetscInt       formDegree, Nk;
1427:   PetscInt       nodeIdxDim, spintdim;
1428:   PetscDualSpace_Lag *lag;
1429:   PetscLagNodeIndices ni, verti;

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

1457:     PetscSectionGetDof(section, p, &dof);
1458:     if (!dof) continue;
1459:     plag = (PetscDualSpace_Lag *) psp->data;
1460:     PetscSectionGetOffset(section, p, &off);
1461:     PetscLagNodeIndicesPushForward(dm, verti, p, plag->vertIndices, plag->intNodeIndices, 0, formDegree, &(ni->nodeIdx[off * nodeIdxDim]), &(ni->nodeVec[off * Nk]));
1462:   }
1463:   lag->allNodeIndices = ni;
1464:   return 0;
1465: }

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

1487:   PetscDualSpaceGetDM(sp, &dm);
1488:   DMGetDimension(dm, &dim);
1489:   PetscDualSpaceGetSection(sp, &section);
1490:   PetscSectionGetStorageSize(section, &nDofs);
1491:   DMPlexGetChart(dm, &pStart, &pEnd);
1492:   PetscDualSpaceGetFormDegree(sp, &k);
1493:   PetscDualSpaceGetNumComponents(sp, &Nc);
1494:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1495:   for (p = pStart, nNodes = 0, maxNzforms = 0; p < pEnd; p++) {
1496:     PetscDualSpace  psp;
1497:     DM              pdm;
1498:     PetscInt        pdim, pNk;
1499:     PetscQuadrature intNodes;
1500:     Mat intMat;

1502:     PetscDualSpaceGetPointSubspace(sp, p, &psp);
1503:     if (!psp) continue;
1504:     PetscDualSpaceGetDM(psp, &pdm);
1505:     DMGetDimension(pdm, &pdim);
1506:     if (pdim < PetscAbsInt(k)) continue;
1507:     PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
1508:     PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
1509:     if (intNodes) {
1510:       PetscInt nNodesp;

1512:       PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, NULL, NULL);
1513:       nNodes += nNodesp;
1514:     }
1515:     if (intMat) {
1516:       PetscInt maxNzsp;
1517:       PetscInt maxNzformsp;

1519:       MatSeqAIJGetMaxRowNonzeros(intMat, &maxNzsp);
1521:       maxNzformsp = maxNzsp / pNk;
1522:       maxNzforms = PetscMax(maxNzforms, maxNzformsp);
1523:     }
1524:   }
1525:   MatCreateSeqAIJ(PETSC_COMM_SELF, nDofs, nNodes * Nc, maxNzforms * Nk, NULL, &allMat);
1526:   MatSetOption(allMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);
1527:   PetscMalloc7(dim, &v0, dim, &pv0, dim * dim, &J, dim * dim, &Jinv, Nk * Nk, &L, maxNzforms * Nk, &work, maxNzforms * Nk, &iwork);
1528:   for (j = 0; j < dim; j++) pv0[j] = -1.;
1529:   PetscMalloc1(dim * nNodes, &nodes);
1530:   for (p = pStart, countNodes = 0; p < pEnd; p++) {
1531:     PetscDualSpace  psp;
1532:     PetscQuadrature intNodes;
1533:     DM pdm;
1534:     PetscInt pdim, pNk;
1535:     PetscInt countNodesIn = countNodes;
1536:     PetscReal detJ;
1537:     Mat intMat;

1539:     PetscDualSpaceGetPointSubspace(sp, p, &psp);
1540:     if (!psp) continue;
1541:     PetscDualSpaceGetDM(psp, &pdm);
1542:     DMGetDimension(pdm, &pdim);
1543:     if (pdim < PetscAbsInt(k)) continue;
1544:     PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
1545:     if (intNodes == NULL && intMat == NULL) continue;
1546:     PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
1547:     if (p) {
1548:       DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, Jinv, &detJ);
1549:     } else { /* identity */
1550:       PetscInt i,j;

1552:       for (i = 0; i < dim; i++) for (j = 0; j < dim; j++) J[i * dim + j] = Jinv[i * dim + j] = 0.;
1553:       for (i = 0; i < dim; i++) J[i * dim + i] = Jinv[i * dim + i] = 1.;
1554:       for (i = 0; i < dim; i++) v0[i] = -1.;
1555:     }
1556:     if (pdim != dim) { /* compactify Jacobian */
1557:       PetscInt i, j;

1559:       for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
1560:     }
1561:     PetscDTAltVPullbackMatrix(pdim, dim, J, k, L);
1562:     if (intNodes) { /* push forward quadrature locations by the affine transformation */
1563:       PetscInt nNodesp;
1564:       const PetscReal *nodesp;
1565:       PetscInt j;

1567:       PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, &nodesp, NULL);
1568:       for (j = 0; j < nNodesp; j++, countNodes++) {
1569:         PetscInt d, e;

1571:         for (d = 0; d < dim; d++) {
1572:           nodes[countNodes * dim + d] = v0[d];
1573:           for (e = 0; e < pdim; e++) {
1574:             nodes[countNodes * dim + d] += J[d * pdim + e] * (nodesp[j * pdim + e] - pv0[e]);
1575:           }
1576:         }
1577:       }
1578:     }
1579:     if (intMat) {
1580:       PetscInt nrows;
1581:       PetscInt off;

1583:       PetscSectionGetDof(section, p, &nrows);
1584:       PetscSectionGetOffset(section, p, &off);
1585:       for (j = 0; j < nrows; j++) {
1586:         PetscInt ncols;
1587:         const PetscInt *cols;
1588:         const PetscScalar *vals;
1589:         PetscInt l, d, e;
1590:         PetscInt row = j + off;

1592:         MatGetRow(intMat, j, &ncols, &cols, &vals);
1594:         for (l = 0; l < ncols / pNk; l++) {
1595:           PetscInt blockcol;

1597:           for (d = 0; d < pNk; d++) {
1599:           }
1600:           blockcol = cols[l * pNk] / pNk;
1601:           for (d = 0; d < Nk; d++) {
1602:             iwork[l * Nk + d] = (blockcol + countNodesIn) * Nk + d;
1603:           }
1604:           for (d = 0; d < Nk; d++) work[l * Nk + d] = 0.;
1605:           for (d = 0; d < Nk; d++) {
1606:             for (e = 0; e < pNk; e++) {
1607:               /* "push forward" dof by pulling back a k-form to be evaluated on the point: multiply on the right by L */
1608:               work[l * Nk + d] += vals[l * pNk + e] * L[e * Nk + d];
1609:             }
1610:           }
1611:         }
1612:         MatSetValues(allMat, 1, &row, (ncols / pNk) * Nk, iwork, work, INSERT_VALUES);
1613:         MatRestoreRow(intMat, j, &ncols, &cols, &vals);
1614:       }
1615:     }
1616:   }
1617:   MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);
1618:   MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);
1619:   PetscQuadratureCreate(PETSC_COMM_SELF, &allNodes);
1620:   PetscQuadratureSetData(allNodes, dim, 0, nNodes, nodes, NULL);
1621:   PetscFree7(v0, pv0, J, Jinv, L, work, iwork);
1622:   MatDestroy(&(sp->allMat));
1623:   sp->allMat = allMat;
1624:   PetscQuadratureDestroy(&(sp->allNodes));
1625:   sp->allNodes = allNodes;
1626:   return 0;
1627: }

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

1647:   PetscDualSpaceGetDM(sp, &dm);
1648:   DMGetDimension(dm, &dim);
1649:   PetscDualSpaceGetNumComponents(sp, &Nc);
1650:   PetscDualSpaceGetFormDegree(sp, &k);
1651:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1652:   PetscDualSpaceGetAllData(sp, &allNodes, &allMat);
1653:   nNodes = 0;
1654:   if (allNodes) {
1655:     PetscQuadratureGetData(allNodes, NULL, NULL, &nNodes, &nodes, NULL);
1656:   }
1657:   MatGetSize(allMat, &nDofs, NULL);
1658:   PetscDualSpaceGetSection(sp, &section);
1659:   PetscSectionGetStorageSize(section, &spdim);
1661:   PetscMalloc1(nDofs, &(sp->functional));
1662:   PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments);
1663:   if (useMoments) {
1664:     Mat              allMat;
1665:     PetscInt         momentOrder, i;
1666:     PetscBool        tensor;
1667:     const PetscReal *weights;
1668:     PetscScalar     *array;

1671:     PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder);
1672:     PetscDualSpaceLagrangeGetTensor(sp, &tensor);
1673:     if (!tensor) PetscDTStroudConicalQuadrature(dim, Nc, PetscMax(momentOrder + 1,1), -1.0, 1.0, &(sp->functional[0]));
1674:     else         PetscDTGaussTensorQuadrature(dim, Nc, PetscMax(momentOrder + 1,1), -1.0, 1.0, &(sp->functional[0]));
1675:     /* Need to replace allNodes and allMat */
1676:     PetscObjectReference((PetscObject) sp->functional[0]);
1677:     PetscQuadratureDestroy(&(sp->allNodes));
1678:     sp->allNodes = sp->functional[0];
1679:     PetscQuadratureGetData(sp->allNodes, NULL, NULL, &nNodes, NULL, &weights);
1680:     MatCreateSeqDense(PETSC_COMM_SELF, nDofs, nNodes * Nc, NULL, &allMat);
1681:     MatDenseGetArrayWrite(allMat, &array);
1682:     for (i = 0; i < nNodes * Nc; ++i) array[i] = weights[i];
1683:     MatDenseRestoreArrayWrite(allMat, &array);
1684:     MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);
1685:     MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);
1686:     MatDestroy(&(sp->allMat));
1687:     sp->allMat = allMat;
1688:     return 0;
1689:   }
1690:   for (f = 0; f < nDofs; f++) {
1691:     PetscInt ncols, c;
1692:     const PetscInt *cols;
1693:     const PetscScalar *vals;
1694:     PetscReal *nodesf;
1695:     PetscReal *weightsf;
1696:     PetscInt nNodesf;
1697:     PetscInt countNodes;

1699:     MatGetRow(allMat, f, &ncols, &cols, &vals);
1701:     for (c = 1, nNodesf = 1; c < ncols; c++) {
1702:       if ((cols[c] / Nc) != (cols[c-1] / Nc)) nNodesf++;
1703:     }
1704:     PetscMalloc1(dim * nNodesf, &nodesf);
1705:     PetscMalloc1(Nc * nNodesf, &weightsf);
1706:     for (c = 0, countNodes = 0; c < ncols; c++) {
1707:       if (!c || ((cols[c] / Nc) != (cols[c-1] / Nc))) {
1708:         PetscInt d;

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

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

1734:   MatGetSize(A, &m, &n);
1736:   PetscMalloc1(m * Ncopies, &nnz);
1737:   for (i = 0, maxnnz = 0; i < m; i++) {
1738:     PetscInt innz;
1739:     MatGetRow(A, i, &innz, NULL, NULL);
1741:     for (j = 0; j < Ncopies; j++) nnz[i * Ncopies + j] = innz;
1742:     maxnnz = PetscMax(maxnnz, innz);
1743:   }
1744:   MatCreateSeqAIJ(PETSC_COMM_SELF, m * Ncopies, n * Ncopies, 0, nnz, &Ac);
1745:   MatSetOption(Ac, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1746:   PetscFree(nnz);
1747:   PetscMalloc1(maxnnz, &iwork);
1748:   for (i = 0; i < m; i++) {
1749:     PetscInt innz;
1750:     const PetscInt    *cols;
1751:     const PetscScalar *vals;

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

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

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

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

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

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

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

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

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

1867:     for (c2 = c + 1; c2 < coneSize; c2++) {
1868:       PetscInt  f2 = cone[c2];
1869:       PetscBool isTensorff2;
1870:       PetscInt f2ConeSize;

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

1875:       DMPlexPointIsTensor_Internal_Given(dm, p, f, f2, &isTensorff2);
1876:       if (isTensorff2) {
1877:         if (isTensor) *isTensor = PETSC_TRUE;
1878:         if (endA) *endA = f;
1879:         if (endB) *endB = f2;
1880:         return 0;
1881:       }
1882:     }
1883:   }
1884:   if (isTensor) *isTensor = PETSC_FALSE;
1885:   if (endA) *endA = -1;
1886:   if (endB) *endB = -1;
1887:   return 0;
1888: }

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

1896:   DMPlexIsInterpolated(dm, &interpolated);
1898:   DMPlexPointIsTensor_Internal(dm, p, isTensor, endA, endB);
1899:   return 0;
1900: }

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

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

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

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

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

2021:   /* step 1: sanitize input */
2022:   PetscObjectGetComm((PetscObject) sp, &comm);
2023:   DMGetDimension(dm, &dim);
2024:   PetscObjectTypeCompare((PetscObject)sp, PETSCDUALSPACEBDM, &isbdm);
2025:   if (isbdm) {
2026:     sp->k = -(dim-1); /* form degree of H-div */
2027:     PetscObjectChangeTypeName((PetscObject)sp, PETSCDUALSPACELAGRANGE);
2028:   }
2029:   PetscDualSpaceGetFormDegree(sp, &formDegree);
2031:   PetscDTBinomialInt(dim,PetscAbsInt(formDegree),&Nk);
2032:   if (sp->Nc <= 0 && lag->numCopies > 0) sp->Nc = Nk * lag->numCopies;
2033:   Nc = sp->Nc;
2035:   if (lag->numCopies <= 0) lag->numCopies = Nc / Nk;
2036:   Ncopies = lag->numCopies;
2038:   if (!dim) sp->order = 0;
2039:   order = sp->order;
2040:   uniform = sp->uniform;
2042:   if (lag->trimmed && !formDegree) lag->trimmed = PETSC_FALSE; /* trimmed spaces are the same as full spaces for 0-forms */
2043:   if (lag->nodeType == PETSCDTNODES_DEFAULT) {
2044:     lag->nodeType = PETSCDTNODES_GAUSSJACOBI;
2045:     lag->nodeExponent = 0.;
2046:     /* trimmed spaces don't include corner vertices, so don't use end nodes by default */
2047:     lag->endNodes = lag->trimmed ? PETSC_FALSE : PETSC_TRUE;
2048:   }
2049:   /* If a trimmed space and the user did choose nodes with endpoints, skip them by default */
2050:   if (lag->numNodeSkip < 0) lag->numNodeSkip = (lag->trimmed && lag->endNodes) ? 1 : 0;
2051:   numNodeSkip = lag->numNodeSkip;
2053:   if (lag->trimmed && PetscAbsInt(formDegree) == dim) { /* convert trimmed n-forms to untrimmed of one polynomial order less */
2054:     sp->order--;
2055:     order--;
2056:     lag->trimmed = PETSC_FALSE;
2057:   }
2058:   trimmed = lag->trimmed;
2059:   if (!order || PetscAbsInt(formDegree) == dim) lag->continuous = PETSC_FALSE;
2060:   continuous = lag->continuous;
2061:   DMPlexGetDepth(dm, &depth);
2062:   DMPlexGetChart(dm, &pStart, &pEnd);
2063:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2066:   DMPlexIsInterpolated(dm, &interpolated);
2067:   if (interpolated != DMPLEX_INTERPOLATED_FULL) {
2068:     DMPlexInterpolate(dm, &dmint);
2069:   } else {
2070:     PetscObjectReference((PetscObject)dm);
2071:     dmint = dm;
2072:   }
2073:   tensorCell = PETSC_FALSE;
2074:   if (dim > 1) {
2075:     DMPlexPointIsTensor(dmint, 0, &tensorCell, &tensorf, &tensorf2);
2076:   }
2077:   lag->tensorCell = tensorCell;
2078:   if (dim < 2 || !lag->tensorCell) lag->tensorSpace = PETSC_FALSE;
2079:   tensorSpace = lag->tensorSpace;
2080:   if (!lag->nodeFamily) {
2081:     Petsc1DNodeFamilyCreate(lag->nodeType, lag->nodeExponent, lag->endNodes, &lag->nodeFamily);
2082:   }
2083:   nodeFamily = lag->nodeFamily;

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

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

2101:       DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, NULL, &detJ);
2102:       DMPlexGetCellType(dm, p, &ptype);

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

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

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

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

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

2162:   if (Ncopies > 1) {
2163:     Mat intMatScalar, allMatScalar;
2164:     PetscDualSpace scalarsp;
2165:     PetscDualSpace_Lag *scalarlag;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2567:         /* compare the oriented permutation indices */
2568:         for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break;
2569:         if (d < nodeIdxDim) break;
2570:       }
2572:     }
2573:     groupSize = nEnd - n;
2574:     /* each pushforward dof vector will be expressed in a basis of the unpermuted dofs */
2575:     for (m = n; m < nEnd; m++) nnz[permOrnt[m]] = groupSize;

2577:     maxGroupSize = PetscMax(maxGroupSize, nEnd - n);
2578:     n = nEnd;
2579:   }
2581:   MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes, nNodes, 0, nnz, &A);
2582:   PetscFree(nnz);
2583:   PetscMalloc3(maxGroupSize * nodeVecDim, &V, maxGroupSize * nodeVecDim, &W, nodeVecDim * 2, &work);
2584:   for (n = 0; n < nNodes;) { /* incremented in the loop */
2585:     PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2586:     PetscInt nEnd;
2587:     PetscInt m;
2588:     PetscInt groupSize;
2589:     for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2590:       PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2591:       PetscInt d;

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

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

2619:       PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&bm,&bn,&bnrhs,V,&blda,W,&bldb,work,&blwork, &info));
2621:       /* repack */
2622:       {
2623:         PetscInt i, j;

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

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

2639:           for (d = 0; d < nodeVecDim; d++) {
2640:             W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
2641:           }
2642:         }
2643:         res = 0.;
2644:         for (PetscInt i = 0; i < groupSize; i++) {
2645:           for (PetscInt j = 0; j < nodeVecDim; j++) {
2646:             for (PetscInt k = 0; k < groupSize; k++) {
2647:               W[i * nodeVecDim + j] -= V[i * groupSize + k] * intNodeIndices->nodeVec[perm[n+k] * nodeVecDim + j];
2648:             }
2649:             res += PetscAbsScalar(W[i * nodeVecDim + j]);
2650:           }
2651:         }
2653:       }
2654:     }
2655:     MatSetValues(A, groupSize, &permOrnt[n], groupSize, &perm[n], V, INSERT_VALUES);
2656:     n = nEnd;
2657:   }
2658:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
2659:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
2660:   *symMat = A;
2661:   PetscFree3(V,W,work);
2662:   PetscLagNodeIndicesDestroy(&ni);
2663:   return 0;
2664: }

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

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

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

2685:   PetscDualSpaceGetOrder(sp,&order);
2686:   PetscDualSpaceGetNumComponents(sp,&Nc);
2687:   DMGetDimension(sp->dm,&dim);
2688:   if (!lag->symComputed) { /* store symmetries */
2689:     PetscInt       pStart, pEnd, p;
2690:     PetscInt       numPoints;
2691:     PetscInt       numFaces;
2692:     PetscInt       spintdim;
2693:     PetscInt       ***symperms;
2694:     PetscScalar    ***symflips;

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

2718:       lag->numSelfSym = 2 * numFaces;
2719:       lag->selfSymOff = numFaces;
2720:       PetscCalloc1(2*numFaces,&cellSymperms);
2721:       PetscCalloc1(2*numFaces,&cellSymflips);
2722:       /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
2723:       symperms[0] = &cellSymperms[numFaces];
2724:       symflips[0] = &cellSymflips[numFaces];
2727:       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 */
2728:         Mat symMat;
2729:         PetscInt *perm;
2730:         PetscScalar *flips;
2731:         PetscInt i;

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

2745:           MatGetRow(symMat, i, &ncols, &cols, &vals);
2746:           for (j = 0; j < ncols; j++) {
2747:             if (PetscAbsScalar(vals[j]) > PETSC_SMALL) {
2749:               nz_seen = PETSC_TRUE;
2753:               for (k = 0; k < nCopies; k++) {
2754:                 perm[cols[j] * nCopies + k] = i * nCopies + k;
2755:               }
2756:               if (PetscRealPart(vals[j]) < 0.) {
2757:                 for (k = 0; k < nCopies; k++) {
2758:                   flips[i * nCopies + k] = -1.;
2759:                 }
2760:               } else {
2761:                 for (k = 0; k < nCopies; k++) {
2762:                   flips[i * nCopies + k] = 1.;
2763:                 }
2764:               }
2765:             }
2766:           }
2767:           MatRestoreRow(symMat, i, &ncols, &cols, &vals);
2768:         }
2769:         MatDestroy(&symMat);
2770:         /* if there were no sign flips, keep NULL */
2771:         for (i = 0; i < spintdim; i++) if (flips[i] != 1.) break;
2772:         if (i == spintdim) {
2773:           PetscFree(flips);
2774:           flips = NULL;
2775:         }
2776:         /* if the permutation is identity, keep NULL */
2777:         for (i = 0; i < spintdim; i++) if (perm[i] != i) break;
2778:         if (i == spintdim) {
2779:           PetscFree(perm);
2780:           perm = NULL;
2781:         }
2782:         symperms[0][ornt] = perm;
2783:         symflips[0][ornt] = flips;
2784:       }
2785:       /* if no orientations produced non-identity permutations, keep NULL */
2786:       for (ornt = -numFaces; ornt < numFaces; ornt++) if (symperms[0][ornt]) break;
2787:       if (ornt == numFaces) {
2788:         PetscFree(cellSymperms);
2789:         symperms[0] = NULL;
2790:       }
2791:       /* if no orientations produced sign flips, keep NULL */
2792:       for (ornt = -numFaces; ornt < numFaces; ornt++) if (symflips[0][ornt]) break;
2793:       if (ornt == numFaces) {
2794:         PetscFree(cellSymflips);
2795:         symflips[0] = NULL;
2796:       }
2797:     }
2798:     { /* get the symmetries of closure points */
2799:       PetscInt closureSize = 0;
2800:       PetscInt *closure = NULL;
2801:       PetscInt r;

2803:       DMPlexGetTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);
2804:       for (r = 0; r < closureSize; r++) {
2805:         PetscDualSpace psp;
2806:         PetscInt point = closure[2 * r];
2807:         PetscInt pspintdim;
2808:         const PetscInt ***psymperms = NULL;
2809:         const PetscScalar ***psymflips = NULL;

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

2841: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2842: {
2843:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2847:   *continuous = lag->continuous;
2848:   return 0;
2849: }

2851: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2852: {
2853:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2856:   lag->continuous = continuous;
2857:   return 0;
2858: }

2860: /*@
2861:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

2863:   Not Collective

2865:   Input Parameter:
2866: . sp         - the PetscDualSpace

2868:   Output Parameter:
2869: . continuous - flag for element continuity

2871:   Level: intermediate

2873: .seealso: PetscDualSpaceLagrangeSetContinuity()
2874: @*/
2875: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2876: {
2879:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2880:   return 0;
2881: }

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

2886:   Logically Collective on sp

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

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

2895:   Level: intermediate

2897: .seealso: PetscDualSpaceLagrangeGetContinuity()
2898: @*/
2899: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2900: {
2903:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2904:   return 0;
2905: }

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

2911:   *tensor = lag->tensorSpace;
2912:   return 0;
2913: }

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

2919:   lag->tensorSpace = tensor;
2920:   return 0;
2921: }

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

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

2931: static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp, PetscBool trimmed)
2932: {
2933:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

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

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

2943:   if (nodeType) *nodeType = lag->nodeType;
2944:   if (boundary) *boundary = lag->endNodes;
2945:   if (exponent) *exponent = lag->nodeExponent;
2946:   return 0;
2947: }

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

2954:   lag->nodeType = nodeType;
2955:   lag->endNodes = boundary;
2956:   lag->nodeExponent = exponent;
2957:   return 0;
2958: }

2960: static PetscErrorCode PetscDualSpaceLagrangeGetUseMoments_Lagrange(PetscDualSpace sp, PetscBool *useMoments)
2961: {
2962:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2964:   *useMoments = lag->useMoments;
2965:   return 0;
2966: }

2968: static PetscErrorCode PetscDualSpaceLagrangeSetUseMoments_Lagrange(PetscDualSpace sp, PetscBool useMoments)
2969: {
2970:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2972:   lag->useMoments = useMoments;
2973:   return 0;
2974: }

2976: static PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt *momentOrder)
2977: {
2978:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2980:   *momentOrder = lag->momentOrder;
2981:   return 0;
2982: }

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

2988:   lag->momentOrder = momentOrder;
2989:   return 0;
2990: }

2992: /*@
2993:   PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space

2995:   Not collective

2997:   Input Parameter:
2998: . sp - The PetscDualSpace

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

3003:   Level: intermediate

3005: .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
3006: @*/
3007: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
3008: {
3011:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));
3012:   return 0;
3013: }

3015: /*@
3016:   PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space

3018:   Not collective

3020:   Input Parameters:
3021: + sp - The PetscDualSpace
3022: - tensor - Whether the dual space has tensor layout (vs. simplicial)

3024:   Level: intermediate

3026: .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
3027: @*/
3028: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
3029: {
3031:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));
3032:   return 0;
3033: }

3035: /*@
3036:   PetscDualSpaceLagrangeGetTrimmed - Get the trimmed nature of the dual space

3038:   Not collective

3040:   Input Parameter:
3041: . sp - The PetscDualSpace

3043:   Output Parameter:
3044: . 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)

3046:   Level: intermediate

3048: .seealso: PetscDualSpaceLagrangeSetTrimmed(), PetscDualSpaceCreate()
3049: @*/
3050: PetscErrorCode PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp, PetscBool *trimmed)
3051: {
3054:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTrimmed_C",(PetscDualSpace,PetscBool *),(sp,trimmed));
3055:   return 0;
3056: }

3058: /*@
3059:   PetscDualSpaceLagrangeSetTrimmed - Set the trimmed nature of the dual space

3061:   Not collective

3063:   Input Parameters:
3064: + sp - The PetscDualSpace
3065: - 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)

3067:   Level: intermediate

3069: .seealso: PetscDualSpaceLagrangeGetTrimmed(), PetscDualSpaceCreate()
3070: @*/
3071: PetscErrorCode PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp, PetscBool trimmed)
3072: {
3074:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTrimmed_C",(PetscDualSpace,PetscBool),(sp,trimmed));
3075:   return 0;
3076: }

3078: /*@
3079:   PetscDualSpaceLagrangeGetNodeType - Get a description of how nodes are laid out for Lagrange polynomials in this
3080:   dual space

3082:   Not collective

3084:   Input Parameter:
3085: . sp - The PetscDualSpace

3087:   Output Parameters:
3088: + nodeType - The type of nodes
3089: . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that
3090:              include the boundary are Gauss-Lobatto-Jacobi nodes)
3091: - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function
3092:              '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type

3094:   Level: advanced

3096: .seealso: PetscDTNodeType, PetscDualSpaceLagrangeSetNodeType()
3097: @*/
3098: PetscErrorCode PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
3099: {
3104:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetNodeType_C",(PetscDualSpace,PetscDTNodeType *,PetscBool *,PetscReal *),(sp,nodeType,boundary,exponent));
3105:   return 0;
3106: }

3108: /*@
3109:   PetscDualSpaceLagrangeSetNodeType - Set a description of how nodes are laid out for Lagrange polynomials in this
3110:   dual space

3112:   Logically collective

3114:   Input Parameters:
3115: + sp - The PetscDualSpace
3116: . nodeType - The type of nodes
3117: . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that
3118:              include the boundary are Gauss-Lobatto-Jacobi nodes)
3119: - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function
3120:              '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type

3122:   Level: advanced

3124: .seealso: PetscDTNodeType, PetscDualSpaceLagrangeGetNodeType()
3125: @*/
3126: PetscErrorCode PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
3127: {
3129:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetNodeType_C",(PetscDualSpace,PetscDTNodeType,PetscBool,PetscReal),(sp,nodeType,boundary,exponent));
3130:   return 0;
3131: }

3133: /*@
3134:   PetscDualSpaceLagrangeGetUseMoments - Get the flag for using moment functionals

3136:   Not collective

3138:   Input Parameter:
3139: . sp - The PetscDualSpace

3141:   Output Parameter:
3142: . useMoments - Moment flag

3144:   Level: advanced

3146: .seealso: PetscDualSpaceLagrangeSetUseMoments()
3147: @*/
3148: PetscErrorCode PetscDualSpaceLagrangeGetUseMoments(PetscDualSpace sp, PetscBool *useMoments)
3149: {
3152:   PetscUseMethod(sp,"PetscDualSpaceLagrangeGetUseMoments_C",(PetscDualSpace,PetscBool *),(sp,useMoments));
3153:   return 0;
3154: }

3156: /*@
3157:   PetscDualSpaceLagrangeSetUseMoments - Set the flag for moment functionals

3159:   Logically collective

3161:   Input Parameters:
3162: + sp - The PetscDualSpace
3163: - useMoments - The flag for moment functionals

3165:   Level: advanced

3167: .seealso: PetscDualSpaceLagrangeGetUseMoments()
3168: @*/
3169: PetscErrorCode PetscDualSpaceLagrangeSetUseMoments(PetscDualSpace sp, PetscBool useMoments)
3170: {
3172:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetUseMoments_C",(PetscDualSpace,PetscBool),(sp,useMoments));
3173:   return 0;
3174: }

3176: /*@
3177:   PetscDualSpaceLagrangeGetMomentOrder - Get the order for moment integration

3179:   Not collective

3181:   Input Parameter:
3182: . sp - The PetscDualSpace

3184:   Output Parameter:
3185: . order - Moment integration order

3187:   Level: advanced

3189: .seealso: PetscDualSpaceLagrangeSetMomentOrder()
3190: @*/
3191: PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder(PetscDualSpace sp, PetscInt *order)
3192: {
3195:   PetscUseMethod(sp,"PetscDualSpaceLagrangeGetMomentOrder_C",(PetscDualSpace,PetscInt *),(sp,order));
3196:   return 0;
3197: }

3199: /*@
3200:   PetscDualSpaceLagrangeSetMomentOrder - Set the order for moment integration

3202:   Logically collective

3204:   Input Parameters:
3205: + sp - The PetscDualSpace
3206: - order - The order for moment integration

3208:   Level: advanced

3210: .seealso: PetscDualSpaceLagrangeGetMomentOrder()
3211: @*/
3212: PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder(PetscDualSpace sp, PetscInt order)
3213: {
3215:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetMomentOrder_C",(PetscDualSpace,PetscInt),(sp,order));
3216:   return 0;
3217: }

3219: static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
3220: {
3221:   sp->ops->destroy              = PetscDualSpaceDestroy_Lagrange;
3222:   sp->ops->view                 = PetscDualSpaceView_Lagrange;
3223:   sp->ops->setfromoptions       = PetscDualSpaceSetFromOptions_Lagrange;
3224:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Lagrange;
3225:   sp->ops->setup                = PetscDualSpaceSetUp_Lagrange;
3226:   sp->ops->createheightsubspace = NULL;
3227:   sp->ops->createpointsubspace  = NULL;
3228:   sp->ops->getsymmetries        = PetscDualSpaceGetSymmetries_Lagrange;
3229:   sp->ops->apply                = PetscDualSpaceApplyDefault;
3230:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
3231:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
3232:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
3233:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;
3234:   return 0;
3235: }

3237: /*MC
3238:   PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals

3240:   Level: intermediate

3242: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
3243: M*/
3244: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
3245: {
3246:   PetscDualSpace_Lag *lag;

3249:   PetscNewLog(sp,&lag);
3250:   sp->data = lag;

3252:   lag->tensorCell  = PETSC_FALSE;
3253:   lag->tensorSpace = PETSC_FALSE;
3254:   lag->continuous  = PETSC_TRUE;
3255:   lag->numCopies   = PETSC_DEFAULT;
3256:   lag->numNodeSkip = PETSC_DEFAULT;
3257:   lag->nodeType    = PETSCDTNODES_DEFAULT;
3258:   lag->useMoments  = PETSC_FALSE;
3259:   lag->momentOrder = 0;

3261:   PetscDualSpaceInitialize_Lagrange(sp);
3262:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
3263:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
3264:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
3265:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
3266:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Lagrange);
3267:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Lagrange);
3268:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Lagrange);
3269:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Lagrange);
3270:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetUseMoments_C", PetscDualSpaceLagrangeGetUseMoments_Lagrange);
3271:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetUseMoments_C", PetscDualSpaceLagrangeSetUseMoments_Lagrange);
3272:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetMomentOrder_C", PetscDualSpaceLagrangeGetMomentOrder_Lagrange);
3273:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetMomentOrder_C", PetscDualSpaceLagrangeSetMomentOrder_Lagrange);
3274:   return 0;
3275: }