Actual source code: dspacelagrange.c

petsc-master 2020-05-31
Report Typos and Errors
  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;

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

 44: static PetscErrorCode Petsc1DNodeFamilyReference(Petsc1DNodeFamily nf)
 45: {
 47:   if (nf) nf->refct++;
 48:   return(0);
 49: }

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

 56:   if (!(*nf)) return(0);
 57:   if (--(*nf)->refct > 0) {
 58:     *nf = NULL;
 59:     return(0);
 60:   }
 61:   nc = (*nf)->nComputed;
 62:   for (i = 0; i < nc; i++) {
 63:     PetscFree((*nf)->nodesets[i]);
 64:   }
 65:   PetscFree((*nf)->nodesets);
 66:   PetscFree(*nf);
 67:   *nf = NULL;
 68:   return(0);
 69: }

 71: static PetscErrorCode Petsc1DNodeFamilyGetNodeSets(Petsc1DNodeFamily f, PetscInt degree, PetscReal ***nodesets)
 72: {
 73:   PetscInt       nc;

 77:   nc = f->nComputed;
 78:   if (degree >= nc) {
 79:     PetscInt    i, j;
 80:     PetscReal **new_nodesets;
 81:     PetscReal  *w;

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

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

129:   w = 0.;
130:   if (dim == 1) {
131:     node[0] = nodesets[degree][tup[0]];
132:     node[1] = nodesets[degree][tup[1]];
133:   } else {
134:     for (i = 0; i < dim + 1; i++) node[i] = 0.;
135:     for (i = 0; i < dim + 1; i++) {
136:       PetscReal wi = nodesets[degree][degree-tup[i]];

138:       for (j = 0; j < dim+1; j++) tup[dim+1+j] = tup[j+(j>=i)];
139:       PetscNodeRecursive_Internal(dim-1,degree-tup[i],nodesets,&tup[dim+1],&node[dim+1]);
140:       for (j = 0; j < dim+1; j++) node[j+(j>=i)] += wi * node[dim+1+j];
141:       w += wi;
142:     }
143:     for (i = 0; i < dim+1; i++) node[i] /= w;
144:   }
145:   return(0);
146: }

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

161:   if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have non-negative dimension\n");
162:   if (degree < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have non-negative degree\n");
163:   if (!dim) return(0);
164:   PetscCalloc1(dim+2, &tup);
165:   k = 0;
166:   PetscDTBinomialInt(degree + dim, dim, &npoints);
167:   Petsc1DNodeFamilyGetNodeSets(f, degree, &nodesets);
168:   worksize = ((dim + 2) * (dim + 3)) / 2;
169:   PetscMalloc2(worksize, &nodework, worksize, &tupwork);
170:   /* loop over the tuples of length dim with sum at most degree */
171:   for (k = 0; k < npoints; k++) {
172:     PetscInt i;

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

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

330: /* this is just here so I can access the values in tests/ex1.c outside the library */
331: PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices ni, PetscInt *nodeIdxDim, PetscInt *nodeVecDim, PetscInt *nNodes, const PetscInt *nodeIdx[], const PetscReal *nodeVec[])
332: {
334:   *nodeIdxDim = ni->nodeIdxDim;
335:   *nodeVecDim = ni->nodeVecDim;
336:   *nNodes = ni->nNodes;
337:   *nodeIdx = ni->nodeIdx;
338:   *nodeVec = ni->nodeVec;
339:   return(0);
340: }

342: static PetscErrorCode PetscLagNodeIndicesReference(PetscLagNodeIndices ni)
343: {
345:   if (ni) ni->refct++;
346:   return(0);
347: }

349: static PetscErrorCode PetscLagNodeIndicesDestroy(PetscLagNodeIndices *ni) {

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

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

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;


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

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


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

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

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

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

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

674:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
675:   PetscDTBinomialInt(dimT, PetscAbsInt(kT), &NkT);
676:   PetscDTBinomialInt(dimF, PetscAbsInt(kF), &NkF);
677:   PetscDTBinomialInt(dim, PetscAbsInt(kT), &MkT);
678:   PetscDTBinomialInt(dim, PetscAbsInt(kF), &MkF);
679:   PetscNew(&ni);
680:   ni->nodeIdxDim = nodeIdxDim = tracei->nodeIdxDim + fiberi->nodeIdxDim;
681:   ni->nodeVecDim = Nk;
682:   ni->nNodes = nNodes = tracei->nNodes * fiberi->nNodes;
683:   ni->refct = 1;
684:   PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
685:   /* first concatenate the indices */
686:   for (l = 0, j = 0; j < fiberi->nNodes; j++) {
687:     for (i = 0; i < tracei->nNodes; i++, l++) {
688:       PetscInt m, n = 0;

690:       for (m = 0; m < tracei->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = tracei->nodeIdx[i * tracei->nodeIdxDim + m];
691:       for (m = 0; m < fiberi->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = fiberi->nodeIdx[j * fiberi->nodeIdxDim + m];
692:     }
693:   }

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

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

757: /* simple union of two sets of nodes */
758: static PetscErrorCode PetscLagNodeIndicesMerge(PetscLagNodeIndices niA, PetscLagNodeIndices niB, PetscLagNodeIndices *nodeIndices)
759: {
760:   PetscLagNodeIndices ni;
761:   PetscInt            nodeIdxDim, nodeVecDim, nNodes;
762:   PetscErrorCode      ierr;

765:   PetscNew(&ni);
766:   ni->nodeIdxDim = nodeIdxDim = niA->nodeIdxDim;
767:   if (niB->nodeIdxDim != nodeIdxDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeIdxDim");
768:   ni->nodeVecDim = nodeVecDim = niA->nodeVecDim;
769:   if (niB->nodeVecDim != nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeVecDim");
770:   ni->nNodes = nNodes = niA->nNodes + niB->nNodes;
771:   ni->refct = 1;
772:   PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
773:   PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
774:   PetscArraycpy(ni->nodeIdx, niA->nodeIdx, niA->nNodes * nodeIdxDim);
775:   PetscArraycpy(ni->nodeVec, niA->nodeVec, niA->nNodes * nodeVecDim);
776:   PetscArraycpy(&(ni->nodeIdx[niA->nNodes * nodeIdxDim]), niB->nodeIdx, niB->nNodes * nodeIdxDim);
777:   PetscArraycpy(&(ni->nodeVec[niA->nNodes * nodeVecDim]), niB->nodeVec, niB->nNodes * nodeVecDim);
778:   *nodeIndices = ni;
779:   return(0);
780: }

782: #define PETSCTUPINTCOMPREVLEX(N)                                   \
783: static int PetscTupIntCompRevlex_##N(const void *a, const void *b) \
784: {                                                                  \
785:   const PetscInt *A = (const PetscInt *) a;                        \
786:   const PetscInt *B = (const PetscInt *) b;                        \
787:   int i;                                                           \
788:   PetscInt diff = 0;                                               \
789:   for (i = 0; i < N; i++) {                                        \
790:     diff = A[N - i] - B[N - i];                                    \
791:     if (diff) break;                                               \
792:   }                                                                \
793:   return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1;                    \
794: }

796: PETSCTUPINTCOMPREVLEX(3)
797: PETSCTUPINTCOMPREVLEX(4)
798: PETSCTUPINTCOMPREVLEX(5)
799: PETSCTUPINTCOMPREVLEX(6)
800: PETSCTUPINTCOMPREVLEX(7)

802: static int PetscTupIntCompRevlex_N(const void *a, const void *b)
803: {
804:   const PetscInt *A = (const PetscInt *) a;
805:   const PetscInt *B = (const PetscInt *) b;
806:   int i;
807:   int N = A[0];
808:   PetscInt diff = 0;
809:   for (i = 0; i < N; i++) {
810:     diff = A[N - i] - B[N - i];
811:     if (diff) break;
812:   }
813:   return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1;
814: }

816: /* The nodes are not necessarily in revlex order wrt nodeIdx: get the permutation
817:  * that puts them in that order */
818: static PetscErrorCode PetscLagNodeIndicesGetPermutation(PetscLagNodeIndices ni, PetscInt *perm[])
819: {

823:   if (!(ni->perm)) {
824:     PetscInt *sorter;
825:     PetscInt m = ni->nNodes;
826:     PetscInt nodeIdxDim = ni->nodeIdxDim;
827:     PetscInt i, j, k, l;
828:     PetscInt *prm;
829:     int (*comp) (const void *, const void *);

831:     PetscMalloc1((nodeIdxDim + 2) * m, &sorter);
832:     for (k = 0, l = 0, i = 0; i < m; i++) {
833:       sorter[k++] = nodeIdxDim + 1;
834:       sorter[k++] = i;
835:       for (j = 0; j < nodeIdxDim; j++) sorter[k++] = ni->nodeIdx[l++];
836:     }
837:     switch (nodeIdxDim) {
838:     case 2:
839:       comp = PetscTupIntCompRevlex_3;
840:       break;
841:     case 3:
842:       comp = PetscTupIntCompRevlex_4;
843:       break;
844:     case 4:
845:       comp = PetscTupIntCompRevlex_5;
846:       break;
847:     case 5:
848:       comp = PetscTupIntCompRevlex_6;
849:       break;
850:     case 6:
851:       comp = PetscTupIntCompRevlex_7;
852:       break;
853:     default:
854:       comp = PetscTupIntCompRevlex_N;
855:       break;
856:     }
857:     qsort(sorter, m, (nodeIdxDim + 2) * sizeof(PetscInt), comp);
858:     PetscMalloc1(m, &prm);
859:     for (i = 0; i < m; i++) prm[i] = sorter[(nodeIdxDim + 2) * i + 1];
860:     ni->perm = prm;
861:     PetscFree(sorter);
862:   }
863:   *perm = ni->perm;
864:   return(0);
865: }

867: static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
868: {
869:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
870:   PetscErrorCode      ierr;

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

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

879:       for (i = 0; i < lag->numSelfSym; i++) {
880:         PetscFree(allocated[i]);
881:       }
882:       PetscFree(allocated);
883:     }
884:     PetscFree(lag->symperms);
885:   }
886:   if (lag->symflips) {
887:     PetscScalar **selfSyms = lag->symflips[0];

889:     if (selfSyms) {
890:       PetscInt i;
891:       PetscScalar **allocated = &selfSyms[-lag->selfSymOff];

893:       for (i = 0; i < lag->numSelfSym; i++) {
894:         PetscFree(allocated[i]);
895:       }
896:       PetscFree(allocated);
897:     }
898:     PetscFree(lag->symflips);
899:   }
900:   Petsc1DNodeFamilyDestroy(&(lag->nodeFamily));
901:   PetscLagNodeIndicesDestroy(&(lag->vertIndices));
902:   PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));
903:   PetscLagNodeIndicesDestroy(&(lag->allNodeIndices));
904:   PetscFree(lag);
905:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
906:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
907:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);
908:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);
909:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL);
910:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL);
911:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL);
912:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL);
913:   return(0);
914: }

916: static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
917: {
918:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
919:   PetscErrorCode      ierr;

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

926: static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
927: {
928:   PetscBool      iascii;

934:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
935:   if (iascii) {PetscDualSpaceLagrangeView_Ascii(sp, viewer);}
936:   return(0);
937: }

939: static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
940: {
941:   PetscBool      continuous, tensor, trimmed, flg, flg2, flg3;
942:   PetscDTNodeType nodeType;
943:   PetscReal      nodeExponent;
944:   PetscBool      nodeEndpoints;

948:   PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
949:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
950:   PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
951:   PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &nodeEndpoints, &nodeExponent);
952:   if (nodeType == PETSCDTNODES_DEFAULT) nodeType = PETSCDTNODES_GAUSSJACOBI;
953:   PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
954:   PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
955:   if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
956:   PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetTensor", tensor, &tensor, &flg);
957:   if (flg) {PetscDualSpaceLagrangeSetTensor(sp, tensor);}
958:   PetscOptionsBool("-petscdualspace_lagrange_trimmed", "Flag for trimmed dual space", "PetscDualSpaceLagrangeSetTrimmed", trimmed, &trimmed, &flg);
959:   if (flg) {PetscDualSpaceLagrangeSetTrimmed(sp, trimmed);}
960:   PetscOptionsEnum("-petscdualspace_lagrange_node_type", "Lagrange node location type", "PetscDualSpaceLagrangeSetNodeType", PetscDTNodeTypes, (PetscEnum)nodeType, (PetscEnum *)&nodeType, &flg);
961:   PetscOptionsBool("-petscdualspace_lagrange_node_endpoints", "Flag for nodes that include endpoints", "PetscDualSpaceLagrangeSetNodeType", nodeEndpoints, &nodeEndpoints, &flg2);
962:   flg3 = PETSC_FALSE;
963:   if (nodeType == PETSCDTNODES_GAUSSJACOBI) {
964:     PetscOptionsReal("-petscdualspace_lagrange_node_exponent", "Gauss-Jacobi weight function exponent", "PetscDualSpaceLagrangeSetNodeType", nodeExponent, &nodeExponent, &flg3);
965:   }
966:   if (flg || flg2 || flg3) {PetscDualSpaceLagrangeSetNodeType(sp, nodeType, nodeEndpoints, nodeExponent);}
967:   PetscOptionsTail();
968:   return(0);
969: }

