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