Actual source code: dspacelagrange.c

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

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

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

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

     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: {
     53:   PetscInt       i, nc;

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

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

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

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

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

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

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

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

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

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

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

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

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

    350: static PetscErrorCode PetscLagNodeIndicesDuplicate(PetscLagNodeIndices ni, PetscLagNodeIndices *niNew)
    351: {

    355:   PetscNew(niNew);
    356:   (*niNew)->refct = 1;
    357:   (*niNew)->nodeIdxDim = ni->nodeIdxDim;
    358:   (*niNew)->nodeVecDim = ni->nodeVecDim;
    359:   (*niNew)->nNodes = ni->nNodes;
    360:   PetscMalloc1(ni->nNodes * ni->nodeIdxDim, &((*niNew)->nodeIdx));
    361:   PetscArraycpy((*niNew)->nodeIdx, ni->nodeIdx, ni->nNodes * ni->nodeIdxDim);
    362:   PetscMalloc1(ni->nNodes * ni->nodeVecDim, &((*niNew)->nodeVec));
    363:   PetscArraycpy((*niNew)->nodeVec, ni->nodeVec, ni->nNodes * ni->nodeVecDim);
    364:   (*niNew)->perm = NULL;
    365:   return(0);
    366: }

    368: static PetscErrorCode PetscLagNodeIndicesDestroy(PetscLagNodeIndices *ni)
    369: {

    373:   if (!(*ni)) return(0);
    374:   if (--(*ni)->refct > 0) {
    375:     *ni = NULL;
    376:     return(0);
    377:   }
    378:   PetscFree((*ni)->nodeIdx);
    379:   PetscFree((*ni)->nodeVec);
    380:   PetscFree((*ni)->perm);
    381:   PetscFree(*ni);
    382:   *ni = NULL;
    383:   return(0);
    384: }

    386: /* The vertices are given nodeIdx coordinates (e.g. the corners of the barycentric triangle).  Those coordinates are
    387:  * in some other order, and to understand the effect of different symmetries, we need them to be in closure order.
    388:  *
    389:  * If sortIdx is PETSC_FALSE, the coordinates are already in revlex order, otherwise we must sort them
    390:  * to that order before we do the real work of this function, which is
    391:  *
    392:  * - mark the vertices in closure order
    393:  * - sort them in revlex order
    394:  * - use the resulting permutation to list the vertex coordinates in closure order
    395:  */
    396: static PetscErrorCode PetscLagNodeIndicesComputeVertexOrder(DM dm, PetscLagNodeIndices ni, PetscBool sortIdx)
    397: {
    398:   PetscInt        v, w, vStart, vEnd, c, d;
    399:   PetscInt        nVerts;
    400:   PetscInt        closureSize = 0;
    401:   PetscInt       *closure = NULL;
    402:   PetscInt       *closureOrder;
    403:   PetscInt       *invClosureOrder;
    404:   PetscInt       *revlexOrder;
    405:   PetscInt       *newNodeIdx;
    406:   PetscInt        dim;
    407:   Vec             coordVec;
    408:   const PetscScalar *coords;
    409:   PetscErrorCode  ierr;

    412:   DMGetDimension(dm, &dim);
    413:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    414:   nVerts = vEnd - vStart;
    415:   PetscMalloc1(nVerts, &closureOrder);
    416:   PetscMalloc1(nVerts, &invClosureOrder);
    417:   PetscMalloc1(nVerts, &revlexOrder);
    418:   if (sortIdx) { /* bubble sort nodeIdx into revlex order */
    419:     PetscInt nodeIdxDim = ni->nodeIdxDim;
    420:     PetscInt *idxOrder;

    422:     PetscMalloc1(nVerts * nodeIdxDim, &newNodeIdx);
    423:     PetscMalloc1(nVerts, &idxOrder);
    424:     for (v = 0; v < nVerts; v++) idxOrder[v] = v;
    425:     for (v = 0; v < nVerts; v++) {
    426:       for (w = v + 1; w < nVerts; w++) {
    427:         const PetscInt *iv = &(ni->nodeIdx[idxOrder[v] * nodeIdxDim]);
    428:         const PetscInt *iw = &(ni->nodeIdx[idxOrder[w] * nodeIdxDim]);
    429:         PetscInt diff = 0;

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

    435:           idxOrder[v] = idxOrder[w];
    436:           idxOrder[w] = swap;
    437:         }
    438:       }
    439:     }
    440:     for (v = 0; v < nVerts; v++) {
    441:       for (d = 0; d < nodeIdxDim; d++) {
    442:         newNodeIdx[v * ni->nodeIdxDim + d] = ni->nodeIdx[idxOrder[v] * nodeIdxDim + d];
    443:       }
    444:     }
    445:     PetscFree(ni->nodeIdx);
    446:     ni->nodeIdx = newNodeIdx;
    447:     newNodeIdx = NULL;
    448:     PetscFree(idxOrder);
    449:   }
    450:   DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);
    451:   c = closureSize - nVerts;
    452:   for (v = 0; v < nVerts; v++) closureOrder[v] = closure[2 * (c + v)] - vStart;
    453:   for (v = 0; v < nVerts; v++) invClosureOrder[closureOrder[v]] = v;
    454:   DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);
    455:   DMGetCoordinatesLocal(dm, &coordVec);
    456:   VecGetArrayRead(coordVec, &coords);
    457:   /* bubble sort closure vertices by coordinates in revlex order */
    458:   for (v = 0; v < nVerts; v++) revlexOrder[v] = v;
    459:   for (v = 0; v < nVerts; v++) {
    460:     for (w = v + 1; w < nVerts; w++) {
    461:       const PetscScalar *cv = &coords[closureOrder[revlexOrder[v]] * dim];
    462:       const PetscScalar *cw = &coords[closureOrder[revlexOrder[w]] * dim];
    463:       PetscReal diff = 0;

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

    469:         revlexOrder[v] = revlexOrder[w];
    470:         revlexOrder[w] = swap;
    471:       }
    472:     }
    473:   }
    474:   VecRestoreArrayRead(coordVec, &coords);
    475:   PetscMalloc1(ni->nodeIdxDim * nVerts, &newNodeIdx);
    476:   /* reorder nodeIdx to be in closure order */
    477:   for (v = 0; v < nVerts; v++) {
    478:     for (d = 0; d < ni->nodeIdxDim; d++) {
    479:       newNodeIdx[revlexOrder[v] * ni->nodeIdxDim + d] = ni->nodeIdx[v * ni->nodeIdxDim + d];
    480:     }
    481:   }
    482:   PetscFree(ni->nodeIdx);
    483:   ni->nodeIdx = newNodeIdx;
    484:   ni->perm = invClosureOrder;
    485:   PetscFree(revlexOrder);
    486:   PetscFree(closureOrder);
    487:   return(0);
    488: }

    490: /* the coordinates of the simplex vertices are the corners of the barycentric simplex.
    491:  * When we stack them on top of each other in revlex order, they look like the identity matrix */
    492: static PetscErrorCode PetscLagNodeIndicesCreateSimplexVertices(DM dm, PetscLagNodeIndices *nodeIndices)
    493: {
    494:   PetscLagNodeIndices ni;
    495:   PetscInt       dim, d;

    500:   PetscNew(&ni);
    501:   DMGetDimension(dm, &dim);
    502:   ni->nodeIdxDim = dim + 1;
    503:   ni->nodeVecDim = 0;
    504:   ni->nNodes = dim + 1;
    505:   ni->refct = 1;
    506:   PetscCalloc1((dim + 1)*(dim + 1), &(ni->nodeIdx));
    507:   for (d = 0; d < dim + 1; d++) ni->nodeIdx[d*(dim + 2)] = 1;
    508:   PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_FALSE);
    509:   *nodeIndices = ni;
    510:   return(0);
    511: }

    513: /* A polytope that is a tensor product of a facet and a segment.
    514:  * We take whatever coordinate system was being used for the facet
    515:  * and we concatenate the barycentric coordinates for the vertices
    516:  * at the end of the segment, (1,0) and (0,1), to get a coordinate
    517:  * system for the tensor product element */
    518: static PetscErrorCode PetscLagNodeIndicesCreateTensorVertices(DM dm, PetscLagNodeIndices facetni, PetscLagNodeIndices *nodeIndices)
    519: {
    520:   PetscLagNodeIndices ni;
    521:   PetscInt       nodeIdxDim, subNodeIdxDim = facetni->nodeIdxDim;
    522:   PetscInt       nVerts, nSubVerts = facetni->nNodes;
    523:   PetscInt       dim, d, e, f, g;

    528:   PetscNew(&ni);
    529:   DMGetDimension(dm, &dim);
    530:   ni->nodeIdxDim = nodeIdxDim = subNodeIdxDim + 2;
    531:   ni->nodeVecDim = 0;
    532:   ni->nNodes = nVerts = 2 * nSubVerts;
    533:   ni->refct = 1;
    534:   PetscCalloc1(nodeIdxDim * nVerts, &(ni->nodeIdx));
    535:   for (f = 0, d = 0; d < 2; d++) {
    536:     for (e = 0; e < nSubVerts; e++, f++) {
    537:       for (g = 0; g < subNodeIdxDim; g++) {
    538:         ni->nodeIdx[f * nodeIdxDim + g] = facetni->nodeIdx[e * subNodeIdxDim + g];
    539:       }
    540:       ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim] = (1 - d);
    541:       ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim + 1] = d;
    542:     }
    543:   }
    544:   PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_TRUE);
    545:   *nodeIndices = ni;
    546:   return(0);
    547: }

    549: /* This helps us compute symmetries, and it also helps us compute coordinates for dofs that are being pushed
    550:  * forward from a boundary mesh point.
    551:  *
    552:  * Input:
    553:  *
    554:  * dm - the target reference cell where we want new coordinates and dof directions to be valid
    555:  * vert - the vertex coordinate system for the target reference cell
    556:  * p - the point in the target reference cell that the dofs are coming from
    557:  * vertp - the vertex coordinate system for p's reference cell
    558:  * ornt - the resulting coordinates and dof vectors will be for p under this orientation
    559:  * nodep - the node coordinates and dof vectors in p's reference cell
    560:  * formDegree - the form degree that the dofs transform as
    561:  *
    562:  * Output:
    563:  *
    564:  * pfNodeIdx - the node coordinates for p's dofs, in the dm reference cell, from the ornt perspective
    565:  * pfNodeVec - the node dof vectors for p's dofs, in the dm reference cell, from the ornt perspective
    566:  */
    567: static PetscErrorCode PetscLagNodeIndicesPushForward(DM dm, PetscLagNodeIndices vert, PetscInt p, PetscLagNodeIndices vertp, PetscLagNodeIndices nodep, PetscInt ornt, PetscInt formDegree, PetscInt pfNodeIdx[], PetscReal pfNodeVec[])
    568: {
    569:   PetscInt       *closureVerts;
    570:   PetscInt        closureSize = 0;
    571:   PetscInt       *closure = NULL;
    572:   PetscInt        dim, pdim, c, i, j, k, n, v, vStart, vEnd;
    573:   PetscInt        nSubVert = vertp->nNodes;
    574:   PetscInt        nodeIdxDim = vert->nodeIdxDim;
    575:   PetscInt        subNodeIdxDim = vertp->nodeIdxDim;
    576:   PetscInt        nNodes = nodep->nNodes;
    577:   const PetscInt  *vertIdx = vert->nodeIdx;
    578:   const PetscInt  *subVertIdx = vertp->nodeIdx;
    579:   const PetscInt  *nodeIdx = nodep->nodeIdx;
    580:   const PetscReal *nodeVec = nodep->nodeVec;
    581:   PetscReal       *J, *Jstar;
    582:   PetscReal       detJ;
    583:   PetscInt        depth, pdepth, Nk, pNk;
    584:   Vec             coordVec;
    585:   PetscScalar      *newCoords = NULL;
    586:   const PetscScalar *oldCoords = NULL;
    587:   PetscErrorCode  ierr;

    590:   DMGetDimension(dm, &dim);
    591:   DMPlexGetDepth(dm, &depth);
    592:   DMGetCoordinatesLocal(dm, &coordVec);
    593:   DMPlexGetPointDepth(dm, p, &pdepth);
    594:   pdim = pdepth != depth ? pdepth != 0 ? pdepth : 0 : dim;
    595:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
    596:   DMGetWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);
    597:   DMPlexGetTransitiveClosure_Internal(dm, p, ornt, PETSC_TRUE, &closureSize, &closure);
    598:   c = closureSize - nSubVert;
    599:   /* we want which cell closure indices the closure of this point corresponds to */
    600:   for (v = 0; v < nSubVert; v++) closureVerts[v] = vert->perm[closure[2 * (c + v)] - vStart];
    601:   DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
    602:   /* push forward indices */
    603:   for (i = 0; i < nodeIdxDim; i++) { /* for every component of the target index space */
    604:     /* check if this is a component that all vertices around this point have in common */
    605:     for (j = 1; j < nSubVert; j++) {
    606:       if (vertIdx[closureVerts[j] * nodeIdxDim + i] != vertIdx[closureVerts[0] * nodeIdxDim + i]) break;
    607:     }
    608:     if (j == nSubVert) { /* all vertices have this component in common, directly copy to output */
    609:       PetscInt val = vertIdx[closureVerts[0] * nodeIdxDim + i];
    610:       for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = val;
    611:     } else {
    612:       PetscInt subi = -1;
    613:       /* there must be a component in vertp that looks the same */
    614:       for (k = 0; k < subNodeIdxDim; k++) {
    615:         for (j = 0; j < nSubVert; j++) {
    616:           if (vertIdx[closureVerts[j] * nodeIdxDim + i] != subVertIdx[j * subNodeIdxDim + k]) break;
    617:         }
    618:         if (j == nSubVert) {
    619:           subi = k;
    620:           break;
    621:         }
    622:       }
    623:       if (subi < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not find matching coordinate\n");
    624:       /* that component in the vertp system becomes component i in the vert system for each dof */
    625:       for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = nodeIdx[n * subNodeIdxDim + subi];
    626:     }
    627:   }
    628:   /* push forward vectors */
    629:   DMGetWorkArray(dm, dim * dim, MPIU_REAL, &J);
    630:   if (ornt != 0) { /* temporarily change the coordinate vector so
    631:                       DMPlexComputeCellGeometryAffineFEM gives us the Jacobian we want */
    632:     PetscInt        closureSize2 = 0;
    633:     PetscInt       *closure2 = NULL;

    635:     DMPlexGetTransitiveClosure_Internal(dm, p, 0, PETSC_TRUE, &closureSize2, &closure2);
    636:     PetscMalloc1(dim * nSubVert, &newCoords);
    637:     VecGetArrayRead(coordVec, &oldCoords);
    638:     for (v = 0; v < nSubVert; v++) {
    639:       PetscInt d;
    640:       for (d = 0; d < dim; d++) {
    641:         newCoords[(closure2[2 * (c + v)] - vStart) * dim + d] = oldCoords[closureVerts[v] * dim + d];
    642:       }
    643:     }
    644:     VecRestoreArrayRead(coordVec, &oldCoords);
    645:     DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize2, &closure2);
    646:     VecPlaceArray(coordVec, newCoords);
    647:   }
    648:   DMPlexComputeCellGeometryAffineFEM(dm, p, NULL, J, NULL, &detJ);
    649:   if (ornt != 0) {
    650:     VecResetArray(coordVec);
    651:     PetscFree(newCoords);
    652:   }
    653:   DMRestoreWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);
    654:   /* compactify */
    655:   for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
    656:   /* We have the Jacobian mapping the point's reference cell to this reference cell:
    657:    * pulling back a function to the point and applying the dof is what we want,
    658:    * so we get the pullback matrix and multiply the dof by that matrix on the right */
    659:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
    660:   PetscDTBinomialInt(pdim, PetscAbsInt(formDegree), &pNk);
    661:   DMGetWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);
    662:   PetscDTAltVPullbackMatrix(pdim, dim, J, formDegree, Jstar);
    663:   for (n = 0; n < nNodes; n++) {
    664:     for (i = 0; i < Nk; i++) {
    665:       PetscReal val = 0.;
    666:       for (j = 0; j < pNk; j++) val += nodeVec[n * pNk + j] * Jstar[j * Nk + i];
    667:       pfNodeVec[n * Nk + i] = val;
    668:     }
    669:   }
    670:   DMRestoreWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);
    671:   DMRestoreWorkArray(dm, dim * dim, MPIU_REAL, &J);
    672:   return(0);
    673: }

    675: /* given to sets of nodes, take the tensor product, where the product of the dof indices is concatenation and the
    676:  * product of the dof vectors is the wedge product */
    677: static PetscErrorCode PetscLagNodeIndicesTensor(PetscLagNodeIndices tracei, PetscInt dimT, PetscInt kT, PetscLagNodeIndices fiberi, PetscInt dimF, PetscInt kF, PetscLagNodeIndices *nodeIndices)
    678: {
    679:   PetscInt       dim = dimT + dimF;
    680:   PetscInt       nodeIdxDim, nNodes;
    681:   PetscInt       formDegree = kT + kF;
    682:   PetscInt       Nk, NkT, NkF;
    683:   PetscInt       MkT, MkF;
    684:   PetscLagNodeIndices ni;
    685:   PetscInt       i, j, l;
    686:   PetscReal      *projF, *projT;
    687:   PetscReal      *projFstar, *projTstar;
    688:   PetscReal      *workF, *workF2, *workT, *workT2, *work, *work2;
    689:   PetscReal      *wedgeMat;
    690:   PetscReal      sign;

    694:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
    695:   PetscDTBinomialInt(dimT, PetscAbsInt(kT), &NkT);
    696:   PetscDTBinomialInt(dimF, PetscAbsInt(kF), &NkF);
    697:   PetscDTBinomialInt(dim, PetscAbsInt(kT), &MkT);
    698:   PetscDTBinomialInt(dim, PetscAbsInt(kF), &MkF);
    699:   PetscNew(&ni);
    700:   ni->nodeIdxDim = nodeIdxDim = tracei->nodeIdxDim + fiberi->nodeIdxDim;
    701:   ni->nodeVecDim = Nk;
    702:   ni->nNodes = nNodes = tracei->nNodes * fiberi->nNodes;
    703:   ni->refct = 1;
    704:   PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
    705:   /* first concatenate the indices */
    706:   for (l = 0, j = 0; j < fiberi->nNodes; j++) {
    707:     for (i = 0; i < tracei->nNodes; i++, l++) {
    708:       PetscInt m, n = 0;

    710:       for (m = 0; m < tracei->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = tracei->nodeIdx[i * tracei->nodeIdxDim + m];
    711:       for (m = 0; m < fiberi->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = fiberi->nodeIdx[j * fiberi->nodeIdxDim + m];
    712:     }
    713:   }

    715:   /* now wedge together the push-forward vectors */
    716:   PetscMalloc1(nNodes * Nk, &(ni->nodeVec));
    717:   PetscCalloc2(dimT*dim, &projT, dimF*dim, &projF);
    718:   for (i = 0; i < dimT; i++) projT[i * (dim + 1)] = 1.;
    719:   for (i = 0; i < dimF; i++) projF[i * (dim + dimT + 1) + dimT] = 1.;
    720:   PetscMalloc2(MkT*NkT, &projTstar, MkF*NkF, &projFstar);
    721:   PetscDTAltVPullbackMatrix(dim, dimT, projT, kT, projTstar);
    722:   PetscDTAltVPullbackMatrix(dim, dimF, projF, kF, projFstar);
    723:   PetscMalloc6(MkT, &workT, MkT, &workT2, MkF, &workF, MkF, &workF2, Nk, &work, Nk, &work2);
    724:   PetscMalloc1(Nk * MkT, &wedgeMat);
    725:   sign = (PetscAbsInt(kT * kF) & 1) ? -1. : 1.;
    726:   for (l = 0, j = 0; j < fiberi->nNodes; j++) {
    727:     PetscInt d, e;

    729:     /* push forward fiber k-form */
    730:     for (d = 0; d < MkF; d++) {
    731:       PetscReal val = 0.;
    732:       for (e = 0; e < NkF; e++) val += projFstar[d * NkF + e] * fiberi->nodeVec[j * NkF + e];
    733:       workF[d] = val;
    734:     }
    735:     /* Hodge star to proper form if necessary */
    736:     if (kF < 0) {
    737:       for (d = 0; d < MkF; d++) workF2[d] = workF[d];
    738:       PetscDTAltVStar(dim, PetscAbsInt(kF), 1, workF2, workF);
    739:     }
    740:     /* Compute the matrix that wedges this form with one of the trace k-form */
    741:     PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kF), PetscAbsInt(kT), workF, wedgeMat);
    742:     for (i = 0; i < tracei->nNodes; i++, l++) {
    743:       /* push forward trace k-form */
    744:       for (d = 0; d < MkT; d++) {
    745:         PetscReal val = 0.;
    746:         for (e = 0; e < NkT; e++) val += projTstar[d * NkT + e] * tracei->nodeVec[i * NkT + e];
    747:         workT[d] = val;
    748:       }
    749:       /* Hodge star to proper form if necessary */
    750:       if (kT < 0) {
    751:         for (d = 0; d < MkT; d++) workT2[d] = workT[d];
    752:         PetscDTAltVStar(dim, PetscAbsInt(kT), 1, workT2, workT);
    753:       }
    754:       /* compute the wedge product of the push-forward trace form and firer forms */
    755:       for (d = 0; d < Nk; d++) {
    756:         PetscReal val = 0.;
    757:         for (e = 0; e < MkT; e++) val += wedgeMat[d * MkT + e] * workT[e];
    758:         work[d] = val;
    759:       }
    760:       /* inverse Hodge star from proper form if necessary */
    761:       if (formDegree < 0) {
    762:         for (d = 0; d < Nk; d++) work2[d] = work[d];
    763:         PetscDTAltVStar(dim, PetscAbsInt(formDegree), -1, work2, work);
    764:       }
    765:       /* insert into the array (adjusting for sign) */
    766:       for (d = 0; d < Nk; d++) ni->nodeVec[l * Nk + d] = sign * work[d];
    767:     }
    768:   }
    769:   PetscFree(wedgeMat);
    770:   PetscFree6(workT, workT2, workF, workF2, work, work2);
    771:   PetscFree2(projTstar, projFstar);
    772:   PetscFree2(projT, projF);
    773:   *nodeIndices = ni;
    774:   return(0);
    775: }

    777: /* simple union of two sets of nodes */
    778: static PetscErrorCode PetscLagNodeIndicesMerge(PetscLagNodeIndices niA, PetscLagNodeIndices niB, PetscLagNodeIndices *nodeIndices)
    779: {
    780:   PetscLagNodeIndices ni;
    781:   PetscInt            nodeIdxDim, nodeVecDim, nNodes;
    782:   PetscErrorCode      ierr;

    785:   PetscNew(&ni);
    786:   ni->nodeIdxDim = nodeIdxDim = niA->nodeIdxDim;
    787:   if (niB->nodeIdxDim != nodeIdxDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeIdxDim");
    788:   ni->nodeVecDim = nodeVecDim = niA->nodeVecDim;
    789:   if (niB->nodeVecDim != nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeVecDim");
    790:   ni->nNodes = nNodes = niA->nNodes + niB->nNodes;
    791:   ni->refct = 1;
    792:   PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
    793:   PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
    794:   PetscArraycpy(ni->nodeIdx, niA->nodeIdx, niA->nNodes * nodeIdxDim);
    795:   PetscArraycpy(ni->nodeVec, niA->nodeVec, niA->nNodes * nodeVecDim);
    796:   PetscArraycpy(&(ni->nodeIdx[niA->nNodes * nodeIdxDim]), niB->nodeIdx, niB->nNodes * nodeIdxDim);
    797:   PetscArraycpy(&(ni->nodeVec[niA->nNodes * nodeVecDim]), niB->nodeVec, niB->nNodes * nodeVecDim);
    798:   *nodeIndices = ni;
    799:   return(0);
    800: }

    802: #define PETSCTUPINTCOMPREVLEX(N)                                   \
    803: static int PetscTupIntCompRevlex_##N(const void *a, const void *b) \
    804: {                                                                  \
    805:   const PetscInt *A = (const PetscInt *) a;                        \
    806:   const PetscInt *B = (const PetscInt *) b;                        \
    807:   int i;                                                           \
    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: PETSCTUPINTCOMPREVLEX(3)
    817: PETSCTUPINTCOMPREVLEX(4)
    818: PETSCTUPINTCOMPREVLEX(5)
    819: PETSCTUPINTCOMPREVLEX(6)
    820: PETSCTUPINTCOMPREVLEX(7)

    822: static int PetscTupIntCompRevlex_N(const void *a, const void *b)
    823: {
    824:   const PetscInt *A = (const PetscInt *) a;
    825:   const PetscInt *B = (const PetscInt *) b;
    826:   int i;
    827:   int N = A[0];
    828:   PetscInt diff = 0;
    829:   for (i = 0; i < N; i++) {
    830:     diff = A[N - i] - B[N - i];
    831:     if (diff) break;
    832:   }
    833:   return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1;
    834: }

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

    843:   if (!(ni->perm)) {
    844:     PetscInt *sorter;
    845:     PetscInt m = ni->nNodes;
    846:     PetscInt nodeIdxDim = ni->nodeIdxDim;
    847:     PetscInt i, j, k, l;
    848:     PetscInt *prm;
    849:     int (*comp) (const void *, const void *);

    851:     PetscMalloc1((nodeIdxDim + 2) * m, &sorter);
    852:     for (k = 0, l = 0, i = 0; i < m; i++) {
    853:       sorter[k++] = nodeIdxDim + 1;
    854:       sorter[k++] = i;
    855:       for (j = 0; j < nodeIdxDim; j++) sorter[k++] = ni->nodeIdx[l++];
    856:     }
    857:     switch (nodeIdxDim) {
    858:     case 2:
    859:       comp = PetscTupIntCompRevlex_3;
    860:       break;
    861:     case 3:
    862:       comp = PetscTupIntCompRevlex_4;
    863:       break;
    864:     case 4:
    865:       comp = PetscTupIntCompRevlex_5;
    866:       break;
    867:     case 5:
    868:       comp = PetscTupIntCompRevlex_6;
    869:       break;
    870:     case 6:
    871:       comp = PetscTupIntCompRevlex_7;
    872:       break;
    873:     default:
    874:       comp = PetscTupIntCompRevlex_N;
    875:       break;
    876:     }
    877:     qsort(sorter, m, (nodeIdxDim + 2) * sizeof(PetscInt), comp);
    878:     PetscMalloc1(m, &prm);
    879:     for (i = 0; i < m; i++) prm[i] = sorter[(nodeIdxDim + 2) * i + 1];
    880:     ni->perm = prm;
    881:     PetscFree(sorter);
    882:   }
    883:   *perm = ni->perm;
    884:   return(0);
    885: }

    887: static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
    888: {
    889:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
    890:   PetscErrorCode      ierr;

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

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

    899:       for (i = 0; i < lag->numSelfSym; i++) {
    900:         PetscFree(allocated[i]);
    901:       }
    902:       PetscFree(allocated);
    903:     }
    904:     PetscFree(lag->symperms);
    905:   }
    906:   if (lag->symflips) {
    907:     PetscScalar **selfSyms = lag->symflips[0];

    909:     if (selfSyms) {
    910:       PetscInt i;
    911:       PetscScalar **allocated = &selfSyms[-lag->selfSymOff];

    913:       for (i = 0; i < lag->numSelfSym; i++) {
    914:         PetscFree(allocated[i]);
    915:       }
    916:       PetscFree(allocated);
    917:     }
    918:     PetscFree(lag->symflips);
    919:   }
    920:   Petsc1DNodeFamilyDestroy(&(lag->nodeFamily));
    921:   PetscLagNodeIndicesDestroy(&(lag->vertIndices));
    922:   PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));
    923:   PetscLagNodeIndicesDestroy(&(lag->allNodeIndices));
    924:   PetscFree(lag);
    925:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
    926:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
    927:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);
    928:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);
    929:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL);
    930:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL);
    931:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL);
    932:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL);
    933:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetUseMoments_C", NULL);
    934:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetUseMoments_C", NULL);
    935:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetMomentOrder_C", NULL);
    936:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetMomentOrder_C", NULL);
    937:   return(0);
    938: }

    940: static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
    941: {
    942:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
    943:   PetscErrorCode      ierr;

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

    950: static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
    951: {
    952:   PetscBool      iascii;

    958:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
    959:   if (iascii) {PetscDualSpaceLagrangeView_Ascii(sp, viewer);}
    960:   return(0);
    961: }

    963: static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
    964: {
    965:   PetscBool      continuous, tensor, trimmed, flg, flg2, flg3;
    966:   PetscDTNodeType nodeType;
    967:   PetscReal      nodeExponent;
    968:   PetscInt       momentOrder;
    969:   PetscBool      nodeEndpoints, useMoments;

    973:   PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
    974:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
    975:   PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
    976:   PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &nodeEndpoints, &nodeExponent);
    977:   if (nodeType == PETSCDTNODES_DEFAULT) nodeType = PETSCDTNODES_GAUSSJACOBI;
    978:   PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments);
    979:   PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder);
    980:   PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
    981:   PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
    982:   if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
    983:   PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetTensor", tensor, &tensor, &flg);
    984:   if (flg) {PetscDualSpaceLagrangeSetTensor(sp, tensor);}
    985:   PetscOptionsBool("-petscdualspace_lagrange_trimmed", "Flag for trimmed dual space", "PetscDualSpaceLagrangeSetTrimmed", trimmed, &trimmed, &flg);
    986:   if (flg) {PetscDualSpaceLagrangeSetTrimmed(sp, trimmed);}
    987:   PetscOptionsEnum("-petscdualspace_lagrange_node_type", "Lagrange node location type", "PetscDualSpaceLagrangeSetNodeType", PetscDTNodeTypes, (PetscEnum)nodeType, (PetscEnum *)&nodeType, &flg);
    988:   PetscOptionsBool("-petscdualspace_lagrange_node_endpoints", "Flag for nodes that include endpoints", "PetscDualSpaceLagrangeSetNodeType", nodeEndpoints, &nodeEndpoints, &flg2);
    989:   flg3 = PETSC_FALSE;
    990:   if (nodeType == PETSCDTNODES_GAUSSJACOBI) {
    991:     PetscOptionsReal("-petscdualspace_lagrange_node_exponent", "Gauss-Jacobi weight function exponent", "PetscDualSpaceLagrangeSetNodeType", nodeExponent, &nodeExponent, &flg3);
    992:   }
    993:   if (flg || flg2 || flg3) {PetscDualSpaceLagrangeSetNodeType(sp, nodeType, nodeEndpoints, nodeExponent);}
    994:   PetscOptionsBool("-petscdualspace_lagrange_use_moments", "Use moments (where appropriate) for functionals", "PetscDualSpaceLagrangeSetUseMoments", useMoments, &useMoments, &flg);
    995:   if (flg) {PetscDualSpaceLagrangeSetUseMoments(sp, useMoments);}
    996:   PetscOptionsInt("-petscdualspace_lagrange_moment_order", "Quadrature order for moment functionals", "PetscDualSpaceLagrangeSetMomentOrder", momentOrder, &momentOrder, &flg);
    997:   if (flg) {PetscDualSpaceLagrangeSetMomentOrder(sp, momentOrder);}
    998:   PetscOptionsTail();
    999:   return(0);
    1000: }

    1002: static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace spNew)
    1003: {
    1004:   PetscBool           cont, tensor, trimmed, boundary;
    1005:   PetscDTNodeType     nodeType;
    1006:   PetscReal           exponent;
    1007:   PetscDualSpace_Lag *lag    = (PetscDualSpace_Lag *) sp->data;
    1008:   PetscErrorCode      ierr;

    1011:   PetscDualSpaceLagrangeGetContinuity(sp, &cont);
    1012:   PetscDualSpaceLagrangeSetContinuity(spNew, cont);
    1013:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
    1014:   PetscDualSpaceLagrangeSetTensor(spNew, tensor);
    1015:   PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
    1016:   PetscDualSpaceLagrangeSetTrimmed(spNew, trimmed);
    1017:   PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &boundary, &exponent);
    1018:   PetscDualSpaceLagrangeSetNodeType(spNew, nodeType, boundary, exponent);
    1019:   if (lag->nodeFamily) {
    1020:     PetscDualSpace_Lag *lagnew = (PetscDualSpace_Lag *) spNew->data;

    1022:     Petsc1DNodeFamilyReference(lag->nodeFamily);
    1023:     lagnew->nodeFamily = lag->nodeFamily;
    1024:   }
    1025:   return(0);
    1026: }

    1028: /* for making tensor product spaces: take a dual space and product a segment space that has all the same
    1029:  * specifications (trimmed, continuous, order, node set), except for the form degree */
    1030: static PetscErrorCode PetscDualSpaceCreateEdgeSubspace_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt k, PetscInt Nc, PetscBool interiorOnly, PetscDualSpace *bdsp)
    1031: {
    1032:   DM                 K;
    1033:   PetscDualSpace_Lag *newlag;
    1034:   PetscErrorCode     ierr;

    1037:   PetscDualSpaceDuplicate(sp,bdsp);
    1038:   PetscDualSpaceSetFormDegree(*bdsp, k);
    1039:   PetscDualSpaceCreateReferenceCell(*bdsp, 1, PETSC_TRUE, &K);
    1040:   PetscDualSpaceSetDM(*bdsp, K);
    1041:   DMDestroy(&K);
    1042:   PetscDualSpaceSetOrder(*bdsp, order);
    1043:   PetscDualSpaceSetNumComponents(*bdsp, Nc);
    1044:   newlag = (PetscDualSpace_Lag *) (*bdsp)->data;
    1045:   newlag->interiorOnly = interiorOnly;
    1046:   PetscDualSpaceSetUp(*bdsp);
    1047:   return(0);
    1048: }

    1050: /* just the points, weights aren't handled */
    1051: static PetscErrorCode PetscQuadratureCreateTensor(PetscQuadrature trace, PetscQuadrature fiber, PetscQuadrature *product)
    1052: {
    1053:   PetscInt         dimTrace, dimFiber;
    1054:   PetscInt         numPointsTrace, numPointsFiber;
    1055:   PetscInt         dim, numPoints;
    1056:   const PetscReal *pointsTrace;
    1057:   const PetscReal *pointsFiber;
    1058:   PetscReal       *points;
    1059:   PetscInt         i, j, k, p;
    1060:   PetscErrorCode   ierr;

    1063:   PetscQuadratureGetData(trace, &dimTrace, NULL, &numPointsTrace, &pointsTrace, NULL);
    1064:   PetscQuadratureGetData(fiber, &dimFiber, NULL, &numPointsFiber, &pointsFiber, NULL);
    1065:   dim = dimTrace + dimFiber;
    1066:   numPoints = numPointsFiber * numPointsTrace;
    1067:   PetscMalloc1(numPoints * dim, &points);
    1068:   for (p = 0, j = 0; j < numPointsFiber; j++) {
    1069:     for (i = 0; i < numPointsTrace; i++, p++) {
    1070:       for (k = 0; k < dimTrace; k++) points[p * dim +            k] = pointsTrace[i * dimTrace + k];
    1071:       for (k = 0; k < dimFiber; k++) points[p * dim + dimTrace + k] = pointsFiber[j * dimFiber + k];
    1072:     }
    1073:   }
    1074:   PetscQuadratureCreate(PETSC_COMM_SELF, product);
    1075:   PetscQuadratureSetData(*product, dim, 0, numPoints, points, NULL);
    1076:   return(0);
    1077: }

    1079: /* Kronecker tensor product where matrix is considered a matrix of k-forms, so that
    1080:  * the entries in the product matrix are wedge products of the entries in the original matrices */
    1081: static PetscErrorCode MatTensorAltV(Mat trace, Mat fiber, PetscInt dimTrace, PetscInt kTrace, PetscInt dimFiber, PetscInt kFiber, Mat *product)
    1082: {
    1083:   PetscInt mTrace, nTrace, mFiber, nFiber, m, n, k, i, j, l;
    1084:   PetscInt dim, NkTrace, NkFiber, Nk;
    1085:   PetscInt dT, dF;
    1086:   PetscInt *nnzTrace, *nnzFiber, *nnz;
    1087:   PetscInt iT, iF, jT, jF, il, jl;
    1088:   PetscReal *workT, *workT2, *workF, *workF2, *work, *workstar;
    1089:   PetscReal *projT, *projF;
    1090:   PetscReal *projTstar, *projFstar;
    1091:   PetscReal *wedgeMat;
    1092:   PetscReal sign;
    1093:   PetscScalar *workS;
    1094:   Mat prod;
    1095:   /* this produces dof groups that look like the identity */

    1099:   MatGetSize(trace, &mTrace, &nTrace);
    1100:   PetscDTBinomialInt(dimTrace, PetscAbsInt(kTrace), &NkTrace);
    1101:   if (nTrace % NkTrace) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point value space of trace matrix is not a multiple of k-form size");
    1102:   MatGetSize(fiber, &mFiber, &nFiber);
    1103:   PetscDTBinomialInt(dimFiber, PetscAbsInt(kFiber), &NkFiber);
    1104:   if (nFiber % NkFiber) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point value space of fiber matrix is not a multiple of k-form size");
    1105:   PetscMalloc2(mTrace, &nnzTrace, mFiber, &nnzFiber);
    1106:   for (i = 0; i < mTrace; i++) {
    1107:     MatGetRow(trace, i, &(nnzTrace[i]), NULL, NULL);
    1108:     if (nnzTrace[i] % NkTrace) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in trace matrix are not in k-form size blocks");
    1109:   }
    1110:   for (i = 0; i < mFiber; i++) {
    1111:     MatGetRow(fiber, i, &(nnzFiber[i]), NULL, NULL);
    1112:     if (nnzFiber[i] % NkFiber) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in fiber matrix are not in k-form size blocks");
    1113:   }
    1114:   dim = dimTrace + dimFiber;
    1115:   k = kFiber + kTrace;
    1116:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
    1117:   m = mTrace * mFiber;
    1118:   PetscMalloc1(m, &nnz);
    1119:   for (l = 0, j = 0; j < mFiber; j++) for (i = 0; i < mTrace; i++, l++) nnz[l] = (nnzTrace[i] / NkTrace) * (nnzFiber[j] / NkFiber) * Nk;
    1120:   n = (nTrace / NkTrace) * (nFiber / NkFiber) * Nk;
    1121:   MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &prod);
    1122:   PetscFree(nnz);
    1123:   PetscFree2(nnzTrace,nnzFiber);
    1124:   /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
    1125:   MatSetOption(prod, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
    1126:   /* compute pullbacks */
    1127:   PetscDTBinomialInt(dim, PetscAbsInt(kTrace), &dT);
    1128:   PetscDTBinomialInt(dim, PetscAbsInt(kFiber), &dF);
    1129:   PetscMalloc4(dimTrace * dim, &projT, dimFiber * dim, &projF, dT * NkTrace, &projTstar, dF * NkFiber, &projFstar);
    1130:   PetscArrayzero(projT, dimTrace * dim);
    1131:   for (i = 0; i < dimTrace; i++) projT[i * (dim + 1)] = 1.;
    1132:   PetscArrayzero(projF, dimFiber * dim);
    1133:   for (i = 0; i < dimFiber; i++) projF[i * (dim + 1) + dimTrace] = 1.;
    1134:   PetscDTAltVPullbackMatrix(dim, dimTrace, projT, kTrace, projTstar);
    1135:   PetscDTAltVPullbackMatrix(dim, dimFiber, projF, kFiber, projFstar);
    1136:   PetscMalloc5(dT, &workT, dF, &workF, Nk, &work, Nk, &workstar, Nk, &workS);
    1137:   PetscMalloc2(dT, &workT2, dF, &workF2);
    1138:   PetscMalloc1(Nk * dT, &wedgeMat);
    1139:   sign = (PetscAbsInt(kTrace * kFiber) & 1) ? -1. : 1.;
    1140:   for (i = 0, iF = 0; iF < mFiber; iF++) {
    1141:     PetscInt           ncolsF, nformsF;
    1142:     const PetscInt    *colsF;
    1143:     const PetscScalar *valsF;

    1145:     MatGetRow(fiber, iF, &ncolsF, &colsF, &valsF);
    1146:     nformsF = ncolsF / NkFiber;
    1147:     for (iT = 0; iT < mTrace; iT++, i++) {
    1148:       PetscInt           ncolsT, nformsT;
    1149:       const PetscInt    *colsT;
    1150:       const PetscScalar *valsT;

    1152:       MatGetRow(trace, iT, &ncolsT, &colsT, &valsT);
    1153:       nformsT = ncolsT / NkTrace;
    1154:       for (j = 0, jF = 0; jF < nformsF; jF++) {
    1155:         PetscInt colF = colsF[jF * NkFiber] / NkFiber;

    1157:         for (il = 0; il < dF; il++) {
    1158:           PetscReal val = 0.;
    1159:           for (jl = 0; jl < NkFiber; jl++) val += projFstar[il * NkFiber + jl] * PetscRealPart(valsF[jF * NkFiber + jl]);
    1160:           workF[il] = val;
    1161:         }
    1162:         if (kFiber < 0) {
    1163:           for (il = 0; il < dF; il++) workF2[il] = workF[il];
    1164:           PetscDTAltVStar(dim, PetscAbsInt(kFiber), 1, workF2, workF);
    1165:         }
    1166:         PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kFiber), PetscAbsInt(kTrace), workF, wedgeMat);
    1167:         for (jT = 0; jT < nformsT; jT++, j++) {
    1168:           PetscInt colT = colsT[jT * NkTrace] / NkTrace;
    1169:           PetscInt col = colF * (nTrace / NkTrace) + colT;
    1170:           const PetscScalar *vals;

    1172:           for (il = 0; il < dT; il++) {
    1173:             PetscReal val = 0.;
    1174:             for (jl = 0; jl < NkTrace; jl++) val += projTstar[il * NkTrace + jl] * PetscRealPart(valsT[jT * NkTrace + jl]);
    1175:             workT[il] = val;
    1176:           }
    1177:           if (kTrace < 0) {
    1178:             for (il = 0; il < dT; il++) workT2[il] = workT[il];
    1179:             PetscDTAltVStar(dim, PetscAbsInt(kTrace), 1, workT2, workT);
    1180:           }

    1182:           for (il = 0; il < Nk; il++) {
    1183:             PetscReal val = 0.;
    1184:             for (jl = 0; jl < dT; jl++) val += sign * wedgeMat[il * dT + jl] * workT[jl];
    1185:             work[il] = val;
    1186:           }
    1187:           if (k < 0) {
    1188:             PetscDTAltVStar(dim, PetscAbsInt(k), -1, work, workstar);
    1189: #if defined(PETSC_USE_COMPLEX)
    1190:             for (l = 0; l < Nk; l++) workS[l] = workstar[l];
    1191:             vals = &workS[0];
    1192: #else
    1193:             vals = &workstar[0];
    1194: #endif
    1195:           } else {
    1196: #if defined(PETSC_USE_COMPLEX)
    1197:             for (l = 0; l < Nk; l++) workS[l] = work[l];
    1198:             vals = &workS[0];
    1199: #else
    1200:             vals = &work[0];
    1201: #endif
    1202:           }
    1203:           for (l = 0; l < Nk; l++) {
    1204:             MatSetValue(prod, i, col * Nk + l, vals[l], INSERT_VALUES);
    1205:           } /* Nk */
    1206:         } /* jT */
    1207:       } /* jF */
    1208:       MatRestoreRow(trace, iT, &ncolsT, &colsT, &valsT);
    1209:     } /* iT */
    1210:     MatRestoreRow(fiber, iF, &ncolsF, &colsF, &valsF);
    1211:   } /* iF */
    1212:   PetscFree(wedgeMat);
    1213:   PetscFree4(projT, projF, projTstar, projFstar);
    1214:   PetscFree2(workT2, workF2);
    1215:   PetscFree5(workT, workF, work, workstar, workS);
    1216:   MatAssemblyBegin(prod, MAT_FINAL_ASSEMBLY);
    1217:   MatAssemblyEnd(prod, MAT_FINAL_ASSEMBLY);
    1218:   *product = prod;
    1219:   return(0);
    1220: }

    1222: /* Union of quadrature points, with an attempt to identify commont points in the two sets */
    1223: static PetscErrorCode PetscQuadraturePointsMerge(PetscQuadrature quadA, PetscQuadrature quadB, PetscQuadrature *quadJoint, PetscInt *aToJoint[], PetscInt *bToJoint[])
    1224: {
    1225:   PetscInt         dimA, dimB;
    1226:   PetscInt         nA, nB, nJoint, i, j, d;
    1227:   const PetscReal *pointsA;
    1228:   const PetscReal *pointsB;
    1229:   PetscReal       *pointsJoint;
    1230:   PetscInt        *aToJ, *bToJ;
    1231:   PetscQuadrature  qJ;
    1232:   PetscErrorCode   ierr;

    1235:   PetscQuadratureGetData(quadA, &dimA, NULL, &nA, &pointsA, NULL);
    1236:   PetscQuadratureGetData(quadB, &dimB, NULL, &nB, &pointsB, NULL);
    1237:   if (dimA != dimB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Quadrature points must be in the same dimension");
    1238:   nJoint = nA;
    1239:   PetscMalloc1(nA, &aToJ);
    1240:   for (i = 0; i < nA; i++) aToJ[i] = i;
    1241:   PetscMalloc1(nB, &bToJ);
    1242:   for (i = 0; i < nB; i++) {
    1243:     for (j = 0; j < nA; j++) {
    1244:       bToJ[i] = -1;
    1245:       for (d = 0; d < dimA; d++) if (PetscAbsReal(pointsB[i * dimA + d] - pointsA[j * dimA + d]) > PETSC_SMALL) break;
    1246:       if (d == dimA) {
    1247:         bToJ[i] = j;
    1248:         break;
    1249:       }
    1250:     }
    1251:     if (bToJ[i] == -1) {
    1252:       bToJ[i] = nJoint++;
    1253:     }
    1254:   }
    1255:   *aToJoint = aToJ;
    1256:   *bToJoint = bToJ;
    1257:   PetscMalloc1(nJoint * dimA, &pointsJoint);
    1258:   PetscArraycpy(pointsJoint, pointsA, nA * dimA);
    1259:   for (i = 0; i < nB; i++) {
    1260:     if (bToJ[i] >= nA) {
    1261:       for (d = 0; d < dimA; d++) pointsJoint[bToJ[i] * dimA + d] = pointsB[i * dimA + d];
    1262:     }
    1263:   }
    1264:   PetscQuadratureCreate(PETSC_COMM_SELF, &qJ);
    1265:   PetscQuadratureSetData(qJ, dimA, 0, nJoint, pointsJoint, NULL);
    1266:   *quadJoint = qJ;
    1267:   return(0);
    1268: }

    1270: /* Matrices matA and matB are both quadrature -> dof matrices: produce a matrix that is joint quadrature -> union of
    1271:  * dofs, where the joint quadrature was produced by PetscQuadraturePointsMerge */
    1272: static PetscErrorCode MatricesMerge(Mat matA, Mat matB, PetscInt dim, PetscInt k, PetscInt numMerged, const PetscInt aToMerged[], const PetscInt bToMerged[], Mat *matMerged)
    1273: {
    1274:   PetscInt m, n, mA, nA, mB, nB, Nk, i, j, l;
    1275:   Mat      M;
    1276:   PetscInt *nnz;
    1277:   PetscInt maxnnz;
    1278:   PetscInt *work;

    1282:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
    1283:   MatGetSize(matA, &mA, &nA);
    1284:   if (nA % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matA column space not a multiple of k-form size");
    1285:   MatGetSize(matB, &mB, &nB);
    1286:   if (nB % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matB column space not a multiple of k-form size");
    1287:   m = mA + mB;
    1288:   n = numMerged * Nk;
    1289:   PetscMalloc1(m, &nnz);
    1290:   maxnnz = 0;
    1291:   for (i = 0; i < mA; i++) {
    1292:     MatGetRow(matA, i, &(nnz[i]), NULL, NULL);
    1293:     if (nnz[i] % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matA are not in k-form size blocks");
    1294:     maxnnz = PetscMax(maxnnz, nnz[i]);
    1295:   }
    1296:   for (i = 0; i < mB; i++) {
    1297:     MatGetRow(matB, i, &(nnz[i+mA]), NULL, NULL);
    1298:     if (nnz[i+mA] % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matB are not in k-form size blocks");
    1299:     maxnnz = PetscMax(maxnnz, nnz[i+mA]);
    1300:   }
    1301:   MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &M);
    1302:   PetscFree(nnz);
    1303:   /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
    1304:   MatSetOption(M, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
    1305:   PetscMalloc1(maxnnz, &work);
    1306:   for (i = 0; i < mA; i++) {
    1307:     const PetscInt *cols;
    1308:     const PetscScalar *vals;
    1309:     PetscInt nCols;
    1310:     MatGetRow(matA, i, &nCols, &cols, &vals);
    1311:     for (j = 0; j < nCols / Nk; j++) {
    1312:       PetscInt newCol = aToMerged[cols[j * Nk] / Nk];
    1313:       for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
    1314:     }
    1315:     MatSetValuesBlocked(M, 1, &i, nCols, work, vals, INSERT_VALUES);
    1316:     MatRestoreRow(matA, i, &nCols, &cols, &vals);
    1317:   }
    1318:   for (i = 0; i < mB; i++) {
    1319:     const PetscInt *cols;
    1320:     const PetscScalar *vals;

    1322:     PetscInt row = i + mA;
    1323:     PetscInt nCols;
    1324:     MatGetRow(matB, i, &nCols, &cols, &vals);
    1325:     for (j = 0; j < nCols / Nk; j++) {
    1326:       PetscInt newCol = bToMerged[cols[j * Nk] / Nk];
    1327:       for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
    1328:     }
    1329:     MatSetValuesBlocked(M, 1, &row, nCols, work, vals, INSERT_VALUES);
    1330:     MatRestoreRow(matB, i, &nCols, &cols, &vals);
    1331:   }
    1332:   PetscFree(work);
    1333:   MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
    1334:   MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
    1335:   *matMerged = M;
    1336:   return(0);
    1337: }

    1339: /* Take a dual space and product a segment space that has all the same specifications (trimmed, continuous, order,
    1340:  * node set), except for the form degree.  For computing boundary dofs and for making tensor product spaces */
    1341: static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Lagrange(PetscDualSpace sp, DM K, PetscInt f, PetscInt k, PetscInt Ncopies, PetscBool interiorOnly, PetscDualSpace *bdsp)
    1342: {
    1343:   PetscInt           Nknew, Ncnew;
    1344:   PetscInt           dim, pointDim = -1;
    1345:   PetscInt           depth;
    1346:   DM                 dm;
    1347:   PetscDualSpace_Lag *newlag;
    1348:   PetscErrorCode     ierr;

    1351:   PetscDualSpaceGetDM(sp,&dm);
    1352:   DMGetDimension(dm,&dim);
    1353:   DMPlexGetDepth(dm,&depth);
    1354:   PetscDualSpaceDuplicate(sp,bdsp);
    1355:   PetscDualSpaceSetFormDegree(*bdsp,k);
    1356:   if (!K) {
    1357:     PetscBool isSimplex;

    1360:     if (depth == dim) {
    1361:       PetscInt coneSize;

    1363:       pointDim = dim - 1;
    1364:       DMPlexGetConeSize(dm,f,&coneSize);
    1365:       isSimplex = (PetscBool) (coneSize == dim);
    1366:       PetscDualSpaceCreateReferenceCell(*bdsp, dim-1, isSimplex, &K);
    1367:     } else if (depth == 1) {
    1368:       pointDim = 0;
    1369:       PetscDualSpaceCreateReferenceCell(*bdsp, 0, PETSC_TRUE, &K);
    1370:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
    1371:   } else {
    1372:     PetscObjectReference((PetscObject)K);
    1373:     DMGetDimension(K, &pointDim);
    1374:   }
    1375:   PetscDualSpaceSetDM(*bdsp, K);
    1376:   DMDestroy(&K);
    1377:   PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew);
    1378:   Ncnew = Nknew * Ncopies;
    1379:   PetscDualSpaceSetNumComponents(*bdsp, Ncnew);
    1380:   newlag = (PetscDualSpace_Lag *) (*bdsp)->data;
    1381:   newlag->interiorOnly = interiorOnly;
    1382:   PetscDualSpaceSetUp(*bdsp);
    1383:   return(0);
    1384: }

    1386: /* Construct simplex nodes from a nodefamily, add Nk dof vectors of length Nk at each node.
    1387:  * Return the (quadrature, matrix) form of the dofs and the nodeIndices form as well.
    1388:  *
    1389:  * Sometimes we want a set of nodes to be contained in the interior of the element,
    1390:  * even when the node scheme puts nodes on the boundaries.  numNodeSkip tells
    1391:  * the routine how many "layers" of nodes need to be skipped.
    1392:  * */
    1393: static PetscErrorCode PetscDualSpaceLagrangeCreateSimplexNodeMat(Petsc1DNodeFamily nodeFamily, PetscInt dim, PetscInt sum, PetscInt Nk, PetscInt numNodeSkip, PetscQuadrature *iNodes, Mat *iMat, PetscLagNodeIndices *nodeIndices)
    1394: {
    1395:   PetscReal *extraNodeCoords, *nodeCoords;
    1396:   PetscInt nNodes, nExtraNodes;
    1397:   PetscInt i, j, k, extraSum = sum + numNodeSkip * (1 + dim);
    1398:   PetscQuadrature intNodes;
    1399:   Mat intMat;
    1400:   PetscLagNodeIndices ni;

    1404:   PetscDTBinomialInt(dim + sum, dim, &nNodes);
    1405:   PetscDTBinomialInt(dim + extraSum, dim, &nExtraNodes);

    1407:   PetscMalloc1(dim * nExtraNodes, &extraNodeCoords);
    1408:   PetscNew(&ni);
    1409:   ni->nodeIdxDim = dim + 1;
    1410:   ni->nodeVecDim = Nk;
    1411:   ni->nNodes = nNodes * Nk;
    1412:   ni->refct = 1;
    1413:   PetscMalloc1(nNodes * Nk * (dim + 1), &(ni->nodeIdx));
    1414:   PetscMalloc1(nNodes * Nk * Nk, &(ni->nodeVec));
    1415:   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.;
    1416:   Petsc1DNodeFamilyComputeSimplexNodes(nodeFamily, dim, extraSum, extraNodeCoords);
    1417:   if (numNodeSkip) {
    1418:     PetscInt k;
    1419:     PetscInt *tup;

    1421:     PetscMalloc1(dim * nNodes, &nodeCoords);
    1422:     PetscMalloc1(dim + 1, &tup);
    1423:     for (k = 0; k < nNodes; k++) {
    1424:       PetscInt j, c;
    1425:       PetscInt index;

    1427:       PetscDTIndexToBary(dim + 1, sum, k, tup);
    1428:       for (j = 0; j < dim + 1; j++) tup[j] += numNodeSkip;
    1429:       for (c = 0; c < Nk; c++) {
    1430:         for (j = 0; j < dim + 1; j++) {
    1431:           ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
    1432:         }
    1433:       }
    1434:       PetscDTBaryToIndex(dim + 1, extraSum, tup, &index);
    1435:       for (j = 0; j < dim; j++) nodeCoords[k * dim + j] = extraNodeCoords[index * dim + j];
    1436:     }
    1437:     PetscFree(tup);
    1438:     PetscFree(extraNodeCoords);
    1439:   } else {
    1440:     PetscInt k;
    1441:     PetscInt *tup;

    1443:     nodeCoords = extraNodeCoords;
    1444:     PetscMalloc1(dim + 1, &tup);
    1445:     for (k = 0; k < nNodes; k++) {
    1446:       PetscInt j, c;

    1448:       PetscDTIndexToBary(dim + 1, sum, k, tup);
    1449:       for (c = 0; c < Nk; c++) {
    1450:         for (j = 0; j < dim + 1; j++) {
    1451:           /* barycentric indices can have zeros, but we don't want to push forward zeros because it makes it harder to
    1452:            * determine which nodes correspond to which under symmetries, so we increase by 1.  This is fine
    1453:            * because the nodeIdx coordinates don't have any meaning other than helping to identify symmetries */
    1454:           ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
    1455:         }
    1456:       }
    1457:     }
    1458:     PetscFree(tup);
    1459:   }
    1460:   PetscQuadratureCreate(PETSC_COMM_SELF, &intNodes);
    1461:   PetscQuadratureSetData(intNodes, dim, 0, nNodes, nodeCoords, NULL);
    1462:   MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes * Nk, nNodes * Nk, Nk, NULL, &intMat);
    1463:   MatSetOption(intMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);
    1464:   for (j = 0; j < nNodes * Nk; j++) {
    1465:     PetscInt rem = j % Nk;
    1466:     PetscInt a, aprev = j - rem;
    1467:     PetscInt anext = aprev + Nk;

    1469:     for (a = aprev; a < anext; a++) {
    1470:       MatSetValue(intMat, j, a, (a == j) ? 1. : 0., INSERT_VALUES);
    1471:     }
    1472:   }
    1473:   MatAssemblyBegin(intMat, MAT_FINAL_ASSEMBLY);
    1474:   MatAssemblyEnd(intMat, MAT_FINAL_ASSEMBLY);
    1475:   *iNodes = intNodes;
    1476:   *iMat = intMat;
    1477:   *nodeIndices = ni;
    1478:   return(0);
    1479: }

    1481: /* once the nodeIndices have been created for the interior of the reference cell, and for all of the boundary cells,
    1482:  * push forward the boudary dofs and concatenate them into the full node indices for the dual space */
    1483: static PetscErrorCode PetscDualSpaceLagrangeCreateAllNodeIdx(PetscDualSpace sp)
    1484: {
    1485:   DM             dm;
    1486:   PetscInt       dim, nDofs;
    1487:   PetscSection   section;
    1488:   PetscInt       pStart, pEnd, p;
    1489:   PetscInt       formDegree, Nk;
    1490:   PetscInt       nodeIdxDim, spintdim;
    1491:   PetscDualSpace_Lag *lag;
    1492:   PetscLagNodeIndices ni, verti;

    1496:   lag = (PetscDualSpace_Lag *) sp->data;
    1497:   verti = lag->vertIndices;
    1498:   PetscDualSpaceGetDM(sp, &dm);
    1499:   DMGetDimension(dm, &dim);
    1500:   PetscDualSpaceGetFormDegree(sp, &formDegree);
    1501:   PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
    1502:   PetscDualSpaceGetSection(sp, &section);
    1503:   PetscSectionGetStorageSize(section, &nDofs);
    1504:   PetscNew(&ni);
    1505:   ni->nodeIdxDim = nodeIdxDim = verti->nodeIdxDim;
    1506:   ni->nodeVecDim = Nk;
    1507:   ni->nNodes = nDofs;
    1508:   ni->refct = 1;
    1509:   PetscMalloc1(nodeIdxDim * nDofs, &(ni->nodeIdx));
    1510:   PetscMalloc1(Nk * nDofs, &(ni->nodeVec));
    1511:   DMPlexGetChart(dm, &pStart, &pEnd);
    1512:   PetscSectionGetDof(section, 0, &spintdim);
    1513:   if (spintdim) {
    1514:     PetscArraycpy(ni->nodeIdx, lag->intNodeIndices->nodeIdx, spintdim * nodeIdxDim);
    1515:     PetscArraycpy(ni->nodeVec, lag->intNodeIndices->nodeVec, spintdim * Nk);
    1516:   }
    1517:   for (p = pStart + 1; p < pEnd; p++) {
    1518:     PetscDualSpace psp = sp->pointSpaces[p];
    1519:     PetscDualSpace_Lag *plag;
    1520:     PetscInt dof, off;

    1522:     PetscSectionGetDof(section, p, &dof);
    1523:     if (!dof) continue;
    1524:     plag = (PetscDualSpace_Lag *) psp->data;
    1525:     PetscSectionGetOffset(section, p, &off);
    1526:     PetscLagNodeIndicesPushForward(dm, verti, p, plag->vertIndices, plag->intNodeIndices, 0, formDegree, &(ni->nodeIdx[off * nodeIdxDim]), &(ni->nodeVec[off * Nk]));
    1527:   }
    1528:   lag->allNodeIndices = ni;
    1529:   return(0);
    1530: }

    1532: /* once the (quadrature, Matrix) forms of the dofs have been created for the interior of the
    1533:  * reference cell and for the boundary cells, jk
    1534:  * push forward the boundary data and concatenate them into the full (quadrature, matrix) data
    1535:  * for the dual space */
    1536: static PetscErrorCode PetscDualSpaceCreateAllDataFromInteriorData(PetscDualSpace sp)
    1537: {
    1538:   DM               dm;
    1539:   PetscSection     section;
    1540:   PetscInt         pStart, pEnd, p, k, Nk, dim, Nc;
    1541:   PetscInt         nNodes;
    1542:   PetscInt         countNodes;
    1543:   Mat              allMat;
    1544:   PetscQuadrature  allNodes;
    1545:   PetscInt         nDofs;
    1546:   PetscInt         maxNzforms, j;
    1547:   PetscScalar      *work;
    1548:   PetscReal        *L, *J, *Jinv, *v0, *pv0;
    1549:   PetscInt         *iwork;
    1550:   PetscReal        *nodes;
    1551:   PetscErrorCode   ierr;

    1554:   PetscDualSpaceGetDM(sp, &dm);
    1555:   DMGetDimension(dm, &dim);
    1556:   PetscDualSpaceGetSection(sp, &section);
    1557:   PetscSectionGetStorageSize(section, &nDofs);
    1558:   DMPlexGetChart(dm, &pStart, &pEnd);
    1559:   PetscDualSpaceGetFormDegree(sp, &k);
    1560:   PetscDualSpaceGetNumComponents(sp, &Nc);
    1561:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
    1562:   for (p = pStart, nNodes = 0, maxNzforms = 0; p < pEnd; p++) {
    1563:     PetscDualSpace  psp;
    1564:     DM              pdm;
    1565:     PetscInt        pdim, pNk;
    1566:     PetscQuadrature intNodes;
    1567:     Mat intMat;

    1569:     PetscDualSpaceGetPointSubspace(sp, p, &psp);
    1570:     if (!psp) continue;
    1571:     PetscDualSpaceGetDM(psp, &pdm);
    1572:     DMGetDimension(pdm, &pdim);
    1573:     if (pdim < PetscAbsInt(k)) continue;
    1574:     PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
    1575:     PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
    1576:     if (intNodes) {
    1577:       PetscInt nNodesp;

    1579:       PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, NULL, NULL);
    1580:       nNodes += nNodesp;
    1581:     }
    1582:     if (intMat) {
    1583:       PetscInt maxNzsp;
    1584:       PetscInt maxNzformsp;

    1586:       MatSeqAIJGetMaxRowNonzeros(intMat, &maxNzsp);
    1587:       if (maxNzsp % pNk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms");
    1588:       maxNzformsp = maxNzsp / pNk;
    1589:       maxNzforms = PetscMax(maxNzforms, maxNzformsp);
    1590:     }
    1591:   }
    1592:   MatCreateSeqAIJ(PETSC_COMM_SELF, nDofs, nNodes * Nc, maxNzforms * Nk, NULL, &allMat);
    1593:   MatSetOption(allMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);
    1594:   PetscMalloc7(dim, &v0, dim, &pv0, dim * dim, &J, dim * dim, &Jinv, Nk * Nk, &L, maxNzforms * Nk, &work, maxNzforms * Nk, &iwork);
    1595:   for (j = 0; j < dim; j++) pv0[j] = -1.;
    1596:   PetscMalloc1(dim * nNodes, &nodes);
    1597:   for (p = pStart, countNodes = 0; p < pEnd; p++) {
    1598:     PetscDualSpace  psp;
    1599:     PetscQuadrature intNodes;
    1600:     DM pdm;
    1601:     PetscInt pdim, pNk;
    1602:     PetscInt countNodesIn = countNodes;
    1603:     PetscReal detJ;
    1604:     Mat intMat;

    1606:     PetscDualSpaceGetPointSubspace(sp, p, &psp);
    1607:     if (!psp) continue;
    1608:     PetscDualSpaceGetDM(psp, &pdm);
    1609:     DMGetDimension(pdm, &pdim);
    1610:     if (pdim < PetscAbsInt(k)) continue;
    1611:     PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
    1612:     if (intNodes == NULL && intMat == NULL) continue;
    1613:     PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
    1614:     if (p) {
    1615:       DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, Jinv, &detJ);
    1616:     } else { /* identity */
    1617:       PetscInt i,j;

    1619:       for (i = 0; i < dim; i++) for (j = 0; j < dim; j++) J[i * dim + j] = Jinv[i * dim + j] = 0.;
    1620:       for (i = 0; i < dim; i++) J[i * dim + i] = Jinv[i * dim + i] = 1.;
    1621:       for (i = 0; i < dim; i++) v0[i] = -1.;
    1622:     }
    1623:     if (pdim != dim) { /* compactify Jacobian */
    1624:       PetscInt i, j;

    1626:       for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
    1627:     }
    1628:     PetscDTAltVPullbackMatrix(pdim, dim, J, k, L);
    1629:     if (intNodes) { /* push forward quadrature locations by the affine transformation */
    1630:       PetscInt nNodesp;
    1631:       const PetscReal *nodesp;
    1632:       PetscInt j;

    1634:       PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, &nodesp, NULL);
    1635:       for (j = 0; j < nNodesp; j++, countNodes++) {
    1636:         PetscInt d, e;

    1638:         for (d = 0; d < dim; d++) {
    1639:           nodes[countNodes * dim + d] = v0[d];
    1640:           for (e = 0; e < pdim; e++) {
    1641:             nodes[countNodes * dim + d] += J[d * pdim + e] * (nodesp[j * pdim + e] - pv0[e]);
    1642:           }
    1643:         }
    1644:       }
    1645:     }
    1646:     if (intMat) {
    1647:       PetscInt nrows;
    1648:       PetscInt off;

    1650:       PetscSectionGetDof(section, p, &nrows);
    1651:       PetscSectionGetOffset(section, p, &off);
    1652:       for (j = 0; j < nrows; j++) {
    1653:         PetscInt ncols;
    1654:         const PetscInt *cols;
    1655:         const PetscScalar *vals;
    1656:         PetscInt l, d, e;
    1657:         PetscInt row = j + off;

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

    1664:           for (d = 0; d < pNk; d++) {
    1665:             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");
    1666:           }
    1667:           blockcol = cols[l * pNk] / pNk;
    1668:           for (d = 0; d < Nk; d++) {
    1669:             iwork[l * Nk + d] = (blockcol + countNodesIn) * Nk + d;
    1670:           }
    1671:           for (d = 0; d < Nk; d++) work[l * Nk + d] = 0.;
    1672:           for (d = 0; d < Nk; d++) {
    1673:             for (e = 0; e < pNk; e++) {
    1674:               /* "push forward" dof by pulling back a k-form to be evaluated on the point: multiply on the right by L */
    1675:               work[l * Nk + d] += vals[l * pNk + e] * L[e * Nk + d];
    1676:             }
    1677:           }
    1678:         }
    1679:         MatSetValues(allMat, 1, &row, (ncols / pNk) * Nk, iwork, work, INSERT_VALUES);
    1680:         MatRestoreRow(intMat, j, &ncols, &cols, &vals);
    1681:       }
    1682:     }
    1683:   }
    1684:   MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);
    1685:   MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);
    1686:   PetscQuadratureCreate(PETSC_COMM_SELF, &allNodes);
    1687:   PetscQuadratureSetData(allNodes, dim, 0, nNodes, nodes, NULL);
    1688:   PetscFree7(v0, pv0, J, Jinv, L, work, iwork);
    1689:   MatDestroy(&(sp->allMat));
    1690:   sp->allMat = allMat;
    1691:   PetscQuadratureDestroy(&(sp->allNodes));
    1692:   sp->allNodes = allNodes;
    1693:   return(0);
    1694: }

    1696: /* rather than trying to get all data from the functionals, we create
    1697:  * the functionals from rows of the quadrature -> dof matrix.
    1698:  *
    1699:  * Ideally most of the uses of PetscDualSpace in PetscFE will switch
    1700:  * to using intMat and allMat, so that the individual functionals
    1701:  * don't need to be constructed at all */
    1702: static PetscErrorCode PetscDualSpaceComputeFunctionalsFromAllData(PetscDualSpace sp)
    1703: {
    1704:   PetscQuadrature allNodes;
    1705:   Mat             allMat;
    1706:   PetscInt        nDofs;
    1707:   PetscInt        dim, k, Nk, Nc, f;
    1708:   DM              dm;
    1709:   PetscInt        nNodes, spdim;
    1710:   const PetscReal *nodes = NULL;
    1711:   PetscSection    section;
    1712:   PetscBool       useMoments;
    1713:   PetscErrorCode  ierr;

    1716:   PetscDualSpaceGetDM(sp, &dm);
    1717:   DMGetDimension(dm, &dim);
    1718:   PetscDualSpaceGetNumComponents(sp, &Nc);
    1719:   PetscDualSpaceGetFormDegree(sp, &k);
    1720:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
    1721:   PetscDualSpaceGetAllData(sp, &allNodes, &allMat);
    1722:   nNodes = 0;
    1723:   if (allNodes) {
    1724:     PetscQuadratureGetData(allNodes, NULL, NULL, &nNodes, &nodes, NULL);
    1725:   }
    1726:   MatGetSize(allMat, &nDofs, NULL);
    1727:   PetscDualSpaceGetSection(sp, &section);
    1728:   PetscSectionGetStorageSize(section, &spdim);
    1729:   if (spdim != nDofs) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "incompatible all matrix size");
    1730:   PetscMalloc1(nDofs, &(sp->functional));
    1731:   PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments);
    1732:   if (useMoments) {
    1733:     Mat              allMat;
    1734:     PetscInt         momentOrder, i;
    1735:     PetscBool        tensor;
    1736:     const PetscReal *weights;
    1737:     PetscScalar     *array;

    1739:     if (nDofs != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "We do not yet support moments beyond P0, nDofs == %D", nDofs);
    1740:     PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder);
    1741:     PetscDualSpaceLagrangeGetTensor(sp, &tensor);
    1742:     if (!tensor) {PetscDTStroudConicalQuadrature(dim, Nc, PetscMax(momentOrder + 1,1), -1.0, 1.0, &(sp->functional[0]));}
    1743:     else         {PetscDTGaussTensorQuadrature(dim, Nc, PetscMax(momentOrder + 1,1), -1.0, 1.0, &(sp->functional[0]));}
    1744:     /* Need to replace allNodes and allMat */
    1745:     PetscObjectReference((PetscObject) sp->functional[0]);
    1746:     PetscQuadratureDestroy(&(sp->allNodes));
    1747:     sp->allNodes = sp->functional[0];
    1748:     PetscQuadratureGetData(sp->allNodes, NULL, NULL, &nNodes, NULL, &weights);
    1749:     MatCreateSeqDense(PETSC_COMM_SELF, nDofs, nNodes * Nc, NULL, &allMat);
    1750:     MatDenseGetArrayWrite(allMat, &array);
    1751:     for (i = 0; i < nNodes * Nc; ++i) array[i] = weights[i];
    1752:     MatDenseRestoreArrayWrite(allMat, &array);
    1753:     MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);
    1754:     MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);
    1755:     MatDestroy(&(sp->allMat));
    1756:     sp->allMat = allMat;
    1757:     return(0);
    1758:   }
    1759:   for (f = 0; f < nDofs; f++) {
    1760:     PetscInt ncols, c;
    1761:     const PetscInt *cols;
    1762:     const PetscScalar *vals;
    1763:     PetscReal *nodesf;
    1764:     PetscReal *weightsf;
    1765:     PetscInt nNodesf;
    1766:     PetscInt countNodes;

    1768:     MatGetRow(allMat, f, &ncols, &cols, &vals);
    1769:     if (ncols % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "all matrix is not laid out as blocks of k-forms");
    1770:     for (c = 1, nNodesf = 1; c < ncols; c++) {
    1771:       if ((cols[c] / Nc) != (cols[c-1] / Nc)) nNodesf++;
    1772:     }
    1773:     PetscMalloc1(dim * nNodesf, &nodesf);
    1774:     PetscMalloc1(Nc * nNodesf, &weightsf);
    1775:     for (c = 0, countNodes = 0; c < ncols; c++) {
    1776:       if (!c || ((cols[c] / Nc) != (cols[c-1] / Nc))) {
    1777:         PetscInt d;

    1779:         for (d = 0; d < Nc; d++) {
    1780:           weightsf[countNodes * Nc + d] = 0.;
    1781:         }
    1782:         for (d = 0; d < dim; d++) {
    1783:           nodesf[countNodes * dim + d] = nodes[(cols[c] / Nc) * dim + d];
    1784:         }
    1785:         countNodes++;
    1786:       }
    1787:       weightsf[(countNodes - 1) * Nc + (cols[c] % Nc)] = PetscRealPart(vals[c]);
    1788:     }
    1789:     PetscQuadratureCreate(PETSC_COMM_SELF, &(sp->functional[f]));
    1790:     PetscQuadratureSetData(sp->functional[f], dim, Nc, nNodesf, nodesf, weightsf);
    1791:     MatRestoreRow(allMat, f, &ncols, &cols, &vals);
    1792:   }
    1793:   return(0);
    1794: }

    1796: /* take a matrix meant for k-forms and expand it to one for Ncopies */
    1797: static PetscErrorCode PetscDualSpaceLagrangeMatrixCreateCopies(Mat A, PetscInt Nk, PetscInt Ncopies, Mat *Abs)
    1798: {
    1799:   PetscInt       m, n, i, j, k;
    1800:   PetscInt       maxnnz, *nnz, *iwork;
    1801:   Mat            Ac;

    1805:   MatGetSize(A, &m, &n);
    1806:   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);
    1807:   PetscMalloc1(m * Ncopies, &nnz);
    1808:   for (i = 0, maxnnz = 0; i < m; i++) {
    1809:     PetscInt innz;
    1810:     MatGetRow(A, i, &innz, NULL, NULL);
    1811:     if (innz % Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "A row %D nnzs is not a multiple of Nk %D", innz, Nk);
    1812:     for (j = 0; j < Ncopies; j++) nnz[i * Ncopies + j] = innz;
    1813:     maxnnz = PetscMax(maxnnz, innz);
    1814:   }
    1815:   MatCreateSeqAIJ(PETSC_COMM_SELF, m * Ncopies, n * Ncopies, 0, nnz, &Ac);
    1816:   MatSetOption(Ac, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
    1817:   PetscFree(nnz);
    1818:   PetscMalloc1(maxnnz, &iwork);
    1819:   for (i = 0; i < m; i++) {
    1820:     PetscInt innz;
    1821:     const PetscInt    *cols;
    1822:     const PetscScalar *vals;

    1824:     MatGetRow(A, i, &innz, &cols, &vals);
    1825:     for (j = 0; j < innz; j++) iwork[j] = (cols[j] / Nk) * (Nk * Ncopies) + (cols[j] % Nk);
    1826:     for (j = 0; j < Ncopies; j++) {
    1827:       PetscInt row = i * Ncopies + j;

    1829:       MatSetValues(Ac, 1, &row, innz, iwork, vals, INSERT_VALUES);
    1830:       for (k = 0; k < innz; k++) iwork[k] += Nk;
    1831:     }
    1832:     MatRestoreRow(A, i, &innz, &cols, &vals);
    1833:   }
    1834:   PetscFree(iwork);
    1835:   MatAssemblyBegin(Ac, MAT_FINAL_ASSEMBLY);
    1836:   MatAssemblyEnd(Ac, MAT_FINAL_ASSEMBLY);
    1837:   *Abs = Ac;
    1838:   return(0);
    1839: }

    1841: /* check if a cell is a tensor product of the segment with a facet,
    1842:  * specifically checking if f and f2 can be the "endpoints" (like the triangles
    1843:  * at either end of a wedge) */
    1844: static PetscErrorCode DMPlexPointIsTensor_Internal_Given(DM dm, PetscInt p, PetscInt f, PetscInt f2, PetscBool *isTensor)
    1845: {
    1846:   PetscInt        coneSize, c;
    1847:   const PetscInt *cone;
    1848:   const PetscInt *fCone;
    1849:   const PetscInt *f2Cone;
    1850:   PetscInt        fs[2];
    1851:   PetscInt        meetSize, nmeet;
    1852:   const PetscInt *meet;
    1853:   PetscErrorCode  ierr;

    1856:   fs[0] = f;
    1857:   fs[1] = f2;
    1858:   DMPlexGetMeet(dm, 2, fs, &meetSize, &meet);
    1859:   nmeet = meetSize;
    1860:   DMPlexRestoreMeet(dm, 2, fs, &meetSize, &meet);
    1861:   /* two points that have a non-empty meet cannot be at opposite ends of a cell */
    1862:   if (nmeet) {
    1863:     *isTensor = PETSC_FALSE;
    1864:     return(0);
    1865:   }
    1866:   DMPlexGetConeSize(dm, p, &coneSize);
    1867:   DMPlexGetCone(dm, p, &cone);
    1868:   DMPlexGetCone(dm, f, &fCone);
    1869:   DMPlexGetCone(dm, f2, &f2Cone);
    1870:   for (c = 0; c < coneSize; c++) {
    1871:     PetscInt e, ef;
    1872:     PetscInt d = -1, d2 = -1;
    1873:     PetscInt dcount, d2count;
    1874:     PetscInt t = cone[c];
    1875:     PetscInt tConeSize;
    1876:     PetscBool tIsTensor;
    1877:     const PetscInt *tCone;

    1879:     if (t == f || t == f2) continue;
    1880:     /* for every other facet in the cone, check that is has
    1881:      * one ridge in common with each end */
    1882:     DMPlexGetConeSize(dm, t, &tConeSize);
    1883:     DMPlexGetCone(dm, t, &tCone);

    1885:     dcount = 0;
    1886:     d2count = 0;
    1887:     for (e = 0; e < tConeSize; e++) {
    1888:       PetscInt q = tCone[e];
    1889:       for (ef = 0; ef < coneSize - 2; ef++) {
    1890:         if (fCone[ef] == q) {
    1891:           if (dcount) {
    1892:             *isTensor = PETSC_FALSE;
    1893:             return(0);
    1894:           }
    1895:           d = q;
    1896:           dcount++;
    1897:         } else if (f2Cone[ef] == q) {
    1898:           if (d2count) {
    1899:             *isTensor = PETSC_FALSE;
    1900:             return(0);
    1901:           }
    1902:           d2 = q;
    1903:           d2count++;
    1904:         }
    1905:       }
    1906:     }
    1907:     /* if the whole cell is a tensor with the segment, then this
    1908:      * facet should be a tensor with the segment */
    1909:     DMPlexPointIsTensor_Internal_Given(dm, t, d, d2, &tIsTensor);
    1910:     if (!tIsTensor) {
    1911:       *isTensor = PETSC_FALSE;
    1912:       return(0);
    1913:     }
    1914:   }
    1915:   *isTensor = PETSC_TRUE;
    1916:   return(0);
    1917: }

    1919: /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair
    1920:  * that could be the opposite ends */
    1921: static PetscErrorCode DMPlexPointIsTensor_Internal(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB)
    1922: {
    1923:   PetscInt        coneSize, c, c2;
    1924:   const PetscInt *cone;
    1925:   PetscErrorCode  ierr;

    1928:   DMPlexGetConeSize(dm, p, &coneSize);
    1929:   if (!coneSize) {
    1930:     if (isTensor) *isTensor = PETSC_FALSE;
    1931:     if (endA) *endA = -1;
    1932:     if (endB) *endB = -1;
    1933:   }
    1934:   DMPlexGetCone(dm, p, &cone);
    1935:   for (c = 0; c < coneSize; c++) {
    1936:     PetscInt f = cone[c];
    1937:     PetscInt fConeSize;

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

    1942:     for (c2 = c + 1; c2 < coneSize; c2++) {
    1943:       PetscInt  f2 = cone[c2];
    1944:       PetscBool isTensorff2;
    1945:       PetscInt f2ConeSize;

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

    1950:       DMPlexPointIsTensor_Internal_Given(dm, p, f, f2, &isTensorff2);
    1951:       if (isTensorff2) {
    1952:         if (isTensor) *isTensor = PETSC_TRUE;
    1953:         if (endA) *endA = f;
    1954:         if (endB) *endB = f2;
    1955:         return(0);
    1956:       }
    1957:     }
    1958:   }
    1959:   if (isTensor) *isTensor = PETSC_FALSE;
    1960:   if (endA) *endA = -1;
    1961:   if (endB) *endB = -1;
    1962:   return(0);
    1963: }

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

    1973:   DMPlexIsInterpolated(dm, &interpolated);
    1974:   if (interpolated != DMPLEX_INTERPOLATED_FULL) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Only for interpolated DMPlex's");
    1975:   DMPlexPointIsTensor_Internal(dm, p, isTensor, endA, endB);
    1976:   return(0);
    1977: }

    1979: /* Let k = formDegree and k' = -sign(k) * dim + k.  Transform a symmetric frame for k-forms on the biunit simplex into
    1980:  * a symmetric frame for k'-forms on the biunit simplex.
    1981:  *
    1982:  * A frame is "symmetric" if the pullback of every symmetry of the biunit simplex is a permutation of the frame.
    1983:  *
    1984:  * forms in the symmetric frame are used as dofs in the untrimmed simplex spaces.  This way, symmetries of the
    1985:  * reference cell result in permutations of dofs grouped by node.
    1986:  *
    1987:  * Use T to transform dof matrices for k'-forms into dof matrices for k-forms as a block diagonal transformation on
    1988:  * the right.
    1989:  */
    1990: static PetscErrorCode BiunitSimplexSymmetricFormTransformation(PetscInt dim, PetscInt formDegree, PetscReal T[])
    1991: {
    1992:   PetscInt       k = formDegree;
    1993:   PetscInt       kd = k < 0 ? dim + k : k - dim;
    1994:   PetscInt       Nk;
    1995:   PetscReal      *biToEq, *eqToBi, *biToEqStar, *eqToBiStar;
    1996:   PetscInt       fact;

    2000:   PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
    2001:   PetscCalloc4(dim * dim, &biToEq, dim * dim, &eqToBi, Nk * Nk, &biToEqStar, Nk * Nk, &eqToBiStar);
    2002:   /* fill in biToEq: Jacobian of the transformation from the biunit simplex to the equilateral simplex */
    2003:   fact = 0;
    2004:   for (PetscInt i = 0; i < dim; i++) {
    2005:     biToEq[i * dim + i] = PetscSqrtReal(((PetscReal)i + 2.) / (2.*((PetscReal)i+1.)));
    2006:     fact += 4*(i+1);
    2007:     for (PetscInt j = i+1; j < dim; j++) {
    2008:       biToEq[i * dim + j] = PetscSqrtReal(1./(PetscReal)fact);
    2009:     }
    2010:   }
    2011:   /* fill in eqToBi: Jacobian of the transformation from the equilateral simplex to the biunit simplex */
    2012:   fact = 0;
    2013:   for (PetscInt j = 0; j < dim; j++) {
    2014:     eqToBi[j * dim + j] = PetscSqrtReal(2.*((PetscReal)j+1.)/((PetscReal)j+2));
    2015:     fact += j+1;
    2016:     for (PetscInt i = 0; i < j; i++) {
    2017:       eqToBi[i * dim + j] = -PetscSqrtReal(1./(PetscReal)fact);
    2018:     }
    2019:   }
    2020:   PetscDTAltVPullbackMatrix(dim, dim, biToEq, kd, biToEqStar);
    2021:   PetscDTAltVPullbackMatrix(dim, dim, eqToBi, k, eqToBiStar);
    2022:   /* product of pullbacks simulates the following steps
    2023:    *
    2024:    * 1. start with frame W = [w_1, w_2, ..., w_m] of k forms that is symmetric on the biunit simplex:
    2025:           if J is the Jacobian of a symmetry of the biunit simplex, then J_k* W = [J_k*w_1, ..., J_k*w_m]
    2026:           is a permutation of W.
    2027:           Even though a k' form --- a (dim - k) form represented by its Hodge star --- has the same geometric
    2028:           content as a k form, W is not a symmetric frame of k' forms on the biunit simplex.  That's because,
    2029:           for general Jacobian J, J_k* != J_k'*.
    2030:    * 2. pullback W to the equilateral triangle using the k pullback, W_eq = eqToBi_k* W.  All symmetries of the
    2031:           equilateral simplex have orthonormal Jacobians.  For an orthonormal Jacobian O, J_k* = J_k'*, so W_eq is
    2032:           also a symmetric frame for k' forms on the equilateral simplex.
    2033:      3. pullback W_eq back to the biunit simplex using the k' pulback, V = biToEq_k'* W_eq = biToEq_k'* eqToBi_k* W.
    2034:           V is a symmetric frame for k' forms on the biunit simplex.
    2035:    */
    2036:   for (PetscInt i = 0; i < Nk; i++) {
    2037:     for (PetscInt j = 0; j < Nk; j++) {
    2038:       PetscReal val = 0.;
    2039:       for (PetscInt k = 0; k < Nk; k++) val += biToEqStar[i * Nk + k] * eqToBiStar[k * Nk + j];
    2040:       T[i * Nk + j] = val;
    2041:     }
    2042:   }
    2043:   PetscFree4(biToEq, eqToBi, biToEqStar, eqToBiStar);
    2044:   return(0);
    2045: }

    2047: /* permute a quadrature -> dof matrix so that its rows are in revlex order by nodeIdx */
    2048: static PetscErrorCode MatPermuteByNodeIdx(Mat A, PetscLagNodeIndices ni, Mat *Aperm)
    2049: {
    2050:   PetscInt       m, n, i, j;
    2051:   PetscInt       nodeIdxDim = ni->nodeIdxDim;
    2052:   PetscInt       nodeVecDim = ni->nodeVecDim;
    2053:   PetscInt       *perm;
    2054:   IS             permIS;
    2055:   IS             id;
    2056:   PetscInt       *nIdxPerm;
    2057:   PetscReal      *nVecPerm;

    2061:   PetscLagNodeIndicesGetPermutation(ni, &perm);
    2062:   MatGetSize(A, &m, &n);
    2063:   PetscMalloc1(nodeIdxDim * m, &nIdxPerm);
    2064:   PetscMalloc1(nodeVecDim * m, &nVecPerm);
    2065:   for (i = 0; i < m; i++) for (j = 0; j < nodeIdxDim; j++) nIdxPerm[i * nodeIdxDim + j] = ni->nodeIdx[perm[i] * nodeIdxDim + j];
    2066:   for (i = 0; i < m; i++) for (j = 0; j < nodeVecDim; j++) nVecPerm[i * nodeVecDim + j] = ni->nodeVec[perm[i] * nodeVecDim + j];
    2067:   ISCreateGeneral(PETSC_COMM_SELF, m, perm, PETSC_USE_POINTER, &permIS);
    2068:   ISSetPermutation(permIS);
    2069:   ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &id);
    2070:   ISSetPermutation(id);
    2071:   MatPermute(A, permIS, id, Aperm);
    2072:   ISDestroy(&permIS);
    2073:   ISDestroy(&id);
    2074:   for (i = 0; i < m; i++) perm[i] = i;
    2075:   PetscFree(ni->nodeIdx);
    2076:   PetscFree(ni->nodeVec);
    2077:   ni->nodeIdx = nIdxPerm;
    2078:   ni->nodeVec = nVecPerm;
    2079:   return(0);
    2080: }

    2082: static PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
    2083: {
    2084:   PetscDualSpace_Lag *lag   = (PetscDualSpace_Lag *) sp->data;
    2085:   DM                  dm    = sp->dm;
    2086:   DM                  dmint = NULL;
    2087:   PetscInt            order;
    2088:   PetscInt            Nc    = sp->Nc;
    2089:   MPI_Comm            comm;
    2090:   PetscBool           continuous;
    2091:   PetscSection        section;
    2092:   PetscInt            depth, dim, pStart, pEnd, cStart, cEnd, p, *pStratStart, *pStratEnd, d;
    2093:   PetscInt            formDegree, Nk, Ncopies;
    2094:   PetscInt            tensorf = -1, tensorf2 = -1;
    2095:   PetscBool           tensorCell, tensorSpace;
    2096:   PetscBool           uniform, trimmed;
    2097:   Petsc1DNodeFamily   nodeFamily;
    2098:   PetscInt            numNodeSkip;
    2099:   DMPlexInterpolatedFlag interpolated;
    2100:   PetscBool           isbdm;
    2101:   PetscErrorCode      ierr;

    2104:   /* step 1: sanitize input */
    2105:   PetscObjectGetComm((PetscObject) sp, &comm);
    2106:   DMGetDimension(dm, &dim);
    2107:   PetscObjectTypeCompare((PetscObject)sp, PETSCDUALSPACEBDM, &isbdm);
    2108:   if (isbdm) {
    2109:     sp->k = -(dim-1); /* form degree of H-div */
    2110:     PetscObjectChangeTypeName((PetscObject)sp, PETSCDUALSPACELAGRANGE);
    2111:   }
    2112:   PetscDualSpaceGetFormDegree(sp, &formDegree);
    2113:   if (PetscAbsInt(formDegree) > dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Form degree must be bounded by dimension");
    2114:   PetscDTBinomialInt(dim,PetscAbsInt(formDegree),&Nk);
    2115:   if (sp->Nc <= 0 && lag->numCopies > 0) sp->Nc = Nk * lag->numCopies;
    2116:   Nc = sp->Nc;
    2117:   if (Nc % Nk) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "Number of components is not a multiple of form degree size");
    2118:   if (lag->numCopies <= 0) lag->numCopies = Nc / Nk;
    2119:   Ncopies = lag->numCopies;
    2120:   if (Nc / Nk != Ncopies) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "Number of copies * (dim choose k) != Nc");
    2121:   if (!dim) sp->order = 0;
    2122:   order = sp->order;
    2123:   uniform = sp->uniform;
    2124:   if (!uniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Variable order not supported yet");
    2125:   if (lag->trimmed && !formDegree) lag->trimmed = PETSC_FALSE; /* trimmed spaces are the same as full spaces for 0-forms */
    2126:   if (lag->nodeType == PETSCDTNODES_DEFAULT) {
    2127:     lag->nodeType = PETSCDTNODES_GAUSSJACOBI;
    2128:     lag->nodeExponent = 0.;
    2129:     /* trimmed spaces don't include corner vertices, so don't use end nodes by default */
    2130:     lag->endNodes = lag->trimmed ? PETSC_FALSE : PETSC_TRUE;
    2131:   }
    2132:   /* If a trimmed space and the user did choose nodes with endpoints, skip them by default */
    2133:   if (lag->numNodeSkip < 0) lag->numNodeSkip = (lag->trimmed && lag->endNodes) ? 1 : 0;
    2134:   numNodeSkip = lag->numNodeSkip;
    2135:   if (lag->trimmed && !order) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zeroth order trimmed elements");
    2136:   if (lag->trimmed && PetscAbsInt(formDegree) == dim) { /* convert trimmed n-forms to untrimmed of one polynomial order less */
    2137:     sp->order--;
    2138:     order--;
    2139:     lag->trimmed = PETSC_FALSE;
    2140:   }
    2141:   trimmed = lag->trimmed;
    2142:   if (!order || PetscAbsInt(formDegree) == dim) lag->continuous = PETSC_FALSE;
    2143:   continuous = lag->continuous;
    2144:   DMPlexGetDepth(dm, &depth);
    2145:   DMPlexGetChart(dm, &pStart, &pEnd);
    2146:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
    2147:   if (pStart != 0 || cStart != 0) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Expect DM with chart starting at zero and cells first");
    2148:   if (cEnd != 1) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Use PETSCDUALSPACEREFINED for multi-cell reference meshes");
    2149:   DMPlexIsInterpolated(dm, &interpolated);
    2150:   if (interpolated != DMPLEX_INTERPOLATED_FULL) {
    2151:     DMPlexInterpolate(dm, &dmint);
    2152:   } else {
    2153:     PetscObjectReference((PetscObject)dm);
    2154:     dmint = dm;
    2155:   }
    2156:   tensorCell = PETSC_FALSE;
    2157:   if (dim > 1) {
    2158:     DMPlexPointIsTensor(dmint, 0, &tensorCell, &tensorf, &tensorf2);
    2159:   }
    2160:   lag->tensorCell = tensorCell;
    2161:   if (dim < 2 || !lag->tensorCell) lag->tensorSpace = PETSC_FALSE;
    2162:   tensorSpace = lag->tensorSpace;
    2163:   if (!lag->nodeFamily) {
    2164:     Petsc1DNodeFamilyCreate(lag->nodeType, lag->nodeExponent, lag->endNodes, &lag->nodeFamily);
    2165:   }
    2166:   nodeFamily = lag->nodeFamily;
    2167:   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");

    2169:   /* step 2: construct the boundary spaces */
    2170:   PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);
    2171:   PetscCalloc1(pEnd,&(sp->pointSpaces));
    2172:   for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);}
    2173:   PetscDualSpaceSectionCreate_Internal(sp, &section);
    2174:   sp->pointSection = section;
    2175:   if (continuous && !(lag->interiorOnly)) {
    2176:     PetscInt h;

    2178:     for (p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */
    2179:       PetscReal v0[3];
    2180:       DMPolytopeType ptype;
    2181:       PetscReal J[9], detJ;
    2182:       PetscInt  q;

    2184:       DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, NULL, &detJ);
    2185:       DMPlexGetCellType(dm, p, &ptype);

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

    2191:         DMPlexGetCellType(dm, q, &qtype);
    2192:         if (qtype == ptype) break;
    2193:       }
    2194:       if (q < p) { /* this facet has the same dual space as that one */
    2195:         PetscObjectReference((PetscObject)sp->pointSpaces[q]);
    2196:         sp->pointSpaces[p] = sp->pointSpaces[q];
    2197:         continue;
    2198:       }
    2199:       /* if not, recursively compute this dual space */
    2200:       PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,p,formDegree,Ncopies,PETSC_FALSE,&sp->pointSpaces[p]);
    2201:     }
    2202:     for (h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */
    2203:       PetscInt hd = depth - h;
    2204:       PetscInt hdim = dim - h;

    2206:       if (hdim < PetscAbsInt(formDegree)) break;
    2207:       for (p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
    2208:         PetscInt suppSize, s;
    2209:         const PetscInt *supp;

    2211:         DMPlexGetSupportSize(dm, p, &suppSize);
    2212:         DMPlexGetSupport(dm, p, &supp);
    2213:         for (s = 0; s < suppSize; s++) {
    2214:           DM             qdm;
    2215:           PetscDualSpace qsp, psp;
    2216:           PetscInt c, coneSize, q;
    2217:           const PetscInt *cone;
    2218:           const PetscInt *refCone;

    2220:           q = supp[0];
    2221:           qsp = sp->pointSpaces[q];
    2222:           DMPlexGetConeSize(dm, q, &coneSize);
    2223:           DMPlexGetCone(dm, q, &cone);
    2224:           for (c = 0; c < coneSize; c++) if (cone[c] == p) break;
    2225:           if (c == coneSize) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "cone/support mismatch");
    2226:           PetscDualSpaceGetDM(qsp, &qdm);
    2227:           DMPlexGetCone(qdm, 0, &refCone);
    2228:           /* get the equivalent dual space from the support dual space */
    2229:           PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp);
    2230:           if (!s) {
    2231:             PetscObjectReference((PetscObject)psp);
    2232:             sp->pointSpaces[p] = psp;
    2233:           }
    2234:         }
    2235:       }
    2236:     }
    2237:     for (p = 1; p < pEnd; p++) {
    2238:       PetscInt pspdim;
    2239:       if (!sp->pointSpaces[p]) continue;
    2240:       PetscDualSpaceGetInteriorDimension(sp->pointSpaces[p], &pspdim);
    2241:       PetscSectionSetDof(section, p, pspdim);
    2242:     }
    2243:   }

    2245:   if (Ncopies > 1) {
    2246:     Mat intMatScalar, allMatScalar;
    2247:     PetscDualSpace scalarsp;
    2248:     PetscDualSpace_Lag *scalarlag;

    2250:     PetscDualSpaceDuplicate(sp, &scalarsp);
    2251:     /* Setting the number of components to Nk is a space with 1 copy of each k-form */
    2252:     PetscDualSpaceSetNumComponents(scalarsp, Nk);
    2253:     PetscDualSpaceSetUp(scalarsp);
    2254:     PetscDualSpaceGetInteriorData(scalarsp, &(sp->intNodes), &intMatScalar);
    2255:     PetscObjectReference((PetscObject)(sp->intNodes));
    2256:     if (intMatScalar) {PetscDualSpaceLagrangeMatrixCreateCopies(intMatScalar, Nk, Ncopies, &(sp->intMat));}
    2257:     PetscDualSpaceGetAllData(scalarsp, &(sp->allNodes), &allMatScalar);
    2258:     PetscObjectReference((PetscObject)(sp->allNodes));
    2259:     PetscDualSpaceLagrangeMatrixCreateCopies(allMatScalar, Nk, Ncopies, &(sp->allMat));
    2260:     sp->spdim = scalarsp->spdim * Ncopies;
    2261:     sp->spintdim = scalarsp->spintdim * Ncopies;
    2262:     scalarlag = (PetscDualSpace_Lag *) scalarsp->data;
    2263:     PetscLagNodeIndicesReference(scalarlag->vertIndices);
    2264:     lag->vertIndices = scalarlag->vertIndices;
    2265:     PetscLagNodeIndicesReference(scalarlag->intNodeIndices);
    2266:     lag->intNodeIndices = scalarlag->intNodeIndices;
    2267:     PetscLagNodeIndicesReference(scalarlag->allNodeIndices);
    2268:     lag->allNodeIndices = scalarlag->allNodeIndices;
    2269:     PetscDualSpaceDestroy(&scalarsp);
    2270:     PetscSectionSetDof(section, 0, sp->spintdim);
    2271:     PetscDualSpaceSectionSetUp_Internal(sp, section);
    2272:     PetscDualSpaceComputeFunctionalsFromAllData(sp);
    2273:     PetscFree2(pStratStart, pStratEnd);
    2274:     DMDestroy(&dmint);
    2275:     return(0);
    2276:   }

    2278:   if (trimmed && !continuous) {
    2279:     /* the dofs of a trimmed space don't have a nice tensor/lattice structure:
    2280:      * just construct the continuous dual space and copy all of the data over,
    2281:      * allocating it all to the cell instead of splitting it up between the boundaries */
    2282:     PetscDualSpace  spcont;
    2283:     PetscInt        spdim, f;
    2284:     PetscQuadrature allNodes;
    2285:     PetscDualSpace_Lag *lagc;
    2286:     Mat             allMat;

    2288:     PetscDualSpaceDuplicate(sp, &spcont);
    2289:     PetscDualSpaceLagrangeSetContinuity(spcont, PETSC_TRUE);
    2290:     PetscDualSpaceSetUp(spcont);
    2291:     PetscDualSpaceGetDimension(spcont, &spdim);
    2292:     sp->spdim = sp->spintdim = spdim;
    2293:     PetscSectionSetDof(section, 0, spdim);
    2294:     PetscDualSpaceSectionSetUp_Internal(sp, section);
    2295:     PetscMalloc1(spdim, &(sp->functional));
    2296:     for (f = 0; f < spdim; f++) {
    2297:       PetscQuadrature fn;

    2299:       PetscDualSpaceGetFunctional(spcont, f, &fn);
    2300:       PetscObjectReference((PetscObject)fn);
    2301:       sp->functional[f] = fn;
    2302:     }
    2303:     PetscDualSpaceGetAllData(spcont, &allNodes, &allMat);
    2304:     PetscObjectReference((PetscObject) allNodes);
    2305:     PetscObjectReference((PetscObject) allNodes);
    2306:     sp->allNodes = sp->intNodes = allNodes;
    2307:     PetscObjectReference((PetscObject) allMat);
    2308:     PetscObjectReference((PetscObject) allMat);
    2309:     sp->allMat = sp->intMat = allMat;
    2310:     lagc = (PetscDualSpace_Lag *) spcont->data;
    2311:     PetscLagNodeIndicesReference(lagc->vertIndices);
    2312:     lag->vertIndices = lagc->vertIndices;
    2313:     PetscLagNodeIndicesReference(lagc->allNodeIndices);
    2314:     PetscLagNodeIndicesReference(lagc->allNodeIndices);
    2315:     lag->intNodeIndices = lagc->allNodeIndices;
    2316:     lag->allNodeIndices = lagc->allNodeIndices;
    2317:     PetscDualSpaceDestroy(&spcont);
    2318:     PetscFree2(pStratStart, pStratEnd);
    2319:     DMDestroy(&dmint);
    2320:     return(0);
    2321:   }

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

    2327:     if (trimmed) {
    2328:       /* there is one dof in the interior of the a trimmed element for each full polynomial of with degree at most
    2329:        * order + k - dim - 1 */
    2330:       if (order + PetscAbsInt(formDegree) > dim) {
    2331:         PetscInt sum = order + PetscAbsInt(formDegree) - dim - 1;
    2332:         PetscInt nDofs;

    2334:         PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));
    2335:         MatGetSize(sp->intMat, &nDofs, NULL);
    2336:         PetscSectionSetDof(section, 0, nDofs);
    2337:       }
    2338:       PetscDualSpaceSectionSetUp_Internal(sp, section);
    2339:       PetscDualSpaceCreateAllDataFromInteriorData(sp);
    2340:       PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
    2341:     } else {
    2342:       if (!continuous) {
    2343:         /* if discontinuous just construct one node for each set of dofs (a set of dofs is a basis for the k-form
    2344:          * space) */
    2345:         PetscInt sum = order;
    2346:         PetscInt nDofs;

    2348:         PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));
    2349:         MatGetSize(sp->intMat, &nDofs, NULL);
    2350:         PetscSectionSetDof(section, 0, nDofs);
    2351:         PetscDualSpaceSectionSetUp_Internal(sp, section);
    2352:         PetscObjectReference((PetscObject)(sp->intNodes));
    2353:         sp->allNodes = sp->intNodes;
    2354:         PetscObjectReference((PetscObject)(sp->intMat));
    2355:         sp->allMat = sp->intMat;
    2356:         PetscLagNodeIndicesReference(lag->intNodeIndices);
    2357:         lag->allNodeIndices = lag->intNodeIndices;
    2358:       } else {
    2359:         /* there is one dof in the interior of the a full element for each trimmed polynomial of with degree at most
    2360:          * order + k - dim, but with complementary form degree */
    2361:         if (order + PetscAbsInt(formDegree) > dim) {
    2362:           PetscDualSpace trimmedsp;
    2363:           PetscDualSpace_Lag *trimmedlag;
    2364:           PetscQuadrature intNodes;
    2365:           PetscInt trFormDegree = formDegree >= 0 ? formDegree - dim : dim - PetscAbsInt(formDegree);
    2366:           PetscInt nDofs;
    2367:           Mat intMat;

    2369:           PetscDualSpaceDuplicate(sp, &trimmedsp);
    2370:           PetscDualSpaceLagrangeSetTrimmed(trimmedsp, PETSC_TRUE);
    2371:           PetscDualSpaceSetOrder(trimmedsp, order + PetscAbsInt(formDegree) - dim);
    2372:           PetscDualSpaceSetFormDegree(trimmedsp, trFormDegree);
    2373:           trimmedlag = (PetscDualSpace_Lag *) trimmedsp->data;
    2374:           trimmedlag->numNodeSkip = numNodeSkip + 1;
    2375:           PetscDualSpaceSetUp(trimmedsp);
    2376:           PetscDualSpaceGetAllData(trimmedsp, &intNodes, &intMat);
    2377:           PetscObjectReference((PetscObject)intNodes);
    2378:           sp->intNodes = intNodes;
    2379:           PetscLagNodeIndicesReference(trimmedlag->allNodeIndices);
    2380:           lag->intNodeIndices = trimmedlag->allNodeIndices;
    2381:           PetscObjectReference((PetscObject)intMat);
    2382:           if (PetscAbsInt(formDegree) > 0 && PetscAbsInt(formDegree) < dim) {
    2383:             PetscReal *T;
    2384:             PetscScalar *work;
    2385:             PetscInt nCols, nRows;
    2386:             Mat intMatT;

    2388:             MatDuplicate(intMat, MAT_COPY_VALUES, &intMatT);
    2389:             MatGetSize(intMat, &nRows, &nCols);
    2390:             PetscMalloc2(Nk * Nk, &T, nCols, &work);
    2391:             BiunitSimplexSymmetricFormTransformation(dim, formDegree, T);
    2392:             for (PetscInt row = 0; row < nRows; row++) {
    2393:               PetscInt nrCols;
    2394:               const PetscInt *rCols;
    2395:               const PetscScalar *rVals;

    2397:               MatGetRow(intMat, row, &nrCols, &rCols, &rVals);
    2398:               if (nrCols % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in intMat matrix are not in k-form size blocks");
    2399:               for (PetscInt b = 0; b < nrCols; b += Nk) {
    2400:                 const PetscScalar *v = &rVals[b];
    2401:                 PetscScalar *w = &work[b];
    2402:                 for (PetscInt j = 0; j < Nk; j++) {
    2403:                   w[j] = 0.;
    2404:                   for (PetscInt i = 0; i < Nk; i++) {
    2405:                     w[j] += v[i] * T[i * Nk + j];
    2406:                   }
    2407:                 }
    2408:               }
    2409:               MatSetValuesBlocked(intMatT, 1, &row, nrCols, rCols, work, INSERT_VALUES);
    2410:               MatRestoreRow(intMat, row, &nrCols, &rCols, &rVals);
    2411:             }
    2412:             MatAssemblyBegin(intMatT, MAT_FINAL_ASSEMBLY);
    2413:             MatAssemblyEnd(intMatT, MAT_FINAL_ASSEMBLY);
    2414:             MatDestroy(&intMat);
    2415:             intMat = intMatT;
    2416:             PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));
    2417:             PetscLagNodeIndicesDuplicate(trimmedlag->allNodeIndices, &(lag->intNodeIndices));
    2418:             {
    2419:               PetscInt nNodes = lag->intNodeIndices->nNodes;
    2420:               PetscReal *newNodeVec = lag->intNodeIndices->nodeVec;
    2421:               const PetscReal *oldNodeVec = trimmedlag->allNodeIndices->nodeVec;

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

    2427:                 for (PetscInt j = 0; j < Nk; j++) {
    2428:                   w[j] = 0.;
    2429:                   for (PetscInt i = 0; i < Nk; i++) {
    2430:                     w[j] += v[i] * T[i * Nk + j];
    2431:                   }
    2432:                 }
    2433:               }
    2434:             }
    2435:             PetscFree2(T, work);
    2436:           }
    2437:           sp->intMat = intMat;
    2438:           MatGetSize(sp->intMat, &nDofs, NULL);
    2439:           PetscDualSpaceDestroy(&trimmedsp);
    2440:           PetscSectionSetDof(section, 0, nDofs);
    2441:         }
    2442:         PetscDualSpaceSectionSetUp_Internal(sp, section);
    2443:         PetscDualSpaceCreateAllDataFromInteriorData(sp);
    2444:         PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
    2445:       }
    2446:     }
    2447:   } else {
    2448:     PetscQuadrature intNodesTrace = NULL;
    2449:     PetscQuadrature intNodesFiber = NULL;
    2450:     PetscQuadrature intNodes = NULL;
    2451:     PetscLagNodeIndices intNodeIndices = NULL;
    2452:     Mat             intMat = NULL;

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

    2460:       if (sp->pointSpaces[tensorf]) {
    2461:         PetscObjectReference((PetscObject)(sp->pointSpaces[tensorf]));
    2462:         trace = sp->pointSpaces[tensorf];
    2463:       } else {
    2464:         PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,formDegree,Ncopies,PETSC_TRUE,&trace);
    2465:       }
    2466:       PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,0,1,PETSC_TRUE,&fiber);
    2467:       tracel = (PetscDualSpace_Lag *) trace->data;
    2468:       fiberl = (PetscDualSpace_Lag *) fiber->data;
    2469:       PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));
    2470:       PetscDualSpaceGetInteriorData(trace, &intNodesTrace, &intMatTrace);
    2471:       PetscDualSpaceGetInteriorData(fiber, &intNodesFiber, &intMatFiber);
    2472:       if (intNodesTrace && intNodesFiber) {
    2473:         PetscQuadratureCreateTensor(intNodesTrace, intNodesFiber, &intNodes);
    2474:         MatTensorAltV(intMatTrace, intMatFiber, dim-1, formDegree, 1, 0, &intMat);
    2475:         PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, formDegree, fiberl->intNodeIndices, 1, 0, &intNodeIndices);
    2476:       }
    2477:       PetscObjectReference((PetscObject) intNodesTrace);
    2478:       PetscObjectReference((PetscObject) intNodesFiber);
    2479:       PetscDualSpaceDestroy(&fiber);
    2480:       PetscDualSpaceDestroy(&trace);
    2481:     }
    2482:     if (PetscAbsInt(formDegree) > 0) { /* get the trace (k-1)-forms on the first facet, and the 1-forms on the edge,
    2483:                                           and wedge them together to create the remaining k-form dofs */
    2484:       PetscDualSpace  trace, fiber;
    2485:       PetscDualSpace_Lag *tracel, *fiberl;
    2486:       PetscQuadrature intNodesTrace2, intNodesFiber2, intNodes2;
    2487:       PetscLagNodeIndices intNodeIndices2;
    2488:       Mat             intMatTrace, intMatFiber, intMat2;
    2489:       PetscInt        traceDegree = formDegree > 0 ? formDegree - 1 : formDegree + 1;
    2490:       PetscInt        fiberDegree = formDegree > 0 ? 1 : -1;

    2492:       PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,traceDegree,Ncopies,PETSC_TRUE,&trace);
    2493:       PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,fiberDegree,1,PETSC_TRUE,&fiber);
    2494:       tracel = (PetscDualSpace_Lag *) trace->data;
    2495:       fiberl = (PetscDualSpace_Lag *) fiber->data;
    2496:       if (!lag->vertIndices) {
    2497:         PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));
    2498:       }
    2499:       PetscDualSpaceGetInteriorData(trace, &intNodesTrace2, &intMatTrace);
    2500:       PetscDualSpaceGetInteriorData(fiber, &intNodesFiber2, &intMatFiber);
    2501:       if (intNodesTrace2 && intNodesFiber2) {
    2502:         PetscQuadratureCreateTensor(intNodesTrace2, intNodesFiber2, &intNodes2);
    2503:         MatTensorAltV(intMatTrace, intMatFiber, dim-1, traceDegree, 1, fiberDegree, &intMat2);
    2504:         PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, traceDegree, fiberl->intNodeIndices, 1, fiberDegree, &intNodeIndices2);
    2505:         if (!intMat) {
    2506:           intMat = intMat2;
    2507:           intNodes = intNodes2;
    2508:           intNodeIndices = intNodeIndices2;
    2509:         } else {
    2510:           /* merge the matrices, quadrature points, and nodes */
    2511:           PetscInt         nM;
    2512:           PetscInt         nDof, nDof2;
    2513:           PetscInt        *toMerged = NULL, *toMerged2 = NULL;
    2514:           PetscQuadrature  merged = NULL;
    2515:           PetscLagNodeIndices intNodeIndicesMerged = NULL;
    2516:           Mat              matMerged = NULL;

    2518:           MatGetSize(intMat, &nDof, NULL);
    2519:           MatGetSize(intMat2, &nDof2, NULL);
    2520:           PetscQuadraturePointsMerge(intNodes, intNodes2, &merged, &toMerged, &toMerged2);
    2521:           PetscQuadratureGetData(merged, NULL, NULL, &nM, NULL, NULL);
    2522:           MatricesMerge(intMat, intMat2, dim, formDegree, nM, toMerged, toMerged2, &matMerged);
    2523:           PetscLagNodeIndicesMerge(intNodeIndices, intNodeIndices2, &intNodeIndicesMerged);
    2524:           PetscFree(toMerged);
    2525:           PetscFree(toMerged2);
    2526:           MatDestroy(&intMat);
    2527:           MatDestroy(&intMat2);
    2528:           PetscQuadratureDestroy(&intNodes);
    2529:           PetscQuadratureDestroy(&intNodes2);
    2530:           PetscLagNodeIndicesDestroy(&intNodeIndices);
    2531:           PetscLagNodeIndicesDestroy(&intNodeIndices2);
    2532:           intNodes = merged;
    2533:           intMat = matMerged;
    2534:           intNodeIndices = intNodeIndicesMerged;
    2535:           if (!trimmed) {
    2536:             /* I think users expect that, when a node has a full basis for the k-forms,
    2537:              * they should be consecutive dofs.  That isn't the case for trimmed spaces,
    2538:              * but is for some of the nodes in untrimmed spaces, so in that case we
    2539:              * sort them to group them by node */
    2540:             Mat intMatPerm;

    2542:             MatPermuteByNodeIdx(intMat, intNodeIndices, &intMatPerm);
    2543:             MatDestroy(&intMat);
    2544:             intMat = intMatPerm;
    2545:           }
    2546:         }
    2547:       }
    2548:       PetscDualSpaceDestroy(&fiber);
    2549:       PetscDualSpaceDestroy(&trace);
    2550:     }
    2551:     PetscQuadratureDestroy(&intNodesTrace);
    2552:     PetscQuadratureDestroy(&intNodesFiber);
    2553:     sp->intNodes = intNodes;
    2554:     sp->intMat = intMat;
    2555:     lag->intNodeIndices = intNodeIndices;
    2556:     {
    2557:       PetscInt nDofs = 0;

    2559:       if (intMat) {
    2560:         MatGetSize(intMat, &nDofs, NULL);
    2561:       }
    2562:       PetscSectionSetDof(section, 0, nDofs);
    2563:     }
    2564:     PetscDualSpaceSectionSetUp_Internal(sp, section);
    2565:     if (continuous) {
    2566:       PetscDualSpaceCreateAllDataFromInteriorData(sp);
    2567:       PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
    2568:     } else {
    2569:       PetscObjectReference((PetscObject) intNodes);
    2570:       sp->allNodes = intNodes;
    2571:       PetscObjectReference((PetscObject) intMat);
    2572:       sp->allMat = intMat;
    2573:       PetscLagNodeIndicesReference(intNodeIndices);
    2574:       lag->allNodeIndices = intNodeIndices;
    2575:     }
    2576:   }
    2577:   PetscSectionGetStorageSize(section, &sp->spdim);
    2578:   PetscSectionGetConstrainedStorageSize(section, &sp->spintdim);
    2579:   PetscDualSpaceComputeFunctionalsFromAllData(sp);
    2580:   PetscFree2(pStratStart, pStratEnd);
    2581:   DMDestroy(&dmint);
    2582:   return(0);
    2583: }

    2585: /* Create a matrix that represents the transformation that DMPlexVecGetClosure() would need
    2586:  * to get the representation of the dofs for a mesh point if the mesh point had this orientation
    2587:  * relative to the cell */
    2588: PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat)
    2589: {
    2590:   PetscDualSpace_Lag *lag;
    2591:   DM dm;
    2592:   PetscLagNodeIndices vertIndices, intNodeIndices;
    2593:   PetscLagNodeIndices ni;
    2594:   PetscInt nodeIdxDim, nodeVecDim, nNodes;
    2595:   PetscInt formDegree;
    2596:   PetscInt *perm, *permOrnt;
    2597:   PetscInt *nnz;
    2598:   PetscInt n;
    2599:   PetscInt maxGroupSize;
    2600:   PetscScalar *V, *W, *work;
    2601:   Mat A;

    2605:   if (!sp->spintdim) {
    2606:     *symMat = NULL;
    2607:     return(0);
    2608:   }
    2609:   lag = (PetscDualSpace_Lag *) sp->data;
    2610:   vertIndices = lag->vertIndices;
    2611:   intNodeIndices = lag->intNodeIndices;
    2612:   PetscDualSpaceGetDM(sp, &dm);
    2613:   PetscDualSpaceGetFormDegree(sp, &formDegree);
    2614:   PetscNew(&ni);
    2615:   ni->refct = 1;
    2616:   ni->nodeIdxDim = nodeIdxDim = intNodeIndices->nodeIdxDim;
    2617:   ni->nodeVecDim = nodeVecDim = intNodeIndices->nodeVecDim;
    2618:   ni->nNodes = nNodes = intNodeIndices->nNodes;
    2619:   PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
    2620:   PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
    2621:   /* push forward the dofs by the symmetry of the reference element induced by ornt */
    2622:   PetscLagNodeIndicesPushForward(dm, vertIndices, 0, vertIndices, intNodeIndices, ornt, formDegree, ni->nodeIdx, ni->nodeVec);
    2623:   /* get the revlex order for both the original and transformed dofs */
    2624:   PetscLagNodeIndicesGetPermutation(intNodeIndices, &perm);
    2625:   PetscLagNodeIndicesGetPermutation(ni, &permOrnt);
    2626:   PetscMalloc1(nNodes, &nnz);
    2627:   for (n = 0, maxGroupSize = 0; n < nNodes;) { /* incremented in the loop */
    2628:     PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
    2629:     PetscInt m, nEnd;
    2630:     PetscInt groupSize;
    2631:     /* for each group of dofs that have the same nodeIdx coordinate */
    2632:     for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
    2633:       PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
    2634:       PetscInt d;

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

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

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

    2652:         /* compare the oriented permutation indices */
    2653:         for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break;
    2654:         if (d < nodeIdxDim) break;
    2655:       }
    2656:       if (m < nEnd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs with same index after symmetry not same block size");
    2657:     }
    2658:     groupSize = nEnd - n;
    2659:     /* each pushforward dof vector will be expressed in a basis of the unpermuted dofs */
    2660:     for (m = n; m < nEnd; m++) nnz[permOrnt[m]] = groupSize;

    2662:     maxGroupSize = PetscMax(maxGroupSize, nEnd - n);
    2663:     n = nEnd;
    2664:   }
    2665:   if (maxGroupSize > nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs are not in blocks that can be solved");
    2666:   MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes, nNodes, 0, nnz, &A);
    2667:   PetscFree(nnz);
    2668:   PetscMalloc3(maxGroupSize * nodeVecDim, &V, maxGroupSize * nodeVecDim, &W, nodeVecDim * 2, &work);
    2669:   for (n = 0; n < nNodes;) { /* incremented in the loop */
    2670:     PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
    2671:     PetscInt nEnd;
    2672:     PetscInt m;
    2673:     PetscInt groupSize;
    2674:     for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
    2675:       PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
    2676:       PetscInt d;

    2678:       /* compare the oriented permutation indices */
    2679:       for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break;
    2680:       if (d < nodeIdxDim) break;
    2681:     }
    2682:     groupSize = nEnd - n;
    2683:     /* get all of the vectors from the original and all of the pushforward vectors */
    2684:     for (m = n; m < nEnd; m++) {
    2685:       PetscInt d;

    2687:       for (d = 0; d < nodeVecDim; d++) {
    2688:         V[(m - n) * nodeVecDim + d] = intNodeIndices->nodeVec[perm[m] * nodeVecDim + d];
    2689:         W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
    2690:       }
    2691:     }
    2692:     /* now we have to solve for W in terms of V: the systems isn't always square, but the span
    2693:      * of V and W should always be the same, so the solution of the normal equations works */
    2694:     {
    2695:       char transpose = 'N';
    2696:       PetscBLASInt bm = nodeVecDim;
    2697:       PetscBLASInt bn = groupSize;
    2698:       PetscBLASInt bnrhs = groupSize;
    2699:       PetscBLASInt blda = bm;
    2700:       PetscBLASInt bldb = bm;
    2701:       PetscBLASInt blwork = 2 * nodeVecDim;
    2702:       PetscBLASInt info;

    2704:       PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&bm,&bn,&bnrhs,V,&blda,W,&bldb,work,&blwork, &info));
    2705:       if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS");
    2706:       /* repack */
    2707:       {
    2708:         PetscInt i, j;

    2710:         for (i = 0; i < groupSize; i++) {
    2711:           for (j = 0; j < groupSize; j++) {
    2712:             /* notice the different leading dimension */
    2713:             V[i * groupSize + j] = W[i * nodeVecDim + j];
    2714:           }
    2715:         }
    2716:       }
    2717:       if (PetscDefined(USE_DEBUG)) {
    2718:         PetscReal res;

    2720:         /* check that the normal error is 0 */
    2721:         for (m = n; m < nEnd; m++) {
    2722:           PetscInt d;

    2724:           for (d = 0; d < nodeVecDim; d++) {
    2725:             W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
    2726:           }
    2727:         }
    2728:         res = 0.;
    2729:         for (PetscInt i = 0; i < groupSize; i++) {
    2730:           for (PetscInt j = 0; j < nodeVecDim; j++) {
    2731:             for (PetscInt k = 0; k < groupSize; k++) {
    2732:               W[i * nodeVecDim + j] -= V[i * groupSize + k] * intNodeIndices->nodeVec[perm[n+k] * nodeVecDim + j];
    2733:             }
    2734:             res += PetscAbsScalar(W[i * nodeVecDim + j]);
    2735:           }
    2736:         }
    2737:         if (res > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Dof block did not solve");
    2738:       }
    2739:     }
    2740:     MatSetValues(A, groupSize, &permOrnt[n], groupSize, &perm[n], V, INSERT_VALUES);
    2741:     n = nEnd;
    2742:   }
    2743:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    2744:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
    2745:   *symMat = A;
    2746:   PetscFree3(V,W,work);
    2747:   PetscLagNodeIndicesDestroy(&ni);
    2748:   return(0);
    2749: }

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

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

    2755: /* the existing interface for symmetries is insufficient for all cases:
    2756:  * - it should be sufficient for form degrees that are scalar (0 and n)
    2757:  * - it should be sufficient for hypercube dofs
    2758:  * - it isn't sufficient for simplex cells with non-scalar form degrees if
    2759:  *   there are any dofs in the interior
    2760:  *
    2761:  * We compute the general transformation matrices, and if they fit, we return them,
    2762:  * otherwise we error (but we should probably change the interface to allow for
    2763:  * these symmetries)
    2764:  */
    2765: static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
    2766: {
    2767:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
    2768:   PetscInt           dim, order, Nc;
    2769:   PetscErrorCode     ierr;

    2772:   PetscDualSpaceGetOrder(sp,&order);
    2773:   PetscDualSpaceGetNumComponents(sp,&Nc);
    2774:   DMGetDimension(sp->dm,&dim);
    2775:   if (!lag->symComputed) { /* store symmetries */
    2776:     PetscInt       pStart, pEnd, p;
    2777:     PetscInt       numPoints;
    2778:     PetscInt       numFaces;
    2779:     PetscInt       spintdim;
    2780:     PetscInt       ***symperms;
    2781:     PetscScalar    ***symflips;

    2783:     DMPlexGetChart(sp->dm, &pStart, &pEnd);
    2784:     numPoints = pEnd - pStart;
    2785:     DMPlexGetConeSize(sp->dm, 0, &numFaces);
    2786:     PetscCalloc1(numPoints,&symperms);
    2787:     PetscCalloc1(numPoints,&symflips);
    2788:     spintdim = sp->spintdim;
    2789:     /* The nodal symmetry behavior is not present when tensorSpace != tensorCell: someone might want this for the "S"
    2790:      * family of FEEC spaces.  Most used in particular are discontinuous polynomial L2 spaces in tensor cells, where
    2791:      * the symmetries are not necessary for FE assembly.  So for now we assume this is the case and don't return
    2792:      * symmetries if tensorSpace != tensorCell */
    2793:     if (spintdim && 0 < dim && dim < 3 && (lag->tensorSpace == lag->tensorCell)) { /* compute self symmetries */
    2794:       PetscInt **cellSymperms;
    2795:       PetscScalar **cellSymflips;
    2796:       PetscInt ornt;
    2797:       PetscInt nCopies = Nc / lag->intNodeIndices->nodeVecDim;
    2798:       PetscInt nNodes = lag->intNodeIndices->nNodes;

    2800:       lag->numSelfSym = 2 * numFaces;
    2801:       lag->selfSymOff = numFaces;
    2802:       PetscCalloc1(2*numFaces,&cellSymperms);
    2803:       PetscCalloc1(2*numFaces,&cellSymflips);
    2804:       /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
    2805:       symperms[0] = &cellSymperms[numFaces];
    2806:       symflips[0] = &cellSymflips[numFaces];
    2807:       if (lag->intNodeIndices->nodeVecDim * nCopies != Nc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs");
    2808:       if (nNodes * nCopies != spintdim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs");
    2809:       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 */
    2810:         Mat symMat;
    2811:         PetscInt *perm;
    2812:         PetscScalar *flips;
    2813:         PetscInt i;

    2815:         if (!ornt) continue;
    2816:         PetscMalloc1(spintdim, &perm);
    2817:         PetscCalloc1(spintdim, &flips);
    2818:         for (i = 0; i < spintdim; i++) perm[i] = -1;
    2819:         PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, ornt, &symMat);
    2820:         for (i = 0; i < nNodes; i++) {
    2821:           PetscInt ncols;
    2822:           PetscInt j, k;
    2823:           const PetscInt *cols;
    2824:           const PetscScalar *vals;
    2825:           PetscBool nz_seen = PETSC_FALSE;

    2827:           MatGetRow(symMat, i, &ncols, &cols, &vals);
    2828:           for (j = 0; j < ncols; j++) {
    2829:             if (PetscAbsScalar(vals[j]) > PETSC_SMALL) {
    2830:               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");
    2831:               nz_seen = PETSC_TRUE;
    2832:               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");
    2833:               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");
    2834:               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");
    2835:               for (k = 0; k < nCopies; k++) {
    2836:                 perm[cols[j] * nCopies + k] = i * nCopies + k;
    2837:               }
    2838:               if (PetscRealPart(vals[j]) < 0.) {
    2839:                 for (k = 0; k < nCopies; k++) {
    2840:                   flips[i * nCopies + k] = -1.;
    2841:                 }
    2842:               } else {
    2843:                 for (k = 0; k < nCopies; k++) {
    2844:                   flips[i * nCopies + k] = 1.;
    2845:                 }
    2846:               }
    2847:             }
    2848:           }
    2849:           MatRestoreRow(symMat, i, &ncols, &cols, &vals);
    2850:         }
    2851:         MatDestroy(&symMat);
    2852:         /* if there were no sign flips, keep NULL */
    2853:         for (i = 0; i < spintdim; i++) if (flips[i] != 1.) break;
    2854:         if (i == spintdim) {
    2855:           PetscFree(flips);
    2856:           flips = NULL;
    2857:         }
    2858:         /* if the permutation is identity, keep NULL */
    2859:         for (i = 0; i < spintdim; i++) if (perm[i] != i) break;
    2860:         if (i == spintdim) {
    2861:           PetscFree(perm);
    2862:           perm = NULL;
    2863:         }
    2864:         symperms[0][ornt] = perm;
    2865:         symflips[0][ornt] = flips;
    2866:       }
    2867:       /* if no orientations produced non-identity permutations, keep NULL */
    2868:       for (ornt = -numFaces; ornt < numFaces; ornt++) if (symperms[0][ornt]) break;
    2869:       if (ornt == numFaces) {
    2870:         PetscFree(cellSymperms);
    2871:         symperms[0] = NULL;
    2872:       }
    2873:       /* if no orientations produced sign flips, keep NULL */
    2874:       for (ornt = -numFaces; ornt < numFaces; ornt++) if (symflips[0][ornt]) break;
    2875:       if (ornt == numFaces) {
    2876:         PetscFree(cellSymflips);
    2877:         symflips[0] = NULL;
    2878:       }
    2879:     }
    2880:     { /* get the symmetries of closure points */
    2881:       PetscInt closureSize = 0;
    2882:       PetscInt *closure = NULL;
    2883:       PetscInt r;

    2885:       DMPlexGetTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);
    2886:       for (r = 0; r < closureSize; r++) {
    2887:         PetscDualSpace psp;
    2888:         PetscInt point = closure[2 * r];
    2889:         PetscInt pspintdim;
    2890:         const PetscInt ***psymperms = NULL;
    2891:         const PetscScalar ***psymflips = NULL;

    2893:         if (!point) continue;
    2894:         PetscDualSpaceGetPointSubspace(sp, point, &psp);
    2895:         if (!psp) continue;
    2896:         PetscDualSpaceGetInteriorDimension(psp, &pspintdim);
    2897:         if (!pspintdim) continue;
    2898:         PetscDualSpaceGetSymmetries(psp,&psymperms,&psymflips);
    2899:         symperms[r] = (PetscInt **) (psymperms ? psymperms[0] : NULL);
    2900:         symflips[r] = (PetscScalar **) (psymflips ? psymflips[0] : NULL);
    2901:       }
    2902:       DMPlexRestoreTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);
    2903:     }
    2904:     for (p = 0; p < pEnd; p++) if (symperms[p]) break;
    2905:     if (p == pEnd) {
    2906:       PetscFree(symperms);
    2907:       symperms = NULL;
    2908:     }
    2909:     for (p = 0; p < pEnd; p++) if (symflips[p]) break;
    2910:     if (p == pEnd) {
    2911:       PetscFree(symflips);
    2912:       symflips = NULL;
    2913:     }
    2914:     lag->symperms = symperms;
    2915:     lag->symflips = symflips;
    2916:     lag->symComputed = PETSC_TRUE;
    2917:   }
    2918:   if (perms) *perms = (const PetscInt ***) lag->symperms;
    2919:   if (flips) *flips = (const PetscScalar ***) lag->symflips;
    2920:   return(0);
    2921: }

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

    2930:   *continuous = lag->continuous;
    2931:   return(0);
    2932: }

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

    2940:   lag->continuous = continuous;
    2941:   return(0);
    2942: }

    2944: /*@
    2945:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

    2947:   Not Collective

    2949:   Input Parameter:
    2950: . sp         - the PetscDualSpace

    2952:   Output Parameter:
    2953: . continuous - flag for element continuity

    2955:   Level: intermediate

    2957: .seealso: PetscDualSpaceLagrangeSetContinuity()
    2958: @*/
    2959: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
    2960: {

    2966:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
    2967:   return(0);
    2968: }

    2970: /*@
    2971:   PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous

    2973:   Logically Collective on sp

    2975:   Input Parameters:
    2976: + sp         - the PetscDualSpace
    2977: - continuous - flag for element continuity

    2979:   Options Database:
    2980: . -petscdualspace_lagrange_continuity <bool>

    2982:   Level: intermediate

    2984: .seealso: PetscDualSpaceLagrangeGetContinuity()
    2985: @*/
    2986: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
    2987: {

    2993:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
    2994:   return(0);
    2995: }

    2997: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
    2998: {
    2999:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

    3002:   *tensor = lag->tensorSpace;
    3003:   return(0);
    3004: }

    3006: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
    3007: {
    3008:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

    3011:   lag->tensorSpace = tensor;
    3012:   return(0);
    3013: }

    3015: static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Lagrange(PetscDualSpace sp, PetscBool *trimmed)
    3016: {
    3017:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

    3020:   *trimmed = lag->trimmed;
    3021:   return(0);
    3022: }

    3024: static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp, PetscBool trimmed)
    3025: {
    3026:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

    3029:   lag->trimmed = trimmed;
    3030:   return(0);
    3031: }

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

    3038:   if (nodeType) *nodeType = lag->nodeType;
    3039:   if (boundary) *boundary = lag->endNodes;
    3040:   if (exponent) *exponent = lag->nodeExponent;
    3041:   return(0);
    3042: }

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

    3049:   if (nodeType == PETSCDTNODES_GAUSSJACOBI && exponent <= -1.) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Exponent must be > -1");
    3050:   lag->nodeType = nodeType;
    3051:   lag->endNodes = boundary;
    3052:   lag->nodeExponent = exponent;
    3053:   return(0);
    3054: }

    3056: static PetscErrorCode PetscDualSpaceLagrangeGetUseMoments_Lagrange(PetscDualSpace sp, PetscBool *useMoments)
    3057: {
    3058:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

    3061:   *useMoments = lag->useMoments;
    3062:   return(0);
    3063: }

    3065: static PetscErrorCode PetscDualSpaceLagrangeSetUseMoments_Lagrange(PetscDualSpace sp, PetscBool useMoments)
    3066: {
    3067:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

    3070:   lag->useMoments = useMoments;
    3071:   return(0);
    3072: }

    3074: static PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt *momentOrder)
    3075: {
    3076:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

    3079:   *momentOrder = lag->momentOrder;
    3080:   return(0);
    3081: }

    3083: static PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt momentOrder)
    3084: {
    3085:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

    3088:   lag->momentOrder = momentOrder;
    3089:   return(0);
    3090: }

    3092: /*@
    3093:   PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space

    3095:   Not collective

    3097:   Input Parameter:
    3098: . sp - The PetscDualSpace

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

    3103:   Level: intermediate

    3105: .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
    3106: @*/
    3107: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
    3108: {

    3114:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));
    3115:   return(0);
    3116: }

    3118: /*@
    3119:   PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space

    3121:   Not collective

    3123:   Input Parameters:
    3124: + sp - The PetscDualSpace
    3125: - tensor - Whether the dual space has tensor layout (vs. simplicial)

    3127:   Level: intermediate

    3129: .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
    3130: @*/
    3131: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
    3132: {

    3137:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));
    3138:   return(0);
    3139: }

    3141: /*@
    3142:   PetscDualSpaceLagrangeGetTrimmed - Get the trimmed nature of the dual space

    3144:   Not collective

    3146:   Input Parameter:
    3147: . sp - The PetscDualSpace

    3149:   Output Parameter:
    3150: . 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)

    3152:   Level: intermediate

    3154: .seealso: PetscDualSpaceLagrangeSetTrimmed(), PetscDualSpaceCreate()
    3155: @*/
    3156: PetscErrorCode PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp, PetscBool *trimmed)
    3157: {

    3163:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTrimmed_C",(PetscDualSpace,PetscBool *),(sp,trimmed));
    3164:   return(0);
    3165: }

    3167: /*@
    3168:   PetscDualSpaceLagrangeSetTrimmed - Set the trimmed nature of the dual space

    3170:   Not collective

    3172:   Input Parameters:
    3173: + sp - The PetscDualSpace
    3174: - 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)

    3176:   Level: intermediate

    3178: .seealso: PetscDualSpaceLagrangeGetTrimmed(), PetscDualSpaceCreate()
    3179: @*/
    3180: PetscErrorCode PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp, PetscBool trimmed)
    3181: {

    3186:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTrimmed_C",(PetscDualSpace,PetscBool),(sp,trimmed));
    3187:   return(0);
    3188: }

    3190: /*@
    3191:   PetscDualSpaceLagrangeGetNodeType - Get a description of how nodes are laid out for Lagrange polynomials in this
    3192:   dual space

    3194:   Not collective

    3196:   Input Parameter:
    3197: . sp - The PetscDualSpace

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

    3206:   Level: advanced

    3208: .seealso: PetscDTNodeType, PetscDualSpaceLagrangeSetNodeType()
    3209: @*/
    3210: PetscErrorCode PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
    3211: {

    3219:   PetscTryMethod(sp,"PetscDualSpaceLagrangeGetNodeType_C",(PetscDualSpace,PetscDTNodeType *,PetscBool *,PetscReal *),(sp,nodeType,boundary,exponent));
    3220:   return(0);
    3221: }

    3223: /*@
    3224:   PetscDualSpaceLagrangeSetNodeType - Set a description of how nodes are laid out for Lagrange polynomials in this
    3225:   dual space

    3227:   Logically collective

    3229:   Input Parameters:
    3230: + sp - The PetscDualSpace
    3231: . nodeType - The type of nodes
    3232: . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that
    3233:              include the boundary are Gauss-Lobatto-Jacobi nodes)
    3234: - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function
    3235:              '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type

    3237:   Level: advanced

    3239: .seealso: PetscDTNodeType, PetscDualSpaceLagrangeGetNodeType()
    3240: @*/
    3241: PetscErrorCode PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
    3242: {

    3247:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetNodeType_C",(PetscDualSpace,PetscDTNodeType,PetscBool,PetscReal),(sp,nodeType,boundary,exponent));
    3248:   return(0);
    3249: }

    3251: /*@
    3252:   PetscDualSpaceLagrangeGetUseMoments - Get the flag for using moment functionals

    3254:   Not collective

    3256:   Input Parameter:
    3257: . sp - The PetscDualSpace

    3259:   Output Parameter:
    3260: . useMoments - Moment flag

    3262:   Level: advanced

    3264: .seealso: PetscDualSpaceLagrangeSetUseMoments()
    3265: @*/
    3266: PetscErrorCode PetscDualSpaceLagrangeGetUseMoments(PetscDualSpace sp, PetscBool *useMoments)
    3267: {

    3273:   PetscUseMethod(sp,"PetscDualSpaceLagrangeGetUseMoments_C",(PetscDualSpace,PetscBool *),(sp,useMoments));
    3274:   return(0);
    3275: }

    3277: /*@
    3278:   PetscDualSpaceLagrangeSetUseMoments - Set the flag for moment functionals

    3280:   Logically collective

    3282:   Input Parameters:
    3283: + sp - The PetscDualSpace
    3284: - useMoments - The flag for moment functionals

    3286:   Level: advanced

    3288: .seealso: PetscDualSpaceLagrangeGetUseMoments()
    3289: @*/
    3290: PetscErrorCode PetscDualSpaceLagrangeSetUseMoments(PetscDualSpace sp, PetscBool useMoments)
    3291: {

    3296:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetUseMoments_C",(PetscDualSpace,PetscBool),(sp,useMoments));
    3297:   return(0);
    3298: }

    3300: /*@
    3301:   PetscDualSpaceLagrangeGetMomentOrder - Get the order for moment integration

    3303:   Not collective

    3305:   Input Parameter:
    3306: . sp - The PetscDualSpace

    3308:   Output Parameter:
    3309: . order - Moment integration order

    3311:   Level: advanced

    3313: .seealso: PetscDualSpaceLagrangeSetMomentOrder()
    3314: @*/
    3315: PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder(PetscDualSpace sp, PetscInt *order)
    3316: {

    3322:   PetscUseMethod(sp,"PetscDualSpaceLagrangeGetMomentOrder_C",(PetscDualSpace,PetscInt *),(sp,order));
    3323:   return(0);
    3324: }

    3326: /*@
    3327:   PetscDualSpaceLagrangeSetMomentOrder - Set the order for moment integration

    3329:   Logically collective

    3331:   Input Parameters:
    3332: + sp - The PetscDualSpace
    3333: - order - The order for moment integration

    3335:   Level: advanced

    3337: .seealso: PetscDualSpaceLagrangeGetMomentOrder()
    3338: @*/
    3339: PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder(PetscDualSpace sp, PetscInt order)
    3340: {

    3345:   PetscTryMethod(sp,"PetscDualSpaceLagrangeSetMomentOrder_C",(PetscDualSpace,PetscInt),(sp,order));
    3346:   return(0);
    3347: }

    3349: static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
    3350: {
    3352:   sp->ops->destroy              = PetscDualSpaceDestroy_Lagrange;
    3353:   sp->ops->view                 = PetscDualSpaceView_Lagrange;
    3354:   sp->ops->setfromoptions       = PetscDualSpaceSetFromOptions_Lagrange;
    3355:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Lagrange;
    3356:   sp->ops->setup                = PetscDualSpaceSetUp_Lagrange;
    3357:   sp->ops->createheightsubspace = NULL;
    3358:   sp->ops->createpointsubspace  = NULL;
    3359:   sp->ops->getsymmetries        = PetscDualSpaceGetSymmetries_Lagrange;
    3360:   sp->ops->apply                = PetscDualSpaceApplyDefault;
    3361:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
    3362:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
    3363:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
    3364:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;
    3365:   return(0);
    3366: }

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

    3371:   Level: intermediate

    3373: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
    3374: M*/
    3375: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
    3376: {
    3377:   PetscDualSpace_Lag *lag;
    3378:   PetscErrorCode      ierr;

    3382:   PetscNewLog(sp,&lag);
    3383:   sp->data = lag;

    3385:   lag->tensorCell  = PETSC_FALSE;
    3386:   lag->tensorSpace = PETSC_FALSE;
    3387:   lag->continuous  = PETSC_TRUE;
    3388:   lag->numCopies   = PETSC_DEFAULT;
    3389:   lag->numNodeSkip = PETSC_DEFAULT;
    3390:   lag->nodeType    = PETSCDTNODES_DEFAULT;
    3391:   lag->useMoments  = PETSC_FALSE;
    3392:   lag->momentOrder = 0;

    3394:   PetscDualSpaceInitialize_Lagrange(sp);
    3395:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
    3396:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
    3397:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
    3398:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
    3399:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Lagrange);
    3400:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Lagrange);
    3401:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Lagrange);
    3402:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Lagrange);
    3403:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetUseMoments_C", PetscDualSpaceLagrangeGetUseMoments_Lagrange);
    3404:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetUseMoments_C", PetscDualSpaceLagrangeSetUseMoments_Lagrange);
    3405:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetMomentOrder_C", PetscDualSpaceLagrangeGetMomentOrder_Lagrange);
    3406:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetMomentOrder_C", PetscDualSpaceLagrangeSetMomentOrder_Lagrange);
    3407:   return(0);
    3408: }