971: static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace spNew)
972: {
973:   PetscBool           cont, tensor, trimmed, boundary;
974:   PetscDTNodeType     nodeType;
975:   PetscReal           exponent;
976:   PetscDualSpace_Lag *lag    = (PetscDualSpace_Lag *) sp->data;
977:   PetscErrorCode      ierr;

980:   PetscDualSpaceLagrangeGetContinuity(sp, &cont);
981:   PetscDualSpaceLagrangeSetContinuity(spNew, cont);
982:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
983:   PetscDualSpaceLagrangeSetTensor(spNew, tensor);
984:   PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
985:   PetscDualSpaceLagrangeSetTrimmed(spNew, trimmed);
986:   PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &boundary, &exponent);
987:   PetscDualSpaceLagrangeSetNodeType(spNew, nodeType, boundary, exponent);
988:   if (lag->nodeFamily) {
989:     PetscDualSpace_Lag *lagnew = (PetscDualSpace_Lag *) spNew->data;

991:     Petsc1DNodeFamilyReference(lag->nodeFamily);
992:     lagnew->nodeFamily = lag->nodeFamily;
993:   }
994:   return(0);
995: }

997: /* for making tensor product spaces: take a dual space and product a segment space that has all the same
998:  * specifications (trimmed, continuous, order, node set), except for the form degree */
999: static PetscErrorCode PetscDualSpaceCreateEdgeSubspace_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt k, PetscInt Nc, PetscBool interiorOnly, PetscDualSpace *bdsp)
1000: {
1001:   DM                 K;
1002:   PetscDualSpace_Lag *newlag;
1003:   PetscErrorCode     ierr;

1006:   PetscDualSpaceDuplicate(sp,bdsp);
1007:   PetscDualSpaceSetFormDegree(*bdsp, k);
1008:   PetscDualSpaceCreateReferenceCell(*bdsp, 1, PETSC_TRUE, &K);
1009:   PetscDualSpaceSetDM(*bdsp, K);
1010:   DMDestroy(&K);
1011:   PetscDualSpaceSetOrder(*bdsp, order);
1012:   PetscDualSpaceSetNumComponents(*bdsp, Nc);
1013:   newlag = (PetscDualSpace_Lag *) (*bdsp)->data;
1014:   newlag->interiorOnly = interiorOnly;
1015:   PetscDualSpaceSetUp(*bdsp);
1016:   return(0);
1017: }

1019: /* just the points, weights aren't handled */
1020: static PetscErrorCode PetscQuadratureCreateTensor(PetscQuadrature trace, PetscQuadrature fiber, PetscQuadrature *product)
1021: {
1022:   PetscInt         dimTrace, dimFiber;
1023:   PetscInt         numPointsTrace, numPointsFiber;
1024:   PetscInt         dim, numPoints;
1025:   const PetscReal *pointsTrace;
1026:   const PetscReal *pointsFiber;
1027:   PetscReal       *points;
1028:   PetscInt         i, j, k, p;
1029:   PetscErrorCode   ierr;

1032:   PetscQuadratureGetData(trace, &dimTrace, NULL, &numPointsTrace, &pointsTrace, NULL);
1033:   PetscQuadratureGetData(fiber, &dimFiber, NULL, &numPointsFiber, &pointsFiber, NULL);
1034:   dim = dimTrace + dimFiber;
1035:   numPoints = numPointsFiber * numPointsTrace;
1036:   PetscMalloc1(numPoints * dim, &points);
1037:   for (p = 0, j = 0; j < numPointsFiber; j++) {
1038:     for (i = 0; i < numPointsTrace; i++, p++) {
1039:       for (k = 0; k < dimTrace; k++) points[p * dim +            k] = pointsTrace[i * dimTrace + k];
1040:       for (k = 0; k < dimFiber; k++) points[p * dim + dimTrace + k] = pointsFiber[j * dimFiber + k];
1041:     }
1042:   }
1043:   PetscQuadratureCreate(PETSC_COMM_SELF, product);
1044:   PetscQuadratureSetData(*product, dim, 0, numPoints, points, NULL);
1045:   return(0);
1046: }

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

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

1114:     MatGetRow(fiber, iF, &ncolsF, &colsF, &valsF);
1115:     nformsF = ncolsF / NkFiber;
1116:     for (iT = 0; iT < mTrace; iT++, i++) {
1117:       PetscInt           ncolsT, nformsT;
1118:       const PetscInt    *colsT;
1119:       const PetscScalar *valsT;

1121:       MatGetRow(trace, iT, &ncolsT, &colsT, &valsT);
1122:       nformsT = ncolsT / NkTrace;
1123:       for (j = 0, jF = 0; jF < nformsF; jF++) {
1124:         PetscInt colF = colsF[jF * NkFiber] / NkFiber;

1126:         for (il = 0; il < dF; il++) {
1127:           PetscReal val = 0.;
1128:           for (jl = 0; jl < NkFiber; jl++) val += projFstar[il * NkFiber + jl] * PetscRealPart(valsF[jF * NkFiber + jl]);
1129:           workF[il] = val;
1130:         }
1131:         if (kFiber < 0) {
1132:           for (il = 0; il < dF; il++) workF2[il] = workF[il];
1133:           PetscDTAltVStar(dim, PetscAbsInt(kFiber), 1, workF2, workF);
1134:         }
1135:         PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kFiber), PetscAbsInt(kTrace), workF, wedgeMat);
1136:         for (jT = 0; jT < nformsT; jT++, j++) {
1137:           PetscInt colT = colsT[jT * NkTrace] / NkTrace;
1138:           PetscInt col = colF * (nTrace / NkTrace) + colT;
1139:           const PetscScalar *vals;

1141:           for (il = 0; il < dT; il++) {
1142:             PetscReal val = 0.;
1143:             for (jl = 0; jl < NkTrace; jl++) val += projTstar[il * NkTrace + jl] * PetscRealPart(valsT[jT * NkTrace + jl]);
1144:             workT[il] = val;
1145:           }
1146:           if (kTrace < 0) {
1147:             for (il = 0; il < dT; il++) workT2[il] = workT[il];
1148:             PetscDTAltVStar(dim, PetscAbsInt(kTrace), 1, workT2, workT);
1149:           }

1151:           for (il = 0; il < Nk; il++) {
1152:             PetscReal val = 0.;
1153:             for (jl = 0; jl < dT; jl++) val += sign * wedgeMat[il * dT + jl] * workT[jl];
1154:             work[il] = val;
1155:           }
1156:           if (k < 0) {
1157:             PetscDTAltVStar(dim, PetscAbsInt(k), -1, work, workstar);
1158: #if defined(PETSC_USE_COMPLEX)
1159:             for (l = 0; l < Nk; l++) workS[l] = workstar[l];
1160:             vals = &workS[0];
1161: #else
1162:             vals = &workstar[0];
1163: #endif
1164:           } else {
1165: #if defined(PETSC_USE_COMPLEX)
1166:             for (l = 0; l < Nk; l++) workS[l] = work[l];
1167:             vals = &workS[0];
1168: #else
1169:             vals = &work[0];
1170: #endif
1171:           }
1172:           for (l = 0; l < Nk; l++) {
1173:             MatSetValue(prod, i, col * Nk + l, vals[l], INSERT_VALUES);
1174:           } /* Nk */
1175:         } /* jT */
1176:       } /* jF */
1177:       MatRestoreRow(trace, iT, &ncolsT, &colsT, &valsT);
1178:     } /* iT */
1179:     MatRestoreRow(fiber, iF, &ncolsF, &colsF, &valsF);
1180:   } /* iF */
1181:   PetscFree(wedgeMat);
1182:   PetscFree4(projT, projF, projTstar, projFstar);
1183:   PetscFree2(workT2, workF2);
1184:   PetscFree5(workT, workF, work, workstar, workS);
1185:   MatAssemblyBegin(prod, MAT_FINAL_ASSEMBLY);
1186:   MatAssemblyEnd(prod, MAT_FINAL_ASSEMBLY);
1187:   *product = prod;
1188:   return(0);
1189: }

