Actual source code: dspacelagrange.c
1: #include <petsc/private/petscfeimpl.h>
2: #include <petscdmplex.h>
3: #include <petscblaslapack.h>
5: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]);
7: struct _n_Petsc1DNodeFamily {
8: PetscInt refct;
9: PetscDTNodeType nodeFamily;
10: PetscReal gaussJacobiExp;
11: PetscInt nComputed;
12: PetscReal **nodesets;
13: PetscBool endpoints;
14: };
16: /* users set node families for PETSCDUALSPACELAGRANGE with just the inputs to this function, but internally we create
17: * an object that can cache the computations across multiple dual spaces */
18: static PetscErrorCode Petsc1DNodeFamilyCreate(PetscDTNodeType family, PetscReal gaussJacobiExp, PetscBool endpoints, Petsc1DNodeFamily *nf)
19: {
20: Petsc1DNodeFamily f;
22: PetscNew(&f);
23: switch (family) {
24: case PETSCDTNODES_GAUSSJACOBI:
25: case PETSCDTNODES_EQUISPACED:
26: f->nodeFamily = family;
27: break;
28: default:
29: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family");
30: }
31: f->endpoints = endpoints;
32: f->gaussJacobiExp = 0.;
33: if (family == PETSCDTNODES_GAUSSJACOBI) {
35: f->gaussJacobiExp = gaussJacobiExp;
36: }
37: f->refct = 1;
38: *nf = f;
39: return 0;
40: }
42: static PetscErrorCode Petsc1DNodeFamilyReference(Petsc1DNodeFamily nf)
43: {
44: if (nf) nf->refct++;
45: return 0;
46: }
48: static PetscErrorCode Petsc1DNodeFamilyDestroy(Petsc1DNodeFamily *nf)
49: {
50: PetscInt i, nc;
52: if (!(*nf)) return 0;
53: if (--(*nf)->refct > 0) {
54: *nf = NULL;
55: return 0;
56: }
57: nc = (*nf)->nComputed;
58: for (i = 0; i < nc; i++) PetscFree((*nf)->nodesets[i]);
59: PetscFree((*nf)->nodesets);
60: PetscFree(*nf);
61: *nf = NULL;
62: return 0;
63: }
65: static PetscErrorCode Petsc1DNodeFamilyGetNodeSets(Petsc1DNodeFamily f, PetscInt degree, PetscReal ***nodesets)
66: {
67: PetscInt nc;
69: nc = f->nComputed;
70: if (degree >= nc) {
71: PetscInt i, j;
72: PetscReal **new_nodesets;
73: PetscReal *w;
75: PetscMalloc1(degree + 1, &new_nodesets);
76: PetscArraycpy(new_nodesets, f->nodesets, nc);
77: PetscFree(f->nodesets);
78: f->nodesets = new_nodesets;
79: PetscMalloc1(degree + 1, &w);
80: for (i = nc; i < degree + 1; i++) {
81: PetscMalloc1(i + 1, &(f->nodesets[i]));
82: if (!i) {
83: f->nodesets[i][0] = 0.5;
84: } else {
85: switch (f->nodeFamily) {
86: case PETSCDTNODES_EQUISPACED:
87: if (f->endpoints) {
88: for (j = 0; j <= i; j++) f->nodesets[i][j] = (PetscReal)j / (PetscReal)i;
89: } else {
90: /* these nodes are at the centroids of the small simplices created by the equispaced nodes that include
91: * the endpoints */
92: for (j = 0; j <= i; j++) f->nodesets[i][j] = ((PetscReal)j + 0.5) / ((PetscReal)i + 1.);
93: }
94: break;
95: case PETSCDTNODES_GAUSSJACOBI:
96: if (f->endpoints) {
97: PetscDTGaussLobattoJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w);
98: } else {
99: PetscDTGaussJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w);
100: }
101: break;
102: default:
103: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family");
104: }
105: }
106: }
107: PetscFree(w);
108: f->nComputed = degree + 1;
109: }
110: *nodesets = f->nodesets;
111: return 0;
112: }
114: /* http://arxiv.org/abs/2002.09421 for details */
115: static PetscErrorCode PetscNodeRecursive_Internal(PetscInt dim, PetscInt degree, PetscReal **nodesets, PetscInt tup[], PetscReal node[])
116: {
117: PetscReal w;
118: PetscInt i, j;
121: w = 0.;
122: if (dim == 1) {
123: node[0] = nodesets[degree][tup[0]];
124: node[1] = nodesets[degree][tup[1]];
125: } else {
126: for (i = 0; i < dim + 1; i++) node[i] = 0.;
127: for (i = 0; i < dim + 1; i++) {
128: PetscReal wi = nodesets[degree][degree - tup[i]];
130: for (j = 0; j < dim + 1; j++) tup[dim + 1 + j] = tup[j + (j >= i)];
131: PetscNodeRecursive_Internal(dim - 1, degree - tup[i], nodesets, &tup[dim + 1], &node[dim + 1]);
132: for (j = 0; j < dim + 1; j++) node[j + (j >= i)] += wi * node[dim + 1 + j];
133: w += wi;
134: }
135: for (i = 0; i < dim + 1; i++) node[i] /= w;
136: }
137: return 0;
138: }
140: /* compute simplex nodes for the biunit simplex from the 1D node family */
141: static PetscErrorCode Petsc1DNodeFamilyComputeSimplexNodes(Petsc1DNodeFamily f, PetscInt dim, PetscInt degree, PetscReal points[])
142: {
143: PetscInt *tup;
144: PetscInt k;
145: PetscInt npoints;
146: PetscReal **nodesets = NULL;
147: PetscInt worksize;
148: PetscReal *nodework;
149: PetscInt *tupwork;
153: if (!dim) return 0;
154: PetscCalloc1(dim + 2, &tup);
155: k = 0;
156: PetscDTBinomialInt(degree + dim, dim, &npoints);
157: Petsc1DNodeFamilyGetNodeSets(f, degree, &nodesets);
158: worksize = ((dim + 2) * (dim + 3)) / 2;
159: PetscMalloc2(worksize, &nodework, worksize, &tupwork);
160: /* loop over the tuples of length dim with sum at most degree */
161: for (k = 0; k < npoints; k++) {
162: PetscInt i;
164: /* turn thm into tuples of length dim + 1 with sum equal to degree (barycentric indice) */
165: tup[0] = degree;
166: for (i = 0; i < dim; i++) tup[0] -= tup[i + 1];
167: switch (f->nodeFamily) {
168: case PETSCDTNODES_EQUISPACED:
169: /* compute equispaces nodes on the unit reference triangle */
170: if (f->endpoints) {
171: for (i = 0; i < dim; i++) points[dim * k + i] = (PetscReal)tup[i + 1] / (PetscReal)degree;
172: } else {
173: for (i = 0; i < dim; i++) {
174: /* these nodes are at the centroids of the small simplices created by the equispaced nodes that include
175: * the endpoints */
176: points[dim * k + i] = ((PetscReal)tup[i + 1] + 1. / (dim + 1.)) / (PetscReal)(degree + 1.);
177: }
178: }
179: break;
180: default:
181: /* compute equispaces nodes on the barycentric reference triangle (the trace on the first dim dimensions are the
182: * unit reference triangle nodes */
183: for (i = 0; i < dim + 1; i++) tupwork[i] = tup[i];
184: PetscNodeRecursive_Internal(dim, degree, nodesets, tupwork, nodework);
185: for (i = 0; i < dim; i++) points[dim * k + i] = nodework[i + 1];
186: break;
187: }
188: PetscDualSpaceLatticePointLexicographic_Internal(dim, degree, &tup[1]);
189: }
190: /* map from unit simplex to biunit simplex */
191: for (k = 0; k < npoints * dim; k++) points[k] = points[k] * 2. - 1.;
192: PetscFree2(nodework, tupwork);
193: PetscFree(tup);
194: return 0;
195: }
197: /* If we need to get the dofs from a mesh point, or add values into dofs at a mesh point, and there is more than one dof
198: * on that mesh point, we have to be careful about getting/adding everything in the right place.
199: *
200: * With nodal dofs like PETSCDUALSPACELAGRANGE makes, the general approach to calculate the value of dofs associate
201: * with a node A is
202: * - transform the node locations x(A) by the map that takes the mesh point to its reorientation, x' = phi(x(A))
203: * - figure out which node was originally at the location of the transformed point, A' = idx(x')
204: * - if the dofs are not scalars, figure out how to represent the transformed dofs in terms of the basis
205: * of dofs at A' (using pushforward/pullback rules)
206: *
207: * The one sticky point with this approach is the "A' = idx(x')" step: trying to go from real valued coordinates
208: * back to indices. I don't want to rely on floating point tolerances. Additionally, PETSCDUALSPACELAGRANGE may
209: * eventually support quasi-Lagrangian dofs, which could involve quadrature at multiple points, so the location "x(A)"
210: * would be ambiguous.
211: *
212: * So each dof gets an integer value coordinate (nodeIdx in the structure below). The choice of integer coordinates
213: * is somewhat arbitrary, as long as all of the relevant symmetries of the mesh point correspond to *permutations* of
214: * the integer coordinates, which do not depend on numerical precision.
215: *
216: * So
217: *
218: * - DMPlexGetTransitiveClosure_Internal() tells me how an orientation turns into a permutation of the vertices of a
219: * mesh point
220: * - The permutation of the vertices, and the nodeIdx values assigned to them, tells what permutation in index space
221: * is associated with the orientation
222: * - I uses that permutation to get xi' = phi(xi(A)), the integer coordinate of the transformed dof
223: * - I can without numerical issues compute A' = idx(xi')
224: *
225: * Here are some examples of how the process works
226: *
227: * - With a triangle:
228: *
229: * The triangle has the following integer coordinates for vertices, taken from the barycentric triangle
230: *
231: * closure order 2
232: * nodeIdx (0,0,1)
233: * \
234: * +
235: * |\
236: * | \
237: * | \
238: * | \ closure order 1
239: * | \ / nodeIdx (0,1,0)
240: * +-----+
241: * \
242: * closure order 0
243: * nodeIdx (1,0,0)
244: *
245: * If I do DMPlexGetTransitiveClosure_Internal() with orientation 1, the vertices would appear
246: * in the order (1, 2, 0)
247: *
248: * If I list the nodeIdx of each vertex in closure order for orientation 0 (0, 1, 2) and orientation 1 (1, 2, 0), I
249: * see
250: *
251: * orientation 0 | orientation 1
252: *
253: * [0] (1,0,0) [1] (0,1,0)
254: * [1] (0,1,0) [2] (0,0,1)
255: * [2] (0,0,1) [0] (1,0,0)
256: * A B
257: *
258: * In other words, B is the result of a row permutation of A. But, there is also
259: * a column permutation that accomplishes the same result, (2,0,1).
260: *
261: * So if a dof has nodeIdx coordinate (a,b,c), after the transformation its nodeIdx coordinate
262: * is (c,a,b), and the transformed degree of freedom will be a linear combination of dofs
263: * that originally had coordinate (c,a,b).
264: *
265: * - With a quadrilateral:
266: *
267: * The quadrilateral has the following integer coordinates for vertices, taken from concatenating barycentric
268: * coordinates for two segments:
269: *
270: * closure order 3 closure order 2
271: * nodeIdx (1,0,0,1) nodeIdx (0,1,0,1)
272: * \ /
273: * +----+
274: * | |
275: * | |
276: * +----+
277: * / \
278: * closure order 0 closure order 1
279: * nodeIdx (1,0,1,0) nodeIdx (0,1,1,0)
280: *
281: * If I do DMPlexGetTransitiveClosure_Internal() with orientation 1, the vertices would appear
282: * in the order (1, 2, 3, 0)
283: *
284: * If I list the nodeIdx of each vertex in closure order for orientation 0 (0, 1, 2, 3) and
285: * orientation 1 (1, 2, 3, 0), I see
286: *
287: * orientation 0 | orientation 1
288: *
289: * [0] (1,0,1,0) [1] (0,1,1,0)
290: * [1] (0,1,1,0) [2] (0,1,0,1)
291: * [2] (0,1,0,1) [3] (1,0,0,1)
292: * [3] (1,0,0,1) [0] (1,0,1,0)
293: * A B
294: *
295: * The column permutation that accomplishes the same result is (3,2,0,1).
296: *
297: * So if a dof has nodeIdx coordinate (a,b,c,d), after the transformation its nodeIdx coordinate
298: * is (d,c,a,b), and the transformed degree of freedom will be a linear combination of dofs
299: * that originally had coordinate (d,c,a,b).
300: *
301: * Previously PETSCDUALSPACELAGRANGE had hardcoded symmetries for the triangle and quadrilateral,
302: * but this approach will work for any polytope, such as the wedge (triangular prism).
303: */
304: struct _n_PetscLagNodeIndices {
305: PetscInt refct;
306: PetscInt nodeIdxDim;
307: PetscInt nodeVecDim;
308: PetscInt nNodes;
309: PetscInt *nodeIdx; /* for each node an index of size nodeIdxDim */
310: PetscReal *nodeVec; /* for each node a vector of size nodeVecDim */
311: PetscInt *perm; /* if these are vertices, perm takes DMPlex point index to closure order;
312: if these are nodes, perm lists nodes in index revlex order */
313: };
315: /* this is just here so I can access the values in tests/ex1.c outside the library */
316: PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices ni, PetscInt *nodeIdxDim, PetscInt *nodeVecDim, PetscInt *nNodes, const PetscInt *nodeIdx[], const PetscReal *nodeVec[])
317: {
318: *nodeIdxDim = ni->nodeIdxDim;
319: *nodeVecDim = ni->nodeVecDim;
320: *nNodes = ni->nNodes;
321: *nodeIdx = ni->nodeIdx;
322: *nodeVec = ni->nodeVec;
323: return 0;
324: }
326: static PetscErrorCode PetscLagNodeIndicesReference(PetscLagNodeIndices ni)
327: {
328: if (ni) ni->refct++;
329: return 0;
330: }
332: static PetscErrorCode PetscLagNodeIndicesDuplicate(PetscLagNodeIndices ni, PetscLagNodeIndices *niNew)
333: {
334: PetscNew(niNew);
335: (*niNew)->refct = 1;
336: (*niNew)->nodeIdxDim = ni->nodeIdxDim;
337: (*niNew)->nodeVecDim = ni->nodeVecDim;
338: (*niNew)->nNodes = ni->nNodes;
339: PetscMalloc1(ni->nNodes * ni->nodeIdxDim, &((*niNew)->nodeIdx));
340: PetscArraycpy((*niNew)->nodeIdx, ni->nodeIdx, ni->nNodes * ni->nodeIdxDim);
341: PetscMalloc1(ni->nNodes * ni->nodeVecDim, &((*niNew)->nodeVec));
342: PetscArraycpy((*niNew)->nodeVec, ni->nodeVec, ni->nNodes * ni->nodeVecDim);
343: (*niNew)->perm = NULL;
344: return 0;
345: }
347: static PetscErrorCode PetscLagNodeIndicesDestroy(PetscLagNodeIndices *ni)
348: {
349: if (!(*ni)) return 0;
350: if (--(*ni)->refct > 0) {
351: *ni = NULL;
352: return 0;
353: }
354: PetscFree((*ni)->nodeIdx);
355: PetscFree((*ni)->nodeVec);
356: PetscFree((*ni)->perm);
357: PetscFree(*ni);
358: *ni = NULL;
359: return 0;
360: }
362: /* The vertices are given nodeIdx coordinates (e.g. the corners of the barycentric triangle). Those coordinates are
363: * in some other order, and to understand the effect of different symmetries, we need them to be in closure order.
364: *
365: * If sortIdx is PETSC_FALSE, the coordinates are already in revlex order, otherwise we must sort them
366: * to that order before we do the real work of this function, which is
367: *
368: * - mark the vertices in closure order
369: * - sort them in revlex order
370: * - use the resulting permutation to list the vertex coordinates in closure order
371: */
372: static PetscErrorCode PetscLagNodeIndicesComputeVertexOrder(DM dm, PetscLagNodeIndices ni, PetscBool sortIdx)
373: {
374: PetscInt v, w, vStart, vEnd, c, d;
375: PetscInt nVerts;
376: PetscInt closureSize = 0;
377: PetscInt *closure = NULL;
378: PetscInt *closureOrder;
379: PetscInt *invClosureOrder;
380: PetscInt *revlexOrder;
381: PetscInt *newNodeIdx;
382: PetscInt dim;
383: Vec coordVec;
384: const PetscScalar *coords;
386: DMGetDimension(dm, &dim);
387: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
388: nVerts = vEnd - vStart;
389: PetscMalloc1(nVerts, &closureOrder);
390: PetscMalloc1(nVerts, &invClosureOrder);
391: PetscMalloc1(nVerts, &revlexOrder);
392: if (sortIdx) { /* bubble sort nodeIdx into revlex order */
393: PetscInt nodeIdxDim = ni->nodeIdxDim;
394: PetscInt *idxOrder;
396: PetscMalloc1(nVerts * nodeIdxDim, &newNodeIdx);
397: PetscMalloc1(nVerts, &idxOrder);
398: for (v = 0; v < nVerts; v++) idxOrder[v] = v;
399: for (v = 0; v < nVerts; v++) {
400: for (w = v + 1; w < nVerts; w++) {
401: const PetscInt *iv = &(ni->nodeIdx[idxOrder[v] * nodeIdxDim]);
402: const PetscInt *iw = &(ni->nodeIdx[idxOrder[w] * nodeIdxDim]);
403: PetscInt diff = 0;
405: for (d = nodeIdxDim - 1; d >= 0; d--)
406: if ((diff = (iv[d] - iw[d]))) break;
407: if (diff > 0) {
408: PetscInt swap = idxOrder[v];
410: idxOrder[v] = idxOrder[w];
411: idxOrder[w] = swap;
412: }
413: }
414: }
415: for (v = 0; v < nVerts; v++) {
416: for (d = 0; d < nodeIdxDim; d++) newNodeIdx[v * ni->nodeIdxDim + d] = ni->nodeIdx[idxOrder[v] * nodeIdxDim + d];
417: }
418: PetscFree(ni->nodeIdx);
419: ni->nodeIdx = newNodeIdx;
420: newNodeIdx = NULL;
421: PetscFree(idxOrder);
422: }
423: DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);
424: c = closureSize - nVerts;
425: for (v = 0; v < nVerts; v++) closureOrder[v] = closure[2 * (c + v)] - vStart;
426: for (v = 0; v < nVerts; v++) invClosureOrder[closureOrder[v]] = v;
427: DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);
428: DMGetCoordinatesLocal(dm, &coordVec);
429: VecGetArrayRead(coordVec, &coords);
430: /* bubble sort closure vertices by coordinates in revlex order */
431: for (v = 0; v < nVerts; v++) revlexOrder[v] = v;
432: for (v = 0; v < nVerts; v++) {
433: for (w = v + 1; w < nVerts; w++) {
434: const PetscScalar *cv = &coords[closureOrder[revlexOrder[v]] * dim];
435: const PetscScalar *cw = &coords[closureOrder[revlexOrder[w]] * dim];
436: PetscReal diff = 0;
438: for (d = dim - 1; d >= 0; d--)
439: if ((diff = PetscRealPart(cv[d] - cw[d])) != 0.) break;
440: if (diff > 0.) {
441: PetscInt swap = revlexOrder[v];
443: revlexOrder[v] = revlexOrder[w];
444: revlexOrder[w] = swap;
445: }
446: }
447: }
448: VecRestoreArrayRead(coordVec, &coords);
449: PetscMalloc1(ni->nodeIdxDim * nVerts, &newNodeIdx);
450: /* reorder nodeIdx to be in closure order */
451: for (v = 0; v < nVerts; v++) {
452: for (d = 0; d < ni->nodeIdxDim; d++) newNodeIdx[revlexOrder[v] * ni->nodeIdxDim + d] = ni->nodeIdx[v * ni->nodeIdxDim + d];
453: }
454: PetscFree(ni->nodeIdx);
455: ni->nodeIdx = newNodeIdx;
456: ni->perm = invClosureOrder;
457: PetscFree(revlexOrder);
458: PetscFree(closureOrder);
459: return 0;
460: }
462: /* the coordinates of the simplex vertices are the corners of the barycentric simplex.
463: * When we stack them on top of each other in revlex order, they look like the identity matrix */
464: static PetscErrorCode PetscLagNodeIndicesCreateSimplexVertices(DM dm, PetscLagNodeIndices *nodeIndices)
465: {
466: PetscLagNodeIndices ni;
467: PetscInt dim, d;
469: PetscNew(&ni);
470: DMGetDimension(dm, &dim);
471: ni->nodeIdxDim = dim + 1;
472: ni->nodeVecDim = 0;
473: ni->nNodes = dim + 1;
474: ni->refct = 1;
475: PetscCalloc1((dim + 1) * (dim + 1), &(ni->nodeIdx));
476: for (d = 0; d < dim + 1; d++) ni->nodeIdx[d * (dim + 2)] = 1;
477: PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_FALSE);
478: *nodeIndices = ni;
479: return 0;
480: }
482: /* A polytope that is a tensor product of a facet and a segment.
483: * We take whatever coordinate system was being used for the facet
484: * and we concatenate the barycentric coordinates for the vertices
485: * at the end of the segment, (1,0) and (0,1), to get a coordinate
486: * system for the tensor product element */
487: static PetscErrorCode PetscLagNodeIndicesCreateTensorVertices(DM dm, PetscLagNodeIndices facetni, PetscLagNodeIndices *nodeIndices)
488: {
489: PetscLagNodeIndices ni;
490: PetscInt nodeIdxDim, subNodeIdxDim = facetni->nodeIdxDim;
491: PetscInt nVerts, nSubVerts = facetni->nNodes;
492: PetscInt dim, d, e, f, g;
494: PetscNew(&ni);
495: DMGetDimension(dm, &dim);
496: ni->nodeIdxDim = nodeIdxDim = subNodeIdxDim + 2;
497: ni->nodeVecDim = 0;
498: ni->nNodes = nVerts = 2 * nSubVerts;
499: ni->refct = 1;
500: PetscCalloc1(nodeIdxDim * nVerts, &(ni->nodeIdx));
501: for (f = 0, d = 0; d < 2; d++) {
502: for (e = 0; e < nSubVerts; e++, f++) {
503: for (g = 0; g < subNodeIdxDim; g++) ni->nodeIdx[f * nodeIdxDim + g] = facetni->nodeIdx[e * subNodeIdxDim + g];
504: ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim] = (1 - d);
505: ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim + 1] = d;
506: }
507: }
508: PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_TRUE);
509: *nodeIndices = ni;
510: return 0;
511: }
513: /* This helps us compute symmetries, and it also helps us compute coordinates for dofs that are being pushed
514: * forward from a boundary mesh point.
515: *
516: * Input:
517: *
518: * dm - the target reference cell where we want new coordinates and dof directions to be valid
519: * vert - the vertex coordinate system for the target reference cell
520: * p - the point in the target reference cell that the dofs are coming from
521: * vertp - the vertex coordinate system for p's reference cell
522: * ornt - the resulting coordinates and dof vectors will be for p under this orientation
523: * nodep - the node coordinates and dof vectors in p's reference cell
524: * formDegree - the form degree that the dofs transform as
525: *
526: * Output:
527: *
528: * pfNodeIdx - the node coordinates for p's dofs, in the dm reference cell, from the ornt perspective
529: * pfNodeVec - the node dof vectors for p's dofs, in the dm reference cell, from the ornt perspective
530: */
531: static PetscErrorCode PetscLagNodeIndicesPushForward(DM dm, PetscLagNodeIndices vert, PetscInt p, PetscLagNodeIndices vertp, PetscLagNodeIndices nodep, PetscInt ornt, PetscInt formDegree, PetscInt pfNodeIdx[], PetscReal pfNodeVec[])
532: {
533: PetscInt *closureVerts;
534: PetscInt closureSize = 0;
535: PetscInt *closure = NULL;
536: PetscInt dim, pdim, c, i, j, k, n, v, vStart, vEnd;
537: PetscInt nSubVert = vertp->nNodes;
538: PetscInt nodeIdxDim = vert->nodeIdxDim;
539: PetscInt subNodeIdxDim = vertp->nodeIdxDim;
540: PetscInt nNodes = nodep->nNodes;
541: const PetscInt *vertIdx = vert->nodeIdx;
542: const PetscInt *subVertIdx = vertp->nodeIdx;
543: const PetscInt *nodeIdx = nodep->nodeIdx;
544: const PetscReal *nodeVec = nodep->nodeVec;
545: PetscReal *J, *Jstar;
546: PetscReal detJ;
547: PetscInt depth, pdepth, Nk, pNk;
548: Vec coordVec;
549: PetscScalar *newCoords = NULL;
550: const PetscScalar *oldCoords = NULL;
552: DMGetDimension(dm, &dim);
553: DMPlexGetDepth(dm, &depth);
554: DMGetCoordinatesLocal(dm, &coordVec);
555: DMPlexGetPointDepth(dm, p, &pdepth);
556: pdim = pdepth != depth ? pdepth != 0 ? pdepth : 0 : dim;
557: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
558: DMGetWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);
559: DMPlexGetTransitiveClosure_Internal(dm, p, ornt, PETSC_TRUE, &closureSize, &closure);
560: c = closureSize - nSubVert;
561: /* we want which cell closure indices the closure of this point corresponds to */
562: for (v = 0; v < nSubVert; v++) closureVerts[v] = vert->perm[closure[2 * (c + v)] - vStart];
563: DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
564: /* push forward indices */
565: for (i = 0; i < nodeIdxDim; i++) { /* for every component of the target index space */
566: /* check if this is a component that all vertices around this point have in common */
567: for (j = 1; j < nSubVert; j++) {
568: if (vertIdx[closureVerts[j] * nodeIdxDim + i] != vertIdx[closureVerts[0] * nodeIdxDim + i]) break;
569: }
570: if (j == nSubVert) { /* all vertices have this component in common, directly copy to output */
571: PetscInt val = vertIdx[closureVerts[0] * nodeIdxDim + i];
572: for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = val;
573: } else {
574: PetscInt subi = -1;
575: /* there must be a component in vertp that looks the same */
576: for (k = 0; k < subNodeIdxDim; k++) {
577: for (j = 0; j < nSubVert; j++) {
578: if (vertIdx[closureVerts[j] * nodeIdxDim + i] != subVertIdx[j * subNodeIdxDim + k]) break;
579: }
580: if (j == nSubVert) {
581: subi = k;
582: break;
583: }
584: }
586: /* that component in the vertp system becomes component i in the vert system for each dof */
587: for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = nodeIdx[n * subNodeIdxDim + subi];
588: }
589: }
590: /* push forward vectors */
591: DMGetWorkArray(dm, dim * dim, MPIU_REAL, &J);
592: if (ornt != 0) { /* temporarily change the coordinate vector so
593: DMPlexComputeCellGeometryAffineFEM gives us the Jacobian we want */
594: PetscInt closureSize2 = 0;
595: PetscInt *closure2 = NULL;
597: DMPlexGetTransitiveClosure_Internal(dm, p, 0, PETSC_TRUE, &closureSize2, &closure2);
598: PetscMalloc1(dim * nSubVert, &newCoords);
599: VecGetArrayRead(coordVec, &oldCoords);
600: for (v = 0; v < nSubVert; v++) {
601: PetscInt d;
602: for (d = 0; d < dim; d++) newCoords[(closure2[2 * (c + v)] - vStart) * dim + d] = oldCoords[closureVerts[v] * dim + d];
603: }
604: VecRestoreArrayRead(coordVec, &oldCoords);
605: DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize2, &closure2);
606: VecPlaceArray(coordVec, newCoords);
607: }
608: DMPlexComputeCellGeometryAffineFEM(dm, p, NULL, J, NULL, &detJ);
609: if (ornt != 0) {
610: VecResetArray(coordVec);
611: PetscFree(newCoords);
612: }
613: DMRestoreWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);
614: /* compactify */
615: for (i = 0; i < dim; i++)
616: for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
617: /* We have the Jacobian mapping the point's reference cell to this reference cell:
618: * pulling back a function to the point and applying the dof is what we want,
619: * so we get the pullback matrix and multiply the dof by that matrix on the right */
620: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
621: PetscDTBinomialInt(pdim, PetscAbsInt(formDegree), &pNk);
622: DMGetWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);
623: PetscDTAltVPullbackMatrix(pdim, dim, J, formDegree, Jstar);
624: for (n = 0; n < nNodes; n++) {
625: for (i = 0; i < Nk; i++) {
626: PetscReal val = 0.;
627: for (j = 0; j < pNk; j++) val += nodeVec[n * pNk + j] * Jstar[j * Nk + i];
628: pfNodeVec[n * Nk + i] = val;
629: }
630: }
631: DMRestoreWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);
632: DMRestoreWorkArray(dm, dim * dim, MPIU_REAL, &J);
633: return 0;
634: }
636: /* given to sets of nodes, take the tensor product, where the product of the dof indices is concatenation and the
637: * product of the dof vectors is the wedge product */
638: static PetscErrorCode PetscLagNodeIndicesTensor(PetscLagNodeIndices tracei, PetscInt dimT, PetscInt kT, PetscLagNodeIndices fiberi, PetscInt dimF, PetscInt kF, PetscLagNodeIndices *nodeIndices)
639: {
640: PetscInt dim = dimT + dimF;
641: PetscInt nodeIdxDim, nNodes;
642: PetscInt formDegree = kT + kF;
643: PetscInt Nk, NkT, NkF;
644: PetscInt MkT, MkF;
645: PetscLagNodeIndices ni;
646: PetscInt i, j, l;
647: PetscReal *projF, *projT;
648: PetscReal *projFstar, *projTstar;
649: PetscReal *workF, *workF2, *workT, *workT2, *work, *work2;
650: PetscReal *wedgeMat;
651: PetscReal sign;
653: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
654: PetscDTBinomialInt(dimT, PetscAbsInt(kT), &NkT);
655: PetscDTBinomialInt(dimF, PetscAbsInt(kF), &NkF);
656: PetscDTBinomialInt(dim, PetscAbsInt(kT), &MkT);
657: PetscDTBinomialInt(dim, PetscAbsInt(kF), &MkF);
658: PetscNew(&ni);
659: ni->nodeIdxDim = nodeIdxDim = tracei->nodeIdxDim + fiberi->nodeIdxDim;
660: ni->nodeVecDim = Nk;
661: ni->nNodes = nNodes = tracei->nNodes * fiberi->nNodes;
662: ni->refct = 1;
663: PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
664: /* first concatenate the indices */
665: for (l = 0, j = 0; j < fiberi->nNodes; j++) {
666: for (i = 0; i < tracei->nNodes; i++, l++) {
667: PetscInt m, n = 0;
669: for (m = 0; m < tracei->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = tracei->nodeIdx[i * tracei->nodeIdxDim + m];
670: for (m = 0; m < fiberi->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = fiberi->nodeIdx[j * fiberi->nodeIdxDim + m];
671: }
672: }
674: /* now wedge together the push-forward vectors */
675: PetscMalloc1(nNodes * Nk, &(ni->nodeVec));
676: PetscCalloc2(dimT * dim, &projT, dimF * dim, &projF);
677: for (i = 0; i < dimT; i++) projT[i * (dim + 1)] = 1.;
678: for (i = 0; i < dimF; i++) projF[i * (dim + dimT + 1) + dimT] = 1.;
679: PetscMalloc2(MkT * NkT, &projTstar, MkF * NkF, &projFstar);
680: PetscDTAltVPullbackMatrix(dim, dimT, projT, kT, projTstar);
681: PetscDTAltVPullbackMatrix(dim, dimF, projF, kF, projFstar);
682: PetscMalloc6(MkT, &workT, MkT, &workT2, MkF, &workF, MkF, &workF2, Nk, &work, Nk, &work2);
683: PetscMalloc1(Nk * MkT, &wedgeMat);
684: sign = (PetscAbsInt(kT * kF) & 1) ? -1. : 1.;
685: for (l = 0, j = 0; j < fiberi->nNodes; j++) {
686: PetscInt d, e;
688: /* push forward fiber k-form */
689: for (d = 0; d < MkF; d++) {
690: PetscReal val = 0.;
691: for (e = 0; e < NkF; e++) val += projFstar[d * NkF + e] * fiberi->nodeVec[j * NkF + e];
692: workF[d] = val;
693: }
694: /* Hodge star to proper form if necessary */
695: if (kF < 0) {
696: for (d = 0; d < MkF; d++) workF2[d] = workF[d];
697: PetscDTAltVStar(dim, PetscAbsInt(kF), 1, workF2, workF);
698: }
699: /* Compute the matrix that wedges this form with one of the trace k-form */
700: PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kF), PetscAbsInt(kT), workF, wedgeMat);
701: for (i = 0; i < tracei->nNodes; i++, l++) {
702: /* push forward trace k-form */
703: for (d = 0; d < MkT; d++) {
704: PetscReal val = 0.;
705: for (e = 0; e < NkT; e++) val += projTstar[d * NkT + e] * tracei->nodeVec[i * NkT + e];
706: workT[d] = val;
707: }
708: /* Hodge star to proper form if necessary */
709: if (kT < 0) {
710: for (d = 0; d < MkT; d++) workT2[d] = workT[d];
711: PetscDTAltVStar(dim, PetscAbsInt(kT), 1, workT2, workT);
712: }
713: /* compute the wedge product of the push-forward trace form and firer forms */
714: for (d = 0; d < Nk; d++) {
715: PetscReal val = 0.;
716: for (e = 0; e < MkT; e++) val += wedgeMat[d * MkT + e] * workT[e];
717: work[d] = val;
718: }
719: /* inverse Hodge star from proper form if necessary */
720: if (formDegree < 0) {
721: for (d = 0; d < Nk; d++) work2[d] = work[d];
722: PetscDTAltVStar(dim, PetscAbsInt(formDegree), -1, work2, work);
723: }
724: /* insert into the array (adjusting for sign) */
725: for (d = 0; d < Nk; d++) ni->nodeVec[l * Nk + d] = sign * work[d];
726: }
727: }
728: PetscFree(wedgeMat);
729: PetscFree6(workT, workT2, workF, workF2, work, work2);
730: PetscFree2(projTstar, projFstar);
731: PetscFree2(projT, projF);
732: *nodeIndices = ni;
733: return 0;
734: }
736: /* simple union of two sets of nodes */
737: static PetscErrorCode PetscLagNodeIndicesMerge(PetscLagNodeIndices niA, PetscLagNodeIndices niB, PetscLagNodeIndices *nodeIndices)
738: {
739: PetscLagNodeIndices ni;
740: PetscInt nodeIdxDim, nodeVecDim, nNodes;
742: PetscNew(&ni);
743: ni->nodeIdxDim = nodeIdxDim = niA->nodeIdxDim;
745: ni->nodeVecDim = nodeVecDim = niA->nodeVecDim;
747: ni->nNodes = nNodes = niA->nNodes + niB->nNodes;
748: ni->refct = 1;
749: PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
750: PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
751: PetscArraycpy(ni->nodeIdx, niA->nodeIdx, niA->nNodes * nodeIdxDim);
752: PetscArraycpy(ni->nodeVec, niA->nodeVec, niA->nNodes * nodeVecDim);
753: PetscArraycpy(&(ni->nodeIdx[niA->nNodes * nodeIdxDim]), niB->nodeIdx, niB->nNodes * nodeIdxDim);
754: PetscArraycpy(&(ni->nodeVec[niA->nNodes * nodeVecDim]), niB->nodeVec, niB->nNodes * nodeVecDim);
755: *nodeIndices = ni;
756: return 0;
757: }
759: #define PETSCTUPINTCOMPREVLEX(N) \
760: static int PetscConcat_(PetscTupIntCompRevlex_, N)(const void *a, const void *b) \
761: { \
762: const PetscInt *A = (const PetscInt *)a; \
763: const PetscInt *B = (const PetscInt *)b; \
764: int i; \
765: PetscInt diff = 0; \
766: for (i = 0; i < N; i++) { \
767: diff = A[N - i] - B[N - i]; \
768: if (diff) break; \
769: } \
770: return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1; \
771: }
773: PETSCTUPINTCOMPREVLEX(3)
774: PETSCTUPINTCOMPREVLEX(4)
775: PETSCTUPINTCOMPREVLEX(5)
776: PETSCTUPINTCOMPREVLEX(6)
777: PETSCTUPINTCOMPREVLEX(7)
779: static int PetscTupIntCompRevlex_N(const void *a, const void *b)
780: {
781: const PetscInt *A = (const PetscInt *)a;
782: const PetscInt *B = (const PetscInt *)b;
783: int i;
784: int N = A[0];
785: PetscInt diff = 0;
786: for (i = 0; i < N; i++) {
787: diff = A[N - i] - B[N - i];
788: if (diff) break;
789: }
790: return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1;
791: }
793: /* The nodes are not necessarily in revlex order wrt nodeIdx: get the permutation
794: * that puts them in that order */
795: static PetscErrorCode PetscLagNodeIndicesGetPermutation(PetscLagNodeIndices ni, PetscInt *perm[])
796: {
797: if (!(ni->perm)) {
798: PetscInt *sorter;
799: PetscInt m = ni->nNodes;
800: PetscInt nodeIdxDim = ni->nodeIdxDim;
801: PetscInt i, j, k, l;
802: PetscInt *prm;
803: int (*comp)(const void *, const void *);
805: PetscMalloc1((nodeIdxDim + 2) * m, &sorter);
806: for (k = 0, l = 0, i = 0; i < m; i++) {
807: sorter[k++] = nodeIdxDim + 1;
808: sorter[k++] = i;
809: for (j = 0; j < nodeIdxDim; j++) sorter[k++] = ni->nodeIdx[l++];
810: }
811: switch (nodeIdxDim) {
812: case 2:
813: comp = PetscTupIntCompRevlex_3;
814: break;
815: case 3:
816: comp = PetscTupIntCompRevlex_4;
817: break;
818: case 4:
819: comp = PetscTupIntCompRevlex_5;
820: break;
821: case 5:
822: comp = PetscTupIntCompRevlex_6;
823: break;
824: case 6:
825: comp = PetscTupIntCompRevlex_7;
826: break;
827: default:
828: comp = PetscTupIntCompRevlex_N;
829: break;
830: }
831: qsort(sorter, m, (nodeIdxDim + 2) * sizeof(PetscInt), comp);
832: PetscMalloc1(m, &prm);
833: for (i = 0; i < m; i++) prm[i] = sorter[(nodeIdxDim + 2) * i + 1];
834: ni->perm = prm;
835: PetscFree(sorter);
836: }
837: *perm = ni->perm;
838: return 0;
839: }
841: static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
842: {
843: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
845: if (lag->symperms) {
846: PetscInt **selfSyms = lag->symperms[0];
848: if (selfSyms) {
849: PetscInt i, **allocated = &selfSyms[-lag->selfSymOff];
851: for (i = 0; i < lag->numSelfSym; i++) PetscFree(allocated[i]);
852: PetscFree(allocated);
853: }
854: PetscFree(lag->symperms);
855: }
856: if (lag->symflips) {
857: PetscScalar **selfSyms = lag->symflips[0];
859: if (selfSyms) {
860: PetscInt i;
861: PetscScalar **allocated = &selfSyms[-lag->selfSymOff];
863: for (i = 0; i < lag->numSelfSym; i++) PetscFree(allocated[i]);
864: PetscFree(allocated);
865: }
866: PetscFree(lag->symflips);
867: }
868: Petsc1DNodeFamilyDestroy(&(lag->nodeFamily));
869: PetscLagNodeIndicesDestroy(&(lag->vertIndices));
870: PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));
871: PetscLagNodeIndicesDestroy(&(lag->allNodeIndices));
872: PetscFree(lag);
873: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
874: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
875: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);
876: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);
877: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL);
878: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL);
879: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL);
880: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL);
881: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", NULL);
882: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", NULL);
883: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", NULL);
884: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", NULL);
885: return 0;
886: }
888: static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
889: {
890: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
892: PetscViewerASCIIPrintf(viewer, "%s %s%sLagrange dual space\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "", lag->trimmed ? "trimmed " : "");
893: return 0;
894: }
896: static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
897: {
898: PetscBool iascii;
902: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
903: if (iascii) PetscDualSpaceLagrangeView_Ascii(sp, viewer);
904: return 0;
905: }
907: static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscDualSpace sp, PetscOptionItems *PetscOptionsObject)
908: {
909: PetscBool continuous, tensor, trimmed, flg, flg2, flg3;
910: PetscDTNodeType nodeType;
911: PetscReal nodeExponent;
912: PetscInt momentOrder;
913: PetscBool nodeEndpoints, useMoments;
915: PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
916: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
917: PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
918: PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &nodeEndpoints, &nodeExponent);
919: if (nodeType == PETSCDTNODES_DEFAULT) nodeType = PETSCDTNODES_GAUSSJACOBI;
920: PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments);
921: PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder);
922: PetscOptionsHeadBegin(PetscOptionsObject, "PetscDualSpace Lagrange Options");
923: PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
924: if (flg) PetscDualSpaceLagrangeSetContinuity(sp, continuous);
925: PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetTensor", tensor, &tensor, &flg);
926: if (flg) PetscDualSpaceLagrangeSetTensor(sp, tensor);
927: PetscOptionsBool("-petscdualspace_lagrange_trimmed", "Flag for trimmed dual space", "PetscDualSpaceLagrangeSetTrimmed", trimmed, &trimmed, &flg);
928: if (flg) PetscDualSpaceLagrangeSetTrimmed(sp, trimmed);
929: PetscOptionsEnum("-petscdualspace_lagrange_node_type", "Lagrange node location type", "PetscDualSpaceLagrangeSetNodeType", PetscDTNodeTypes, (PetscEnum)nodeType, (PetscEnum *)&nodeType, &flg);
930: PetscOptionsBool("-petscdualspace_lagrange_node_endpoints", "Flag for nodes that include endpoints", "PetscDualSpaceLagrangeSetNodeType", nodeEndpoints, &nodeEndpoints, &flg2);
931: flg3 = PETSC_FALSE;
932: if (nodeType == PETSCDTNODES_GAUSSJACOBI) PetscOptionsReal("-petscdualspace_lagrange_node_exponent", "Gauss-Jacobi weight function exponent", "PetscDualSpaceLagrangeSetNodeType", nodeExponent, &nodeExponent, &flg3);
933: if (flg || flg2 || flg3) PetscDualSpaceLagrangeSetNodeType(sp, nodeType, nodeEndpoints, nodeExponent);
934: PetscOptionsBool("-petscdualspace_lagrange_use_moments", "Use moments (where appropriate) for functionals", "PetscDualSpaceLagrangeSetUseMoments", useMoments, &useMoments, &flg);
935: if (flg) PetscDualSpaceLagrangeSetUseMoments(sp, useMoments);
936: PetscOptionsInt("-petscdualspace_lagrange_moment_order", "Quadrature order for moment functionals", "PetscDualSpaceLagrangeSetMomentOrder", momentOrder, &momentOrder, &flg);
937: if (flg) PetscDualSpaceLagrangeSetMomentOrder(sp, momentOrder);
938: PetscOptionsHeadEnd();
939: return 0;
940: }
942: static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace spNew)
943: {
944: PetscBool cont, tensor, trimmed, boundary;
945: PetscDTNodeType nodeType;
946: PetscReal exponent;
947: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
949: PetscDualSpaceLagrangeGetContinuity(sp, &cont);
950: PetscDualSpaceLagrangeSetContinuity(spNew, cont);
951: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
952: PetscDualSpaceLagrangeSetTensor(spNew, tensor);
953: PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
954: PetscDualSpaceLagrangeSetTrimmed(spNew, trimmed);
955: PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &boundary, &exponent);
956: PetscDualSpaceLagrangeSetNodeType(spNew, nodeType, boundary, exponent);
957: if (lag->nodeFamily) {
958: PetscDualSpace_Lag *lagnew = (PetscDualSpace_Lag *)spNew->data;
960: Petsc1DNodeFamilyReference(lag->nodeFamily);
961: lagnew->nodeFamily = lag->nodeFamily;
962: }
963: return 0;
964: }
966: /* for making tensor product spaces: take a dual space and product a segment space that has all the same
967: * specifications (trimmed, continuous, order, node set), except for the form degree */
968: static PetscErrorCode PetscDualSpaceCreateEdgeSubspace_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt k, PetscInt Nc, PetscBool interiorOnly, PetscDualSpace *bdsp)
969: {
970: DM K;
971: PetscDualSpace_Lag *newlag;
973: PetscDualSpaceDuplicate(sp, bdsp);
974: PetscDualSpaceSetFormDegree(*bdsp, k);
975: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(1, PETSC_TRUE), &K);
976: PetscDualSpaceSetDM(*bdsp, K);
977: DMDestroy(&K);
978: PetscDualSpaceSetOrder(*bdsp, order);
979: PetscDualSpaceSetNumComponents(*bdsp, Nc);
980: newlag = (PetscDualSpace_Lag *)(*bdsp)->data;
981: newlag->interiorOnly = interiorOnly;
982: PetscDualSpaceSetUp(*bdsp);
983: return 0;
984: }
986: /* just the points, weights aren't handled */
987: static PetscErrorCode PetscQuadratureCreateTensor(PetscQuadrature trace, PetscQuadrature fiber, PetscQuadrature *product)
988: {
989: PetscInt dimTrace, dimFiber;
990: PetscInt numPointsTrace, numPointsFiber;
991: PetscInt dim, numPoints;
992: const PetscReal *pointsTrace;
993: const PetscReal *pointsFiber;
994: PetscReal *points;
995: PetscInt i, j, k, p;
997: PetscQuadratureGetData(trace, &dimTrace, NULL, &numPointsTrace, &pointsTrace, NULL);
998: PetscQuadratureGetData(fiber, &dimFiber, NULL, &numPointsFiber, &pointsFiber, NULL);
999: dim = dimTrace + dimFiber;
1000: numPoints = numPointsFiber * numPointsTrace;
1001: PetscMalloc1(numPoints * dim, &points);
1002: for (p = 0, j = 0; j < numPointsFiber; j++) {
1003: for (i = 0; i < numPointsTrace; i++, p++) {
1004: for (k = 0; k < dimTrace; k++) points[p * dim + k] = pointsTrace[i * dimTrace + k];
1005: for (k = 0; k < dimFiber; k++) points[p * dim + dimTrace + k] = pointsFiber[j * dimFiber + k];
1006: }
1007: }
1008: PetscQuadratureCreate(PETSC_COMM_SELF, product);
1009: PetscQuadratureSetData(*product, dim, 0, numPoints, points, NULL);
1010: return 0;
1011: }
1013: /* Kronecker tensor product where matrix is considered a matrix of k-forms, so that
1014: * the entries in the product matrix are wedge products of the entries in the original matrices */
1015: static PetscErrorCode MatTensorAltV(Mat trace, Mat fiber, PetscInt dimTrace, PetscInt kTrace, PetscInt dimFiber, PetscInt kFiber, Mat *product)
1016: {
1017: PetscInt mTrace, nTrace, mFiber, nFiber, m, n, k, i, j, l;
1018: PetscInt dim, NkTrace, NkFiber, Nk;
1019: PetscInt dT, dF;
1020: PetscInt *nnzTrace, *nnzFiber, *nnz;
1021: PetscInt iT, iF, jT, jF, il, jl;
1022: PetscReal *workT, *workT2, *workF, *workF2, *work, *workstar;
1023: PetscReal *projT, *projF;
1024: PetscReal *projTstar, *projFstar;
1025: PetscReal *wedgeMat;
1026: PetscReal sign;
1027: PetscScalar *workS;
1028: Mat prod;
1029: /* this produces dof groups that look like the identity */
1031: MatGetSize(trace, &mTrace, &nTrace);
1032: PetscDTBinomialInt(dimTrace, PetscAbsInt(kTrace), &NkTrace);
1034: MatGetSize(fiber, &mFiber, &nFiber);
1035: PetscDTBinomialInt(dimFiber, PetscAbsInt(kFiber), &NkFiber);
1037: PetscMalloc2(mTrace, &nnzTrace, mFiber, &nnzFiber);
1038: for (i = 0; i < mTrace; i++) {
1039: MatGetRow(trace, i, &(nnzTrace[i]), NULL, NULL);
1041: }
1042: for (i = 0; i < mFiber; i++) {
1043: MatGetRow(fiber, i, &(nnzFiber[i]), NULL, NULL);
1045: }
1046: dim = dimTrace + dimFiber;
1047: k = kFiber + kTrace;
1048: PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1049: m = mTrace * mFiber;
1050: PetscMalloc1(m, &nnz);
1051: for (l = 0, j = 0; j < mFiber; j++)
1052: for (i = 0; i < mTrace; i++, l++) nnz[l] = (nnzTrace[i] / NkTrace) * (nnzFiber[j] / NkFiber) * Nk;
1053: n = (nTrace / NkTrace) * (nFiber / NkFiber) * Nk;
1054: MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &prod);
1055: PetscFree(nnz);
1056: PetscFree2(nnzTrace, nnzFiber);
1057: /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1058: MatSetOption(prod, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1059: /* compute pullbacks */
1060: PetscDTBinomialInt(dim, PetscAbsInt(kTrace), &dT);
1061: PetscDTBinomialInt(dim, PetscAbsInt(kFiber), &dF);
1062: PetscMalloc4(dimTrace * dim, &projT, dimFiber * dim, &projF, dT * NkTrace, &projTstar, dF * NkFiber, &projFstar);
1063: PetscArrayzero(projT, dimTrace * dim);
1064: for (i = 0; i < dimTrace; i++) projT[i * (dim + 1)] = 1.;
1065: PetscArrayzero(projF, dimFiber * dim);
1066: for (i = 0; i < dimFiber; i++) projF[i * (dim + 1) + dimTrace] = 1.;
1067: PetscDTAltVPullbackMatrix(dim, dimTrace, projT, kTrace, projTstar);
1068: PetscDTAltVPullbackMatrix(dim, dimFiber, projF, kFiber, projFstar);
1069: PetscMalloc5(dT, &workT, dF, &workF, Nk, &work, Nk, &workstar, Nk, &workS);
1070: PetscMalloc2(dT, &workT2, dF, &workF2);
1071: PetscMalloc1(Nk * dT, &wedgeMat);
1072: sign = (PetscAbsInt(kTrace * kFiber) & 1) ? -1. : 1.;
1073: for (i = 0, iF = 0; iF < mFiber; iF++) {
1074: PetscInt ncolsF, nformsF;
1075: const PetscInt *colsF;
1076: const PetscScalar *valsF;
1078: MatGetRow(fiber, iF, &ncolsF, &colsF, &valsF);
1079: nformsF = ncolsF / NkFiber;
1080: for (iT = 0; iT < mTrace; iT++, i++) {
1081: PetscInt ncolsT, nformsT;
1082: const PetscInt *colsT;
1083: const PetscScalar *valsT;
1085: MatGetRow(trace, iT, &ncolsT, &colsT, &valsT);
1086: nformsT = ncolsT / NkTrace;
1087: for (j = 0, jF = 0; jF < nformsF; jF++) {
1088: PetscInt colF = colsF[jF * NkFiber] / NkFiber;
1090: for (il = 0; il < dF; il++) {
1091: PetscReal val = 0.;
1092: for (jl = 0; jl < NkFiber; jl++) val += projFstar[il * NkFiber + jl] * PetscRealPart(valsF[jF * NkFiber + jl]);
1093: workF[il] = val;
1094: }
1095: if (kFiber < 0) {
1096: for (il = 0; il < dF; il++) workF2[il] = workF[il];
1097: PetscDTAltVStar(dim, PetscAbsInt(kFiber), 1, workF2, workF);
1098: }
1099: PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kFiber), PetscAbsInt(kTrace), workF, wedgeMat);
1100: for (jT = 0; jT < nformsT; jT++, j++) {
1101: PetscInt colT = colsT[jT * NkTrace] / NkTrace;
1102: PetscInt col = colF * (nTrace / NkTrace) + colT;
1103: const PetscScalar *vals;
1105: for (il = 0; il < dT; il++) {
1106: PetscReal val = 0.;
1107: for (jl = 0; jl < NkTrace; jl++) val += projTstar[il * NkTrace + jl] * PetscRealPart(valsT[jT * NkTrace + jl]);
1108: workT[il] = val;
1109: }
1110: if (kTrace < 0) {
1111: for (il = 0; il < dT; il++) workT2[il] = workT[il];
1112: PetscDTAltVStar(dim, PetscAbsInt(kTrace), 1, workT2, workT);
1113: }
1115: for (il = 0; il < Nk; il++) {
1116: PetscReal val = 0.;
1117: for (jl = 0; jl < dT; jl++) val += sign * wedgeMat[il * dT + jl] * workT[jl];
1118: work[il] = val;
1119: }
1120: if (k < 0) {
1121: PetscDTAltVStar(dim, PetscAbsInt(k), -1, work, workstar);
1122: #if defined(PETSC_USE_COMPLEX)
1123: for (l = 0; l < Nk; l++) workS[l] = workstar[l];
1124: vals = &workS[0];
1125: #else
1126: vals = &workstar[0];
1127: #endif
1128: } else {
1129: #if defined(PETSC_USE_COMPLEX)
1130: for (l = 0; l < Nk; l++) workS[l] = work[l];
1131: vals = &workS[0];
1132: #else
1133: vals = &work[0];
1134: #endif
1135: }
1136: for (l = 0; l < Nk; l++) MatSetValue(prod, i, col * Nk + l, vals[l], INSERT_VALUES); /* Nk */
1137: } /* jT */
1138: } /* jF */
1139: MatRestoreRow(trace, iT, &ncolsT, &colsT, &valsT);
1140: } /* iT */
1141: MatRestoreRow(fiber, iF, &ncolsF, &colsF, &valsF);
1142: } /* iF */
1143: PetscFree(wedgeMat);
1144: PetscFree4(projT, projF, projTstar, projFstar);
1145: PetscFree2(workT2, workF2);
1146: PetscFree5(workT, workF, work, workstar, workS);
1147: MatAssemblyBegin(prod, MAT_FINAL_ASSEMBLY);
1148: MatAssemblyEnd(prod, MAT_FINAL_ASSEMBLY);
1149: *product = prod;
1150: return 0;
1151: }
1153: /* Union of quadrature points, with an attempt to identify commont points in the two sets */
1154: static PetscErrorCode PetscQuadraturePointsMerge(PetscQuadrature quadA, PetscQuadrature quadB, PetscQuadrature *quadJoint, PetscInt *aToJoint[], PetscInt *bToJoint[])
1155: {
1156: PetscInt dimA, dimB;
1157: PetscInt nA, nB, nJoint, i, j, d;
1158: const PetscReal *pointsA;
1159: const PetscReal *pointsB;
1160: PetscReal *pointsJoint;
1161: PetscInt *aToJ, *bToJ;
1162: PetscQuadrature qJ;
1164: PetscQuadratureGetData(quadA, &dimA, NULL, &nA, &pointsA, NULL);
1165: PetscQuadratureGetData(quadB, &dimB, NULL, &nB, &pointsB, NULL);
1167: nJoint = nA;
1168: PetscMalloc1(nA, &aToJ);
1169: for (i = 0; i < nA; i++) aToJ[i] = i;
1170: PetscMalloc1(nB, &bToJ);
1171: for (i = 0; i < nB; i++) {
1172: for (j = 0; j < nA; j++) {
1173: bToJ[i] = -1;
1174: for (d = 0; d < dimA; d++)
1175: if (PetscAbsReal(pointsB[i * dimA + d] - pointsA[j * dimA + d]) > PETSC_SMALL) break;
1176: if (d == dimA) {
1177: bToJ[i] = j;
1178: break;
1179: }
1180: }
1181: if (bToJ[i] == -1) bToJ[i] = nJoint++;
1182: }
1183: *aToJoint = aToJ;
1184: *bToJoint = bToJ;
1185: PetscMalloc1(nJoint * dimA, &pointsJoint);
1186: PetscArraycpy(pointsJoint, pointsA, nA * dimA);
1187: for (i = 0; i < nB; i++) {
1188: if (bToJ[i] >= nA) {
1189: for (d = 0; d < dimA; d++) pointsJoint[bToJ[i] * dimA + d] = pointsB[i * dimA + d];
1190: }
1191: }
1192: PetscQuadratureCreate(PETSC_COMM_SELF, &qJ);
1193: PetscQuadratureSetData(qJ, dimA, 0, nJoint, pointsJoint, NULL);
1194: *quadJoint = qJ;
1195: return 0;
1196: }
1198: /* Matrices matA and matB are both quadrature -> dof matrices: produce a matrix that is joint quadrature -> union of
1199: * dofs, where the joint quadrature was produced by PetscQuadraturePointsMerge */
1200: static PetscErrorCode MatricesMerge(Mat matA, Mat matB, PetscInt dim, PetscInt k, PetscInt numMerged, const PetscInt aToMerged[], const PetscInt bToMerged[], Mat *matMerged)
1201: {
1202: PetscInt m, n, mA, nA, mB, nB, Nk, i, j, l;
1203: Mat M;
1204: PetscInt *nnz;
1205: PetscInt maxnnz;
1206: PetscInt *work;
1208: PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1209: MatGetSize(matA, &mA, &nA);
1211: MatGetSize(matB, &mB, &nB);
1213: m = mA + mB;
1214: n = numMerged * Nk;
1215: PetscMalloc1(m, &nnz);
1216: maxnnz = 0;
1217: for (i = 0; i < mA; i++) {
1218: MatGetRow(matA, i, &(nnz[i]), NULL, NULL);
1220: maxnnz = PetscMax(maxnnz, nnz[i]);
1221: }
1222: for (i = 0; i < mB; i++) {
1223: MatGetRow(matB, i, &(nnz[i + mA]), NULL, NULL);
1225: maxnnz = PetscMax(maxnnz, nnz[i + mA]);
1226: }
1227: MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &M);
1228: PetscFree(nnz);
1229: /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1230: MatSetOption(M, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1231: PetscMalloc1(maxnnz, &work);
1232: for (i = 0; i < mA; i++) {
1233: const PetscInt *cols;
1234: const PetscScalar *vals;
1235: PetscInt nCols;
1236: MatGetRow(matA, i, &nCols, &cols, &vals);
1237: for (j = 0; j < nCols / Nk; j++) {
1238: PetscInt newCol = aToMerged[cols[j * Nk] / Nk];
1239: for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1240: }
1241: MatSetValuesBlocked(M, 1, &i, nCols, work, vals, INSERT_VALUES);
1242: MatRestoreRow(matA, i, &nCols, &cols, &vals);
1243: }
1244: for (i = 0; i < mB; i++) {
1245: const PetscInt *cols;
1246: const PetscScalar *vals;
1248: PetscInt row = i + mA;
1249: PetscInt nCols;
1250: MatGetRow(matB, i, &nCols, &cols, &vals);
1251: for (j = 0; j < nCols / Nk; j++) {
1252: PetscInt newCol = bToMerged[cols[j * Nk] / Nk];
1253: for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1254: }
1255: MatSetValuesBlocked(M, 1, &row, nCols, work, vals, INSERT_VALUES);
1256: MatRestoreRow(matB, i, &nCols, &cols, &vals);
1257: }
1258: PetscFree(work);
1259: MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
1260: MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
1261: *matMerged = M;
1262: return 0;
1263: }
1265: /* Take a dual space and product a segment space that has all the same specifications (trimmed, continuous, order,
1266: * node set), except for the form degree. For computing boundary dofs and for making tensor product spaces */
1267: static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Lagrange(PetscDualSpace sp, DM K, PetscInt f, PetscInt k, PetscInt Ncopies, PetscBool interiorOnly, PetscDualSpace *bdsp)
1268: {
1269: PetscInt Nknew, Ncnew;
1270: PetscInt dim, pointDim = -1;
1271: PetscInt depth;
1272: DM dm;
1273: PetscDualSpace_Lag *newlag;
1275: PetscDualSpaceGetDM(sp, &dm);
1276: DMGetDimension(dm, &dim);
1277: DMPlexGetDepth(dm, &depth);
1278: PetscDualSpaceDuplicate(sp, bdsp);
1279: PetscDualSpaceSetFormDegree(*bdsp, k);
1280: if (!K) {
1281: if (depth == dim) {
1282: DMPolytopeType ct;
1284: pointDim = dim - 1;
1285: DMPlexGetCellType(dm, f, &ct);
1286: DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K);
1287: } else if (depth == 1) {
1288: pointDim = 0;
1289: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DM_POLYTOPE_POINT, &K);
1290: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
1291: } else {
1292: PetscObjectReference((PetscObject)K);
1293: DMGetDimension(K, &pointDim);
1294: }
1295: PetscDualSpaceSetDM(*bdsp, K);
1296: DMDestroy(&K);
1297: PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew);
1298: Ncnew = Nknew * Ncopies;
1299: PetscDualSpaceSetNumComponents(*bdsp, Ncnew);
1300: newlag = (PetscDualSpace_Lag *)(*bdsp)->data;
1301: newlag->interiorOnly = interiorOnly;
1302: PetscDualSpaceSetUp(*bdsp);
1303: return 0;
1304: }
1306: /* Construct simplex nodes from a nodefamily, add Nk dof vectors of length Nk at each node.
1307: * Return the (quadrature, matrix) form of the dofs and the nodeIndices form as well.
1308: *
1309: * Sometimes we want a set of nodes to be contained in the interior of the element,
1310: * even when the node scheme puts nodes on the boundaries. numNodeSkip tells
1311: * the routine how many "layers" of nodes need to be skipped.
1312: * */
1313: static PetscErrorCode PetscDualSpaceLagrangeCreateSimplexNodeMat(Petsc1DNodeFamily nodeFamily, PetscInt dim, PetscInt sum, PetscInt Nk, PetscInt numNodeSkip, PetscQuadrature *iNodes, Mat *iMat, PetscLagNodeIndices *nodeIndices)
1314: {
1315: PetscReal *extraNodeCoords, *nodeCoords;
1316: PetscInt nNodes, nExtraNodes;
1317: PetscInt i, j, k, extraSum = sum + numNodeSkip * (1 + dim);
1318: PetscQuadrature intNodes;
1319: Mat intMat;
1320: PetscLagNodeIndices ni;
1322: PetscDTBinomialInt(dim + sum, dim, &nNodes);
1323: PetscDTBinomialInt(dim + extraSum, dim, &nExtraNodes);
1325: PetscMalloc1(dim * nExtraNodes, &extraNodeCoords);
1326: PetscNew(&ni);
1327: ni->nodeIdxDim = dim + 1;
1328: ni->nodeVecDim = Nk;
1329: ni->nNodes = nNodes * Nk;
1330: ni->refct = 1;
1331: PetscMalloc1(nNodes * Nk * (dim + 1), &(ni->nodeIdx));
1332: PetscMalloc1(nNodes * Nk * Nk, &(ni->nodeVec));
1333: for (i = 0; i < nNodes; i++)
1334: for (j = 0; j < Nk; j++)
1335: for (k = 0; k < Nk; k++) ni->nodeVec[(i * Nk + j) * Nk + k] = (j == k) ? 1. : 0.;
1336: Petsc1DNodeFamilyComputeSimplexNodes(nodeFamily, dim, extraSum, extraNodeCoords);
1337: if (numNodeSkip) {
1338: PetscInt k;
1339: PetscInt *tup;
1341: PetscMalloc1(dim * nNodes, &nodeCoords);
1342: PetscMalloc1(dim + 1, &tup);
1343: for (k = 0; k < nNodes; k++) {
1344: PetscInt j, c;
1345: PetscInt index;
1347: PetscDTIndexToBary(dim + 1, sum, k, tup);
1348: for (j = 0; j < dim + 1; j++) tup[j] += numNodeSkip;
1349: for (c = 0; c < Nk; c++) {
1350: for (j = 0; j < dim + 1; j++) ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1351: }
1352: PetscDTBaryToIndex(dim + 1, extraSum, tup, &index);
1353: for (j = 0; j < dim; j++) nodeCoords[k * dim + j] = extraNodeCoords[index * dim + j];
1354: }
1355: PetscFree(tup);
1356: PetscFree(extraNodeCoords);
1357: } else {
1358: PetscInt k;
1359: PetscInt *tup;
1361: nodeCoords = extraNodeCoords;
1362: PetscMalloc1(dim + 1, &tup);
1363: for (k = 0; k < nNodes; k++) {
1364: PetscInt j, c;
1366: PetscDTIndexToBary(dim + 1, sum, k, tup);
1367: for (c = 0; c < Nk; c++) {
1368: for (j = 0; j < dim + 1; j++) {
1369: /* barycentric indices can have zeros, but we don't want to push forward zeros because it makes it harder to
1370: * determine which nodes correspond to which under symmetries, so we increase by 1. This is fine
1371: * because the nodeIdx coordinates don't have any meaning other than helping to identify symmetries */
1372: ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1373: }
1374: }
1375: }
1376: PetscFree(tup);
1377: }
1378: PetscQuadratureCreate(PETSC_COMM_SELF, &intNodes);
1379: PetscQuadratureSetData(intNodes, dim, 0, nNodes, nodeCoords, NULL);
1380: MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes * Nk, nNodes * Nk, Nk, NULL, &intMat);
1381: MatSetOption(intMat, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1382: for (j = 0; j < nNodes * Nk; j++) {
1383: PetscInt rem = j % Nk;
1384: PetscInt a, aprev = j - rem;
1385: PetscInt anext = aprev + Nk;
1387: for (a = aprev; a < anext; a++) MatSetValue(intMat, j, a, (a == j) ? 1. : 0., INSERT_VALUES);
1388: }
1389: MatAssemblyBegin(intMat, MAT_FINAL_ASSEMBLY);
1390: MatAssemblyEnd(intMat, MAT_FINAL_ASSEMBLY);
1391: *iNodes = intNodes;
1392: *iMat = intMat;
1393: *nodeIndices = ni;
1394: return 0;
1395: }
1397: /* once the nodeIndices have been created for the interior of the reference cell, and for all of the boundary cells,
1398: * push forward the boundary dofs and concatenate them into the full node indices for the dual space */
1399: static PetscErrorCode PetscDualSpaceLagrangeCreateAllNodeIdx(PetscDualSpace sp)
1400: {
1401: DM dm;
1402: PetscInt dim, nDofs;
1403: PetscSection section;
1404: PetscInt pStart, pEnd, p;
1405: PetscInt formDegree, Nk;
1406: PetscInt nodeIdxDim, spintdim;
1407: PetscDualSpace_Lag *lag;
1408: PetscLagNodeIndices ni, verti;
1410: lag = (PetscDualSpace_Lag *)sp->data;
1411: verti = lag->vertIndices;
1412: PetscDualSpaceGetDM(sp, &dm);
1413: DMGetDimension(dm, &dim);
1414: PetscDualSpaceGetFormDegree(sp, &formDegree);
1415: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
1416: PetscDualSpaceGetSection(sp, §ion);
1417: PetscSectionGetStorageSize(section, &nDofs);
1418: PetscNew(&ni);
1419: ni->nodeIdxDim = nodeIdxDim = verti->nodeIdxDim;
1420: ni->nodeVecDim = Nk;
1421: ni->nNodes = nDofs;
1422: ni->refct = 1;
1423: PetscMalloc1(nodeIdxDim * nDofs, &(ni->nodeIdx));
1424: PetscMalloc1(Nk * nDofs, &(ni->nodeVec));
1425: DMPlexGetChart(dm, &pStart, &pEnd);
1426: PetscSectionGetDof(section, 0, &spintdim);
1427: if (spintdim) {
1428: PetscArraycpy(ni->nodeIdx, lag->intNodeIndices->nodeIdx, spintdim * nodeIdxDim);
1429: PetscArraycpy(ni->nodeVec, lag->intNodeIndices->nodeVec, spintdim * Nk);
1430: }
1431: for (p = pStart + 1; p < pEnd; p++) {
1432: PetscDualSpace psp = sp->pointSpaces[p];
1433: PetscDualSpace_Lag *plag;
1434: PetscInt dof, off;
1436: PetscSectionGetDof(section, p, &dof);
1437: if (!dof) continue;
1438: plag = (PetscDualSpace_Lag *)psp->data;
1439: PetscSectionGetOffset(section, p, &off);
1440: PetscLagNodeIndicesPushForward(dm, verti, p, plag->vertIndices, plag->intNodeIndices, 0, formDegree, &(ni->nodeIdx[off * nodeIdxDim]), &(ni->nodeVec[off * Nk]));
1441: }
1442: lag->allNodeIndices = ni;
1443: return 0;
1444: }
1446: /* once the (quadrature, Matrix) forms of the dofs have been created for the interior of the
1447: * reference cell and for the boundary cells, jk
1448: * push forward the boundary data and concatenate them into the full (quadrature, matrix) data
1449: * for the dual space */
1450: static PetscErrorCode PetscDualSpaceCreateAllDataFromInteriorData(PetscDualSpace sp)
1451: {
1452: DM dm;
1453: PetscSection section;
1454: PetscInt pStart, pEnd, p, k, Nk, dim, Nc;
1455: PetscInt nNodes;
1456: PetscInt countNodes;
1457: Mat allMat;
1458: PetscQuadrature allNodes;
1459: PetscInt nDofs;
1460: PetscInt maxNzforms, j;
1461: PetscScalar *work;
1462: PetscReal *L, *J, *Jinv, *v0, *pv0;
1463: PetscInt *iwork;
1464: PetscReal *nodes;
1466: PetscDualSpaceGetDM(sp, &dm);
1467: DMGetDimension(dm, &dim);
1468: PetscDualSpaceGetSection(sp, §ion);
1469: PetscSectionGetStorageSize(section, &nDofs);
1470: DMPlexGetChart(dm, &pStart, &pEnd);
1471: PetscDualSpaceGetFormDegree(sp, &k);
1472: PetscDualSpaceGetNumComponents(sp, &Nc);
1473: PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1474: for (p = pStart, nNodes = 0, maxNzforms = 0; p < pEnd; p++) {
1475: PetscDualSpace psp;
1476: DM pdm;
1477: PetscInt pdim, pNk;
1478: PetscQuadrature intNodes;
1479: Mat intMat;
1481: PetscDualSpaceGetPointSubspace(sp, p, &psp);
1482: if (!psp) continue;
1483: PetscDualSpaceGetDM(psp, &pdm);
1484: DMGetDimension(pdm, &pdim);
1485: if (pdim < PetscAbsInt(k)) continue;
1486: PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
1487: PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
1488: if (intNodes) {
1489: PetscInt nNodesp;
1491: PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, NULL, NULL);
1492: nNodes += nNodesp;
1493: }
1494: if (intMat) {
1495: PetscInt maxNzsp;
1496: PetscInt maxNzformsp;
1498: MatSeqAIJGetMaxRowNonzeros(intMat, &maxNzsp);
1500: maxNzformsp = maxNzsp / pNk;
1501: maxNzforms = PetscMax(maxNzforms, maxNzformsp);
1502: }
1503: }
1504: MatCreateSeqAIJ(PETSC_COMM_SELF, nDofs, nNodes * Nc, maxNzforms * Nk, NULL, &allMat);
1505: MatSetOption(allMat, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1506: PetscMalloc7(dim, &v0, dim, &pv0, dim * dim, &J, dim * dim, &Jinv, Nk * Nk, &L, maxNzforms * Nk, &work, maxNzforms * Nk, &iwork);
1507: for (j = 0; j < dim; j++) pv0[j] = -1.;
1508: PetscMalloc1(dim * nNodes, &nodes);
1509: for (p = pStart, countNodes = 0; p < pEnd; p++) {
1510: PetscDualSpace psp;
1511: PetscQuadrature intNodes;
1512: DM pdm;
1513: PetscInt pdim, pNk;
1514: PetscInt countNodesIn = countNodes;
1515: PetscReal detJ;
1516: Mat intMat;
1518: PetscDualSpaceGetPointSubspace(sp, p, &psp);
1519: if (!psp) continue;
1520: PetscDualSpaceGetDM(psp, &pdm);
1521: DMGetDimension(pdm, &pdim);
1522: if (pdim < PetscAbsInt(k)) continue;
1523: PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
1524: if (intNodes == NULL && intMat == NULL) continue;
1525: PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
1526: if (p) {
1527: DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, Jinv, &detJ);
1528: } else { /* identity */
1529: PetscInt i, j;
1531: for (i = 0; i < dim; i++)
1532: for (j = 0; j < dim; j++) J[i * dim + j] = Jinv[i * dim + j] = 0.;
1533: for (i = 0; i < dim; i++) J[i * dim + i] = Jinv[i * dim + i] = 1.;
1534: for (i = 0; i < dim; i++) v0[i] = -1.;
1535: }
1536: if (pdim != dim) { /* compactify Jacobian */
1537: PetscInt i, j;
1539: for (i = 0; i < dim; i++)
1540: for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
1541: }
1542: PetscDTAltVPullbackMatrix(pdim, dim, J, k, L);
1543: if (intNodes) { /* push forward quadrature locations by the affine transformation */
1544: PetscInt nNodesp;
1545: const PetscReal *nodesp;
1546: PetscInt j;
1548: PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, &nodesp, NULL);
1549: for (j = 0; j < nNodesp; j++, countNodes++) {
1550: PetscInt d, e;
1552: for (d = 0; d < dim; d++) {
1553: nodes[countNodes * dim + d] = v0[d];
1554: for (e = 0; e < pdim; e++) nodes[countNodes * dim + d] += J[d * pdim + e] * (nodesp[j * pdim + e] - pv0[e]);
1555: }
1556: }
1557: }
1558: if (intMat) {
1559: PetscInt nrows;
1560: PetscInt off;
1562: PetscSectionGetDof(section, p, &nrows);
1563: PetscSectionGetOffset(section, p, &off);
1564: for (j = 0; j < nrows; j++) {
1565: PetscInt ncols;
1566: const PetscInt *cols;
1567: const PetscScalar *vals;
1568: PetscInt l, d, e;
1569: PetscInt row = j + off;
1571: MatGetRow(intMat, j, &ncols, &cols, &vals);
1573: for (l = 0; l < ncols / pNk; l++) {
1574: PetscInt blockcol;
1577: blockcol = cols[l * pNk] / pNk;
1578: for (d = 0; d < Nk; d++) iwork[l * Nk + d] = (blockcol + countNodesIn) * Nk + d;
1579: for (d = 0; d < Nk; d++) work[l * Nk + d] = 0.;
1580: for (d = 0; d < Nk; d++) {
1581: for (e = 0; e < pNk; e++) {
1582: /* "push forward" dof by pulling back a k-form to be evaluated on the point: multiply on the right by L */
1583: work[l * Nk + d] += vals[l * pNk + e] * L[e * Nk + d];
1584: }
1585: }
1586: }
1587: MatSetValues(allMat, 1, &row, (ncols / pNk) * Nk, iwork, work, INSERT_VALUES);
1588: MatRestoreRow(intMat, j, &ncols, &cols, &vals);
1589: }
1590: }
1591: }
1592: MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);
1593: MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);
1594: PetscQuadratureCreate(PETSC_COMM_SELF, &allNodes);
1595: PetscQuadratureSetData(allNodes, dim, 0, nNodes, nodes, NULL);
1596: PetscFree7(v0, pv0, J, Jinv, L, work, iwork);
1597: MatDestroy(&(sp->allMat));
1598: sp->allMat = allMat;
1599: PetscQuadratureDestroy(&(sp->allNodes));
1600: sp->allNodes = allNodes;
1601: return 0;
1602: }
1604: /* rather than trying to get all data from the functionals, we create
1605: * the functionals from rows of the quadrature -> dof matrix.
1606: *
1607: * Ideally most of the uses of PetscDualSpace in PetscFE will switch
1608: * to using intMat and allMat, so that the individual functionals
1609: * don't need to be constructed at all */
1610: static PetscErrorCode PetscDualSpaceComputeFunctionalsFromAllData(PetscDualSpace sp)
1611: {
1612: PetscQuadrature allNodes;
1613: Mat allMat;
1614: PetscInt nDofs;
1615: PetscInt dim, k, Nk, Nc, f;
1616: DM dm;
1617: PetscInt nNodes, spdim;
1618: const PetscReal *nodes = NULL;
1619: PetscSection section;
1620: PetscBool useMoments;
1622: PetscDualSpaceGetDM(sp, &dm);
1623: DMGetDimension(dm, &dim);
1624: PetscDualSpaceGetNumComponents(sp, &Nc);
1625: PetscDualSpaceGetFormDegree(sp, &k);
1626: PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1627: PetscDualSpaceGetAllData(sp, &allNodes, &allMat);
1628: nNodes = 0;
1629: if (allNodes) PetscQuadratureGetData(allNodes, NULL, NULL, &nNodes, &nodes, NULL);
1630: MatGetSize(allMat, &nDofs, NULL);
1631: PetscDualSpaceGetSection(sp, §ion);
1632: PetscSectionGetStorageSize(section, &spdim);
1634: PetscMalloc1(nDofs, &(sp->functional));
1635: PetscDualSpaceLagrangeGetUseMoments(sp, &useMoments);
1636: if (useMoments) {
1637: Mat allMat;
1638: PetscInt momentOrder, i;
1639: PetscBool tensor;
1640: const PetscReal *weights;
1641: PetscScalar *array;
1644: PetscDualSpaceLagrangeGetMomentOrder(sp, &momentOrder);
1645: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
1646: if (!tensor) PetscDTStroudConicalQuadrature(dim, Nc, PetscMax(momentOrder + 1, 1), -1.0, 1.0, &(sp->functional[0]));
1647: else PetscDTGaussTensorQuadrature(dim, Nc, PetscMax(momentOrder + 1, 1), -1.0, 1.0, &(sp->functional[0]));
1648: /* Need to replace allNodes and allMat */
1649: PetscObjectReference((PetscObject)sp->functional[0]);
1650: PetscQuadratureDestroy(&(sp->allNodes));
1651: sp->allNodes = sp->functional[0];
1652: PetscQuadratureGetData(sp->allNodes, NULL, NULL, &nNodes, NULL, &weights);
1653: MatCreateSeqDense(PETSC_COMM_SELF, nDofs, nNodes * Nc, NULL, &allMat);
1654: MatDenseGetArrayWrite(allMat, &array);
1655: for (i = 0; i < nNodes * Nc; ++i) array[i] = weights[i];
1656: MatDenseRestoreArrayWrite(allMat, &array);
1657: MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);
1658: MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);
1659: MatDestroy(&(sp->allMat));
1660: sp->allMat = allMat;
1661: return 0;
1662: }
1663: for (f = 0; f < nDofs; f++) {
1664: PetscInt ncols, c;
1665: const PetscInt *cols;
1666: const PetscScalar *vals;
1667: PetscReal *nodesf;
1668: PetscReal *weightsf;
1669: PetscInt nNodesf;
1670: PetscInt countNodes;
1672: MatGetRow(allMat, f, &ncols, &cols, &vals);
1674: for (c = 1, nNodesf = 1; c < ncols; c++) {
1675: if ((cols[c] / Nc) != (cols[c - 1] / Nc)) nNodesf++;
1676: }
1677: PetscMalloc1(dim * nNodesf, &nodesf);
1678: PetscMalloc1(Nc * nNodesf, &weightsf);
1679: for (c = 0, countNodes = 0; c < ncols; c++) {
1680: if (!c || ((cols[c] / Nc) != (cols[c - 1] / Nc))) {
1681: PetscInt d;
1683: for (d = 0; d < Nc; d++) weightsf[countNodes * Nc + d] = 0.;
1684: for (d = 0; d < dim; d++) nodesf[countNodes * dim + d] = nodes[(cols[c] / Nc) * dim + d];
1685: countNodes++;
1686: }
1687: weightsf[(countNodes - 1) * Nc + (cols[c] % Nc)] = PetscRealPart(vals[c]);
1688: }
1689: PetscQuadratureCreate(PETSC_COMM_SELF, &(sp->functional[f]));
1690: PetscQuadratureSetData(sp->functional[f], dim, Nc, nNodesf, nodesf, weightsf);
1691: MatRestoreRow(allMat, f, &ncols, &cols, &vals);
1692: }
1693: return 0;
1694: }
1696: /* take a matrix meant for k-forms and expand it to one for Ncopies */
1697: static PetscErrorCode PetscDualSpaceLagrangeMatrixCreateCopies(Mat A, PetscInt Nk, PetscInt Ncopies, Mat *Abs)
1698: {
1699: PetscInt m, n, i, j, k;
1700: PetscInt maxnnz, *nnz, *iwork;
1701: Mat Ac;
1703: MatGetSize(A, &m, &n);
1705: PetscMalloc1(m * Ncopies, &nnz);
1706: for (i = 0, maxnnz = 0; i < m; i++) {
1707: PetscInt innz;
1708: MatGetRow(A, i, &innz, NULL, NULL);
1710: for (j = 0; j < Ncopies; j++) nnz[i * Ncopies + j] = innz;
1711: maxnnz = PetscMax(maxnnz, innz);
1712: }
1713: MatCreateSeqAIJ(PETSC_COMM_SELF, m * Ncopies, n * Ncopies, 0, nnz, &Ac);
1714: MatSetOption(Ac, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1715: PetscFree(nnz);
1716: PetscMalloc1(maxnnz, &iwork);
1717: for (i = 0; i < m; i++) {
1718: PetscInt innz;
1719: const PetscInt *cols;
1720: const PetscScalar *vals;
1722: MatGetRow(A, i, &innz, &cols, &vals);
1723: for (j = 0; j < innz; j++) iwork[j] = (cols[j] / Nk) * (Nk * Ncopies) + (cols[j] % Nk);
1724: for (j = 0; j < Ncopies; j++) {
1725: PetscInt row = i * Ncopies + j;
1727: MatSetValues(Ac, 1, &row, innz, iwork, vals, INSERT_VALUES);
1728: for (k = 0; k < innz; k++) iwork[k] += Nk;
1729: }
1730: MatRestoreRow(A, i, &innz, &cols, &vals);
1731: }
1732: PetscFree(iwork);
1733: MatAssemblyBegin(Ac, MAT_FINAL_ASSEMBLY);
1734: MatAssemblyEnd(Ac, MAT_FINAL_ASSEMBLY);
1735: *Abs = Ac;
1736: return 0;
1737: }
1739: /* check if a cell is a tensor product of the segment with a facet,
1740: * specifically checking if f and f2 can be the "endpoints" (like the triangles
1741: * at either end of a wedge) */
1742: static PetscErrorCode DMPlexPointIsTensor_Internal_Given(DM dm, PetscInt p, PetscInt f, PetscInt f2, PetscBool *isTensor)
1743: {
1744: PetscInt coneSize, c;
1745: const PetscInt *cone;
1746: const PetscInt *fCone;
1747: const PetscInt *f2Cone;
1748: PetscInt fs[2];
1749: PetscInt meetSize, nmeet;
1750: const PetscInt *meet;
1752: fs[0] = f;
1753: fs[1] = f2;
1754: DMPlexGetMeet(dm, 2, fs, &meetSize, &meet);
1755: nmeet = meetSize;
1756: DMPlexRestoreMeet(dm, 2, fs, &meetSize, &meet);
1757: /* two points that have a non-empty meet cannot be at opposite ends of a cell */
1758: if (nmeet) {
1759: *isTensor = PETSC_FALSE;
1760: return 0;
1761: }
1762: DMPlexGetConeSize(dm, p, &coneSize);
1763: DMPlexGetCone(dm, p, &cone);
1764: DMPlexGetCone(dm, f, &fCone);
1765: DMPlexGetCone(dm, f2, &f2Cone);
1766: for (c = 0; c < coneSize; c++) {
1767: PetscInt e, ef;
1768: PetscInt d = -1, d2 = -1;
1769: PetscInt dcount, d2count;
1770: PetscInt t = cone[c];
1771: PetscInt tConeSize;
1772: PetscBool tIsTensor;
1773: const PetscInt *tCone;
1775: if (t == f || t == f2) continue;
1776: /* for every other facet in the cone, check that is has
1777: * one ridge in common with each end */
1778: DMPlexGetConeSize(dm, t, &tConeSize);
1779: DMPlexGetCone(dm, t, &tCone);
1781: dcount = 0;
1782: d2count = 0;
1783: for (e = 0; e < tConeSize; e++) {
1784: PetscInt q = tCone[e];
1785: for (ef = 0; ef < coneSize - 2; ef++) {
1786: if (fCone[ef] == q) {
1787: if (dcount) {
1788: *isTensor = PETSC_FALSE;
1789: return 0;
1790: }
1791: d = q;
1792: dcount++;
1793: } else if (f2Cone[ef] == q) {
1794: if (d2count) {
1795: *isTensor = PETSC_FALSE;
1796: return 0;
1797: }
1798: d2 = q;
1799: d2count++;
1800: }
1801: }
1802: }
1803: /* if the whole cell is a tensor with the segment, then this
1804: * facet should be a tensor with the segment */
1805: DMPlexPointIsTensor_Internal_Given(dm, t, d, d2, &tIsTensor);
1806: if (!tIsTensor) {
1807: *isTensor = PETSC_FALSE;
1808: return 0;
1809: }
1810: }
1811: *isTensor = PETSC_TRUE;
1812: return 0;
1813: }
1815: /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair
1816: * that could be the opposite ends */
1817: static PetscErrorCode DMPlexPointIsTensor_Internal(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB)
1818: {
1819: PetscInt coneSize, c, c2;
1820: const PetscInt *cone;
1822: DMPlexGetConeSize(dm, p, &coneSize);
1823: if (!coneSize) {
1824: if (isTensor) *isTensor = PETSC_FALSE;
1825: if (endA) *endA = -1;
1826: if (endB) *endB = -1;
1827: }
1828: DMPlexGetCone(dm, p, &cone);
1829: for (c = 0; c < coneSize; c++) {
1830: PetscInt f = cone[c];
1831: PetscInt fConeSize;
1833: DMPlexGetConeSize(dm, f, &fConeSize);
1834: if (fConeSize != coneSize - 2) continue;
1836: for (c2 = c + 1; c2 < coneSize; c2++) {
1837: PetscInt f2 = cone[c2];
1838: PetscBool isTensorff2;
1839: PetscInt f2ConeSize;
1841: DMPlexGetConeSize(dm, f2, &f2ConeSize);
1842: if (f2ConeSize != coneSize - 2) continue;
1844: DMPlexPointIsTensor_Internal_Given(dm, p, f, f2, &isTensorff2);
1845: if (isTensorff2) {
1846: if (isTensor) *isTensor = PETSC_TRUE;
1847: if (endA) *endA = f;
1848: if (endB) *endB = f2;
1849: return 0;
1850: }
1851: }
1852: }
1853: if (isTensor) *isTensor = PETSC_FALSE;
1854: if (endA) *endA = -1;
1855: if (endB) *endB = -1;
1856: return 0;
1857: }
1859: /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair
1860: * that could be the opposite ends */
1861: static PetscErrorCode DMPlexPointIsTensor(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB)
1862: {
1863: DMPlexInterpolatedFlag interpolated;
1865: DMPlexIsInterpolated(dm, &interpolated);
1867: DMPlexPointIsTensor_Internal(dm, p, isTensor, endA, endB);
1868: return 0;
1869: }
1871: /* Let k = formDegree and k' = -sign(k) * dim + k. Transform a symmetric frame for k-forms on the biunit simplex into
1872: * a symmetric frame for k'-forms on the biunit simplex.
1873: *
1874: * A frame is "symmetric" if the pullback of every symmetry of the biunit simplex is a permutation of the frame.
1875: *
1876: * forms in the symmetric frame are used as dofs in the untrimmed simplex spaces. This way, symmetries of the
1877: * reference cell result in permutations of dofs grouped by node.
1878: *
1879: * Use T to transform dof matrices for k'-forms into dof matrices for k-forms as a block diagonal transformation on
1880: * the right.
1881: */
1882: static PetscErrorCode BiunitSimplexSymmetricFormTransformation(PetscInt dim, PetscInt formDegree, PetscReal T[])
1883: {
1884: PetscInt k = formDegree;
1885: PetscInt kd = k < 0 ? dim + k : k - dim;
1886: PetscInt Nk;
1887: PetscReal *biToEq, *eqToBi, *biToEqStar, *eqToBiStar;
1888: PetscInt fact;
1890: PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1891: PetscCalloc4(dim * dim, &biToEq, dim * dim, &eqToBi, Nk * Nk, &biToEqStar, Nk * Nk, &eqToBiStar);
1892: /* fill in biToEq: Jacobian of the transformation from the biunit simplex to the equilateral simplex */
1893: fact = 0;
1894: for (PetscInt i = 0; i < dim; i++) {
1895: biToEq[i * dim + i] = PetscSqrtReal(((PetscReal)i + 2.) / (2. * ((PetscReal)i + 1.)));
1896: fact += 4 * (i + 1);
1897: for (PetscInt j = i + 1; j < dim; j++) biToEq[i * dim + j] = PetscSqrtReal(1. / (PetscReal)fact);
1898: }
1899: /* fill in eqToBi: Jacobian of the transformation from the equilateral simplex to the biunit simplex */
1900: fact = 0;
1901: for (PetscInt j = 0; j < dim; j++) {
1902: eqToBi[j * dim + j] = PetscSqrtReal(2. * ((PetscReal)j + 1.) / ((PetscReal)j + 2));
1903: fact += j + 1;
1904: for (PetscInt i = 0; i < j; i++) eqToBi[i * dim + j] = -PetscSqrtReal(1. / (PetscReal)fact);
1905: }
1906: PetscDTAltVPullbackMatrix(dim, dim, biToEq, kd, biToEqStar);
1907: PetscDTAltVPullbackMatrix(dim, dim, eqToBi, k, eqToBiStar);
1908: /* product of pullbacks simulates the following steps
1909: *
1910: * 1. start with frame W = [w_1, w_2, ..., w_m] of k forms that is symmetric on the biunit simplex:
1911: if J is the Jacobian of a symmetry of the biunit simplex, then J_k* W = [J_k*w_1, ..., J_k*w_m]
1912: is a permutation of W.
1913: Even though a k' form --- a (dim - k) form represented by its Hodge star --- has the same geometric
1914: content as a k form, W is not a symmetric frame of k' forms on the biunit simplex. That's because,
1915: for general Jacobian J, J_k* != J_k'*.
1916: * 2. pullback W to the equilateral triangle using the k pullback, W_eq = eqToBi_k* W. All symmetries of the
1917: equilateral simplex have orthonormal Jacobians. For an orthonormal Jacobian O, J_k* = J_k'*, so W_eq is
1918: also a symmetric frame for k' forms on the equilateral simplex.
1919: 3. pullback W_eq back to the biunit simplex using the k' pulback, V = biToEq_k'* W_eq = biToEq_k'* eqToBi_k* W.
1920: V is a symmetric frame for k' forms on the biunit simplex.
1921: */
1922: for (PetscInt i = 0; i < Nk; i++) {
1923: for (PetscInt j = 0; j < Nk; j++) {
1924: PetscReal val = 0.;
1925: for (PetscInt k = 0; k < Nk; k++) val += biToEqStar[i * Nk + k] * eqToBiStar[k * Nk + j];
1926: T[i * Nk + j] = val;
1927: }
1928: }
1929: PetscFree4(biToEq, eqToBi, biToEqStar, eqToBiStar);
1930: return 0;
1931: }
1933: /* permute a quadrature -> dof matrix so that its rows are in revlex order by nodeIdx */
1934: static PetscErrorCode MatPermuteByNodeIdx(Mat A, PetscLagNodeIndices ni, Mat *Aperm)
1935: {
1936: PetscInt m, n, i, j;
1937: PetscInt nodeIdxDim = ni->nodeIdxDim;
1938: PetscInt nodeVecDim = ni->nodeVecDim;
1939: PetscInt *perm;
1940: IS permIS;
1941: IS id;
1942: PetscInt *nIdxPerm;
1943: PetscReal *nVecPerm;
1945: PetscLagNodeIndicesGetPermutation(ni, &perm);
1946: MatGetSize(A, &m, &n);
1947: PetscMalloc1(nodeIdxDim * m, &nIdxPerm);
1948: PetscMalloc1(nodeVecDim * m, &nVecPerm);
1949: for (i = 0; i < m; i++)
1950: for (j = 0; j < nodeIdxDim; j++) nIdxPerm[i * nodeIdxDim + j] = ni->nodeIdx[perm[i] * nodeIdxDim + j];
1951: for (i = 0; i < m; i++)
1952: for (j = 0; j < nodeVecDim; j++) nVecPerm[i * nodeVecDim + j] = ni->nodeVec[perm[i] * nodeVecDim + j];
1953: ISCreateGeneral(PETSC_COMM_SELF, m, perm, PETSC_USE_POINTER, &permIS);
1954: ISSetPermutation(permIS);
1955: ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &id);
1956: ISSetPermutation(id);
1957: MatPermute(A, permIS, id, Aperm);
1958: ISDestroy(&permIS);
1959: ISDestroy(&id);
1960: for (i = 0; i < m; i++) perm[i] = i;
1961: PetscFree(ni->nodeIdx);
1962: PetscFree(ni->nodeVec);
1963: ni->nodeIdx = nIdxPerm;
1964: ni->nodeVec = nVecPerm;
1965: return 0;
1966: }
1968: static PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1969: {
1970: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
1971: DM dm = sp->dm;
1972: DM dmint = NULL;
1973: PetscInt order;
1974: PetscInt Nc = sp->Nc;
1975: MPI_Comm comm;
1976: PetscBool continuous;
1977: PetscSection section;
1978: PetscInt depth, dim, pStart, pEnd, cStart, cEnd, p, *pStratStart, *pStratEnd, d;
1979: PetscInt formDegree, Nk, Ncopies;
1980: PetscInt tensorf = -1, tensorf2 = -1;
1981: PetscBool tensorCell, tensorSpace;
1982: PetscBool uniform, trimmed;
1983: Petsc1DNodeFamily nodeFamily;
1984: PetscInt numNodeSkip;
1985: DMPlexInterpolatedFlag interpolated;
1986: PetscBool isbdm;
1988: /* step 1: sanitize input */
1989: PetscObjectGetComm((PetscObject)sp, &comm);
1990: DMGetDimension(dm, &dim);
1991: PetscObjectTypeCompare((PetscObject)sp, PETSCDUALSPACEBDM, &isbdm);
1992: if (isbdm) {
1993: sp->k = -(dim - 1); /* form degree of H-div */
1994: PetscObjectChangeTypeName((PetscObject)sp, PETSCDUALSPACELAGRANGE);
1995: }
1996: PetscDualSpaceGetFormDegree(sp, &formDegree);
1998: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
1999: if (sp->Nc <= 0 && lag->numCopies > 0) sp->Nc = Nk * lag->numCopies;
2000: Nc = sp->Nc;
2002: if (lag->numCopies <= 0) lag->numCopies = Nc / Nk;
2003: Ncopies = lag->numCopies;
2005: if (!dim) sp->order = 0;
2006: order = sp->order;
2007: uniform = sp->uniform;
2009: if (lag->trimmed && !formDegree) lag->trimmed = PETSC_FALSE; /* trimmed spaces are the same as full spaces for 0-forms */
2010: if (lag->nodeType == PETSCDTNODES_DEFAULT) {
2011: lag->nodeType = PETSCDTNODES_GAUSSJACOBI;
2012: lag->nodeExponent = 0.;
2013: /* trimmed spaces don't include corner vertices, so don't use end nodes by default */
2014: lag->endNodes = lag->trimmed ? PETSC_FALSE : PETSC_TRUE;
2015: }
2016: /* If a trimmed space and the user did choose nodes with endpoints, skip them by default */
2017: if (lag->numNodeSkip < 0) lag->numNodeSkip = (lag->trimmed && lag->endNodes) ? 1 : 0;
2018: numNodeSkip = lag->numNodeSkip;
2020: if (lag->trimmed && PetscAbsInt(formDegree) == dim) { /* convert trimmed n-forms to untrimmed of one polynomial order less */
2021: sp->order--;
2022: order--;
2023: lag->trimmed = PETSC_FALSE;
2024: }
2025: trimmed = lag->trimmed;
2026: if (!order || PetscAbsInt(formDegree) == dim) lag->continuous = PETSC_FALSE;
2027: continuous = lag->continuous;
2028: DMPlexGetDepth(dm, &depth);
2029: DMPlexGetChart(dm, &pStart, &pEnd);
2030: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2033: DMPlexIsInterpolated(dm, &interpolated);
2034: if (interpolated != DMPLEX_INTERPOLATED_FULL) {
2035: DMPlexInterpolate(dm, &dmint);
2036: } else {
2037: PetscObjectReference((PetscObject)dm);
2038: dmint = dm;
2039: }
2040: tensorCell = PETSC_FALSE;
2041: if (dim > 1) DMPlexPointIsTensor(dmint, 0, &tensorCell, &tensorf, &tensorf2);
2042: lag->tensorCell = tensorCell;
2043: if (dim < 2 || !lag->tensorCell) lag->tensorSpace = PETSC_FALSE;
2044: tensorSpace = lag->tensorSpace;
2045: if (!lag->nodeFamily) Petsc1DNodeFamilyCreate(lag->nodeType, lag->nodeExponent, lag->endNodes, &lag->nodeFamily);
2046: nodeFamily = lag->nodeFamily;
2049: /* step 2: construct the boundary spaces */
2050: PetscMalloc2(depth + 1, &pStratStart, depth + 1, &pStratEnd);
2051: PetscCalloc1(pEnd, &(sp->pointSpaces));
2052: for (d = 0; d <= depth; ++d) DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);
2053: PetscDualSpaceSectionCreate_Internal(sp, §ion);
2054: sp->pointSection = section;
2055: if (continuous && !(lag->interiorOnly)) {
2056: PetscInt h;
2058: for (p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */
2059: PetscReal v0[3];
2060: DMPolytopeType ptype;
2061: PetscReal J[9], detJ;
2062: PetscInt q;
2064: DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, NULL, &detJ);
2065: DMPlexGetCellType(dm, p, &ptype);
2067: /* compare to previous facets: if computed, reference that dualspace */
2068: for (q = pStratStart[depth - 1]; q < p; q++) {
2069: DMPolytopeType qtype;
2071: DMPlexGetCellType(dm, q, &qtype);
2072: if (qtype == ptype) break;
2073: }
2074: if (q < p) { /* this facet has the same dual space as that one */
2075: PetscObjectReference((PetscObject)sp->pointSpaces[q]);
2076: sp->pointSpaces[p] = sp->pointSpaces[q];
2077: continue;
2078: }
2079: /* if not, recursively compute this dual space */
2080: PetscDualSpaceCreateFacetSubspace_Lagrange(sp, NULL, p, formDegree, Ncopies, PETSC_FALSE, &sp->pointSpaces[p]);
2081: }
2082: for (h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */
2083: PetscInt hd = depth - h;
2084: PetscInt hdim = dim - h;
2086: if (hdim < PetscAbsInt(formDegree)) break;
2087: for (p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
2088: PetscInt suppSize, s;
2089: const PetscInt *supp;
2091: DMPlexGetSupportSize(dm, p, &suppSize);
2092: DMPlexGetSupport(dm, p, &supp);
2093: for (s = 0; s < suppSize; s++) {
2094: DM qdm;
2095: PetscDualSpace qsp, psp;
2096: PetscInt c, coneSize, q;
2097: const PetscInt *cone;
2098: const PetscInt *refCone;
2100: q = supp[0];
2101: qsp = sp->pointSpaces[q];
2102: DMPlexGetConeSize(dm, q, &coneSize);
2103: DMPlexGetCone(dm, q, &cone);
2104: for (c = 0; c < coneSize; c++)
2105: if (cone[c] == p) break;
2107: PetscDualSpaceGetDM(qsp, &qdm);
2108: DMPlexGetCone(qdm, 0, &refCone);
2109: /* get the equivalent dual space from the support dual space */
2110: PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp);
2111: if (!s) {
2112: PetscObjectReference((PetscObject)psp);
2113: sp->pointSpaces[p] = psp;
2114: }
2115: }
2116: }
2117: }
2118: for (p = 1; p < pEnd; p++) {
2119: PetscInt pspdim;
2120: if (!sp->pointSpaces[p]) continue;
2121: PetscDualSpaceGetInteriorDimension(sp->pointSpaces[p], &pspdim);
2122: PetscSectionSetDof(section, p, pspdim);
2123: }
2124: }
2126: if (Ncopies > 1) {
2127: Mat intMatScalar, allMatScalar;
2128: PetscDualSpace scalarsp;
2129: PetscDualSpace_Lag *scalarlag;
2131: PetscDualSpaceDuplicate(sp, &scalarsp);
2132: /* Setting the number of components to Nk is a space with 1 copy of each k-form */
2133: PetscDualSpaceSetNumComponents(scalarsp, Nk);
2134: PetscDualSpaceSetUp(scalarsp);
2135: PetscDualSpaceGetInteriorData(scalarsp, &(sp->intNodes), &intMatScalar);
2136: PetscObjectReference((PetscObject)(sp->intNodes));
2137: if (intMatScalar) PetscDualSpaceLagrangeMatrixCreateCopies(intMatScalar, Nk, Ncopies, &(sp->intMat));
2138: PetscDualSpaceGetAllData(scalarsp, &(sp->allNodes), &allMatScalar);
2139: PetscObjectReference((PetscObject)(sp->allNodes));
2140: PetscDualSpaceLagrangeMatrixCreateCopies(allMatScalar, Nk, Ncopies, &(sp->allMat));
2141: sp->spdim = scalarsp->spdim * Ncopies;
2142: sp->spintdim = scalarsp->spintdim * Ncopies;
2143: scalarlag = (PetscDualSpace_Lag *)scalarsp->data;
2144: PetscLagNodeIndicesReference(scalarlag->vertIndices);
2145: lag->vertIndices = scalarlag->vertIndices;
2146: PetscLagNodeIndicesReference(scalarlag->intNodeIndices);
2147: lag->intNodeIndices = scalarlag->intNodeIndices;
2148: PetscLagNodeIndicesReference(scalarlag->allNodeIndices);
2149: lag->allNodeIndices = scalarlag->allNodeIndices;
2150: PetscDualSpaceDestroy(&scalarsp);
2151: PetscSectionSetDof(section, 0, sp->spintdim);
2152: PetscDualSpaceSectionSetUp_Internal(sp, section);
2153: PetscDualSpaceComputeFunctionalsFromAllData(sp);
2154: PetscFree2(pStratStart, pStratEnd);
2155: DMDestroy(&dmint);
2156: return 0;
2157: }
2159: if (trimmed && !continuous) {
2160: /* the dofs of a trimmed space don't have a nice tensor/lattice structure:
2161: * just construct the continuous dual space and copy all of the data over,
2162: * allocating it all to the cell instead of splitting it up between the boundaries */
2163: PetscDualSpace spcont;
2164: PetscInt spdim, f;
2165: PetscQuadrature allNodes;
2166: PetscDualSpace_Lag *lagc;
2167: Mat allMat;
2169: PetscDualSpaceDuplicate(sp, &spcont);
2170: PetscDualSpaceLagrangeSetContinuity(spcont, PETSC_TRUE);
2171: PetscDualSpaceSetUp(spcont);
2172: PetscDualSpaceGetDimension(spcont, &spdim);
2173: sp->spdim = sp->spintdim = spdim;
2174: PetscSectionSetDof(section, 0, spdim);
2175: PetscDualSpaceSectionSetUp_Internal(sp, section);
2176: PetscMalloc1(spdim, &(sp->functional));
2177: for (f = 0; f < spdim; f++) {
2178: PetscQuadrature fn;
2180: PetscDualSpaceGetFunctional(spcont, f, &fn);
2181: PetscObjectReference((PetscObject)fn);
2182: sp->functional[f] = fn;
2183: }
2184: PetscDualSpaceGetAllData(spcont, &allNodes, &allMat);
2185: PetscObjectReference((PetscObject)allNodes);
2186: PetscObjectReference((PetscObject)allNodes);
2187: sp->allNodes = sp->intNodes = allNodes;
2188: PetscObjectReference((PetscObject)allMat);
2189: PetscObjectReference((PetscObject)allMat);
2190: sp->allMat = sp->intMat = allMat;
2191: lagc = (PetscDualSpace_Lag *)spcont->data;
2192: PetscLagNodeIndicesReference(lagc->vertIndices);
2193: lag->vertIndices = lagc->vertIndices;
2194: PetscLagNodeIndicesReference(lagc->allNodeIndices);
2195: PetscLagNodeIndicesReference(lagc->allNodeIndices);
2196: lag->intNodeIndices = lagc->allNodeIndices;
2197: lag->allNodeIndices = lagc->allNodeIndices;
2198: PetscDualSpaceDestroy(&spcont);
2199: PetscFree2(pStratStart, pStratEnd);
2200: DMDestroy(&dmint);
2201: return 0;
2202: }
2204: /* step 3: construct intNodes, and intMat, and combine it with boundray data to make allNodes and allMat */
2205: if (!tensorSpace) {
2206: if (!tensorCell) PetscLagNodeIndicesCreateSimplexVertices(dm, &(lag->vertIndices));
2208: if (trimmed) {
2209: /* there is one dof in the interior of the a trimmed element for each full polynomial of with degree at most
2210: * order + k - dim - 1 */
2211: if (order + PetscAbsInt(formDegree) > dim) {
2212: PetscInt sum = order + PetscAbsInt(formDegree) - dim - 1;
2213: PetscInt nDofs;
2215: PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));
2216: MatGetSize(sp->intMat, &nDofs, NULL);
2217: PetscSectionSetDof(section, 0, nDofs);
2218: }
2219: PetscDualSpaceSectionSetUp_Internal(sp, section);
2220: PetscDualSpaceCreateAllDataFromInteriorData(sp);
2221: PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2222: } else {
2223: if (!continuous) {
2224: /* if discontinuous just construct one node for each set of dofs (a set of dofs is a basis for the k-form
2225: * space) */
2226: PetscInt sum = order;
2227: PetscInt nDofs;
2229: PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));
2230: MatGetSize(sp->intMat, &nDofs, NULL);
2231: PetscSectionSetDof(section, 0, nDofs);
2232: PetscDualSpaceSectionSetUp_Internal(sp, section);
2233: PetscObjectReference((PetscObject)(sp->intNodes));
2234: sp->allNodes = sp->intNodes;
2235: PetscObjectReference((PetscObject)(sp->intMat));
2236: sp->allMat = sp->intMat;
2237: PetscLagNodeIndicesReference(lag->intNodeIndices);
2238: lag->allNodeIndices = lag->intNodeIndices;
2239: } else {
2240: /* there is one dof in the interior of the a full element for each trimmed polynomial of with degree at most
2241: * order + k - dim, but with complementary form degree */
2242: if (order + PetscAbsInt(formDegree) > dim) {
2243: PetscDualSpace trimmedsp;
2244: PetscDualSpace_Lag *trimmedlag;
2245: PetscQuadrature intNodes;
2246: PetscInt trFormDegree = formDegree >= 0 ? formDegree - dim : dim - PetscAbsInt(formDegree);
2247: PetscInt nDofs;
2248: Mat intMat;
2250: PetscDualSpaceDuplicate(sp, &trimmedsp);
2251: PetscDualSpaceLagrangeSetTrimmed(trimmedsp, PETSC_TRUE);
2252: PetscDualSpaceSetOrder(trimmedsp, order + PetscAbsInt(formDegree) - dim);
2253: PetscDualSpaceSetFormDegree(trimmedsp, trFormDegree);
2254: trimmedlag = (PetscDualSpace_Lag *)trimmedsp->data;
2255: trimmedlag->numNodeSkip = numNodeSkip + 1;
2256: PetscDualSpaceSetUp(trimmedsp);
2257: PetscDualSpaceGetAllData(trimmedsp, &intNodes, &intMat);
2258: PetscObjectReference((PetscObject)intNodes);
2259: sp->intNodes = intNodes;
2260: PetscLagNodeIndicesReference(trimmedlag->allNodeIndices);
2261: lag->intNodeIndices = trimmedlag->allNodeIndices;
2262: PetscObjectReference((PetscObject)intMat);
2263: if (PetscAbsInt(formDegree) > 0 && PetscAbsInt(formDegree) < dim) {
2264: PetscReal *T;
2265: PetscScalar *work;
2266: PetscInt nCols, nRows;
2267: Mat intMatT;
2269: MatDuplicate(intMat, MAT_COPY_VALUES, &intMatT);
2270: MatGetSize(intMat, &nRows, &nCols);
2271: PetscMalloc2(Nk * Nk, &T, nCols, &work);
2272: BiunitSimplexSymmetricFormTransformation(dim, formDegree, T);
2273: for (PetscInt row = 0; row < nRows; row++) {
2274: PetscInt nrCols;
2275: const PetscInt *rCols;
2276: const PetscScalar *rVals;
2278: MatGetRow(intMat, row, &nrCols, &rCols, &rVals);
2280: for (PetscInt b = 0; b < nrCols; b += Nk) {
2281: const PetscScalar *v = &rVals[b];
2282: PetscScalar *w = &work[b];
2283: for (PetscInt j = 0; j < Nk; j++) {
2284: w[j] = 0.;
2285: for (PetscInt i = 0; i < Nk; i++) w[j] += v[i] * T[i * Nk + j];
2286: }
2287: }
2288: MatSetValuesBlocked(intMatT, 1, &row, nrCols, rCols, work, INSERT_VALUES);
2289: MatRestoreRow(intMat, row, &nrCols, &rCols, &rVals);
2290: }
2291: MatAssemblyBegin(intMatT, MAT_FINAL_ASSEMBLY);
2292: MatAssemblyEnd(intMatT, MAT_FINAL_ASSEMBLY);
2293: MatDestroy(&intMat);
2294: intMat = intMatT;
2295: PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));
2296: PetscLagNodeIndicesDuplicate(trimmedlag->allNodeIndices, &(lag->intNodeIndices));
2297: {
2298: PetscInt nNodes = lag->intNodeIndices->nNodes;
2299: PetscReal *newNodeVec = lag->intNodeIndices->nodeVec;
2300: const PetscReal *oldNodeVec = trimmedlag->allNodeIndices->nodeVec;
2302: for (PetscInt n = 0; n < nNodes; n++) {
2303: PetscReal *w = &newNodeVec[n * Nk];
2304: const PetscReal *v = &oldNodeVec[n * Nk];
2306: for (PetscInt j = 0; j < Nk; j++) {
2307: w[j] = 0.;
2308: for (PetscInt i = 0; i < Nk; i++) w[j] += v[i] * T[i * Nk + j];
2309: }
2310: }
2311: }
2312: PetscFree2(T, work);
2313: }
2314: sp->intMat = intMat;
2315: MatGetSize(sp->intMat, &nDofs, NULL);
2316: PetscDualSpaceDestroy(&trimmedsp);
2317: PetscSectionSetDof(section, 0, nDofs);
2318: }
2319: PetscDualSpaceSectionSetUp_Internal(sp, section);
2320: PetscDualSpaceCreateAllDataFromInteriorData(sp);
2321: PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2322: }
2323: }
2324: } else {
2325: PetscQuadrature intNodesTrace = NULL;
2326: PetscQuadrature intNodesFiber = NULL;
2327: PetscQuadrature intNodes = NULL;
2328: PetscLagNodeIndices intNodeIndices = NULL;
2329: Mat intMat = NULL;
2331: if (PetscAbsInt(formDegree) < dim) { /* get the trace k-forms on the first facet, and the 0-forms on the edge,
2332: and wedge them together to create some of the k-form dofs */
2333: PetscDualSpace trace, fiber;
2334: PetscDualSpace_Lag *tracel, *fiberl;
2335: Mat intMatTrace, intMatFiber;
2337: if (sp->pointSpaces[tensorf]) {
2338: PetscObjectReference((PetscObject)(sp->pointSpaces[tensorf]));
2339: trace = sp->pointSpaces[tensorf];
2340: } else {
2341: PetscDualSpaceCreateFacetSubspace_Lagrange(sp, NULL, tensorf, formDegree, Ncopies, PETSC_TRUE, &trace);
2342: }
2343: PetscDualSpaceCreateEdgeSubspace_Lagrange(sp, order, 0, 1, PETSC_TRUE, &fiber);
2344: tracel = (PetscDualSpace_Lag *)trace->data;
2345: fiberl = (PetscDualSpace_Lag *)fiber->data;
2346: PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));
2347: PetscDualSpaceGetInteriorData(trace, &intNodesTrace, &intMatTrace);
2348: PetscDualSpaceGetInteriorData(fiber, &intNodesFiber, &intMatFiber);
2349: if (intNodesTrace && intNodesFiber) {
2350: PetscQuadratureCreateTensor(intNodesTrace, intNodesFiber, &intNodes);
2351: MatTensorAltV(intMatTrace, intMatFiber, dim - 1, formDegree, 1, 0, &intMat);
2352: PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, formDegree, fiberl->intNodeIndices, 1, 0, &intNodeIndices);
2353: }
2354: PetscObjectReference((PetscObject)intNodesTrace);
2355: PetscObjectReference((PetscObject)intNodesFiber);
2356: PetscDualSpaceDestroy(&fiber);
2357: PetscDualSpaceDestroy(&trace);
2358: }
2359: if (PetscAbsInt(formDegree) > 0) { /* get the trace (k-1)-forms on the first facet, and the 1-forms on the edge,
2360: and wedge them together to create the remaining k-form dofs */
2361: PetscDualSpace trace, fiber;
2362: PetscDualSpace_Lag *tracel, *fiberl;
2363: PetscQuadrature intNodesTrace2, intNodesFiber2, intNodes2;
2364: PetscLagNodeIndices intNodeIndices2;
2365: Mat intMatTrace, intMatFiber, intMat2;
2366: PetscInt traceDegree = formDegree > 0 ? formDegree - 1 : formDegree + 1;
2367: PetscInt fiberDegree = formDegree > 0 ? 1 : -1;
2369: PetscDualSpaceCreateFacetSubspace_Lagrange(sp, NULL, tensorf, traceDegree, Ncopies, PETSC_TRUE, &trace);
2370: PetscDualSpaceCreateEdgeSubspace_Lagrange(sp, order, fiberDegree, 1, PETSC_TRUE, &fiber);
2371: tracel = (PetscDualSpace_Lag *)trace->data;
2372: fiberl = (PetscDualSpace_Lag *)fiber->data;
2373: if (!lag->vertIndices) PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));
2374: PetscDualSpaceGetInteriorData(trace, &intNodesTrace2, &intMatTrace);
2375: PetscDualSpaceGetInteriorData(fiber, &intNodesFiber2, &intMatFiber);
2376: if (intNodesTrace2 && intNodesFiber2) {
2377: PetscQuadratureCreateTensor(intNodesTrace2, intNodesFiber2, &intNodes2);
2378: MatTensorAltV(intMatTrace, intMatFiber, dim - 1, traceDegree, 1, fiberDegree, &intMat2);
2379: PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, traceDegree, fiberl->intNodeIndices, 1, fiberDegree, &intNodeIndices2);
2380: if (!intMat) {
2381: intMat = intMat2;
2382: intNodes = intNodes2;
2383: intNodeIndices = intNodeIndices2;
2384: } else {
2385: /* merge the matrices, quadrature points, and nodes */
2386: PetscInt nM;
2387: PetscInt nDof, nDof2;
2388: PetscInt *toMerged = NULL, *toMerged2 = NULL;
2389: PetscQuadrature merged = NULL;
2390: PetscLagNodeIndices intNodeIndicesMerged = NULL;
2391: Mat matMerged = NULL;
2393: MatGetSize(intMat, &nDof, NULL);
2394: MatGetSize(intMat2, &nDof2, NULL);
2395: PetscQuadraturePointsMerge(intNodes, intNodes2, &merged, &toMerged, &toMerged2);
2396: PetscQuadratureGetData(merged, NULL, NULL, &nM, NULL, NULL);
2397: MatricesMerge(intMat, intMat2, dim, formDegree, nM, toMerged, toMerged2, &matMerged);
2398: PetscLagNodeIndicesMerge(intNodeIndices, intNodeIndices2, &intNodeIndicesMerged);
2399: PetscFree(toMerged);
2400: PetscFree(toMerged2);
2401: MatDestroy(&intMat);
2402: MatDestroy(&intMat2);
2403: PetscQuadratureDestroy(&intNodes);
2404: PetscQuadratureDestroy(&intNodes2);
2405: PetscLagNodeIndicesDestroy(&intNodeIndices);
2406: PetscLagNodeIndicesDestroy(&intNodeIndices2);
2407: intNodes = merged;
2408: intMat = matMerged;
2409: intNodeIndices = intNodeIndicesMerged;
2410: if (!trimmed) {
2411: /* I think users expect that, when a node has a full basis for the k-forms,
2412: * they should be consecutive dofs. That isn't the case for trimmed spaces,
2413: * but is for some of the nodes in untrimmed spaces, so in that case we
2414: * sort them to group them by node */
2415: Mat intMatPerm;
2417: MatPermuteByNodeIdx(intMat, intNodeIndices, &intMatPerm);
2418: MatDestroy(&intMat);
2419: intMat = intMatPerm;
2420: }
2421: }
2422: }
2423: PetscDualSpaceDestroy(&fiber);
2424: PetscDualSpaceDestroy(&trace);
2425: }
2426: PetscQuadratureDestroy(&intNodesTrace);
2427: PetscQuadratureDestroy(&intNodesFiber);
2428: sp->intNodes = intNodes;
2429: sp->intMat = intMat;
2430: lag->intNodeIndices = intNodeIndices;
2431: {
2432: PetscInt nDofs = 0;
2434: if (intMat) MatGetSize(intMat, &nDofs, NULL);
2435: PetscSectionSetDof(section, 0, nDofs);
2436: }
2437: PetscDualSpaceSectionSetUp_Internal(sp, section);
2438: if (continuous) {
2439: PetscDualSpaceCreateAllDataFromInteriorData(sp);
2440: PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2441: } else {
2442: PetscObjectReference((PetscObject)intNodes);
2443: sp->allNodes = intNodes;
2444: PetscObjectReference((PetscObject)intMat);
2445: sp->allMat = intMat;
2446: PetscLagNodeIndicesReference(intNodeIndices);
2447: lag->allNodeIndices = intNodeIndices;
2448: }
2449: }
2450: PetscSectionGetStorageSize(section, &sp->spdim);
2451: PetscSectionGetConstrainedStorageSize(section, &sp->spintdim);
2452: PetscDualSpaceComputeFunctionalsFromAllData(sp);
2453: PetscFree2(pStratStart, pStratEnd);
2454: DMDestroy(&dmint);
2455: return 0;
2456: }
2458: /* Create a matrix that represents the transformation that DMPlexVecGetClosure() would need
2459: * to get the representation of the dofs for a mesh point if the mesh point had this orientation
2460: * relative to the cell */
2461: PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat)
2462: {
2463: PetscDualSpace_Lag *lag;
2464: DM dm;
2465: PetscLagNodeIndices vertIndices, intNodeIndices;
2466: PetscLagNodeIndices ni;
2467: PetscInt nodeIdxDim, nodeVecDim, nNodes;
2468: PetscInt formDegree;
2469: PetscInt *perm, *permOrnt;
2470: PetscInt *nnz;
2471: PetscInt n;
2472: PetscInt maxGroupSize;
2473: PetscScalar *V, *W, *work;
2474: Mat A;
2476: if (!sp->spintdim) {
2477: *symMat = NULL;
2478: return 0;
2479: }
2480: lag = (PetscDualSpace_Lag *)sp->data;
2481: vertIndices = lag->vertIndices;
2482: intNodeIndices = lag->intNodeIndices;
2483: PetscDualSpaceGetDM(sp, &dm);
2484: PetscDualSpaceGetFormDegree(sp, &formDegree);
2485: PetscNew(&ni);
2486: ni->refct = 1;
2487: ni->nodeIdxDim = nodeIdxDim = intNodeIndices->nodeIdxDim;
2488: ni->nodeVecDim = nodeVecDim = intNodeIndices->nodeVecDim;
2489: ni->nNodes = nNodes = intNodeIndices->nNodes;
2490: PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
2491: PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
2492: /* push forward the dofs by the symmetry of the reference element induced by ornt */
2493: PetscLagNodeIndicesPushForward(dm, vertIndices, 0, vertIndices, intNodeIndices, ornt, formDegree, ni->nodeIdx, ni->nodeVec);
2494: /* get the revlex order for both the original and transformed dofs */
2495: PetscLagNodeIndicesGetPermutation(intNodeIndices, &perm);
2496: PetscLagNodeIndicesGetPermutation(ni, &permOrnt);
2497: PetscMalloc1(nNodes, &nnz);
2498: for (n = 0, maxGroupSize = 0; n < nNodes;) { /* incremented in the loop */
2499: PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2500: PetscInt m, nEnd;
2501: PetscInt groupSize;
2502: /* for each group of dofs that have the same nodeIdx coordinate */
2503: for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2504: PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2505: PetscInt d;
2507: /* compare the oriented permutation indices */
2508: for (d = 0; d < nodeIdxDim; d++)
2509: if (mind[d] != nind[d]) break;
2510: if (d < nodeIdxDim) break;
2511: }
2512: /* permOrnt[[n, nEnd)] is a group of dofs that, under the symmetry are at the same location */
2514: /* the symmetry had better map the group of dofs with the same permuted nodeIdx
2515: * to a group of dofs with the same size, otherwise we messed up */
2516: if (PetscDefined(USE_DEBUG)) {
2517: PetscInt m;
2518: PetscInt *nind = &(intNodeIndices->nodeIdx[perm[n] * nodeIdxDim]);
2520: for (m = n + 1; m < nEnd; m++) {
2521: PetscInt *mind = &(intNodeIndices->nodeIdx[perm[m] * nodeIdxDim]);
2522: PetscInt d;
2524: /* compare the oriented permutation indices */
2525: for (d = 0; d < nodeIdxDim; d++)
2526: if (mind[d] != nind[d]) break;
2527: if (d < nodeIdxDim) break;
2528: }
2530: }
2531: groupSize = nEnd - n;
2532: /* each pushforward dof vector will be expressed in a basis of the unpermuted dofs */
2533: for (m = n; m < nEnd; m++) nnz[permOrnt[m]] = groupSize;
2535: maxGroupSize = PetscMax(maxGroupSize, nEnd - n);
2536: n = nEnd;
2537: }
2539: MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes, nNodes, 0, nnz, &A);
2540: PetscFree(nnz);
2541: PetscMalloc3(maxGroupSize * nodeVecDim, &V, maxGroupSize * nodeVecDim, &W, nodeVecDim * 2, &work);
2542: for (n = 0; n < nNodes;) { /* incremented in the loop */
2543: PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2544: PetscInt nEnd;
2545: PetscInt m;
2546: PetscInt groupSize;
2547: for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2548: PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2549: PetscInt d;
2551: /* compare the oriented permutation indices */
2552: for (d = 0; d < nodeIdxDim; d++)
2553: if (mind[d] != nind[d]) break;
2554: if (d < nodeIdxDim) break;
2555: }
2556: groupSize = nEnd - n;
2557: /* get all of the vectors from the original and all of the pushforward vectors */
2558: for (m = n; m < nEnd; m++) {
2559: PetscInt d;
2561: for (d = 0; d < nodeVecDim; d++) {
2562: V[(m - n) * nodeVecDim + d] = intNodeIndices->nodeVec[perm[m] * nodeVecDim + d];
2563: W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
2564: }
2565: }
2566: /* now we have to solve for W in terms of V: the systems isn't always square, but the span
2567: * of V and W should always be the same, so the solution of the normal equations works */
2568: {
2569: char transpose = 'N';
2570: PetscBLASInt bm = nodeVecDim;
2571: PetscBLASInt bn = groupSize;
2572: PetscBLASInt bnrhs = groupSize;
2573: PetscBLASInt blda = bm;
2574: PetscBLASInt bldb = bm;
2575: PetscBLASInt blwork = 2 * nodeVecDim;
2576: PetscBLASInt info;
2578: PetscCallBLAS("LAPACKgels", LAPACKgels_(&transpose, &bm, &bn, &bnrhs, V, &blda, W, &bldb, work, &blwork, &info));
2580: /* repack */
2581: {
2582: PetscInt i, j;
2584: for (i = 0; i < groupSize; i++) {
2585: for (j = 0; j < groupSize; j++) {
2586: /* notice the different leading dimension */
2587: V[i * groupSize + j] = W[i * nodeVecDim + j];
2588: }
2589: }
2590: }
2591: if (PetscDefined(USE_DEBUG)) {
2592: PetscReal res;
2594: /* check that the normal error is 0 */
2595: for (m = n; m < nEnd; m++) {
2596: PetscInt d;
2598: for (d = 0; d < nodeVecDim; d++) W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
2599: }
2600: res = 0.;
2601: for (PetscInt i = 0; i < groupSize; i++) {
2602: for (PetscInt j = 0; j < nodeVecDim; j++) {
2603: for (PetscInt k = 0; k < groupSize; k++) W[i * nodeVecDim + j] -= V[i * groupSize + k] * intNodeIndices->nodeVec[perm[n + k] * nodeVecDim + j];
2604: res += PetscAbsScalar(W[i * nodeVecDim + j]);
2605: }
2606: }
2608: }
2609: }
2610: MatSetValues(A, groupSize, &permOrnt[n], groupSize, &perm[n], V, INSERT_VALUES);
2611: n = nEnd;
2612: }
2613: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
2614: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
2615: *symMat = A;
2616: PetscFree3(V, W, work);
2617: PetscLagNodeIndicesDestroy(&ni);
2618: return 0;
2619: }
2621: #define BaryIndex(perEdge, a, b, c) (((b) * (2 * perEdge + 1 - (b))) / 2) + (c)
2623: #define CartIndex(perEdge, a, b) (perEdge * (a) + b)
2625: /* the existing interface for symmetries is insufficient for all cases:
2626: * - it should be sufficient for form degrees that are scalar (0 and n)
2627: * - it should be sufficient for hypercube dofs
2628: * - it isn't sufficient for simplex cells with non-scalar form degrees if
2629: * there are any dofs in the interior
2630: *
2631: * We compute the general transformation matrices, and if they fit, we return them,
2632: * otherwise we error (but we should probably change the interface to allow for
2633: * these symmetries)
2634: */
2635: static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
2636: {
2637: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2638: PetscInt dim, order, Nc;
2640: PetscDualSpaceGetOrder(sp, &order);
2641: PetscDualSpaceGetNumComponents(sp, &Nc);
2642: DMGetDimension(sp->dm, &dim);
2643: if (!lag->symComputed) { /* store symmetries */
2644: PetscInt pStart, pEnd, p;
2645: PetscInt numPoints;
2646: PetscInt numFaces;
2647: PetscInt spintdim;
2648: PetscInt ***symperms;
2649: PetscScalar ***symflips;
2651: DMPlexGetChart(sp->dm, &pStart, &pEnd);
2652: numPoints = pEnd - pStart;
2653: {
2654: DMPolytopeType ct;
2655: /* The number of arrangements is no longer based on the number of faces */
2656: DMPlexGetCellType(sp->dm, 0, &ct);
2657: numFaces = DMPolytopeTypeGetNumArrangments(ct) / 2;
2658: }
2659: PetscCalloc1(numPoints, &symperms);
2660: PetscCalloc1(numPoints, &symflips);
2661: spintdim = sp->spintdim;
2662: /* The nodal symmetry behavior is not present when tensorSpace != tensorCell: someone might want this for the "S"
2663: * family of FEEC spaces. Most used in particular are discontinuous polynomial L2 spaces in tensor cells, where
2664: * the symmetries are not necessary for FE assembly. So for now we assume this is the case and don't return
2665: * symmetries if tensorSpace != tensorCell */
2666: if (spintdim && 0 < dim && dim < 3 && (lag->tensorSpace == lag->tensorCell)) { /* compute self symmetries */
2667: PetscInt **cellSymperms;
2668: PetscScalar **cellSymflips;
2669: PetscInt ornt;
2670: PetscInt nCopies = Nc / lag->intNodeIndices->nodeVecDim;
2671: PetscInt nNodes = lag->intNodeIndices->nNodes;
2673: lag->numSelfSym = 2 * numFaces;
2674: lag->selfSymOff = numFaces;
2675: PetscCalloc1(2 * numFaces, &cellSymperms);
2676: PetscCalloc1(2 * numFaces, &cellSymflips);
2677: /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
2678: symperms[0] = &cellSymperms[numFaces];
2679: symflips[0] = &cellSymflips[numFaces];
2682: for (ornt = -numFaces; ornt < numFaces; ornt++) { /* for every symmetry, compute the symmetry matrix, and extract rows to see if it fits in the perm + flip framework */
2683: Mat symMat;
2684: PetscInt *perm;
2685: PetscScalar *flips;
2686: PetscInt i;
2688: if (!ornt) continue;
2689: PetscMalloc1(spintdim, &perm);
2690: PetscCalloc1(spintdim, &flips);
2691: for (i = 0; i < spintdim; i++) perm[i] = -1;
2692: PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, ornt, &symMat);
2693: for (i = 0; i < nNodes; i++) {
2694: PetscInt ncols;
2695: PetscInt j, k;
2696: const PetscInt *cols;
2697: const PetscScalar *vals;
2698: PetscBool nz_seen = PETSC_FALSE;
2700: MatGetRow(symMat, i, &ncols, &cols, &vals);
2701: for (j = 0; j < ncols; j++) {
2702: if (PetscAbsScalar(vals[j]) > PETSC_SMALL) {
2704: nz_seen = PETSC_TRUE;
2708: for (k = 0; k < nCopies; k++) perm[cols[j] * nCopies + k] = i * nCopies + k;
2709: if (PetscRealPart(vals[j]) < 0.) {
2710: for (k = 0; k < nCopies; k++) flips[i * nCopies + k] = -1.;
2711: } else {
2712: for (k = 0; k < nCopies; k++) flips[i * nCopies + k] = 1.;
2713: }
2714: }
2715: }
2716: MatRestoreRow(symMat, i, &ncols, &cols, &vals);
2717: }
2718: MatDestroy(&symMat);
2719: /* if there were no sign flips, keep NULL */
2720: for (i = 0; i < spintdim; i++)
2721: if (flips[i] != 1.) break;
2722: if (i == spintdim) {
2723: PetscFree(flips);
2724: flips = NULL;
2725: }
2726: /* if the permutation is identity, keep NULL */
2727: for (i = 0; i < spintdim; i++)
2728: if (perm[i] != i) break;
2729: if (i == spintdim) {
2730: PetscFree(perm);
2731: perm = NULL;
2732: }
2733: symperms[0][ornt] = perm;
2734: symflips[0][ornt] = flips;
2735: }
2736: /* if no orientations produced non-identity permutations, keep NULL */
2737: for (ornt = -numFaces; ornt < numFaces; ornt++)
2738: if (symperms[0][ornt]) break;
2739: if (ornt == numFaces) {
2740: PetscFree(cellSymperms);
2741: symperms[0] = NULL;
2742: }
2743: /* if no orientations produced sign flips, keep NULL */
2744: for (ornt = -numFaces; ornt < numFaces; ornt++)
2745: if (symflips[0][ornt]) break;
2746: if (ornt == numFaces) {
2747: PetscFree(cellSymflips);
2748: symflips[0] = NULL;
2749: }
2750: }
2751: { /* get the symmetries of closure points */
2752: PetscInt closureSize = 0;
2753: PetscInt *closure = NULL;
2754: PetscInt r;
2756: DMPlexGetTransitiveClosure(sp->dm, 0, PETSC_TRUE, &closureSize, &closure);
2757: for (r = 0; r < closureSize; r++) {
2758: PetscDualSpace psp;
2759: PetscInt point = closure[2 * r];
2760: PetscInt pspintdim;
2761: const PetscInt ***psymperms = NULL;
2762: const PetscScalar ***psymflips = NULL;
2764: if (!point) continue;
2765: PetscDualSpaceGetPointSubspace(sp, point, &psp);
2766: if (!psp) continue;
2767: PetscDualSpaceGetInteriorDimension(psp, &pspintdim);
2768: if (!pspintdim) continue;
2769: PetscDualSpaceGetSymmetries(psp, &psymperms, &psymflips);
2770: symperms[r] = (PetscInt **)(psymperms ? psymperms[0] : NULL);
2771: symflips[r] = (PetscScalar **)(psymflips ? psymflips[0] : NULL);
2772: }
2773: DMPlexRestoreTransitiveClosure(sp->dm, 0, PETSC_TRUE, &closureSize, &closure);
2774: }
2775: for (p = 0; p < pEnd; p++)
2776: if (symperms[p]) break;
2777: if (p == pEnd) {
2778: PetscFree(symperms);
2779: symperms = NULL;
2780: }
2781: for (p = 0; p < pEnd; p++)
2782: if (symflips[p]) break;
2783: if (p == pEnd) {
2784: PetscFree(symflips);
2785: symflips = NULL;
2786: }
2787: lag->symperms = symperms;
2788: lag->symflips = symflips;
2789: lag->symComputed = PETSC_TRUE;
2790: }
2791: if (perms) *perms = (const PetscInt ***)lag->symperms;
2792: if (flips) *flips = (const PetscScalar ***)lag->symflips;
2793: return 0;
2794: }
2796: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2797: {
2798: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2802: *continuous = lag->continuous;
2803: return 0;
2804: }
2806: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2807: {
2808: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2811: lag->continuous = continuous;
2812: return 0;
2813: }
2815: /*@
2816: PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity
2818: Not Collective
2820: Input Parameter:
2821: . sp - the `PetscDualSpace`
2823: Output Parameter:
2824: . continuous - flag for element continuity
2826: Level: intermediate
2828: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeSetContinuity()`
2829: @*/
2830: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2831: {
2834: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace, PetscBool *), (sp, continuous));
2835: return 0;
2836: }
2838: /*@
2839: PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous
2841: Logically Collective on sp
2843: Input Parameters:
2844: + sp - the `PetscDualSpace`
2845: - continuous - flag for element continuity
2847: Options Database:
2848: . -petscdualspace_lagrange_continuity <bool> - use a continuous element
2850: Level: intermediate
2852: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeGetContinuity()`
2853: @*/
2854: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2855: {
2858: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace, PetscBool), (sp, continuous));
2859: return 0;
2860: }
2862: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
2863: {
2864: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2866: *tensor = lag->tensorSpace;
2867: return 0;
2868: }
2870: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
2871: {
2872: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2874: lag->tensorSpace = tensor;
2875: return 0;
2876: }
2878: static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Lagrange(PetscDualSpace sp, PetscBool *trimmed)
2879: {
2880: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2882: *trimmed = lag->trimmed;
2883: return 0;
2884: }
2886: static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp, PetscBool trimmed)
2887: {
2888: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2890: lag->trimmed = trimmed;
2891: return 0;
2892: }
2894: static PetscErrorCode PetscDualSpaceLagrangeGetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
2895: {
2896: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2898: if (nodeType) *nodeType = lag->nodeType;
2899: if (boundary) *boundary = lag->endNodes;
2900: if (exponent) *exponent = lag->nodeExponent;
2901: return 0;
2902: }
2904: static PetscErrorCode PetscDualSpaceLagrangeSetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
2905: {
2906: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2909: lag->nodeType = nodeType;
2910: lag->endNodes = boundary;
2911: lag->nodeExponent = exponent;
2912: return 0;
2913: }
2915: static PetscErrorCode PetscDualSpaceLagrangeGetUseMoments_Lagrange(PetscDualSpace sp, PetscBool *useMoments)
2916: {
2917: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2919: *useMoments = lag->useMoments;
2920: return 0;
2921: }
2923: static PetscErrorCode PetscDualSpaceLagrangeSetUseMoments_Lagrange(PetscDualSpace sp, PetscBool useMoments)
2924: {
2925: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2927: lag->useMoments = useMoments;
2928: return 0;
2929: }
2931: static PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt *momentOrder)
2932: {
2933: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2935: *momentOrder = lag->momentOrder;
2936: return 0;
2937: }
2939: static PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder_Lagrange(PetscDualSpace sp, PetscInt momentOrder)
2940: {
2941: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2943: lag->momentOrder = momentOrder;
2944: return 0;
2945: }
2947: /*@
2948: PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space
2950: Not collective
2952: Input Parameter:
2953: . sp - The `PetscDualSpace`
2955: Output Parameter:
2956: . tensor - Whether the dual space has tensor layout (vs. simplicial)
2958: Level: intermediate
2960: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeSetTensor()`, `PetscDualSpaceCreate()`
2961: @*/
2962: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
2963: {
2966: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetTensor_C", (PetscDualSpace, PetscBool *), (sp, tensor));
2967: return 0;
2968: }
2970: /*@
2971: PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space
2973: Not collective
2975: Input Parameters:
2976: + sp - The `PetscDualSpace`
2977: - tensor - Whether the dual space has tensor layout (vs. simplicial)
2979: Level: intermediate
2981: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeGetTensor()`, `PetscDualSpaceCreate()`
2982: @*/
2983: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
2984: {
2986: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetTensor_C", (PetscDualSpace, PetscBool), (sp, tensor));
2987: return 0;
2988: }
2990: /*@
2991: PetscDualSpaceLagrangeGetTrimmed - Get the trimmed nature of the dual space
2993: Not collective
2995: Input Parameter:
2996: . sp - The `PetscDualSpace`
2998: Output Parameter:
2999: . trimmed - Whether the dual space represents to dual basis of a trimmed polynomial space (e.g. Raviart-Thomas and higher order / other form degree variants)
3001: Level: intermediate
3003: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeSetTrimmed()`, `PetscDualSpaceCreate()`
3004: @*/
3005: PetscErrorCode PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp, PetscBool *trimmed)
3006: {
3009: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetTrimmed_C", (PetscDualSpace, PetscBool *), (sp, trimmed));
3010: return 0;
3011: }
3013: /*@
3014: PetscDualSpaceLagrangeSetTrimmed - Set the trimmed nature of the dual space
3016: Not collective
3018: Input Parameters:
3019: + sp - The `PetscDualSpace`
3020: - trimmed - Whether the dual space represents to dual basis of a trimmed polynomial space (e.g. Raviart-Thomas and higher order / other form degree variants)
3022: Level: intermediate
3024: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeGetTrimmed()`, `PetscDualSpaceCreate()`
3025: @*/
3026: PetscErrorCode PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp, PetscBool trimmed)
3027: {
3029: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetTrimmed_C", (PetscDualSpace, PetscBool), (sp, trimmed));
3030: return 0;
3031: }
3033: /*@
3034: PetscDualSpaceLagrangeGetNodeType - Get a description of how nodes are laid out for Lagrange polynomials in this
3035: dual space
3037: Not collective
3039: Input Parameter:
3040: . sp - The `PetscDualSpace`
3042: Output Parameters:
3043: + nodeType - The type of nodes
3044: . boundary - Whether the node type is one that includes endpoints (if nodeType is `PETSCDTNODES_GAUSSJACOBI`, nodes that
3045: include the boundary are Gauss-Lobatto-Jacobi nodes)
3046: - exponent - If nodeType is `PETSCDTNODES_GAUSJACOBI`, indicates the exponent used for both ends of the 1D Jacobi weight function
3047: '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type
3049: Level: advanced
3051: .seealso: `PetscDualSpace`, `PetscDTNodeType`, `PetscDualSpaceLagrangeSetNodeType()`
3052: @*/
3053: PetscErrorCode PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
3054: {
3059: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetNodeType_C", (PetscDualSpace, PetscDTNodeType *, PetscBool *, PetscReal *), (sp, nodeType, boundary, exponent));
3060: return 0;
3061: }
3063: /*@
3064: PetscDualSpaceLagrangeSetNodeType - Set a description of how nodes are laid out for Lagrange polynomials in this
3065: dual space
3067: Logically collective
3069: Input Parameters:
3070: + sp - The `PetscDualSpace`
3071: . nodeType - The type of nodes
3072: . boundary - Whether the node type is one that includes endpoints (if nodeType is `PETSCDTNODES_GAUSSJACOBI`, nodes that
3073: include the boundary are Gauss-Lobatto-Jacobi nodes)
3074: - exponent - If nodeType is `PETSCDTNODES_GAUSJACOBI`, indicates the exponent used for both ends of the 1D Jacobi weight function
3075: '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type
3077: Level: advanced
3079: .seealso: `PetscDualSpace`, `PetscDTNodeType`, `PetscDualSpaceLagrangeGetNodeType()`
3080: @*/
3081: PetscErrorCode PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
3082: {
3084: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetNodeType_C", (PetscDualSpace, PetscDTNodeType, PetscBool, PetscReal), (sp, nodeType, boundary, exponent));
3085: return 0;
3086: }
3088: /*@
3089: PetscDualSpaceLagrangeGetUseMoments - Get the flag for using moment functionals
3091: Not collective
3093: Input Parameter:
3094: . sp - The `PetscDualSpace`
3096: Output Parameter:
3097: . useMoments - Moment flag
3099: Level: advanced
3101: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeSetUseMoments()`
3102: @*/
3103: PetscErrorCode PetscDualSpaceLagrangeGetUseMoments(PetscDualSpace sp, PetscBool *useMoments)
3104: {
3107: PetscUseMethod(sp, "PetscDualSpaceLagrangeGetUseMoments_C", (PetscDualSpace, PetscBool *), (sp, useMoments));
3108: return 0;
3109: }
3111: /*@
3112: PetscDualSpaceLagrangeSetUseMoments - Set the flag for moment functionals
3114: Logically collective
3116: Input Parameters:
3117: + sp - The `PetscDualSpace`
3118: - useMoments - The flag for moment functionals
3120: Level: advanced
3122: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeGetUseMoments()`
3123: @*/
3124: PetscErrorCode PetscDualSpaceLagrangeSetUseMoments(PetscDualSpace sp, PetscBool useMoments)
3125: {
3127: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetUseMoments_C", (PetscDualSpace, PetscBool), (sp, useMoments));
3128: return 0;
3129: }
3131: /*@
3132: PetscDualSpaceLagrangeGetMomentOrder - Get the order for moment integration
3134: Not collective
3136: Input Parameter:
3137: . sp - The `PetscDualSpace`
3139: Output Parameter:
3140: . order - Moment integration order
3142: Level: advanced
3144: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeSetMomentOrder()`
3145: @*/
3146: PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder(PetscDualSpace sp, PetscInt *order)
3147: {
3150: PetscUseMethod(sp, "PetscDualSpaceLagrangeGetMomentOrder_C", (PetscDualSpace, PetscInt *), (sp, order));
3151: return 0;
3152: }
3154: /*@
3155: PetscDualSpaceLagrangeSetMomentOrder - Set the order for moment integration
3157: Logically collective
3159: Input Parameters:
3160: + sp - The `PetscDualSpace`
3161: - order - The order for moment integration
3163: Level: advanced
3165: .seealso: `PetscDualSpace`, `PetscDualSpaceLagrangeGetMomentOrder()`
3166: @*/
3167: PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder(PetscDualSpace sp, PetscInt order)
3168: {
3170: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetMomentOrder_C", (PetscDualSpace, PetscInt), (sp, order));
3171: return 0;
3172: }
3174: static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
3175: {
3176: sp->ops->destroy = PetscDualSpaceDestroy_Lagrange;
3177: sp->ops->view = PetscDualSpaceView_Lagrange;
3178: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
3179: sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange;
3180: sp->ops->setup = PetscDualSpaceSetUp_Lagrange;
3181: sp->ops->createheightsubspace = NULL;
3182: sp->ops->createpointsubspace = NULL;
3183: sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Lagrange;
3184: sp->ops->apply = PetscDualSpaceApplyDefault;
3185: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
3186: sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
3187: sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
3188: sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
3189: return 0;
3190: }
3192: /*MC
3193: PETSCDUALSPACELAGRANGE = "lagrange" - A `PetscDualSpaceType` that encapsulates a dual space of pointwise evaluation functionals
3195: Level: intermediate
3197: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`
3198: M*/
3199: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
3200: {
3201: PetscDualSpace_Lag *lag;
3204: PetscNew(&lag);
3205: sp->data = lag;
3207: lag->tensorCell = PETSC_FALSE;
3208: lag->tensorSpace = PETSC_FALSE;
3209: lag->continuous = PETSC_TRUE;
3210: lag->numCopies = PETSC_DEFAULT;
3211: lag->numNodeSkip = PETSC_DEFAULT;
3212: lag->nodeType = PETSCDTNODES_DEFAULT;
3213: lag->useMoments = PETSC_FALSE;
3214: lag->momentOrder = 0;
3216: PetscDualSpaceInitialize_Lagrange(sp);
3217: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
3218: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
3219: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
3220: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
3221: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Lagrange);
3222: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Lagrange);
3223: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Lagrange);
3224: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Lagrange);
3225: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", PetscDualSpaceLagrangeGetUseMoments_Lagrange);
3226: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", PetscDualSpaceLagrangeSetUseMoments_Lagrange);
3227: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", PetscDualSpaceLagrangeGetMomentOrder_Lagrange);
3228: PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", PetscDualSpaceLagrangeSetMomentOrder_Lagrange);
3229: return 0;
3230: }