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