1191: /* Union of quadrature points, with an attempt to identify commont points in the two sets */
1192: static PetscErrorCode PetscQuadraturePointsMerge(PetscQuadrature quadA, PetscQuadrature quadB, PetscQuadrature *quadJoint, PetscInt *aToJoint[], PetscInt *bToJoint[])
1193: {
1194:   PetscInt         dimA, dimB;
1195:   PetscInt         nA, nB, nJoint, i, j, d;
1196:   const PetscReal *pointsA;
1197:   const PetscReal *pointsB;
1198:   PetscReal       *pointsJoint;
1199:   PetscInt        *aToJ, *bToJ;
1200:   PetscQuadrature  qJ;
1201:   PetscErrorCode   ierr;

1204:   PetscQuadratureGetData(quadA, &dimA, NULL, &nA, &pointsA, NULL);
1205:   PetscQuadratureGetData(quadB, &dimB, NULL, &nB, &pointsB, NULL);
1206:   if (dimA != dimB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Quadrature points must be in the same dimension");
1207:   nJoint = nA;
1208:   PetscMalloc1(nA, &aToJ);
1209:   for (i = 0; i < nA; i++) aToJ[i] = i;
1210:   PetscMalloc1(nB, &bToJ);
1211:   for (i = 0; i < nB; i++) {
1212:     for (j = 0; j < nA; j++) {
1213:       bToJ[i] = -1;
1214:       for (d = 0; d < dimA; d++) if (PetscAbsReal(pointsB[i * dimA + d] - pointsA[j * dimA + d]) > PETSC_SMALL) break;
1215:       if (d == dimA) {
1216:         bToJ[i] = j;
1217:         break;
1218:       }
1219:     }
1220:     if (bToJ[i] == -1) {
1221:       bToJ[i] = nJoint++;
1222:     }
1223:   }
1224:   *aToJoint = aToJ;
1225:   *bToJoint = bToJ;
1226:   PetscMalloc1(nJoint * dimA, &pointsJoint);
1227:   PetscArraycpy(pointsJoint, pointsA, nA * dimA);
1228:   for (i = 0; i < nB; i++) {
1229:     if (bToJ[i] >= nA) {
1230:       for (d = 0; d < dimA; d++) pointsJoint[bToJ[i] * dimA + d] = pointsB[i * dimA + d];
1231:     }
1232:   }
1233:   PetscQuadratureCreate(PETSC_COMM_SELF, &qJ);
1234:   PetscQuadratureSetData(qJ, dimA, 0, nJoint, pointsJoint, NULL);
1235:   *quadJoint = qJ;
1236:   return(0);
1237: }

1239: /* Matrices matA and matB are both quadrature -> dof matrices: produce a matrix that is joint quadrature -> union of
1240:  * dofs, where the joint quadrature was produced by PetscQuadraturePointsMerge */
1241: static PetscErrorCode MatricesMerge(Mat matA, Mat matB, PetscInt dim, PetscInt k, PetscInt numMerged, const PetscInt aToMerged[], const PetscInt bToMerged[], Mat *matMerged)
1242: {
1243:   PetscInt m, n, mA, nA, mB, nB, Nk, i, j, l;
1244:   Mat      M;
1245:   PetscInt *nnz;
1246:   PetscInt maxnnz;
1247:   PetscInt *work;

1251:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1252:   MatGetSize(matA, &mA, &nA);
1253:   if (nA % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matA column space not a multiple of k-form size");
1254:   MatGetSize(matB, &mB, &nB);
1255:   if (nB % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matB column space not a multiple of k-form size");
1256:   m = mA + mB;
1257:   n = numMerged * Nk;
1258:   PetscMalloc1(m, &nnz);
1259:   maxnnz = 0;
1260:   for (i = 0; i < mA; i++) {
1261:     MatGetRow(matA, i, &(nnz[i]), NULL, NULL);
1262:     if (nnz[i] % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matA are not in k-form size blocks");
1263:     maxnnz = PetscMax(maxnnz, nnz[i]);
1264:   }
1265:   for (i = 0; i < mB; i++) {
1266:     MatGetRow(matB, i, &(nnz[i+mA]), NULL, NULL);
1267:     if (nnz[i+mA] % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matB are not in k-form size blocks");
1268:     maxnnz = PetscMax(maxnnz, nnz[i+mA]);
1269:   }
1270:   MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &M);
1271:   PetscFree(nnz);
1272:   /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1273:   MatSetOption(M, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1274:   PetscMalloc1(maxnnz, &work);
1275:   for (i = 0; i < mA; i++) {
1276:     const PetscInt *cols;
1277:     const PetscScalar *vals;
1278:     PetscInt nCols;
1279:     MatGetRow(matA, i, &nCols, &cols, &vals);
1280:     for (j = 0; j < nCols / Nk; j++) {
1281:       PetscInt newCol = aToMerged[cols[j * Nk] / Nk];
1282:       for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1283:     }
1284:     MatSetValuesBlocked(M, 1, &i, nCols, work, vals, INSERT_VALUES);
1285:     MatRestoreRow(matA, i, &nCols, &cols, &vals);
1286:   }
1287:   for (i = 0; i < mB; i++) {
1288:     const PetscInt *cols;
1289:     const PetscScalar *vals;

1291:     PetscInt row = i + mA;
1292:     PetscInt nCols;
1293:     MatGetRow(matB, i, &nCols, &cols, &vals);
1294:     for (j = 0; j < nCols / Nk; j++) {
1295:       PetscInt newCol = bToMerged[cols[j * Nk] / Nk];
1296:       for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1297:     }
1298:     MatSetValuesBlocked(M, 1, &row, nCols, work, vals, INSERT_VALUES);
1299:     MatRestoreRow(matB, i, &nCols, &cols, &vals);
1300:   }
1301:   PetscFree(work);
1302:   MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
1303:   MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
1304:   *matMerged = M;
1305:   return(0);
1306: }

1308: /* Take a dual space and product a segment space that has all the same specifications (trimmed, continuous, order,
1309:  * node set), except for the form degree.  For computing boundary dofs and for making tensor product spaces */
1310: static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Lagrange(PetscDualSpace sp, DM K, PetscInt f, PetscInt k, PetscInt Ncopies, PetscBool interiorOnly, PetscDualSpace *bdsp)
1311: {
1312:   PetscInt           Nknew, Ncnew;
1313:   PetscInt           dim, pointDim = -1;
1314:   PetscInt           depth;
1315:   DM                 dm;
1316:   PetscDualSpace_Lag *newlag;
1317:   PetscErrorCode     ierr;

1320:   PetscDualSpaceGetDM(sp,&dm);
1321:   DMGetDimension(dm,&dim);
1322:   DMPlexGetDepth(dm,&depth);
1323:   PetscDualSpaceDuplicate(sp,bdsp);
1324:   PetscDualSpaceSetFormDegree(*bdsp,k);
1325:   if (!K) {
1326:     PetscBool isSimplex;


1329:     if (depth == dim) {
1330:       PetscInt coneSize;

1332:       pointDim = dim - 1;
1333:       DMPlexGetConeSize(dm,f,&coneSize);
1334:       isSimplex = (PetscBool) (coneSize == dim);
1335:       PetscDualSpaceCreateReferenceCell(*bdsp, dim-1, isSimplex, &K);
1336:     } else if (depth == 1) {
1337:       pointDim = 0;
1338:       PetscDualSpaceCreateReferenceCell(*bdsp, 0, PETSC_TRUE, &K);
1339:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
1340:   } else {
1341:     PetscObjectReference((PetscObject)K);
1342:     DMGetDimension(K, &pointDim);
1343:   }
1344:   PetscDualSpaceSetDM(*bdsp, K);
1345:   DMDestroy(&K);
1346:   PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew);
1347:   Ncnew = Nknew * Ncopies;
1348:   PetscDualSpaceSetNumComponents(*bdsp, Ncnew);
1349:   newlag = (PetscDualSpace_Lag *) (*bdsp)->data;
1350:   newlag->interiorOnly = interiorOnly;
1351:   PetscDualSpaceSetUp(*bdsp);
1352:   return(0);
1353: }

1355: /* Construct simplex nodes from a nodefamily, add Nk dof vectors of length Nk at each node.
1356:  * Return the (quadrature, matrix) form of the dofs and the nodeIndices form as well.
1357:  *
1358:  * Sometimes we want a set of nodes to be contained in the interior of the element,
1359:  * even when the node scheme puts nodes on the boundaries.  numNodeSkip tells
1360:  * the routine how many "layers" of nodes need to be skipped.
1361:  * */
1362: static PetscErrorCode PetscDualSpaceLagrangeCreateSimplexNodeMat(Petsc1DNodeFamily nodeFamily, PetscInt dim, PetscInt sum, PetscInt Nk, PetscInt numNodeSkip, PetscQuadrature *iNodes, Mat *iMat, PetscLagNodeIndices *nodeIndices)
1363: {
1364:   PetscReal *extraNodeCoords, *nodeCoords;
1365:   PetscInt nNodes, nExtraNodes;
1366:   PetscInt i, j, k, extraSum = sum + numNodeSkip * (1 + dim);
1367:   PetscQuadrature intNodes;
1368:   Mat intMat;
1369:   PetscLagNodeIndices ni;

1373:   PetscDTBinomialInt(dim + sum, dim, &nNodes);
1374:   PetscDTBinomialInt(dim + extraSum, dim, &nExtraNodes);

1376:   PetscMalloc1(dim * nExtraNodes, &extraNodeCoords);
1377:   PetscNew(&ni);
1378:   ni->nodeIdxDim = dim + 1;
1379:   ni->nodeVecDim = Nk;
1380:   ni->nNodes = nNodes * Nk;
1381:   ni->refct = 1;
1382:   PetscMalloc1(nNodes * Nk * (dim + 1), &(ni->nodeIdx));
1383:   PetscMalloc1(nNodes * Nk * Nk, &(ni->nodeVec));
1384:   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.;
1385:   Petsc1DNodeFamilyComputeSimplexNodes(nodeFamily, dim, extraSum, extraNodeCoords);
1386:   if (numNodeSkip) {
1387:     PetscInt k;
1388:     PetscInt *tup;

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

1396:       PetscDTIndexToBary(dim + 1, sum, k, tup);
1397:       for (j = 0; j < dim + 1; j++) tup[j] += numNodeSkip;
1398:       for (c = 0; c < Nk; c++) {
1399:         for (j = 0; j < dim + 1; j++) {
1400:           ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1401:         }
1402:       }
1403:       PetscDTBaryToIndex(dim + 1, extraSum, tup, &index);
1404:       for (j = 0; j < dim; j++) nodeCoords[k * dim + j] = extraNodeCoords[index * dim + j];
1405:     }
1406:     PetscFree(tup);
1407:     PetscFree(extraNodeCoords);
1408:   } else {
1409:     PetscInt k;
1410:     PetscInt *tup;

1412:     nodeCoords = extraNodeCoords;
1413:     PetscMalloc1(dim + 1, &tup);
1414:     for (k = 0; k < nNodes; k++) {
1415:       PetscInt j, c;

1417:       PetscDTIndexToBary(dim + 1, sum, k, tup);
1418:       for (c = 0; c < Nk; c++) {
1419:         for (j = 0; j < dim + 1; j++) {
1420:           /* barycentric indices can have zeros, but we don't want to push forward zeros because it makes it harder to
1421:            * determine which nodes correspond to which under symmetries, so we increase by 1.  This is fine
1422:            * because the nodeIdx coordinates don't have any meaning other than helping to identify symmetries */
1423:           ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1424:         }
1425:       }
1426:     }
1427:     PetscFree(tup);
1428:   }
1429:   PetscQuadratureCreate(PETSC_COMM_SELF, &intNodes);
1430:   PetscQuadratureSetData(intNodes, dim, 0, nNodes, nodeCoords, NULL);
1431:   MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes * Nk, nNodes * Nk, Nk, NULL, &intMat);
1432:   MatSetOption(intMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);
1433:   for (j = 0; j < nNodes * Nk; j++) {
1434:     PetscInt rem = j % Nk;
1435:     PetscInt a, aprev = j - rem;
1436:     PetscInt anext = aprev + Nk;

1438:     for (a = aprev; a < anext; a++) {
1439:       MatSetValue(intMat, j, a, (a == j) ? 1. : 0., INSERT_VALUES);
1440:     }
1441:   }
1442:   MatAssemblyBegin(intMat, MAT_FINAL_ASSEMBLY);
1443:   MatAssemblyEnd(intMat, MAT_FINAL_ASSEMBLY);
1444:   *iNodes = intNodes;
1445:   *iMat = intMat;
1446:   *nodeIndices = ni;
1447:   return(0);
1448: }

1450: /* once the nodeIndices have been created for the interior of the reference cell, and for all of the boundary cells,
1451:  * push forward the boudary dofs and concatenate them into the full node indices for the dual space */
1452: static PetscErrorCode PetscDualSpaceLagrangeCreateAllNodeIdx(PetscDualSpace sp)
1453: {
1454:   DM             dm;
1455:   PetscInt       dim, nDofs;
1456:   PetscSection   section;
1457:   PetscInt       pStart, pEnd, p;
1458:   PetscInt       formDegree, Nk;
1459:   PetscInt       nodeIdxDim, spintdim;
1460:   PetscDualSpace_Lag *lag;
1461:   PetscLagNodeIndices ni, verti;

1465:   lag = (PetscDualSpace_Lag *) sp->data;
1466:   verti = lag->vertIndices;
1467:   PetscDualSpaceGetDM(sp, &dm);
1468:   DMGetDimension(dm, &dim);
1469:   PetscDualSpaceGetFormDegree(sp, &formDegree);
1470:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
1471:   PetscDualSpaceGetSection(sp, &section);
1472:   PetscSectionGetStorageSize(section, &nDofs);
1473:   PetscNew(&ni);
1474:   ni->nodeIdxDim = nodeIdxDim = verti->nodeIdxDim;
1475:   ni->nodeVecDim = Nk;
1476:   ni->nNodes = nDofs;
1477:   ni->refct = 1;
1478:   PetscMalloc1(nodeIdxDim * nDofs, &(ni->nodeIdx));
1479:   PetscMalloc1(Nk * nDofs, &(ni->nodeVec));
1480:   DMPlexGetChart(dm, &pStart, &pEnd);
1481:   PetscSectionGetDof(section, 0, &spintdim);
1482:   if (spintdim) {
1483:     PetscArraycpy(ni->nodeIdx, lag->intNodeIndices->nodeIdx, spintdim * nodeIdxDim);
1484:     PetscArraycpy(ni->nodeVec, lag->intNodeIndices->nodeVec, spintdim * Nk);
1485:   }
1486:   for (p = pStart + 1; p < pEnd; p++) {
1487:     PetscDualSpace psp = sp->pointSpaces[p];
1488:     PetscDualSpace_Lag *plag;
1489:     PetscInt dof, off;

1491:     PetscSectionGetDof(section, p, &dof);
1492:     if (!dof) continue;
1493:     plag = (PetscDualSpace_Lag *) psp->data;
1494:     PetscSectionGetOffset(section, p, &off);
1495:     PetscLagNodeIndicesPushForward(dm, verti, p, plag->vertIndices, plag->intNodeIndices, 0, formDegree, &(ni->nodeIdx[off * nodeIdxDim]), &(ni->nodeVec[off * Nk]));
1496:   }
1497:   lag->allNodeIndices = ni;
1498:   return(0);
1499: }

1501: /* once the (quadrature, Matrix) forms of the dofs have been created for the interior of the
1502:  * reference cell and for the boundary cells, jk
1503:  * push forward the boundary data and concatenate them into the full (quadrature, matrix) data
1504:  * for the dual space */
1505: static PetscErrorCode PetscDualSpaceCreateAllDataFromInteriorData(PetscDualSpace sp)
1506: {
1507:   DM               dm;
1508:   PetscSection     section;
1509:   PetscInt         pStart, pEnd, p, k, Nk, dim, Nc;
1510:   PetscInt         nNodes;
1511:   PetscInt         countNodes;
1512:   Mat              allMat;
1513:   PetscQuadrature  allNodes;
1514:   PetscInt         nDofs;
1515:   PetscInt         maxNzforms, j;
1516:   PetscScalar      *work;
1517:   PetscReal        *L, *J, *Jinv, *v0, *pv0;
1518:   PetscInt         *iwork;
1519:   PetscReal        *nodes;
1520:   PetscErrorCode   ierr;

1523:   PetscDualSpaceGetDM(sp, &dm);
1524:   DMGetDimension(dm, &dim);
1525:   PetscDualSpaceGetSection(sp, &section);
1526:   PetscSectionGetStorageSize(section, &nDofs);
1527:   DMPlexGetChart(dm, &pStart, &pEnd);
1528:   PetscDualSpaceGetFormDegree(sp, &k);
1529:   PetscDualSpaceGetNumComponents(sp, &Nc);
1530:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1531:   for (p = pStart, nNodes = 0, maxNzforms = 0; p < pEnd; p++) {
1532:     PetscDualSpace  psp;
1533:     DM              pdm;
1534:     PetscInt        pdim, pNk;
1535:     PetscQuadrature intNodes;
1536:     Mat intMat;

1538:     PetscDualSpaceGetPointSubspace(sp, p, &psp);
1539:     if (!psp) continue;
1540:     PetscDualSpaceGetDM(psp, &pdm);
1541:     DMGetDimension(pdm, &pdim);
1542:     if (pdim < PetscAbsInt(k)) continue;
1543:     PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
1544:     PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
1545:     if (intNodes) {
1546:       PetscInt nNodesp;

1548:       PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, NULL, NULL);
1549:       nNodes += nNodesp;
1550:     }
1551:     if (intMat) {
1552:       PetscInt maxNzsp;
1553:       PetscInt maxNzformsp;

1555:       MatSeqAIJGetMaxRowNonzeros(intMat, &maxNzsp);
1556:       if (maxNzsp % pNk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms");
1557:       maxNzformsp = maxNzsp / pNk;
1558:       maxNzforms = PetscMax(maxNzforms, maxNzformsp);
1559:     }
1560:   }
1561:   MatCreateSeqAIJ(PETSC_COMM_SELF, nDofs, nNodes * Nc, maxNzforms * Nk, NULL, &allMat);
1562:   MatSetOption(allMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);
1563:   PetscMalloc7(dim, &v0, dim, &pv0, dim * dim, &J, dim * dim, &Jinv, Nk * Nk, &L, maxNzforms * Nk, &work, maxNzforms * Nk, &iwork);
1564:   for (j = 0; j < dim; j++) pv0[j] = -1.;
1565:   PetscMalloc1(dim * nNodes, &nodes);
1566:   for (p = pStart, countNodes = 0; p < pEnd; p++) {
1567:     PetscDualSpace  psp;
1568:     PetscQuadrature intNodes;
1569:     DM pdm;
1570:     PetscInt pdim, pNk;
1571:     PetscInt countNodesIn = countNodes;
1572:     PetscReal detJ;
1573:     Mat intMat;

1575:     PetscDualSpaceGetPointSubspace(sp, p, &psp);
1576:     if (!psp) continue;
1577:     PetscDualSpaceGetDM(psp, &pdm);
1578:     DMGetDimension(pdm, &pdim);
1579:     if (pdim < PetscAbsInt(k)) continue;
1580:     PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
1581:     if (intNodes == NULL && intMat == NULL) continue;
1582:     PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
1583:     if (p) {
1584:       DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, Jinv, &detJ);
1585:     } else { /* identity */
1586:       PetscInt i,j;

1588:       for (i = 0; i < dim; i++) for (j = 0; j < dim; j++) J[i * dim + j] = Jinv[i * dim + j] = 0.;
1589:       for (i = 0; i < dim; i++) J[i * dim + i] = Jinv[i * dim + i] = 1.;
1590:       for (i = 0; i < dim; i++) v0[i] = -1.;
1591:     }
1592:     if (pdim != dim) { /* compactify Jacobian */
1593:       PetscInt i, j;

1595:       for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
1596:     }
1597:     PetscDTAltVPullbackMatrix(pdim, dim, J, k, L);
1598:     if (intNodes) { /* push forward quadrature locations by the affine transformation */
1599:       PetscInt nNodesp;
1600:       const PetscReal *nodesp;
1601:       PetscInt j;

1603:       PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, &nodesp, NULL);
1604:       for (j = 0; j < nNodesp; j++, countNodes++) {
1605:         PetscInt d, e;

1607:         for (d = 0; d < dim; d++) {
1608:           nodes[countNodes * dim + d] = v0[d];
1609:           for (e = 0; e < pdim; e++) {
1610:             nodes[countNodes * dim + d] += J[d * pdim + e] * (nodesp[j * pdim + e] - pv0[e]);
1611:           }
1612:         }
1613:       }
1614:     }
1615:     if (intMat) {
1616:       PetscInt nrows;
1617:       PetscInt off;

1619:       PetscSectionGetDof(section, p, &nrows);
1620:       PetscSectionGetOffset(section, p, &off);
1621:       for (j = 0; j < nrows; j++) {
1622:         PetscInt ncols;
1623:         const PetscInt *cols;
1624:         const PetscScalar *vals;
1625:         PetscInt l, d, e;
1626:         PetscInt row = j + off;

1628:         MatGetRow(intMat, j, &ncols, &cols, &vals);
1629:         if (ncols % pNk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms");
1630:         for (l = 0; l < ncols / pNk; l++) {
1631:           PetscInt blockcol;

1633:           for (d = 0; d < pNk; d++) {
1634:             if ((cols[l * pNk + d] % pNk) != d) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms");
1635:           }
1636:           blockcol = cols[l * pNk] / pNk;
1637:           for (d = 0; d < Nk; d++) {
1638:             iwork[l * Nk + d] = (blockcol + countNodesIn) * Nk + d;
1639:           }
1640:           for (d = 0; d < Nk; d++) work[l * Nk + d] = 0.;
1641:           for (d = 0; d < Nk; d++) {
1642:             for (e = 0; e < pNk; e++) {
1643:               /* "push forward" dof by pulling back a k-form to be evaluated on the point: multiply on the right by L */
1644:               work[l * Nk + d] += vals[l * pNk + e] * L[e * pNk + d];
1645:             }
1646:           }
1647:         }
1648:         MatSetValues(allMat, 1, &row, (ncols / pNk) * Nk, iwork, work, INSERT_VALUES);
1649:         MatRestoreRow(intMat, j, &ncols, &cols, &vals);
1650:       }
1651:     }
1652:   }
1653:   MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);
1654:   MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);
1655:   PetscQuadratureCreate(PETSC_COMM_SELF, &allNodes);
1656:   PetscQuadratureSetData(allNodes, dim, 0, nNodes, nodes, NULL);
1657:   PetscFree7(v0, pv0, J, Jinv, L, work, iwork);
1658:   MatDestroy(&(sp->allMat));
1659:   sp->allMat = allMat;
1660:   PetscQuadratureDestroy(&(sp->allNodes));
1661:   sp->allNodes = allNodes;
1662:   return(0);
1663: }

