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