1665: /* rather than trying to get all data from the functionals, we create
1666:  * the functionals from rows of the quadrature -> dof matrix.
1667:  *
1668:  * Ideally most of the uses of PetscDualSpace in PetscFE will switch
1669:  * to using intMat and allMat, so that the individual functionals
1670:  * don't need to be constructed at all */
1671: static PetscErrorCode PetscDualSpaceComputeFunctionalsFromAllData(PetscDualSpace sp)
1672: {
1673:   PetscQuadrature allNodes;
1674:   Mat             allMat;
1675:   PetscInt        nDofs;
1676:   PetscInt        dim, k, Nk, Nc, f;
1677:   DM              dm;
1678:   PetscInt        nNodes, spdim;
1679:   const PetscReal *nodes = NULL;
1680:   PetscSection    section;
1681:   PetscErrorCode  ierr;

1684:   PetscDualSpaceGetDM(sp, &dm);
1685:   DMGetDimension(dm, &dim);
1686:   PetscDualSpaceGetNumComponents(sp, &Nc);
1687:   PetscDualSpaceGetFormDegree(sp, &k);
1688:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1689:   PetscDualSpaceGetAllData(sp, &allNodes, &allMat);
1690:   nNodes = 0;
1691:   if (allNodes) {
1692:     PetscQuadratureGetData(allNodes, NULL, NULL, &nNodes, &nodes, NULL);
1693:   }
1694:   MatGetSize(allMat, &nDofs, NULL);
1695:   PetscDualSpaceGetSection(sp, &section);
1696:   PetscSectionGetStorageSize(section, &spdim);
1697:   if (spdim != nDofs) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "incompatible all matrix size");
1698:   PetscMalloc1(nDofs, &(sp->functional));
1699:   for (f = 0; f < nDofs; f++) {
1700:     PetscInt ncols, c;
1701:     const PetscInt *cols;
1702:     const PetscScalar *vals;
1703:     PetscReal *nodesf;
1704:     PetscReal *weightsf;
1705:     PetscInt nNodesf;
1706:     PetscInt countNodes;

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

1719:         for (d = 0; d < Nc; d++) {
1720:           weightsf[countNodes * Nc + d] = 0.;
1721:         }
1722:         for (d = 0; d < dim; d++) {
1723:           nodesf[countNodes * dim + d] = nodes[(cols[c] / Nc) * dim + d];
1724:         }
1725:         countNodes++;
1726:       }
1727:       weightsf[(countNodes - 1) * Nc + (cols[c] % Nc)] = PetscRealPart(vals[c]);
1728:     }
1729:     PetscQuadratureCreate(PETSC_COMM_SELF, &(sp->functional[f]));
1730:     PetscQuadratureSetData(sp->functional[f], dim, Nc, nNodesf, nodesf, weightsf);
1731:     MatRestoreRow(allMat, f, &ncols, &cols, &vals);
1732:   }
1733:   return(0);
1734: }

1736: /* take a matrix meant for k-forms and expand it to one for Ncopies */
1737: static PetscErrorCode PetscDualSpaceLagrangeMatrixCreateCopies(Mat A, PetscInt Nk, PetscInt Ncopies, Mat *Abs)
1738: {
1739:   PetscInt       m, n, i, j, k;
1740:   PetscInt       maxnnz, *nnz, *iwork;
1741:   Mat            Ac;

1745:   MatGetSize(A, &m, &n);
1746:   if (n % Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of columns in A %D is not a multiple of Nk %D", n, Nk);
1747:   PetscMalloc1(m * Ncopies, &nnz);
1748:   for (i = 0, maxnnz = 0; i < m; i++) {
1749:     PetscInt innz;
1750:     MatGetRow(A, i, &innz, NULL, NULL);
1751:     if (innz % Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "A row %D nnzs is not a multiple of Nk %D", innz, Nk);
1752:     for (j = 0; j < Ncopies; j++) nnz[i * Ncopies + j] = innz;
1753:     maxnnz = PetscMax(maxnnz, innz);
1754:   }
1755:   MatCreateSeqAIJ(PETSC_COMM_SELF, m * Ncopies, n * Ncopies, 0, nnz, &Ac);
1756:   MatSetOption(Ac, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1757:   PetscFree(nnz);
1758:   PetscMalloc1(maxnnz, &iwork);
1759:   for (i = 0; i < m; i++) {
1760:     PetscInt innz;
1761:     const PetscInt    *cols;
1762:     const PetscScalar *vals;

1764:     MatGetRow(A, i, &innz, &cols, &vals);
1765:     for (j = 0; j < innz; j++) iwork[j] = (cols[j] / Nk) * (Nk * Ncopies) + (cols[j] % Nk);
1766:     for (j = 0; j < Ncopies; j++) {
1767:       PetscInt row = i * Ncopies + j;

1769:       MatSetValues(Ac, 1, &row, innz, iwork, vals, INSERT_VALUES);
1770:       for (k = 0; k < innz; k++) iwork[k] += Nk;
1771:     }
1772:     MatRestoreRow(A, i, &innz, &cols, &vals);
1773:   }
1774:   PetscFree(iwork);
1775:   MatAssemblyBegin(Ac, MAT_FINAL_ASSEMBLY);
1776:   MatAssemblyEnd(Ac, MAT_FINAL_ASSEMBLY);
1777:   *Abs = Ac;
1778:   return(0);
1779: }

1781: /* check if a cell is a tensor product of the segment with a facet,
1782:  * specifically checking if f and f2 can be the "endpoints" (like the triangles
1783:  * at either end of a wedge) */
1784: static PetscErrorCode DMPlexPointIsTensor_Internal_Given(DM dm, PetscInt p, PetscInt f, PetscInt f2, PetscBool *isTensor)
1785: {
1786:   PetscInt        coneSize, c;
1787:   const PetscInt *cone;
1788:   const PetscInt *fCone;
1789:   const PetscInt *f2Cone;
1790:   PetscInt        fs[2];
1791:   PetscInt        meetSize, nmeet;
1792:   const PetscInt *meet;
1793:   PetscErrorCode  ierr;

1796:   fs[0] = f;
1797:   fs[1] = f2;
1798:   DMPlexGetMeet(dm, 2, fs, &meetSize, &meet);
1799:   nmeet = meetSize;
1800:   DMPlexRestoreMeet(dm, 2, fs, &meetSize, &meet);
1801:   /* two points that have a non-empty meet cannot be at opposite ends of a cell */
1802:   if (nmeet) {
1803:     *isTensor = PETSC_FALSE;
1804:     return(0);
1805:   }
1806:   DMPlexGetConeSize(dm, p, &coneSize);
1807:   DMPlexGetCone(dm, p, &cone);
1808:   DMPlexGetCone(dm, f, &fCone);
1809:   DMPlexGetCone(dm, f2, &f2Cone);
1810:   for (c = 0; c < coneSize; c++) {
1811:     PetscInt e, ef;
1812:     PetscInt d = -1, d2 = -1;
1813:     PetscInt dcount, d2count;
1814:     PetscInt t = cone[c];
1815:     PetscInt tConeSize;
1816:     PetscBool tIsTensor;
1817:     const PetscInt *tCone;

1819:     if (t == f || t == f2) continue;
1820:     /* for every other facet in the cone, check that is has
1821:      * one ridge in common with each end */
1822:     DMPlexGetConeSize(dm, t, &tConeSize);
1823:     DMPlexGetCone(dm, t, &tCone);

1825:     dcount = 0;
1826:     d2count = 0;
1827:     for (e = 0; e < tConeSize; e++) {
1828:       PetscInt q = tCone[e];
1829:       for (ef = 0; ef < coneSize - 2; ef++) {
1830:         if (fCone[ef] == q) {
1831:           if (dcount) {
1832:             *isTensor = PETSC_FALSE;
1833:             return(0);
1834:           }
1835:           d = q;
1836:           dcount++;
1837:         } else if (f2Cone[ef] == q) {
1838:           if (d2count) {
1839:             *isTensor = PETSC_FALSE;
1840:             return(0);
1841:           }
1842:           d2 = q;
1843:           d2count++;
1844:         }
1845:       }
1846:     }
1847:     /* if the whole cell is a tensor with the segment, then this
1848:      * facet should be a tensor with the segment */
1849:     DMPlexPointIsTensor_Internal_Given(dm, t, d, d2, &tIsTensor);
1850:     if (!tIsTensor) {
1851:       *isTensor = PETSC_FALSE;
1852:       return(0);
1853:     }
1854:   }
1855:   *isTensor = PETSC_TRUE;
1856:   return(0);
1857: }

1859: /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair
1860:  * that could be the opposite ends */
1861: static PetscErrorCode DMPlexPointIsTensor_Internal(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB)
1862: {
1863:   PetscInt        coneSize, c, c2;
1864:   const PetscInt *cone;
1865:   PetscErrorCode  ierr;

1868:   DMPlexGetConeSize(dm, p, &coneSize);
1869:   if (!coneSize) {
1870:     if (isTensor) *isTensor = PETSC_FALSE;
1871:     if (endA) *endA = -1;
1872:     if (endB) *endB = -1;
1873:   }
1874:   DMPlexGetCone(dm, p, &cone);
1875:   for (c = 0; c < coneSize; c++) {
1876:     PetscInt f = cone[c];
1877:     PetscInt fConeSize;

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

1882:     for (c2 = c + 1; c2 < coneSize; c2++) {
1883:       PetscInt  f2 = cone[c2];
1884:       PetscBool isTensorff2;
1885:       PetscInt f2ConeSize;

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

1890:       DMPlexPointIsTensor_Internal_Given(dm, p, f, f2, &isTensorff2);
1891:       if (isTensorff2) {
1892:         if (isTensor) *isTensor = PETSC_TRUE;
1893:         if (endA) *endA = f;
1894:         if (endB) *endB = f2;
1895:         return(0);
1896:       }
1897:     }
1898:   }
1899:   if (isTensor) *isTensor = PETSC_FALSE;
1900:   if (endA) *endA = -1;
1901:   if (endB) *endB = -1;
1902:   return(0);
1903: }

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

1913:   DMPlexIsInterpolated(dm, &interpolated);
1914:   if (interpolated != DMPLEX_INTERPOLATED_FULL) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Only for interpolated DMPlex's");
1915:   DMPlexPointIsTensor_Internal(dm, p, isTensor, endA, endB);
1916:   return(0);
1917: }

1919: /* permute a quadrature -> dof matrix so that its rows are in revlex order by nodeIdx */
1920: static PetscErrorCode MatPermuteByNodeIdx(Mat A, PetscLagNodeIndices ni, Mat *Aperm)
1921: {
1922:   PetscInt       m, n, i, j;
1923:   PetscInt       nodeIdxDim = ni->nodeIdxDim;
1924:   PetscInt       nodeVecDim = ni->nodeVecDim;
1925:   PetscInt       *perm;
1926:   IS             permIS;
1927:   IS             id;
1928:   PetscInt       *nIdxPerm;
1929:   PetscReal      *nVecPerm;

1933:   PetscLagNodeIndicesGetPermutation(ni, &perm);
1934:   MatGetSize(A, &m, &n);
1935:   PetscMalloc1(nodeIdxDim * m, &nIdxPerm);
1936:   PetscMalloc1(nodeVecDim * m, &nVecPerm);
1937:   for (i = 0; i < m; i++) for (j = 0; j < nodeIdxDim; j++) nIdxPerm[i * nodeIdxDim + j] = ni->nodeIdx[perm[i] * nodeIdxDim + j];
1938:   for (i = 0; i < m; i++) for (j = 0; j < nodeVecDim; j++) nVecPerm[i * nodeVecDim + j] = ni->nodeVec[perm[i] * nodeVecDim + j];
1939:   ISCreateGeneral(PETSC_COMM_SELF, m, perm, PETSC_USE_POINTER, &permIS);
1940:   ISSetPermutation(permIS);
1941:   ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &id);
1942:   ISSetPermutation(id);
1943:   MatPermute(A, permIS, id, Aperm);
1944:   ISDestroy(&permIS);
1945:   ISDestroy(&id);
1946:   for (i = 0; i < m; i++) perm[i] = i;
1947:   PetscFree(ni->nodeIdx);
1948:   PetscFree(ni->nodeVec);
1949:   ni->nodeIdx = nIdxPerm;
1950:   ni->nodeVec = nVecPerm;
1951:   return(0);
1952: }

1954: static PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1955: {
1956:   PetscDualSpace_Lag *lag   = (PetscDualSpace_Lag *) sp->data;
1957:   DM                  dm    = sp->dm;
1958:   DM                  dmint = NULL;
1959:   PetscInt            order;
1960:   PetscInt            Nc    = sp->Nc;
1961:   MPI_Comm            comm;
1962:   PetscBool           continuous;
1963:   PetscSection        section;
1964:   PetscInt            depth, dim, pStart, pEnd, cStart, cEnd, p, *pStratStart, *pStratEnd, d;
1965:   PetscInt            formDegree, Nk, Ncopies;
1966:   PetscInt            tensorf = -1, tensorf2 = -1;
1967:   PetscBool           tensorCell, tensorSpace;
1968:   PetscBool           uniform, trimmed;
1969:   Petsc1DNodeFamily   nodeFamily;
1970:   PetscInt            numNodeSkip;
1971:   DMPlexInterpolatedFlag interpolated;
1972:   PetscBool           isbdm;
1973:   PetscErrorCode      ierr;

1976:   /* step 1: sanitize input */
1977:   PetscObjectGetComm((PetscObject) sp, &comm);
1978:   DMGetDimension(dm, &dim);
1979:   PetscObjectTypeCompare((PetscObject)sp, "bdm", &isbdm);
1980:   if (isbdm) {
1981:     sp->k = -(dim-1); /* form degree of H-div */
1982:     PetscObjectChangeTypeName((PetscObject)sp, PETSCDUALSPACELAGRANGE);
1983:   }
1984:   PetscDualSpaceGetFormDegree(sp, &formDegree);
1985:   if (PetscAbsInt(formDegree) > dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Form degree must be bounded by dimension");
1986:   PetscDTBinomialInt(dim,PetscAbsInt(formDegree),&Nk);
1987:   if (sp->Nc <= 0 && lag->numCopies > 0) sp->Nc = Nk * lag->numCopies;
1988:   Nc = sp->Nc;
1989:   if (Nc % Nk) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "Number of components is not a multiple of form degree size");
1990:   if (lag->numCopies <= 0) lag->numCopies = Nc / Nk;
1991:   Ncopies = lag->numCopies;
1992:   if (Nc / Nk != Ncopies) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "Number of copies * (dim choose k) != Nc");
1993:   if (!dim) sp->order = 0;
1994:   order = sp->order;
1995:   uniform = sp->uniform;
1996:   if (!uniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Variable order not supported yet");
1997:   if (lag->trimmed && !formDegree) lag->trimmed = PETSC_FALSE; /* trimmed spaces are the same as full spaces for 0-forms */
1998:   if (lag->nodeType == PETSCDTNODES_DEFAULT) {
1999:     lag->nodeType = PETSCDTNODES_GAUSSJACOBI;
2000:     lag->nodeExponent = 0.;
2001:     /* trimmed spaces don't include corner vertices, so don't use end nodes by default */
2002:     lag->endNodes = lag->trimmed ? PETSC_FALSE : PETSC_TRUE;
2003:   }
2004:   /* If a trimmed space and the user did choose nodes with endpoints, skip them by default */
2005:   if (lag->numNodeSkip < 0) lag->numNodeSkip = (lag->trimmed && lag->endNodes) ? 1 : 0;
2006:   numNodeSkip = lag->numNodeSkip;
2007:   if (lag->trimmed && !order) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zeroth order trimmed elements");
2008:   if (lag->trimmed && PetscAbsInt(formDegree) == dim) { /* convert trimmed n-forms to untrimmed of one polynomial order less */
2009:     sp->order--;
2010:     order--;
2011:     lag->trimmed = PETSC_FALSE;
2012:   }
2013:   trimmed = lag->trimmed;
2014:   if (!order || PetscAbsInt(formDegree) == dim) lag->continuous = PETSC_FALSE;
2015:   continuous = lag->continuous;
2016:   DMPlexGetDepth(dm, &depth);
2017:   DMPlexGetChart(dm, &pStart, &pEnd);
2018:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2019:   if (pStart != 0 || cStart != 0) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Expect DM with chart starting at zero and cells first");
2020:   if (cEnd != 1) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Use PETSCDUALSPACEREFINED for multi-cell reference meshes");
2021:   DMPlexIsInterpolated(dm, &interpolated);
2022:   if (interpolated != DMPLEX_INTERPOLATED_FULL) {
2023:     DMPlexInterpolate(dm, &dmint);
2024:   } else {
2025:     PetscObjectReference((PetscObject)dm);
2026:     dmint = dm;
2027:   }
2028:   tensorCell = PETSC_FALSE;
2029:   if (dim > 1) {
2030:     DMPlexPointIsTensor(dmint, 0, &tensorCell, &tensorf, &tensorf2);
2031:   }
2032:   lag->tensorCell = tensorCell;
2033:   if (dim < 2 || !lag->tensorCell) lag->tensorSpace = PETSC_FALSE;
2034:   tensorSpace = lag->tensorSpace;
2035:   if (!lag->nodeFamily) {
2036:     Petsc1DNodeFamilyCreate(lag->nodeType, lag->nodeExponent, lag->endNodes, &lag->nodeFamily);
2037:   }
2038:   nodeFamily = lag->nodeFamily;
2039:   if (interpolated != DMPLEX_INTERPOLATED_FULL && continuous && (PetscAbsInt(formDegree) > 0 || order > 1)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Reference element won't support all boundary nodes");

2041:   /* step 2: construct the boundary spaces */
2042:   PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);
2043:   PetscCalloc1(pEnd,&(sp->pointSpaces));
2044:   for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);}
2045:   PetscDualSpaceSectionCreate_Internal(sp, &section);
2046:   sp->pointSection = section;
2047:   if (continuous && !(lag->interiorOnly)) {
2048:     PetscInt h;

2050:     for (p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */
2051:       PetscReal v0[3];
2052:       DMPolytopeType ptype;
2053:       PetscReal J[9], detJ;
2054:       PetscInt  q;

2056:       DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, NULL, &detJ);
2057:       DMPlexGetCellType(dm, p, &ptype);

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

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

2078:       if (hdim < PetscAbsInt(formDegree)) break;
2079:       for (p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
2080:         PetscInt suppSize, s;
2081:         const PetscInt *supp;

2083:         DMPlexGetSupportSize(dm, p, &suppSize);
2084:         DMPlexGetSupport(dm, p, &supp);
2085:         for (s = 0; s < suppSize; s++) {
2086:           DM             qdm;
2087:           PetscDualSpace qsp, psp;
2088:           PetscInt c, coneSize, q;
2089:           const PetscInt *cone;
2090:           const PetscInt *refCone;

2092:           q = supp[0];
2093:           qsp = sp->pointSpaces[q];
2094:           DMPlexGetConeSize(dm, q, &coneSize);
2095:           DMPlexGetCone(dm, q, &cone);
2096:           for (c = 0; c < coneSize; c++) if (cone[c] == p) break;
2097:           if (c == coneSize) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "cone/suppport mismatch");
2098:           PetscDualSpaceGetDM(qsp, &qdm);
2099:           DMPlexGetCone(qdm, 0, &refCone);
2100:           /* get the equivalent dual space from the support dual space */
2101:           PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp);
2102:           if (!s) {
2103:             PetscObjectReference((PetscObject)psp);
2104:             sp->pointSpaces[p] = psp;
2105:           }
2106:         }
2107:       }
2108:     }
2109:     for (p = 1; p < pEnd; p++) {
2110:       PetscInt pspdim;
2111:       if (!sp->pointSpaces[p]) continue;
2112:       PetscDualSpaceGetInteriorDimension(sp->pointSpaces[p], &pspdim);
2113:       PetscSectionSetDof(section, p, pspdim);
2114:     }
2115:   }

2117:   if (Ncopies > 1) {
2118:     Mat intMatScalar, allMatScalar;
2119:     PetscDualSpace scalarsp;
2120:     PetscDualSpace_Lag *scalarlag;

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

2150:   if (trimmed && !continuous) {
2151:     /* the dofs of a trimmed space don't have a nice tensor/lattice structure:
2152:      * just construct the continuous dual space and copy all of the data over,
2153:      * allocating it all to the cell instead of splitting it up between the boundaries */
2154:     PetscDualSpace  spcont;
2155:     PetscInt        spdim, f;
2156:     PetscQuadrature allNodes;
2157:     PetscDualSpace_Lag *lagc;
2158:     Mat             allMat;

2160:     PetscDualSpaceDuplicate(sp, &spcont);
2161:     PetscDualSpaceLagrangeSetContinuity(spcont, PETSC_TRUE);
2162:     PetscDualSpaceSetUp(spcont);
2163:     PetscDualSpaceGetDimension(spcont, &spdim);
2164:     sp->spdim = sp->spintdim = spdim;
2165:     PetscSectionSetDof(section, 0, spdim);
2166:     PetscDualSpaceSectionSetUp_Internal(sp, section);
2167:     PetscMalloc1(spdim, &(sp->functional));
2168:     for (f = 0; f < spdim; f++) {
2169:       PetscQuadrature fn;

2171:       PetscDualSpaceGetFunctional(spcont, f, &fn);
2172:       PetscObjectReference((PetscObject)fn);
2173:       sp->functional[f] = fn;
2174:     }
2175:     PetscDualSpaceGetAllData(spcont, &allNodes, &allMat);
2176:     PetscObjectReference((PetscObject) allNodes);
2177:     PetscObjectReference((PetscObject) allNodes);
2178:     sp->allNodes = sp->intNodes = allNodes;
2179:     PetscObjectReference((PetscObject) allMat);
2180:     PetscObjectReference((PetscObject) allMat);
2181:     sp->allMat = sp->intMat = allMat;
2182:     lagc = (PetscDualSpace_Lag *) spcont->data;
2183:     PetscLagNodeIndicesReference(lagc->vertIndices);
2184:     lag->vertIndices = lagc->vertIndices;
2185:     PetscLagNodeIndicesReference(lagc->allNodeIndices);
2186:     PetscLagNodeIndicesReference(lagc->allNodeIndices);
2187:     lag->intNodeIndices = lagc->allNodeIndices;
2188:     lag->allNodeIndices = lagc->allNodeIndices;
2189:     PetscDualSpaceDestroy(&spcont);
2190:     PetscFree2(pStratStart, pStratEnd);
2191:     DMDestroy(&dmint);
2192:     return(0);
2193:   }

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

2199:     if (trimmed) {
2200:       /* there is one dof in the interior of the a trimmed element for each full polynomial of with degree at most
2201:        * order + k - dim - 1 */
2202:       if (order + PetscAbsInt(formDegree) > dim) {
2203:         PetscInt sum = order + PetscAbsInt(formDegree) - dim - 1;
2204:         PetscInt nDofs;

2206:         PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));
2207:         MatGetSize(sp->intMat, &nDofs, NULL);
2208:         PetscSectionSetDof(section, 0, nDofs);
2209:       }
2210:       PetscDualSpaceSectionSetUp_Internal(sp, section);
2211:       PetscDualSpaceCreateAllDataFromInteriorData(sp);
2212:       PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2213:     } else {
2214:       if (!continuous) {
2215:         /* if discontinuous just construct one node for each set of dofs (a set of dofs is a basis for the k-form
2216:          * space) */
2217:         PetscInt sum = order;
2218:         PetscInt nDofs;

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

2241:           PetscDualSpaceDuplicate(sp, &trimmedsp);
2242:           PetscDualSpaceLagrangeSetTrimmed(trimmedsp, PETSC_TRUE);
2243:           PetscDualSpaceSetOrder(trimmedsp, order + PetscAbsInt(formDegree) - dim);
2244:           PetscDualSpaceSetFormDegree(trimmedsp, trFormDegree);
2245:           trimmedlag = (PetscDualSpace_Lag *) trimmedsp->data;
2246:           trimmedlag->numNodeSkip = numNodeSkip + 1;
2247:           PetscDualSpaceSetUp(trimmedsp);
2248:           PetscDualSpaceGetAllData(trimmedsp, &intNodes, &intMat);
2249:           PetscObjectReference((PetscObject)intNodes);
2250:           sp->intNodes = intNodes;
2251:           PetscObjectReference((PetscObject)intMat);
2252:           sp->intMat = intMat;
2253:           MatGetSize(sp->intMat, &nDofs, NULL);
2254:           PetscLagNodeIndicesReference(trimmedlag->allNodeIndices);
2255:           lag->intNodeIndices = trimmedlag->allNodeIndices;
2256:           PetscDualSpaceDestroy(&trimmedsp);
2257:           PetscSectionSetDof(section, 0, nDofs);
2258:         }
2259:         PetscDualSpaceSectionSetUp_Internal(sp, section);
2260:         PetscDualSpaceCreateAllDataFromInteriorData(sp);
2261:         PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2262:       }
2263:     }
2264:   } else {
2265:     PetscQuadrature intNodesTrace = NULL;
2266:     PetscQuadrature intNodesFiber = NULL;
2267:     PetscQuadrature intNodes = NULL;
2268:     PetscLagNodeIndices intNodeIndices = NULL;
2269:     Mat             intMat = NULL;

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

2277:       if (sp->pointSpaces[tensorf]) {
2278:         PetscObjectReference((PetscObject)(sp->pointSpaces[tensorf]));
2279:         trace = sp->pointSpaces[tensorf];
2280:       } else {
2281:         PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,formDegree,Ncopies,PETSC_TRUE,&trace);
2282:       }
2283:       PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,0,1,PETSC_TRUE,&fiber);
2284:       tracel = (PetscDualSpace_Lag *) trace->data;
2285:       fiberl = (PetscDualSpace_Lag *) fiber->data;
2286:       PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));
2287:       PetscDualSpaceGetInteriorData(trace, &intNodesTrace, &intMatTrace);
2288:       PetscDualSpaceGetInteriorData(fiber, &intNodesFiber, &intMatFiber);
2289:       if (intNodesTrace && intNodesFiber) {
2290:         PetscQuadratureCreateTensor(intNodesTrace, intNodesFiber, &intNodes);
2291:         MatTensorAltV(intMatTrace, intMatFiber, dim-1, formDegree, 1, 0, &intMat);
2292:         PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, formDegree, fiberl->intNodeIndices, 1, 0, &intNodeIndices);
2293:       }
2294:       PetscObjectReference((PetscObject) intNodesTrace);
2295:       PetscObjectReference((PetscObject) intNodesFiber);
2296:       PetscDualSpaceDestroy(&fiber);
2297:       PetscDualSpaceDestroy(&trace);
2298:     }
2299:     if (PetscAbsInt(formDegree) > 0) { /* get the trace (k-1)-forms on the first facet, and the 1-forms on the edge,
2300:                                           and wedge them together to create the remaining k-form dofs */
2301:       PetscDualSpace  trace, fiber;
2302:       PetscDualSpace_Lag *tracel, *fiberl;
2303:       PetscQuadrature intNodesTrace2, intNodesFiber2, intNodes2;
2304:       PetscLagNodeIndices intNodeIndices2;
2305:       Mat             intMatTrace, intMatFiber, intMat2;
2306:       PetscInt        traceDegree = formDegree > 0 ? formDegree - 1 : formDegree + 1;
2307:       PetscInt        fiberDegree = formDegree > 0 ? 1 : -1;

2309:       PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,traceDegree,Ncopies,PETSC_TRUE,&trace);
2310:       PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,fiberDegree,1,PETSC_TRUE,&fiber);
2311:       tracel = (PetscDualSpace_Lag *) trace->data;
2312:       fiberl = (PetscDualSpace_Lag *) fiber->data;
2313:       if (!lag->vertIndices) {
2314:         PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));
2315:       }
2316:       PetscDualSpaceGetInteriorData(trace, &intNodesTrace2, &intMatTrace);
2317:       PetscDualSpaceGetInteriorData(fiber, &intNodesFiber2, &intMatFiber);
2318:       if (intNodesTrace2 && intNodesFiber2) {
2319:         PetscQuadratureCreateTensor(intNodesTrace2, intNodesFiber2, &intNodes2);
2320:         MatTensorAltV(intMatTrace, intMatFiber, dim-1, traceDegree, 1, fiberDegree, &intMat2);
2321:         PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, traceDegree, fiberl->intNodeIndices, 1, fiberDegree, &intNodeIndices2);
2322:         if (!intMat) {
2323:           intMat = intMat2;
2324:           intNodes = intNodes2;
2325:           intNodeIndices = intNodeIndices2;
2326:         } else {
2327:           /* merge the matrices, quadrature points, and nodes */
2328:           PetscInt         nM;
2329:           PetscInt         nDof, nDof2;
2330:           PetscInt        *toMerged = NULL, *toMerged2 = NULL;
2331:           PetscQuadrature  merged = NULL;
2332:           PetscLagNodeIndices intNodeIndicesMerged = NULL;
2333:           Mat              matMerged = NULL;

2335:           MatGetSize(intMat, &nDof, 0);
2336:           MatGetSize(intMat2, &nDof2, 0);
2337:           PetscQuadraturePointsMerge(intNodes, intNodes2, &merged, &toMerged, &toMerged2);
2338:           PetscQuadratureGetData(merged, NULL, NULL, &nM, NULL, NULL);
2339:           MatricesMerge(intMat, intMat2, dim, formDegree, nM, toMerged, toMerged2, &matMerged);
2340:           PetscLagNodeIndicesMerge(intNodeIndices, intNodeIndices2, &intNodeIndicesMerged);
2341:           PetscFree(toMerged);
2342:           PetscFree(toMerged2);
2343:           MatDestroy(&intMat);
2344:           MatDestroy(&intMat2);
2345:           PetscQuadratureDestroy(&intNodes);
2346:           PetscQuadratureDestroy(&intNodes2);
2347:           PetscLagNodeIndicesDestroy(&intNodeIndices);
2348:           PetscLagNodeIndicesDestroy(&intNodeIndices2);
2349:           intNodes = merged;
2350:           intMat = matMerged;
2351:           intNodeIndices = intNodeIndicesMerged;
2352:           if (!trimmed) {
2353:             /* I think users expect that, when a node has a full basis for the k-forms,
2354:              * they should be consecutive dofs.  That isn't the case for trimmed spaces,
2355:              * but is for some of the nodes in untrimmed spaces, so in that case we
2356:              * sort them to group them by node */
2357:             Mat intMatPerm;

2359:             MatPermuteByNodeIdx(intMat, intNodeIndices, &intMatPerm);
2360:             MatDestroy(&intMat);
2361:             intMat = intMatPerm;
2362:           }
2363:         }
2364:       }
2365:       PetscDualSpaceDestroy(&fiber);
2366:       PetscDualSpaceDestroy(&trace);
2367:     }
2368:     PetscQuadratureDestroy(&intNodesTrace);
2369:     PetscQuadratureDestroy(&intNodesFiber);
2370:     sp->intNodes = intNodes;
2371:     sp->intMat = intMat;
2372:     lag->intNodeIndices = intNodeIndices;
2373:     {
2374:       PetscInt nDofs = 0;

2376:       if (intMat) {
2377:         MatGetSize(intMat, &nDofs, NULL);
2378:       }
2379:       PetscSectionSetDof(section, 0, nDofs);
2380:     }
2381:     PetscDualSpaceSectionSetUp_Internal(sp, section);
2382:     if (continuous) {
2383:       PetscDualSpaceCreateAllDataFromInteriorData(sp);
2384:       PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2385:     } else {
2386:       PetscObjectReference((PetscObject) intNodes);
2387:       sp->allNodes = intNodes;
2388:       PetscObjectReference((PetscObject) intMat);
2389:       sp->allMat = intMat;
2390:       PetscLagNodeIndicesReference(intNodeIndices);
2391:       lag->allNodeIndices = intNodeIndices;
2392:     }
2393:   }
2394:   PetscSectionGetStorageSize(section, &sp->spdim);
2395:   PetscSectionGetConstrainedStorageSize(section, &sp->spintdim);
2396:   PetscDualSpaceComputeFunctionalsFromAllData(sp);
2397:   PetscFree2(pStratStart, pStratEnd);
2398:   DMDestroy(&dmint);
2399:   return(0);
2400: }

2402: /* Create a matrix that represents the transformation that DMPlexVecGetClosure() would need
2403:  * to get the representation of the dofs for a mesh point if the mesh point had this orientation
2404:  * relative to the cell */
2405: PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat)
2406: {
2407:   PetscDualSpace_Lag *lag;
2408:   DM dm;
2409:   PetscLagNodeIndices vertIndices, intNodeIndices;
2410:   PetscLagNodeIndices ni;
2411:   PetscInt nodeIdxDim, nodeVecDim, nNodes;
2412:   PetscInt formDegree;
2413:   PetscInt *perm, *permOrnt;
2414:   PetscInt *nnz;
2415:   PetscInt n;
2416:   PetscInt maxGroupSize;
2417:   PetscScalar *V, *W, *work;
2418:   Mat A;

2422:   if (!sp->spintdim) {
2423:     *symMat = NULL;
2424:     return(0);
2425:   }
2426:   lag = (PetscDualSpace_Lag *) sp->data;
2427:   vertIndices = lag->vertIndices;
2428:   intNodeIndices = lag->intNodeIndices;
2429:   PetscDualSpaceGetDM(sp, &dm);
2430:   PetscDualSpaceGetFormDegree(sp, &formDegree);
2431:   PetscNew(&ni);
2432:   ni->refct = 1;
2433:   ni->nodeIdxDim = nodeIdxDim = intNodeIndices->nodeIdxDim;
2434:   ni->nodeVecDim = nodeVecDim = intNodeIndices->nodeVecDim;
2435:   ni->nNodes = nNodes = intNodeIndices->nNodes;
2436:   PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
2437:   PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
2438:   /* push forward the dofs by the symmetry of the reference element induced by ornt */
2439:   PetscLagNodeIndicesPushForward(dm, vertIndices, 0, vertIndices, intNodeIndices, ornt, formDegree, ni->nodeIdx, ni->nodeVec);
2440:   /* get the revlex order for both the original and transformed dofs */
2441:   PetscLagNodeIndicesGetPermutation(intNodeIndices, &perm);
2442:   PetscLagNodeIndicesGetPermutation(ni, &permOrnt);
2443:   PetscMalloc1(nNodes, &nnz);
2444:   for (n = 0, maxGroupSize = 0; n < nNodes;) { /* incremented in the loop */
2445:     PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2446:     PetscInt m, nEnd;
2447:     PetscInt groupSize;
2448:     /* for each group of dofs that have the same nodeIdx coordinate */
2449:     for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2450:       PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2451:       PetscInt d;

2453:       /* compare the oriented permutation indices */
2454:       for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break;
2455:       if (d < nodeIdxDim) break;
2456:     }
2457:     /* permOrnt[[n, nEnd)] is a group of dofs that, under the symmetry are at the same location */

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

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

2469:         /* compare the oriented permutation indices */
2470:         for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break;
2471:         if (d < nodeIdxDim) break;
2472:       }
2473:       if (m < nEnd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs with same index after symmetry not same block size");
2474:     }
2475:     groupSize = nEnd - n;
2476:     /* each pushforward dof vector will be expressed in a basis of the unpermuted dofs */
2477:     for (m = n; m < nEnd; m++) nnz[permOrnt[m]] = groupSize;

2479:     maxGroupSize = PetscMax(maxGroupSize, nEnd - n);
2480:     n = nEnd;
2481:   }
2482:   if (maxGroupSize > nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs are not in blocks that can be solved");
2483:   MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes, nNodes, 0, nnz, &A);
2484:   PetscFree(nnz);
2485:   PetscMalloc3(maxGroupSize * nodeVecDim, &V, maxGroupSize * nodeVecDim, &W, nodeVecDim * 2, &work);
2486:   for (n = 0; n < nNodes;) { /* incremented in the loop */
2487:     PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2488:     PetscInt nEnd;
2489:     PetscInt m;
2490:     PetscInt groupSize;
2491:     for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2492:       PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2493:       PetscInt d;

2495:       /* compare the oriented permutation indices */
2496:       for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break;
2497:       if (d < nodeIdxDim) break;
2498:     }
2499:     groupSize = nEnd - n;
2500:     /* get all of the vectors from the original and all of the pushforward vectors */
2501:     for (m = n; m < nEnd; m++) {
2502:       PetscInt d;

2504:       for (d = 0; d < nodeVecDim; d++) {
2505:         V[(m - n) * nodeVecDim + d] = intNodeIndices->nodeVec[perm[m] * nodeVecDim + d];
2506:         W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
2507:       }
2508:     }
2509:     /* now we have to solve for W in terms of V: the systems isn't always square, but the span
2510:      * of V and W should always be the same, so the solution of the normal equations works */
2511:     {
2512:       char transpose = 'N';
2513:       PetscBLASInt bm = nodeVecDim;
2514:       PetscBLASInt bn = groupSize;
2515:       PetscBLASInt bnrhs = groupSize;
2516:       PetscBLASInt blda = bm;
2517:       PetscBLASInt bldb = bm;
2518:       PetscBLASInt blwork = 2 * nodeVecDim;
2519:       PetscBLASInt info;

2521:       PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&bm,&bn,&bnrhs,V,&blda,W,&bldb,work,&blwork, &info));
2522:       if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS");
2523:       /* repack */
2524:       {
2525:         PetscInt i, j;

2527:         for (i = 0; i < groupSize; i++) {
2528:           for (j = 0; j < groupSize; j++) {
2529:             /* notice the different leading dimension */
2530:             V[i * groupSize + j] = W[i * nodeVecDim + j];
2531:           }
2532:         }
2533:       }
2534:     }
2535:     MatSetValues(A, groupSize, &permOrnt[n], groupSize, &perm[n], V, INSERT_VALUES);
2536:     n = nEnd;
2537:   }
2538:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
2539:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
2540:   *symMat = A;
2541:   PetscFree3(V,W,work);
2542:   PetscLagNodeIndicesDestroy(&ni);
2543:   return(0);
2544: }

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

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

2550: /* the existing interface for symmetries is insufficient for all cases:
2551:  * - it should be sufficient for form degrees that are scalar (0 and n)
2552:  * - it should be sufficient for hypercube dofs
2553:  * - it isn't sufficient for simplex cells with non-scalar form degrees if
2554:  *   there are any dofs in the interior
2555:  *
2556:  * We compute the general transformation matrices, and if they fit, we return them,
2557:  * otherwise we error (but we should probably change the interface to allow for
2558:  * these symmetries)
2559:  */
2560: static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
2561: {
2562:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2563:   PetscInt           dim, order, Nc;
2564:   PetscErrorCode     ierr;

2567:   PetscDualSpaceGetOrder(sp,&order);
2568:   PetscDualSpaceGetNumComponents(sp,&Nc);
2569:   DMGetDimension(sp->dm,&dim);
2570:   if (!lag->symComputed) { /* store symmetries */
2571:     PetscInt       pStart, pEnd, p;
2572:     PetscInt       numPoints;
2573:     PetscInt       numFaces;
2574:     PetscInt       spintdim;
2575:     PetscInt       ***symperms;
2576:     PetscScalar    ***symflips;

2578:     DMPlexGetChart(sp->dm, &pStart, &pEnd);
2579:     numPoints = pEnd - pStart;
2580:     DMPlexGetConeSize(sp->dm, 0, &numFaces);
2581:     PetscCalloc1(numPoints,&symperms);
2582:     PetscCalloc1(numPoints,&symflips);
2583:     spintdim = sp->spintdim;
2584:     /* The nodal symmetry behavior is not present when tensorSpace != tensorCell: someone might want this for the "S"
2585:      * family of FEEC spaces.  Most used in particular are discontinuous polynomial L2 spaces in tensor cells, where
2586:      * the symmetries are not necessary for FE assembly.  So for now we assume this is the case and don't return
2587:      * symmetries if tensorSpace != tensorCell */
2588:     if (spintdim && 0 < dim && dim < 3 && (lag->tensorSpace == lag->tensorCell)) { /* compute self symmetries */
2589:       PetscInt **cellSymperms;
2590:       PetscScalar **cellSymflips;
2591:       PetscInt ornt;
2592:       PetscInt nCopies = Nc / lag->intNodeIndices->nodeVecDim;
2593:       PetscInt nNodes = lag->intNodeIndices->nNodes;

2595:       lag->numSelfSym = 2 * numFaces;
2596:       lag->selfSymOff = numFaces;
2597:       PetscCalloc1(2*numFaces,&cellSymperms);
2598:       PetscCalloc1(2*numFaces,&cellSymflips);
2599:       /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
2600:       symperms[0] = &cellSymperms[numFaces];
2601:       symflips[0] = &cellSymflips[numFaces];
2602:       if (lag->intNodeIndices->nodeVecDim * nCopies != Nc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs");
2603:       if (nNodes * nCopies != spintdim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs");
2604:       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 */
2605:         Mat symMat;
2606:         PetscInt *perm;
2607:         PetscScalar *flips;
2608:         PetscInt i;

2610:         if (!ornt) continue;
2611:         PetscMalloc1(spintdim, &perm);
2612:         PetscCalloc1(spintdim, &flips);
2613:         for (i = 0; i < spintdim; i++) perm[i] = -1;
2614:         PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, ornt, &symMat);
2615:         for (i = 0; i < nNodes; i++) {
2616:           PetscInt ncols;
2617:           PetscInt j, k;
2618:           const PetscInt *cols;
2619:           const PetscScalar *vals;
2620:           PetscBool nz_seen = PETSC_FALSE;

2622:           MatGetRow(symMat, i, &ncols, &cols, &vals);
2623:           for (j = 0; j < ncols; j++) {
2624:             if (PetscAbsScalar(vals[j]) > PETSC_SMALL) {
2625:               if (nz_seen) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2626:               nz_seen = PETSC_TRUE;
2627:               if (PetscAbsReal(PetscAbsScalar(vals[j]) - PetscRealConstant(1.)) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2628:               if (PetscAbsReal(PetscImaginaryPart(vals[j])) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2629:               if (perm[cols[j] * nCopies] >= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2630:               for (k = 0; k < nCopies; k++) {
2631:                 perm[cols[j] * nCopies + k] = i * nCopies + k;
2632:               }
2633:               if (PetscRealPart(vals[j]) < 0.) {
2634:                 for (k = 0; k < nCopies; k++) {
2635:                   flips[i * nCopies + k] = -1.;
2636:                 }
2637:               } else {
2638:                 for (k = 0; k < nCopies; k++) {
2639:                   flips[i * nCopies + k] = 1.;
2640:                 }
2641:               }
2642:             }
2643:           }
2644:           MatRestoreRow(symMat, i, &ncols, &cols, &vals);
2645:         }
2646:         MatDestroy(&symMat);
2647:         /* if there were no sign flips, keep NULL */
2648:         for (i = 0; i < spintdim; i++) if (flips[i] != 1.) break;
2649:         if (i == spintdim) {
2650:           PetscFree(flips);
2651:           flips = NULL;
2652:         }
2653:         /* if the permutation is identity, keep NULL */
2654:         for (i = 0; i < spintdim; i++) if (perm[i] != i) break;
2655:         if (i == spintdim) {
2656:           PetscFree(perm);
2657:           perm = NULL;
2658:         }
2659:         symperms[0][ornt] = perm;
2660:         symflips[0][ornt] = flips;
2661:       }
2662:       /* if no orientations produced non-identity permutations, keep NULL */
2663:       for (ornt = -numFaces; ornt < numFaces; ornt++) if (symperms[0][ornt]) break;
2664:       if (ornt == numFaces) {
2665:         PetscFree(cellSymperms);
2666:         symperms[0] = NULL;
2667:       }
2668:       /* if no orientations produced sign flips, keep NULL */
2669:       for (ornt = -numFaces; ornt < numFaces; ornt++) if (symflips[0][ornt]) break;
2670:       if (ornt == numFaces) {
2671:         PetscFree(cellSymflips);
2672:         symflips[0] = NULL;
2673:       }
2674:     }
2675:     { /* get the symmetries of closure points */
2676:       PetscInt closureSize = 0;
2677:       PetscInt *closure = NULL;
2678:       PetscInt r;

2680:       DMPlexGetTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);
2681:       for (r = 0; r < closureSize; r++) {
2682:         PetscDualSpace psp;
2683:         PetscInt point = closure[2 * r];
2684:         PetscInt pspintdim;
2685:         const PetscInt ***psymperms = NULL;
2686:         const PetscScalar ***psymflips = NULL;

2688:         if (!point) continue;
2689:         PetscDualSpaceGetPointSubspace(sp, point, &psp);
2690:         if (!psp) continue;
2691:         PetscDualSpaceGetInteriorDimension(psp, &pspintdim);
2692:         if (!pspintdim) continue;
2693:         PetscDualSpaceGetSymmetries(psp,&psymperms,&psymflips);
2694:         symperms[r] = (PetscInt **) (psymperms ? psymperms[0] : NULL);
2695:         symflips[r] = (PetscScalar **) (psymflips ? psymflips[0] : NULL);
2696:       }
2697:       DMPlexRestoreTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);
2698:     }
2699:     for (p = 0; p < pEnd; p++) if (symperms[p]) break;
2700:     if (p == pEnd) {
2701:       PetscFree(symperms);
2702:       symperms = NULL;
2703:     }
2704:     for (p = 0; p < pEnd; p++) if (symflips[p]) break;
2705:     if (p == pEnd) {
2706:       PetscFree(symflips);
2707:       symflips = NULL;
2708:     }
2709:     lag->symperms = symperms;
2710:     lag->symflips = symflips;
2711:     lag->symComputed = PETSC_TRUE;
2712:   }
2713:   if (perms) *perms = (const PetscInt ***) lag->symperms;
2714:   if (flips) *flips = (const PetscScalar ***) lag->symflips;
2715:   return(0);
2716: }

2718: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2719: {
2720:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2725:   *continuous = lag->continuous;
2726:   return(0);
2727: }

2729: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2730: {
2731:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

2735:   lag->continuous = continuous;
2736:   return(0);
2737: }

2739: /*@
2740:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

2742:   Not Collective

2744:   Input Parameter:
2745: . sp         - the PetscDualSpace

2747:   Output Parameter:
2748: . continuous - flag for element continuity

2750:   Level: intermediate

2752: .seealso: PetscDualSpaceLagrangeSetContinuity()
2753: @*/
2754: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2755: {

2761:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2762:   return(0);
2763: }

2765: /*@
2766:   PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous

2768:   Logically Collective on sp

2770:   Input Parameters:
2771: + sp         - the PetscDualSpace
2772: - continuous - flag for element continuity

2774:   Options Database:
2775: . -petscdualspace_lagrange_continuity <bool>

2777:   Level: intermediate

2779: .seealso: PetscDualSpaceLagrangeGetContinuity()
2780: @*/
2781: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2782: {

2788:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2789:   return(0);
2790: }

2792: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
2793: {
2794:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2797:   *tensor = lag->tensorSpace;
2798:   return(0);
2799: }

2801: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
2802: {
2803:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2806:   lag->tensorSpace = tensor;
2807:   return(0);
2808: }

2810: static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Lagrange(PetscDualSpace sp, PetscBool *trimmed)
2811: {
2812:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2815:   *trimmed = lag->trimmed;
2816:   return(0);
2817: }

2819: static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp, PetscBool trimmed)
2820: {
2821:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

2824:   lag->trimmed = trimmed;
2825:   return(0);
2826: }

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

2833:   if (nodeType) *nodeType = lag->nodeType;
2834:   if (boundary) *boundary = lag->endNodes;
2835:   if (exponent) *exponent = lag->nodeExponent;
2836:   return(0);
2837: }

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

2844:   if (nodeType == PETSCDTNODES_GAUSSJACOBI && exponent <= -1.) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Exponent must be > -1");
2845:   lag->nodeType = nodeType;
2846:   lag->endNodes = boundary;
2847:   lag->nodeExponent = exponent;
2848:   return(0);
2849: }

2851: /*@
2852:   PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space

2854:   Not collective

2856:   Input Parameter:
2857: . sp - The PetscDualSpace

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

2862:   Level: intermediate

2864: .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
2865: @*/
2866: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
2867: {

2873:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));
2874:   return(0);
2875: }

2877: /*@
2878:   PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space

2880:   Not collective

2882:   Input Parameters:
2883: + sp - The PetscDualSpace
2884: - tensor - Whether the dual space has tensor layout (vs. simplicial)

2886:   Level: intermediate

2888: .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
2889: @*/
2890: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
2891: {

2896:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));
2897:   return(0);
2898: }

2900: /*@
2901:   PetscDualSpaceLagrangeGetTrimmed - Get the trimmed nature of the dual space

2903:   Not collective

2905:   Input Parameter:
2906: . sp - The PetscDualSpace

2908:   Output Parameter:
2909: . 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)

2911:   Level: intermediate

2913: .seealso: PetscDualSpaceLagrangeSetTrimmed(), PetscDualSpaceCreate()
2914: @*/
2915: PetscErrorCode PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp, PetscBool *trimmed)
2916: {

2922:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTrimmed_C",(PetscDualSpace,PetscBool *),(sp,trimmed));
2923:   return(0);
2924: }

2926: /*@
2927:   PetscDualSpaceLagrangeSetTrimmed - Set the trimmed nature of the dual space

2929:   Not collective

2931:   Input Parameters:
2932: + sp - The PetscDualSpace
2933: - 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)

2935:   Level: intermediate

2937: .seealso: PetscDualSpaceLagrangeGetTrimmed(), PetscDualSpaceCreate()
2938: @*/
2939: PetscErrorCode PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp, PetscBool trimmed)
2940: {

2945:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTrimmed_C",(PetscDualSpace,PetscBool),(sp,trimmed));
2946:   return(0);
2947: }

2949: /*@
2950:   PetscDualSpaceLagrangeGetNodeType - Get a description of how nodes are laid out for Lagrange polynomials in this
2951:   dual space

2953:   Not collective

2955:   Input Parameter:
2956: . sp - The PetscDualSpace

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

2965:   Level: advanced

2967: .seealso: PetscDTNodeType, PetscDualSpaceLagrangeSetNodeType()
2968: @*/
2969: PetscErrorCode PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
2970: {

2978:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetNodeType_C",(PetscDualSpace,PetscDTNodeType *,PetscBool *,PetscReal *),(sp,nodeType,boundary,exponent));
2979:   return(0);
2980: }

2982: /*@
2983:   PetscDualSpaceLagrangeSetNodeType - Set a description of how nodes are laid out for Lagrange polynomials in this
2984:   dual space

2986:   Logically collective

2988:   Input Parameters:
2989: + sp - The PetscDualSpace
2990: . nodeType - The type of nodes
2991: . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that
2992:              include the boundary are Gauss-Lobatto-Jacobi nodes)
2993: - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function
2994:              '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type

2996:   Level: advanced

2998: .seealso: PetscDTNodeType, PetscDualSpaceLagrangeGetNodeType()
2999: @*/
3000: PetscErrorCode PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
3001: {

3006:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetNodeType_C",(PetscDualSpace,PetscDTNodeType,PetscBool,PetscReal),(sp,nodeType,boundary,exponent));
3007:   return(0);
3008: }


3011: static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
3012: {
3014:   sp->ops->destroy              = PetscDualSpaceDestroy_Lagrange;
3015:   sp->ops->view                 = PetscDualSpaceView_Lagrange;
3016:   sp->ops->setfromoptions       = PetscDualSpaceSetFromOptions_Lagrange;
3017:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Lagrange;
3018:   sp->ops->setup                = PetscDualSpaceSetUp_Lagrange;
3019:   sp->ops->createheightsubspace = NULL;
3020:   sp->ops->createpointsubspace  = NULL;
3021:   sp->ops->getsymmetries        = PetscDualSpaceGetSymmetries_Lagrange;
3022:   sp->ops->apply                = PetscDualSpaceApplyDefault;
3023:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
3024:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
3025:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
3026:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;
3027:   return(0);
3028: }

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

3033:   Level: intermediate

3035: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
3036: M*/
3037: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
3038: {
3039:   PetscDualSpace_Lag *lag;
3040:   PetscErrorCode      ierr;

3044:   PetscNewLog(sp,&lag);
3045:   sp->data = lag;

3047:   lag->tensorCell  = PETSC_FALSE;
3048:   lag->tensorSpace = PETSC_FALSE;
3049:   lag->continuous  = PETSC_TRUE;
3050:   lag->numCopies   = PETSC_DEFAULT;
3051:   lag->numNodeSkip = PETSC_DEFAULT;
3052:   lag->nodeType    = PETSCDTNODES_DEFAULT;

3054:   PetscDualSpaceInitialize_Lagrange(sp);
3055:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
3056:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
3057:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
3058:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
3059:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Lagrange);
3060:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Lagrange);
3061:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Lagrange);
3062:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Lagrange);
3063:   return(0);
3064: }