Actual source code: dspacelagrange.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/petscfeimpl.h>
2: #include <petscdmplex.h>
3: #include <petscblaslapack.h>
5: PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]);
7: struct _n_Petsc1DNodeFamily
8: {
9: PetscInt refct;
10: PetscDTNodeType nodeFamily;
11: PetscReal gaussJacobiExp;
12: PetscInt nComputed;
13: PetscReal **nodesets;
14: PetscBool endpoints;
15: };
17: /* users set node families for PETSCDUALSPACELAGRANGE with just the inputs to this function, but internally we create
18: * an object that can cache the computations across multiple dual spaces */
19: static PetscErrorCode Petsc1DNodeFamilyCreate(PetscDTNodeType family, PetscReal gaussJacobiExp, PetscBool endpoints, Petsc1DNodeFamily *nf)
20: {
21: Petsc1DNodeFamily f;
25: PetscNew(&f);
26: switch (family) {
27: case PETSCDTNODES_GAUSSJACOBI:
28: case PETSCDTNODES_EQUISPACED:
29: f->nodeFamily = family;
30: break;
31: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family");
32: }
33: f->endpoints = endpoints;
34: f->gaussJacobiExp = 0.;
35: if (family == PETSCDTNODES_GAUSSJACOBI) {
36: if (gaussJacobiExp <= -1.) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Gauss-Jacobi exponent must be > -1.\n");
37: f->gaussJacobiExp = gaussJacobiExp;
38: }
39: f->refct = 1;
40: *nf = f;
41: return(0);
42: }
44: static PetscErrorCode Petsc1DNodeFamilyReference(Petsc1DNodeFamily nf)
45: {
47: if (nf) nf->refct++;
48: return(0);
49: }
51: static PetscErrorCode Petsc1DNodeFamilyDestroy(Petsc1DNodeFamily *nf) {
52: PetscInt i, nc;
56: if (!(*nf)) return(0);
57: if (--(*nf)->refct > 0) {
58: *nf = NULL;
59: return(0);
60: }
61: nc = (*nf)->nComputed;
62: for (i = 0; i < nc; i++) {
63: PetscFree((*nf)->nodesets[i]);
64: }
65: PetscFree((*nf)->nodesets);
66: PetscFree(*nf);
67: *nf = NULL;
68: return(0);
69: }
71: static PetscErrorCode Petsc1DNodeFamilyGetNodeSets(Petsc1DNodeFamily f, PetscInt degree, PetscReal ***nodesets)
72: {
73: PetscInt nc;
77: nc = f->nComputed;
78: if (degree >= nc) {
79: PetscInt i, j;
80: PetscReal **new_nodesets;
81: PetscReal *w;
83: PetscMalloc1(degree + 1, &new_nodesets);
84: PetscArraycpy(new_nodesets, f->nodesets, nc);
85: PetscFree(f->nodesets);
86: f->nodesets = new_nodesets;
87: PetscMalloc1(degree + 1, &w);
88: for (i = nc; i < degree + 1; i++) {
89: PetscMalloc1(i + 1, &(f->nodesets[i]));
90: if (!i) {
91: f->nodesets[i][0] = 0.5;
92: } else {
93: switch (f->nodeFamily) {
94: case PETSCDTNODES_EQUISPACED:
95: if (f->endpoints) {
96: for (j = 0; j <= i; j++) f->nodesets[i][j] = (PetscReal) j / (PetscReal) i;
97: } else {
98: /* these nodes are at the centroids of the small simplices created by the equispaced nodes that include
99: * the endpoints */
100: for (j = 0; j <= i; j++) f->nodesets[i][j] = ((PetscReal) j + 0.5) / ((PetscReal) i + 1.);
101: }
102: break;
103: case PETSCDTNODES_GAUSSJACOBI:
104: if (f->endpoints) {
105: PetscDTGaussLobattoJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w);
106: } else {
107: PetscDTGaussJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w);
108: }
109: break;
110: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family");
111: }
112: }
113: }
114: PetscFree(w);
115: f->nComputed = degree + 1;
116: }
117: *nodesets = f->nodesets;
118: return(0);
119: }
121: /* http://arxiv.org/abs/2002.09421 for details */
122: static PetscErrorCode PetscNodeRecursive_Internal(PetscInt dim, PetscInt degree, PetscReal **nodesets, PetscInt tup[], PetscReal node[])
123: {
124: PetscReal w;
125: PetscInt i, j;
129: w = 0.;
130: if (dim == 1) {
131: node[0] = nodesets[degree][tup[0]];
132: node[1] = nodesets[degree][tup[1]];
133: } else {
134: for (i = 0; i < dim + 1; i++) node[i] = 0.;
135: for (i = 0; i < dim + 1; i++) {
136: PetscReal wi = nodesets[degree][degree-tup[i]];
138: for (j = 0; j < dim+1; j++) tup[dim+1+j] = tup[j+(j>=i)];
139: PetscNodeRecursive_Internal(dim-1,degree-tup[i],nodesets,&tup[dim+1],&node[dim+1]);
140: for (j = 0; j < dim+1; j++) node[j+(j>=i)] += wi * node[dim+1+j];
141: w += wi;
142: }
143: for (i = 0; i < dim+1; i++) node[i] /= w;
144: }
145: return(0);
146: }
148: /* compute simplex nodes for the biunit simplex from the 1D node family */
149: static PetscErrorCode Petsc1DNodeFamilyComputeSimplexNodes(Petsc1DNodeFamily f, PetscInt dim, PetscInt degree, PetscReal points[])
150: {
151: PetscInt *tup;
152: PetscInt k;
153: PetscInt npoints;
154: PetscReal **nodesets = NULL;
155: PetscInt worksize;
156: PetscReal *nodework;
157: PetscInt *tupwork;
161: if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have non-negative dimension\n");
162: if (degree < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have non-negative degree\n");
163: if (!dim) return(0);
164: PetscCalloc1(dim+2, &tup);
165: k = 0;
166: PetscDTBinomialInt(degree + dim, dim, &npoints);
167: Petsc1DNodeFamilyGetNodeSets(f, degree, &nodesets);
168: worksize = ((dim + 2) * (dim + 3)) / 2;
169: PetscMalloc2(worksize, &nodework, worksize, &tupwork);
170: /* loop over the tuples of length dim with sum at most degree */
171: for (k = 0; k < npoints; k++) {
172: PetscInt i;
174: /* turn thm into tuples of length dim + 1 with sum equal to degree (barycentric indice) */
175: tup[0] = degree;
176: for (i = 0; i < dim; i++) {
177: tup[0] -= tup[i+1];
178: }
179: switch(f->nodeFamily) {
180: case PETSCDTNODES_EQUISPACED:
181: /* compute equispaces nodes on the unit reference triangle */
182: if (f->endpoints) {
183: for (i = 0; i < dim; i++) {
184: points[dim*k + i] = (PetscReal) tup[i+1] / (PetscReal) degree;
185: }
186: } else {
187: for (i = 0; i < dim; i++) {
188: /* these nodes are at the centroids of the small simplices created by the equispaced nodes that include
189: * the endpoints */
190: points[dim*k + i] = ((PetscReal) tup[i+1] + 1./(dim+1.)) / (PetscReal) (degree + 1.);
191: }
192: }
193: break;
194: default:
195: /* compute equispaces nodes on the barycentric reference triangle (the trace on the first dim dimensions are the
196: * unit reference triangle nodes */
197: for (i = 0; i < dim + 1; i++) tupwork[i] = tup[i];
198: PetscNodeRecursive_Internal(dim, degree, nodesets, tupwork, nodework);
199: for (i = 0; i < dim; i++) points[dim*k + i] = nodework[i + 1];
200: break;
201: }
202: PetscDualSpaceLatticePointLexicographic_Internal(dim, degree, &tup[1]);
203: }
204: /* map from unit simplex to biunit simplex */
205: for (k = 0; k < npoints * dim; k++) points[k] = points[k] * 2. - 1.;
206: PetscFree2(nodework, tupwork);
207: PetscFree(tup);
208: return(0);
209: }
211: /* 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
212: * on that mesh point, we have to be careful about getting/adding everything in the right place.
213: *
214: * With nodal dofs like PETSCDUALSPACELAGRANGE makes, the general approach to calculate the value of dofs associate
215: * with a node A is
216: * - transform the node locations x(A) by the map that takes the mesh point to its reorientation, x' = phi(x(A))
217: * - figure out which node was originally at the location of the transformed point, A' = idx(x')
218: * - if the dofs are not scalars, figure out how to represent the transformed dofs in terms of the basis
219: * of dofs at A' (using pushforward/pullback rules)
220: *
221: * The one sticky point with this approach is the "A' = idx(x')" step: trying to go from real valued coordinates
222: * back to indices. I don't want to rely on floating point tolerances. Additionally, PETSCDUALSPACELAGRANGE may
223: * eventually support quasi-Lagrangian dofs, which could involve quadrature at multiple points, so the location "x(A)"
224: * would be ambiguous.
225: *
226: * So each dof gets an integer value coordinate (nodeIdx in the structure below). The choice of integer coordinates
227: * is somewhat arbitrary, as long as all of the relevant symmetries of the mesh point correspond to *permutations* of
228: * the integer coordinates, which do not depend on numerical precision.
229: *
230: * So
231: *
232: * - DMPlexGetTransitiveClosure_Internal() tells me how an orientation turns into a permutation of the vertices of a
233: * mesh point
234: * - The permutation of the vertices, and the nodeIdx values assigned to them, tells what permutation in index space
235: * is associated with the orientation
236: * - I uses that permutation to get xi' = phi(xi(A)), the integer coordinate of the transformed dof
237: * - I can without numerical issues compute A' = idx(xi')
238: *
239: * Here are some examples of how the process works
240: *
241: * - With a triangle:
242: *
243: * The triangle has the following integer coordinates for vertices, taken from the barycentric triangle
244: *
245: * closure order 2
246: * nodeIdx (0,0,1)
247: * \
248: * +
249: * |\
250: * | \
251: * | \
252: * | \ closure order 1
253: * | \ / nodeIdx (0,1,0)
254: * +-----+
255: * \
256: * closure order 0
257: * nodeIdx (1,0,0)
258: *
259: * If I do DMPlexGetTransitiveClosure_Internal() with orientation 1, the vertices would appear
260: * in the order (1, 2, 0)
261: *
262: * If I list the nodeIdx of each vertex in closure order for orientation 0 (0, 1, 2) and orientation 1 (1, 2, 0), I
263: * see
264: *
265: * orientation 0 | orientation 1
266: *
267: * [0] (1,0,0) [1] (0,1,0)
268: * [1] (0,1,0) [2] (0,0,1)
269: * [2] (0,0,1) [0] (1,0,0)
270: * A B
271: *
272: * In other words, B is the result of a row permutation of A. But, there is also
273: * a column permutation that accomplishes the same result, (2,0,1).
274: *
275: * So if a dof has nodeIdx coordinate (a,b,c), after the transformation its nodeIdx coordinate
276: * is (c,a,b), and the transformed degree of freedom will be a linear combination of dofs
277: * that originally had coordinate (c,a,b).
278: *
279: * - With a quadrilateral:
280: *
281: * The quadrilateral has the following integer coordinates for vertices, taken from concatenating barycentric
282: * coordinates for two segments:
283: *
284: * closure order 3 closure order 2
285: * nodeIdx (1,0,0,1) nodeIdx (0,1,0,1)
286: * \ /
287: * +----+
288: * | |
289: * | |
290: * +----+
291: * / \
292: * closure order 0 closure order 1
293: * nodeIdx (1,0,1,0) nodeIdx (0,1,1,0)
294: *
295: * If I do DMPlexGetTransitiveClosure_Internal() with orientation 1, the vertices would appear
296: * in the order (1, 2, 3, 0)
297: *
298: * If I list the nodeIdx of each vertex in closure order for orientation 0 (0, 1, 2, 3) and
299: * orientation 1 (1, 2, 3, 0), I see
300: *
301: * orientation 0 | orientation 1
302: *
303: * [0] (1,0,1,0) [1] (0,1,1,0)
304: * [1] (0,1,1,0) [2] (0,1,0,1)
305: * [2] (0,1,0,1) [3] (1,0,0,1)
306: * [3] (1,0,0,1) [0] (1,0,1,0)
307: * A B
308: *
309: * The column permutation that accomplishes the same result is (3,2,0,1).
310: *
311: * So if a dof has nodeIdx coordinate (a,b,c,d), after the transformation its nodeIdx coordinate
312: * is (d,c,a,b), and the transformed degree of freedom will be a linear combination of dofs
313: * that originally had coordinate (d,c,a,b).
314: *
315: * Previously PETSCDUALSPACELAGRANGE had hardcoded symmetries for the triangle and quadrilateral,
316: * but this approach will work for any polytope, such as the wedge (triangular prism).
317: */
318: struct _n_PetscLagNodeIndices
319: {
320: PetscInt refct;
321: PetscInt nodeIdxDim;
322: PetscInt nodeVecDim;
323: PetscInt nNodes;
324: PetscInt *nodeIdx; /* for each node an index of size nodeIdxDim */
325: PetscReal *nodeVec; /* for each node a vector of size nodeVecDim */
326: PetscInt *perm; /* if these are vertices, perm takes DMPlex point index to closure order;
327: if these are nodes, perm lists nodes in index revlex order */
328: };
330: /* this is just here so I can access the values in tests/ex1.c outside the library */
331: PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices ni, PetscInt *nodeIdxDim, PetscInt *nodeVecDim, PetscInt *nNodes, const PetscInt *nodeIdx[], const PetscReal *nodeVec[])
332: {
334: *nodeIdxDim = ni->nodeIdxDim;
335: *nodeVecDim = ni->nodeVecDim;
336: *nNodes = ni->nNodes;
337: *nodeIdx = ni->nodeIdx;
338: *nodeVec = ni->nodeVec;
339: return(0);
340: }
342: static PetscErrorCode PetscLagNodeIndicesReference(PetscLagNodeIndices ni)
343: {
345: if (ni) ni->refct++;
346: return(0);
347: }
349: static PetscErrorCode PetscLagNodeIndicesDestroy(PetscLagNodeIndices *ni) {
353: if (!(*ni)) return(0);
354: if (--(*ni)->refct > 0) {
355: *ni = NULL;
356: return(0);
357: }
358: PetscFree((*ni)->nodeIdx);
359: PetscFree((*ni)->nodeVec);
360: PetscFree((*ni)->perm);
361: PetscFree(*ni);
362: *ni = NULL;
363: return(0);
364: }
366: /* The vertices are given nodeIdx coordinates (e.g. the corners of the barycentric triangle). Those coordinates are
367: * in some other order, and to understand the effect of different symmetries, we need them to be in closure order.
368: *
369: * If sortIdx is PETSC_FALSE, the coordinates are already in revlex order, otherwise we must sort them
370: * to that order before we do the real work of this function, which is
371: *
372: * - mark the vertices in closure order
373: * - sort them in revlex order
374: * - use the resulting permutation to list the vertex coordinates in closure order
375: */
376: static PetscErrorCode PetscLagNodeIndicesComputeVertexOrder(DM dm, PetscLagNodeIndices ni, PetscBool sortIdx)
377: {
378: PetscInt v, w, vStart, vEnd, c, d;
379: PetscInt nVerts;
380: PetscInt closureSize = 0;
381: PetscInt *closure = NULL;
382: PetscInt *closureOrder;
383: PetscInt *invClosureOrder;
384: PetscInt *revlexOrder;
385: PetscInt *newNodeIdx;
386: PetscInt dim;
387: Vec coordVec;
388: const PetscScalar *coords;
389: PetscErrorCode ierr;
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;
480: PetscNew(&ni);
481: DMGetDimension(dm, &dim);
482: ni->nodeIdxDim = dim + 1;
483: ni->nodeVecDim = 0;
484: ni->nNodes = dim + 1;
485: ni->refct = 1;
486: PetscCalloc1((dim + 1)*(dim + 1), &(ni->nodeIdx));
487: for (d = 0; d < dim + 1; d++) ni->nodeIdx[d*(dim + 2)] = 1;
488: PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_FALSE);
489: *nodeIndices = ni;
490: return(0);
491: }
493: /* A polytope that is a tensor product of a facet and a segment.
494: * We take whatever coordinate system was being used for the facet
495: * and we concatenaty the barycentric coordinates for the vertices
496: * at the end of the segment, (1,0) and (0,1), to get a coordinate
497: * system for the tensor product element */
498: static PetscErrorCode PetscLagNodeIndicesCreateTensorVertices(DM dm, PetscLagNodeIndices facetni, PetscLagNodeIndices *nodeIndices)
499: {
500: PetscLagNodeIndices ni;
501: PetscInt nodeIdxDim, subNodeIdxDim = facetni->nodeIdxDim;
502: PetscInt nVerts, nSubVerts = facetni->nNodes;
503: PetscInt dim, d, e, f, g;
508: PetscNew(&ni);
509: DMGetDimension(dm, &dim);
510: ni->nodeIdxDim = nodeIdxDim = subNodeIdxDim + 2;
511: ni->nodeVecDim = 0;
512: ni->nNodes = nVerts = 2 * nSubVerts;
513: ni->refct = 1;
514: PetscCalloc1(nodeIdxDim * nVerts, &(ni->nodeIdx));
515: for (f = 0, d = 0; d < 2; d++) {
516: for (e = 0; e < nSubVerts; e++, f++) {
517: for (g = 0; g < subNodeIdxDim; g++) {
518: ni->nodeIdx[f * nodeIdxDim + g] = facetni->nodeIdx[e * subNodeIdxDim + g];
519: }
520: ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim] = (1 - d);
521: ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim + 1] = d;
522: }
523: }
524: PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_TRUE);
525: *nodeIndices = ni;
526: return(0);
527: }
529: /* This helps us compute symmetries, and it also helps us compute coordinates for dofs that are being pushed
530: * forward from a boundary mesh point.
531: *
532: * Input:
533: *
534: * dm - the target reference cell where we want new coordinates and dof directions to be valid
535: * vert - the vertex coordinate system for the target reference cell
536: * p - the point in the target reference cell that the dofs are coming from
537: * vertp - the vertex coordinate system for p's reference cell
538: * ornt - the resulting coordinates and dof vectors will be for p under this orientation
539: * nodep - the node coordinates and dof vectors in p's reference cell
540: * formDegree - the form degree that the dofs transform as
541: *
542: * Output:
543: *
544: * pfNodeIdx - the node coordinates for p's dofs, in the dm reference cell, from the ornt perspective
545: * pfNodeVec - the node dof vectors for p's dofs, in the dm reference cell, from the ornt perspective
546: */
547: static PetscErrorCode PetscLagNodeIndicesPushForward(DM dm, PetscLagNodeIndices vert, PetscInt p, PetscLagNodeIndices vertp, PetscLagNodeIndices nodep, PetscInt ornt, PetscInt formDegree, PetscInt pfNodeIdx[], PetscReal pfNodeVec[])
548: {
549: PetscInt *closureVerts;
550: PetscInt closureSize = 0;
551: PetscInt *closure = NULL;
552: PetscInt dim, pdim, c, i, j, k, n, v, vStart, vEnd;
553: PetscInt nSubVert = vertp->nNodes;
554: PetscInt nodeIdxDim = vert->nodeIdxDim;
555: PetscInt subNodeIdxDim = vertp->nodeIdxDim;
556: PetscInt nNodes = nodep->nNodes;
557: const PetscInt *vertIdx = vert->nodeIdx;
558: const PetscInt *subVertIdx = vertp->nodeIdx;
559: const PetscInt *nodeIdx = nodep->nodeIdx;
560: const PetscReal *nodeVec = nodep->nodeVec;
561: PetscReal *J, *Jstar;
562: PetscReal detJ;
563: PetscInt depth, pdepth, Nk, pNk;
564: Vec coordVec;
565: PetscScalar *newCoords = NULL;
566: const PetscScalar *oldCoords = NULL;
567: PetscErrorCode ierr;
570: DMGetDimension(dm, &dim);
571: DMPlexGetDepth(dm, &depth);
572: DMGetCoordinatesLocal(dm, &coordVec);
573: DMPlexGetPointDepth(dm, p, &pdepth);
574: pdim = pdepth != depth ? pdepth != 0 ? pdepth : 0 : dim;
575: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
576: DMGetWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);
577: DMPlexGetTransitiveClosure_Internal(dm, p, ornt, PETSC_TRUE, &closureSize, &closure);
578: c = closureSize - nSubVert;
579: /* we want which cell closure indices the closure of this point corresponds to */
580: for (v = 0; v < nSubVert; v++) closureVerts[v] = vert->perm[closure[2 * (c + v)] - vStart];
581: DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
582: /* push forward indices */
583: for (i = 0; i < nodeIdxDim; i++) { /* for every component of the target index space */
584: /* check if this is a component that all vertices around this point have in common */
585: for (j = 1; j < nSubVert; j++) {
586: if (vertIdx[closureVerts[j] * nodeIdxDim + i] != vertIdx[closureVerts[0] * nodeIdxDim + i]) break;
587: }
588: if (j == nSubVert) { /* all vertices have this component in common, directly copy to output */
589: PetscInt val = vertIdx[closureVerts[0] * nodeIdxDim + i];
590: for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = val;
591: } else {
592: PetscInt subi = -1;
593: /* there must be a component in vertp that looks the same */
594: for (k = 0; k < subNodeIdxDim; k++) {
595: for (j = 0; j < nSubVert; j++) {
596: if (vertIdx[closureVerts[j] * nodeIdxDim + i] != subVertIdx[j * subNodeIdxDim + k]) break;
597: }
598: if (j == nSubVert) {
599: subi = k;
600: break;
601: }
602: }
603: if (subi < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not find matching coordinate\n");
604: /* that component in the vertp system becomes component i in the vert system for each dof */
605: for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = nodeIdx[n * subNodeIdxDim + subi];
606: }
607: }
608: /* push forward vectors */
609: DMGetWorkArray(dm, dim * dim, MPIU_REAL, &J);
610: if (ornt != 0) { /* temporarily change the coordinate vector so
611: DMPlexComputeCellGeometryAffineFEM gives us the Jacobian we want */
612: PetscInt closureSize2 = 0;
613: PetscInt *closure2 = NULL;
615: DMPlexGetTransitiveClosure_Internal(dm, p, 0, PETSC_TRUE, &closureSize2, &closure2);
616: PetscMalloc1(dim * nSubVert, &newCoords);
617: VecGetArrayRead(coordVec, &oldCoords);
618: for (v = 0; v < nSubVert; v++) {
619: PetscInt d;
620: for (d = 0; d < dim; d++) {
621: newCoords[(closure2[2 * (c + v)] - vStart) * dim + d] = oldCoords[closureVerts[v] * dim + d];
622: }
623: }
624: VecRestoreArrayRead(coordVec, &oldCoords);
625: DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize2, &closure2);
626: VecPlaceArray(coordVec, newCoords);
627: }
628: DMPlexComputeCellGeometryAffineFEM(dm, p, NULL, J, NULL, &detJ);
629: if (ornt != 0) {
630: VecResetArray(coordVec);
631: PetscFree(newCoords);
632: }
633: DMRestoreWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);
634: /* compactify */
635: for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
636: /* We have the Jacobian mapping the point's reference cell to this reference cell:
637: * pulling back a function to the point and applying the dof is what we want,
638: * so we get the pullback matrix and multiply the dof by that matrix on the right */
639: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
640: PetscDTBinomialInt(pdim, PetscAbsInt(formDegree), &pNk);
641: DMGetWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);
642: PetscDTAltVPullbackMatrix(pdim, dim, J, formDegree, Jstar);
643: for (n = 0; n < nNodes; n++) {
644: for (i = 0; i < Nk; i++) {
645: PetscReal val = 0.;
646: for (j = 0; j < pNk; j++) val += nodeVec[n * pNk + j] * Jstar[j * pNk + i];
647: pfNodeVec[n * Nk + i] = val;
648: }
649: }
650: DMRestoreWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);
651: DMRestoreWorkArray(dm, dim * dim, MPIU_REAL, &J);
652: return(0);
653: }
655: /* given to sets of nodes, take the tensor product, where the product of the dof indices is concatenation and the
656: * product of the dof vectors is the wedge product */
657: static PetscErrorCode PetscLagNodeIndicesTensor(PetscLagNodeIndices tracei, PetscInt dimT, PetscInt kT, PetscLagNodeIndices fiberi, PetscInt dimF, PetscInt kF, PetscLagNodeIndices *nodeIndices)
658: {
659: PetscInt dim = dimT + dimF;
660: PetscInt nodeIdxDim, nNodes;
661: PetscInt formDegree = kT + kF;
662: PetscInt Nk, NkT, NkF;
663: PetscInt MkT, MkF;
664: PetscLagNodeIndices ni;
665: PetscInt i, j, l;
666: PetscReal *projF, *projT;
667: PetscReal *projFstar, *projTstar;
668: PetscReal *workF, *workF2, *workT, *workT2, *work, *work2;
669: PetscReal *wedgeMat;
670: PetscReal sign;
674: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
675: PetscDTBinomialInt(dimT, PetscAbsInt(kT), &NkT);
676: PetscDTBinomialInt(dimF, PetscAbsInt(kF), &NkF);
677: PetscDTBinomialInt(dim, PetscAbsInt(kT), &MkT);
678: PetscDTBinomialInt(dim, PetscAbsInt(kF), &MkF);
679: PetscNew(&ni);
680: ni->nodeIdxDim = nodeIdxDim = tracei->nodeIdxDim + fiberi->nodeIdxDim;
681: ni->nodeVecDim = Nk;
682: ni->nNodes = nNodes = tracei->nNodes * fiberi->nNodes;
683: ni->refct = 1;
684: PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
685: /* first concatenate the indices */
686: for (l = 0, j = 0; j < fiberi->nNodes; j++) {
687: for (i = 0; i < tracei->nNodes; i++, l++) {
688: PetscInt m, n = 0;
690: for (m = 0; m < tracei->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = tracei->nodeIdx[i * tracei->nodeIdxDim + m];
691: for (m = 0; m < fiberi->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = fiberi->nodeIdx[j * fiberi->nodeIdxDim + m];
692: }
693: }
695: /* now wedge together the push-forward vectors */
696: PetscMalloc1(nNodes * Nk, &(ni->nodeVec));
697: PetscCalloc2(dimT*dim, &projT, dimF*dim, &projF);
698: for (i = 0; i < dimT; i++) projT[i * (dim + 1)] = 1.;
699: for (i = 0; i < dimF; i++) projF[i * (dim + dimT + 1) + dimT] = 1.;
700: PetscMalloc2(MkT*NkT, &projTstar, MkF*NkF, &projFstar);
701: PetscDTAltVPullbackMatrix(dim, dimT, projT, kT, projTstar);
702: PetscDTAltVPullbackMatrix(dim, dimF, projF, kF, projFstar);
703: PetscMalloc6(MkT, &workT, MkT, &workT2, MkF, &workF, MkF, &workF2, Nk, &work, Nk, &work2);
704: PetscMalloc1(Nk * MkT, &wedgeMat);
705: sign = (PetscAbsInt(kT * kF) & 1) ? -1. : 1.;
706: for (l = 0, j = 0; j < fiberi->nNodes; j++) {
707: PetscInt d, e;
709: /* push forward fiber k-form */
710: for (d = 0; d < MkF; d++) {
711: PetscReal val = 0.;
712: for (e = 0; e < NkF; e++) val += projFstar[d * NkF + e] * fiberi->nodeVec[j * NkF + e];
713: workF[d] = val;
714: }
715: /* Hodge star to proper form if necessary */
716: if (kF < 0) {
717: for (d = 0; d < MkF; d++) workF2[d] = workF[d];
718: PetscDTAltVStar(dim, PetscAbsInt(kF), 1, workF2, workF);
719: }
720: /* Compute the matrix that wedges this form with one of the trace k-form */
721: PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kF), PetscAbsInt(kT), workF, wedgeMat);
722: for (i = 0; i < tracei->nNodes; i++, l++) {
723: /* push forward trace k-form */
724: for (d = 0; d < MkT; d++) {
725: PetscReal val = 0.;
726: for (e = 0; e < NkT; e++) val += projTstar[d * NkT + e] * tracei->nodeVec[i * NkT + e];
727: workT[d] = val;
728: }
729: /* Hodge star to proper form if necessary */
730: if (kT < 0) {
731: for (d = 0; d < MkT; d++) workT2[d] = workT[d];
732: PetscDTAltVStar(dim, PetscAbsInt(kT), 1, workT2, workT);
733: }
734: /* compute the wedge product of the push-forward trace form and firer forms */
735: for (d = 0; d < Nk; d++) {
736: PetscReal val = 0.;
737: for (e = 0; e < MkT; e++) val += wedgeMat[d * MkT + e] * workT[e];
738: work[d] = val;
739: }
740: /* inverse Hodge star from proper form if necessary */
741: if (formDegree < 0) {
742: for (d = 0; d < Nk; d++) work2[d] = work[d];
743: PetscDTAltVStar(dim, PetscAbsInt(formDegree), -1, work2, work);
744: }
745: /* insert into the array (adjusting for sign) */
746: for (d = 0; d < Nk; d++) ni->nodeVec[l * Nk + d] = sign * work[d];
747: }
748: }
749: PetscFree(wedgeMat);
750: PetscFree6(workT, workT2, workF, workF2, work, work2);
751: PetscFree2(projTstar, projFstar);
752: PetscFree2(projT, projF);
753: *nodeIndices = ni;
754: return(0);
755: }
757: /* simple union of two sets of nodes */
758: static PetscErrorCode PetscLagNodeIndicesMerge(PetscLagNodeIndices niA, PetscLagNodeIndices niB, PetscLagNodeIndices *nodeIndices)
759: {
760: PetscLagNodeIndices ni;
761: PetscInt nodeIdxDim, nodeVecDim, nNodes;
762: PetscErrorCode ierr;
765: PetscNew(&ni);
766: ni->nodeIdxDim = nodeIdxDim = niA->nodeIdxDim;
767: if (niB->nodeIdxDim != nodeIdxDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeIdxDim");
768: ni->nodeVecDim = nodeVecDim = niA->nodeVecDim;
769: if (niB->nodeVecDim != nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeVecDim");
770: ni->nNodes = nNodes = niA->nNodes + niB->nNodes;
771: ni->refct = 1;
772: PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
773: PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
774: PetscArraycpy(ni->nodeIdx, niA->nodeIdx, niA->nNodes * nodeIdxDim);
775: PetscArraycpy(ni->nodeVec, niA->nodeVec, niA->nNodes * nodeVecDim);
776: PetscArraycpy(&(ni->nodeIdx[niA->nNodes * nodeIdxDim]), niB->nodeIdx, niB->nNodes * nodeIdxDim);
777: PetscArraycpy(&(ni->nodeVec[niA->nNodes * nodeVecDim]), niB->nodeVec, niB->nNodes * nodeVecDim);
778: *nodeIndices = ni;
779: return(0);
780: }
782: #define PETSCTUPINTCOMPREVLEX(N) \
783: static int PetscTupIntCompRevlex_##N(const void *a, const void *b) \
784: { \
785: const PetscInt *A = (const PetscInt *) a; \
786: const PetscInt *B = (const PetscInt *) b; \
787: int i; \
788: PetscInt diff = 0; \
789: for (i = 0; i < N; i++) { \
790: diff = A[N - i] - B[N - i]; \
791: if (diff) break; \
792: } \
793: return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1; \
794: }
796: PETSCTUPINTCOMPREVLEX(3)
797: PETSCTUPINTCOMPREVLEX(4)
798: PETSCTUPINTCOMPREVLEX(5)
799: PETSCTUPINTCOMPREVLEX(6)
800: PETSCTUPINTCOMPREVLEX(7)
802: static int PetscTupIntCompRevlex_N(const void *a, const void *b)
803: {
804: const PetscInt *A = (const PetscInt *) a;
805: const PetscInt *B = (const PetscInt *) b;
806: int i;
807: int N = A[0];
808: PetscInt diff = 0;
809: for (i = 0; i < N; i++) {
810: diff = A[N - i] - B[N - i];
811: if (diff) break;
812: }
813: return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1;
814: }
816: /* The nodes are not necessarily in revlex order wrt nodeIdx: get the permutation
817: * that puts them in that order */
818: static PetscErrorCode PetscLagNodeIndicesGetPermutation(PetscLagNodeIndices ni, PetscInt *perm[])
819: {
823: if (!(ni->perm)) {
824: PetscInt *sorter;
825: PetscInt m = ni->nNodes;
826: PetscInt nodeIdxDim = ni->nodeIdxDim;
827: PetscInt i, j, k, l;
828: PetscInt *prm;
829: int (*comp) (const void *, const void *);
831: PetscMalloc1((nodeIdxDim + 2) * m, &sorter);
832: for (k = 0, l = 0, i = 0; i < m; i++) {
833: sorter[k++] = nodeIdxDim + 1;
834: sorter[k++] = i;
835: for (j = 0; j < nodeIdxDim; j++) sorter[k++] = ni->nodeIdx[l++];
836: }
837: switch (nodeIdxDim) {
838: case 2:
839: comp = PetscTupIntCompRevlex_3;
840: break;
841: case 3:
842: comp = PetscTupIntCompRevlex_4;
843: break;
844: case 4:
845: comp = PetscTupIntCompRevlex_5;
846: break;
847: case 5:
848: comp = PetscTupIntCompRevlex_6;
849: break;
850: case 6:
851: comp = PetscTupIntCompRevlex_7;
852: break;
853: default:
854: comp = PetscTupIntCompRevlex_N;
855: break;
856: }
857: qsort(sorter, m, (nodeIdxDim + 2) * sizeof(PetscInt), comp);
858: PetscMalloc1(m, &prm);
859: for (i = 0; i < m; i++) prm[i] = sorter[(nodeIdxDim + 2) * i + 1];
860: ni->perm = prm;
861: PetscFree(sorter);
862: }
863: *perm = ni->perm;
864: return(0);
865: }
867: static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
868: {
869: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
870: PetscErrorCode ierr;
873: if (lag->symperms) {
874: PetscInt **selfSyms = lag->symperms[0];
876: if (selfSyms) {
877: PetscInt i, **allocated = &selfSyms[-lag->selfSymOff];
879: for (i = 0; i < lag->numSelfSym; i++) {
880: PetscFree(allocated[i]);
881: }
882: PetscFree(allocated);
883: }
884: PetscFree(lag->symperms);
885: }
886: if (lag->symflips) {
887: PetscScalar **selfSyms = lag->symflips[0];
889: if (selfSyms) {
890: PetscInt i;
891: PetscScalar **allocated = &selfSyms[-lag->selfSymOff];
893: for (i = 0; i < lag->numSelfSym; i++) {
894: PetscFree(allocated[i]);
895: }
896: PetscFree(allocated);
897: }
898: PetscFree(lag->symflips);
899: }
900: Petsc1DNodeFamilyDestroy(&(lag->nodeFamily));
901: PetscLagNodeIndicesDestroy(&(lag->vertIndices));
902: PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));
903: PetscLagNodeIndicesDestroy(&(lag->allNodeIndices));
904: PetscFree(lag);
905: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
906: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
907: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);
908: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);
909: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL);
910: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL);
911: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL);
912: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL);
913: return(0);
914: }
916: static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
917: {
918: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
919: PetscErrorCode ierr;
922: PetscViewerASCIIPrintf(viewer, "%s %s%sLagrange dual space\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "", lag->trimmed ? "trimmed " : "");
923: return(0);
924: }
926: static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
927: {
928: PetscBool iascii;
934: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
935: if (iascii) {PetscDualSpaceLagrangeView_Ascii(sp, viewer);}
936: return(0);
937: }
939: static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
940: {
941: PetscBool continuous, tensor, trimmed, flg, flg2, flg3;
942: PetscDTNodeType nodeType;
943: PetscReal nodeExponent;
944: PetscBool nodeEndpoints;
948: PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
949: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
950: PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
951: PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &nodeEndpoints, &nodeExponent);
952: if (nodeType == PETSCDTNODES_DEFAULT) nodeType = PETSCDTNODES_GAUSSJACOBI;
953: PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
954: PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
955: if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
956: PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetTensor", tensor, &tensor, &flg);
957: if (flg) {PetscDualSpaceLagrangeSetTensor(sp, tensor);}
958: PetscOptionsBool("-petscdualspace_lagrange_trimmed", "Flag for trimmed dual space", "PetscDualSpaceLagrangeSetTrimmed", trimmed, &trimmed, &flg);
959: if (flg) {PetscDualSpaceLagrangeSetTrimmed(sp, trimmed);}
960: PetscOptionsEnum("-petscdualspace_lagrange_node_type", "Lagrange node location type", "PetscDualSpaceLagrangeSetNodeType", PetscDTNodeTypes, (PetscEnum)nodeType, (PetscEnum *)&nodeType, &flg);
961: PetscOptionsBool("-petscdualspace_lagrange_node_endpoints", "Flag for nodes that include endpoints", "PetscDualSpaceLagrangeSetNodeType", nodeEndpoints, &nodeEndpoints, &flg2);
962: flg3 = PETSC_FALSE;
963: if (nodeType == PETSCDTNODES_GAUSSJACOBI) {
964: PetscOptionsReal("-petscdualspace_lagrange_node_exponent", "Gauss-Jacobi weight function exponent", "PetscDualSpaceLagrangeSetNodeType", nodeExponent, &nodeExponent, &flg3);
965: }
966: if (flg || flg2 || flg3) {PetscDualSpaceLagrangeSetNodeType(sp, nodeType, nodeEndpoints, nodeExponent);}
967: PetscOptionsTail();
968: return(0);
969: }
971: static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace spNew)
972: {
973: PetscBool cont, tensor, trimmed, boundary;
974: PetscDTNodeType nodeType;
975: PetscReal exponent;
976: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
977: PetscErrorCode ierr;
980: PetscDualSpaceLagrangeGetContinuity(sp, &cont);
981: PetscDualSpaceLagrangeSetContinuity(spNew, cont);
982: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
983: PetscDualSpaceLagrangeSetTensor(spNew, tensor);
984: PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);
985: PetscDualSpaceLagrangeSetTrimmed(spNew, trimmed);
986: PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &boundary, &exponent);
987: PetscDualSpaceLagrangeSetNodeType(spNew, nodeType, boundary, exponent);
988: if (lag->nodeFamily) {
989: PetscDualSpace_Lag *lagnew = (PetscDualSpace_Lag *) spNew->data;
991: Petsc1DNodeFamilyReference(lag->nodeFamily);
992: lagnew->nodeFamily = lag->nodeFamily;
993: }
994: return(0);
995: }
997: /* for making tensor product spaces: take a dual space and product a segment space that has all the same
998: * specifications (trimmed, continuous, order, node set), except for the form degree */
999: static PetscErrorCode PetscDualSpaceCreateEdgeSubspace_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt k, PetscInt Nc, PetscBool interiorOnly, PetscDualSpace *bdsp)
1000: {
1001: DM K;
1002: PetscDualSpace_Lag *newlag;
1003: PetscErrorCode ierr;
1006: PetscDualSpaceDuplicate(sp,bdsp);
1007: PetscDualSpaceSetFormDegree(*bdsp, k);
1008: PetscDualSpaceCreateReferenceCell(*bdsp, 1, PETSC_TRUE, &K);
1009: PetscDualSpaceSetDM(*bdsp, K);
1010: DMDestroy(&K);
1011: PetscDualSpaceSetOrder(*bdsp, order);
1012: PetscDualSpaceSetNumComponents(*bdsp, Nc);
1013: newlag = (PetscDualSpace_Lag *) (*bdsp)->data;
1014: newlag->interiorOnly = interiorOnly;
1015: PetscDualSpaceSetUp(*bdsp);
1016: return(0);
1017: }
1019: /* just the points, weights aren't handled */
1020: static PetscErrorCode PetscQuadratureCreateTensor(PetscQuadrature trace, PetscQuadrature fiber, PetscQuadrature *product)
1021: {
1022: PetscInt dimTrace, dimFiber;
1023: PetscInt numPointsTrace, numPointsFiber;
1024: PetscInt dim, numPoints;
1025: const PetscReal *pointsTrace;
1026: const PetscReal *pointsFiber;
1027: PetscReal *points;
1028: PetscInt i, j, k, p;
1029: PetscErrorCode ierr;
1032: PetscQuadratureGetData(trace, &dimTrace, NULL, &numPointsTrace, &pointsTrace, NULL);
1033: PetscQuadratureGetData(fiber, &dimFiber, NULL, &numPointsFiber, &pointsFiber, NULL);
1034: dim = dimTrace + dimFiber;
1035: numPoints = numPointsFiber * numPointsTrace;
1036: PetscMalloc1(numPoints * dim, &points);
1037: for (p = 0, j = 0; j < numPointsFiber; j++) {
1038: for (i = 0; i < numPointsTrace; i++, p++) {
1039: for (k = 0; k < dimTrace; k++) points[p * dim + k] = pointsTrace[i * dimTrace + k];
1040: for (k = 0; k < dimFiber; k++) points[p * dim + dimTrace + k] = pointsFiber[j * dimFiber + k];
1041: }
1042: }
1043: PetscQuadratureCreate(PETSC_COMM_SELF, product);
1044: PetscQuadratureSetData(*product, dim, 0, numPoints, points, NULL);
1045: return(0);
1046: }
1048: /* Kronecker tensor product where matrix is considered a matrix of k-forms, so that
1049: * the entries in the product matrix are wedge products of the entries in the original matrices */
1050: static PetscErrorCode MatTensorAltV(Mat trace, Mat fiber, PetscInt dimTrace, PetscInt kTrace, PetscInt dimFiber, PetscInt kFiber, Mat *product)
1051: {
1052: PetscInt mTrace, nTrace, mFiber, nFiber, m, n, k, i, j, l;
1053: PetscInt dim, NkTrace, NkFiber, Nk;
1054: PetscInt dT, dF;
1055: PetscInt *nnzTrace, *nnzFiber, *nnz;
1056: PetscInt iT, iF, jT, jF, il, jl;
1057: PetscReal *workT, *workT2, *workF, *workF2, *work, *workstar;
1058: PetscReal *projT, *projF;
1059: PetscReal *projTstar, *projFstar;
1060: PetscReal *wedgeMat;
1061: PetscReal sign;
1062: PetscScalar *workS;
1063: Mat prod;
1064: /* this produces dof groups that look like the identity */
1068: MatGetSize(trace, &mTrace, &nTrace);
1069: PetscDTBinomialInt(dimTrace, PetscAbsInt(kTrace), &NkTrace);
1070: if (nTrace % NkTrace) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point value space of trace matrix is not a multiple of k-form size");
1071: MatGetSize(fiber, &mFiber, &nFiber);
1072: PetscDTBinomialInt(dimFiber, PetscAbsInt(kFiber), &NkFiber);
1073: if (nFiber % NkFiber) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point value space of fiber matrix is not a multiple of k-form size");
1074: PetscMalloc2(mTrace, &nnzTrace, mFiber, &nnzFiber);
1075: for (i = 0; i < mTrace; i++) {
1076: MatGetRow(trace, i, &(nnzTrace[i]), NULL, NULL);
1077: if (nnzTrace[i] % NkTrace) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in trace matrix are not in k-form size blocks");
1078: }
1079: for (i = 0; i < mFiber; i++) {
1080: MatGetRow(fiber, i, &(nnzFiber[i]), NULL, NULL);
1081: if (nnzFiber[i] % NkFiber) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in fiber matrix are not in k-form size blocks");
1082: }
1083: dim = dimTrace + dimFiber;
1084: k = kFiber + kTrace;
1085: PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1086: m = mTrace * mFiber;
1087: PetscMalloc1(m, &nnz);
1088: for (l = 0, j = 0; j < mFiber; j++) for (i = 0; i < mTrace; i++, l++) nnz[l] = (nnzTrace[i] / NkTrace) * (nnzFiber[j] / NkFiber) * Nk;
1089: n = (nTrace / NkTrace) * (nFiber / NkFiber) * Nk;
1090: MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &prod);
1091: PetscFree(nnz);
1092: PetscFree2(nnzTrace,nnzFiber);
1093: /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1094: MatSetOption(prod, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1095: /* compute pullbacks */
1096: PetscDTBinomialInt(dim, PetscAbsInt(kTrace), &dT);
1097: PetscDTBinomialInt(dim, PetscAbsInt(kFiber), &dF);
1098: PetscMalloc4(dimTrace * dim, &projT, dimFiber * dim, &projF, dT * NkTrace, &projTstar, dF * NkFiber, &projFstar);
1099: PetscArrayzero(projT, dimTrace * dim);
1100: for (i = 0; i < dimTrace; i++) projT[i * (dim + 1)] = 1.;
1101: PetscArrayzero(projF, dimFiber * dim);
1102: for (i = 0; i < dimFiber; i++) projF[i * (dim + 1) + dimTrace] = 1.;
1103: PetscDTAltVPullbackMatrix(dim, dimTrace, projT, kTrace, projTstar);
1104: PetscDTAltVPullbackMatrix(dim, dimFiber, projF, kFiber, projFstar);
1105: PetscMalloc5(dT, &workT, dF, &workF, Nk, &work, Nk, &workstar, Nk, &workS);
1106: PetscMalloc2(dT, &workT2, dF, &workF2);
1107: PetscMalloc1(Nk * dT, &wedgeMat);
1108: sign = (PetscAbsInt(kTrace * kFiber) & 1) ? -1. : 1.;
1109: for (i = 0, iF = 0; iF < mFiber; iF++) {
1110: PetscInt ncolsF, nformsF;
1111: const PetscInt *colsF;
1112: const PetscScalar *valsF;
1114: MatGetRow(fiber, iF, &ncolsF, &colsF, &valsF);
1115: nformsF = ncolsF / NkFiber;
1116: for (iT = 0; iT < mTrace; iT++, i++) {
1117: PetscInt ncolsT, nformsT;
1118: const PetscInt *colsT;
1119: const PetscScalar *valsT;
1121: MatGetRow(trace, iT, &ncolsT, &colsT, &valsT);
1122: nformsT = ncolsT / NkTrace;
1123: for (j = 0, jF = 0; jF < nformsF; jF++) {
1124: PetscInt colF = colsF[jF * NkFiber] / NkFiber;
1126: for (il = 0; il < dF; il++) {
1127: PetscReal val = 0.;
1128: for (jl = 0; jl < NkFiber; jl++) val += projFstar[il * NkFiber + jl] * PetscRealPart(valsF[jF * NkFiber + jl]);
1129: workF[il] = val;
1130: }
1131: if (kFiber < 0) {
1132: for (il = 0; il < dF; il++) workF2[il] = workF[il];
1133: PetscDTAltVStar(dim, PetscAbsInt(kFiber), 1, workF2, workF);
1134: }
1135: PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kFiber), PetscAbsInt(kTrace), workF, wedgeMat);
1136: for (jT = 0; jT < nformsT; jT++, j++) {
1137: PetscInt colT = colsT[jT * NkTrace] / NkTrace;
1138: PetscInt col = colF * (nTrace / NkTrace) + colT;
1139: const PetscScalar *vals;
1141: for (il = 0; il < dT; il++) {
1142: PetscReal val = 0.;
1143: for (jl = 0; jl < NkTrace; jl++) val += projTstar[il * NkTrace + jl] * PetscRealPart(valsT[jT * NkTrace + jl]);
1144: workT[il] = val;
1145: }
1146: if (kTrace < 0) {
1147: for (il = 0; il < dT; il++) workT2[il] = workT[il];
1148: PetscDTAltVStar(dim, PetscAbsInt(kTrace), 1, workT2, workT);
1149: }
1151: for (il = 0; il < Nk; il++) {
1152: PetscReal val = 0.;
1153: for (jl = 0; jl < dT; jl++) val += sign * wedgeMat[il * dT + jl] * workT[jl];
1154: work[il] = val;
1155: }
1156: if (k < 0) {
1157: PetscDTAltVStar(dim, PetscAbsInt(k), -1, work, workstar);
1158: #if defined(PETSC_USE_COMPLEX)
1159: for (l = 0; l < Nk; l++) workS[l] = workstar[l];
1160: vals = &workS[0];
1161: #else
1162: vals = &workstar[0];
1163: #endif
1164: } else {
1165: #if defined(PETSC_USE_COMPLEX)
1166: for (l = 0; l < Nk; l++) workS[l] = work[l];
1167: vals = &workS[0];
1168: #else
1169: vals = &work[0];
1170: #endif
1171: }
1172: for (l = 0; l < Nk; l++) {
1173: MatSetValue(prod, i, col * Nk + l, vals[l], INSERT_VALUES);
1174: } /* Nk */
1175: } /* jT */
1176: } /* jF */
1177: MatRestoreRow(trace, iT, &ncolsT, &colsT, &valsT);
1178: } /* iT */
1179: MatRestoreRow(fiber, iF, &ncolsF, &colsF, &valsF);
1180: } /* iF */
1181: PetscFree(wedgeMat);
1182: PetscFree4(projT, projF, projTstar, projFstar);
1183: PetscFree2(workT2, workF2);
1184: PetscFree5(workT, workF, work, workstar, workS);
1185: MatAssemblyBegin(prod, MAT_FINAL_ASSEMBLY);
1186: MatAssemblyEnd(prod, MAT_FINAL_ASSEMBLY);
1187: *product = prod;
1188: return(0);
1189: }
1191: /* Union of quadrature points, with an attempt to identify commont points in the two sets */
1192: static PetscErrorCode PetscQuadraturePointsMerge(PetscQuadrature quadA, PetscQuadrature quadB, PetscQuadrature *quadJoint, PetscInt *aToJoint[], PetscInt *bToJoint[])
1193: {
1194: PetscInt dimA, dimB;
1195: PetscInt nA, nB, nJoint, i, j, d;
1196: const PetscReal *pointsA;
1197: const PetscReal *pointsB;
1198: PetscReal *pointsJoint;
1199: PetscInt *aToJ, *bToJ;
1200: PetscQuadrature qJ;
1201: PetscErrorCode ierr;
1204: PetscQuadratureGetData(quadA, &dimA, NULL, &nA, &pointsA, NULL);
1205: PetscQuadratureGetData(quadB, &dimB, NULL, &nB, &pointsB, NULL);
1206: if (dimA != dimB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Quadrature points must be in the same dimension");
1207: nJoint = nA;
1208: PetscMalloc1(nA, &aToJ);
1209: for (i = 0; i < nA; i++) aToJ[i] = i;
1210: PetscMalloc1(nB, &bToJ);
1211: for (i = 0; i < nB; i++) {
1212: for (j = 0; j < nA; j++) {
1213: bToJ[i] = -1;
1214: for (d = 0; d < dimA; d++) if (PetscAbsReal(pointsB[i * dimA + d] - pointsA[j * dimA + d]) > PETSC_SMALL) break;
1215: if (d == dimA) {
1216: bToJ[i] = j;
1217: break;
1218: }
1219: }
1220: if (bToJ[i] == -1) {
1221: bToJ[i] = nJoint++;
1222: }
1223: }
1224: *aToJoint = aToJ;
1225: *bToJoint = bToJ;
1226: PetscMalloc1(nJoint * dimA, &pointsJoint);
1227: PetscArraycpy(pointsJoint, pointsA, nA * dimA);
1228: for (i = 0; i < nB; i++) {
1229: if (bToJ[i] >= nA) {
1230: for (d = 0; d < dimA; d++) pointsJoint[bToJ[i] * dimA + d] = pointsB[i * dimA + d];
1231: }
1232: }
1233: PetscQuadratureCreate(PETSC_COMM_SELF, &qJ);
1234: PetscQuadratureSetData(qJ, dimA, 0, nJoint, pointsJoint, NULL);
1235: *quadJoint = qJ;
1236: return(0);
1237: }
1239: /* Matrices matA and matB are both quadrature -> dof matrices: produce a matrix that is joint quadrature -> union of
1240: * dofs, where the joint quadrature was produced by PetscQuadraturePointsMerge */
1241: static PetscErrorCode MatricesMerge(Mat matA, Mat matB, PetscInt dim, PetscInt k, PetscInt numMerged, const PetscInt aToMerged[], const PetscInt bToMerged[], Mat *matMerged)
1242: {
1243: PetscInt m, n, mA, nA, mB, nB, Nk, i, j, l;
1244: Mat M;
1245: PetscInt *nnz;
1246: PetscInt maxnnz;
1247: PetscInt *work;
1251: PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1252: MatGetSize(matA, &mA, &nA);
1253: if (nA % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matA column space not a multiple of k-form size");
1254: MatGetSize(matB, &mB, &nB);
1255: if (nB % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matB column space not a multiple of k-form size");
1256: m = mA + mB;
1257: n = numMerged * Nk;
1258: PetscMalloc1(m, &nnz);
1259: maxnnz = 0;
1260: for (i = 0; i < mA; i++) {
1261: MatGetRow(matA, i, &(nnz[i]), NULL, NULL);
1262: if (nnz[i] % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matA are not in k-form size blocks");
1263: maxnnz = PetscMax(maxnnz, nnz[i]);
1264: }
1265: for (i = 0; i < mB; i++) {
1266: MatGetRow(matB, i, &(nnz[i+mA]), NULL, NULL);
1267: if (nnz[i+mA] % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matB are not in k-form size blocks");
1268: maxnnz = PetscMax(maxnnz, nnz[i+mA]);
1269: }
1270: MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &M);
1271: PetscFree(nnz);
1272: /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1273: MatSetOption(M, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1274: PetscMalloc1(maxnnz, &work);
1275: for (i = 0; i < mA; i++) {
1276: const PetscInt *cols;
1277: const PetscScalar *vals;
1278: PetscInt nCols;
1279: MatGetRow(matA, i, &nCols, &cols, &vals);
1280: for (j = 0; j < nCols / Nk; j++) {
1281: PetscInt newCol = aToMerged[cols[j * Nk] / Nk];
1282: for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1283: }
1284: MatSetValuesBlocked(M, 1, &i, nCols, work, vals, INSERT_VALUES);
1285: MatRestoreRow(matA, i, &nCols, &cols, &vals);
1286: }
1287: for (i = 0; i < mB; i++) {
1288: const PetscInt *cols;
1289: const PetscScalar *vals;
1291: PetscInt row = i + mA;
1292: PetscInt nCols;
1293: MatGetRow(matB, i, &nCols, &cols, &vals);
1294: for (j = 0; j < nCols / Nk; j++) {
1295: PetscInt newCol = bToMerged[cols[j * Nk] / Nk];
1296: for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l;
1297: }
1298: MatSetValuesBlocked(M, 1, &row, nCols, work, vals, INSERT_VALUES);
1299: MatRestoreRow(matB, i, &nCols, &cols, &vals);
1300: }
1301: PetscFree(work);
1302: MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
1303: MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
1304: *matMerged = M;
1305: return(0);
1306: }
1308: /* Take a dual space and product a segment space that has all the same specifications (trimmed, continuous, order,
1309: * node set), except for the form degree. For computing boundary dofs and for making tensor product spaces */
1310: static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Lagrange(PetscDualSpace sp, DM K, PetscInt f, PetscInt k, PetscInt Ncopies, PetscBool interiorOnly, PetscDualSpace *bdsp)
1311: {
1312: PetscInt Nknew, Ncnew;
1313: PetscInt dim, pointDim = -1;
1314: PetscInt depth;
1315: DM dm;
1316: PetscDualSpace_Lag *newlag;
1317: PetscErrorCode ierr;
1320: PetscDualSpaceGetDM(sp,&dm);
1321: DMGetDimension(dm,&dim);
1322: DMPlexGetDepth(dm,&depth);
1323: PetscDualSpaceDuplicate(sp,bdsp);
1324: PetscDualSpaceSetFormDegree(*bdsp,k);
1325: if (!K) {
1326: PetscBool isSimplex;
1329: if (depth == dim) {
1330: PetscInt coneSize;
1332: pointDim = dim - 1;
1333: DMPlexGetConeSize(dm,f,&coneSize);
1334: isSimplex = (PetscBool) (coneSize == dim);
1335: PetscDualSpaceCreateReferenceCell(*bdsp, dim-1, isSimplex, &K);
1336: } else if (depth == 1) {
1337: pointDim = 0;
1338: PetscDualSpaceCreateReferenceCell(*bdsp, 0, PETSC_TRUE, &K);
1339: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
1340: } else {
1341: PetscObjectReference((PetscObject)K);
1342: DMGetDimension(K, &pointDim);
1343: }
1344: PetscDualSpaceSetDM(*bdsp, K);
1345: DMDestroy(&K);
1346: PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew);
1347: Ncnew = Nknew * Ncopies;
1348: PetscDualSpaceSetNumComponents(*bdsp, Ncnew);
1349: newlag = (PetscDualSpace_Lag *) (*bdsp)->data;
1350: newlag->interiorOnly = interiorOnly;
1351: PetscDualSpaceSetUp(*bdsp);
1352: return(0);
1353: }
1355: /* Construct simplex nodes from a nodefamily, add Nk dof vectors of length Nk at each node.
1356: * Return the (quadrature, matrix) form of the dofs and the nodeIndices form as well.
1357: *
1358: * Sometimes we want a set of nodes to be contained in the interior of the element,
1359: * even when the node scheme puts nodes on the boundaries. numNodeSkip tells
1360: * the routine how many "layers" of nodes need to be skipped.
1361: * */
1362: static PetscErrorCode PetscDualSpaceLagrangeCreateSimplexNodeMat(Petsc1DNodeFamily nodeFamily, PetscInt dim, PetscInt sum, PetscInt Nk, PetscInt numNodeSkip, PetscQuadrature *iNodes, Mat *iMat, PetscLagNodeIndices *nodeIndices)
1363: {
1364: PetscReal *extraNodeCoords, *nodeCoords;
1365: PetscInt nNodes, nExtraNodes;
1366: PetscInt i, j, k, extraSum = sum + numNodeSkip * (1 + dim);
1367: PetscQuadrature intNodes;
1368: Mat intMat;
1369: PetscLagNodeIndices ni;
1373: PetscDTBinomialInt(dim + sum, dim, &nNodes);
1374: PetscDTBinomialInt(dim + extraSum, dim, &nExtraNodes);
1376: PetscMalloc1(dim * nExtraNodes, &extraNodeCoords);
1377: PetscNew(&ni);
1378: ni->nodeIdxDim = dim + 1;
1379: ni->nodeVecDim = Nk;
1380: ni->nNodes = nNodes * Nk;
1381: ni->refct = 1;
1382: PetscMalloc1(nNodes * Nk * (dim + 1), &(ni->nodeIdx));
1383: PetscMalloc1(nNodes * Nk * Nk, &(ni->nodeVec));
1384: 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.;
1385: Petsc1DNodeFamilyComputeSimplexNodes(nodeFamily, dim, extraSum, extraNodeCoords);
1386: if (numNodeSkip) {
1387: PetscInt k;
1388: PetscInt *tup;
1390: PetscMalloc1(dim * nNodes, &nodeCoords);
1391: PetscMalloc1(dim + 1, &tup);
1392: for (k = 0; k < nNodes; k++) {
1393: PetscInt j, c;
1394: PetscInt index;
1396: PetscDTIndexToBary(dim + 1, sum, k, tup);
1397: for (j = 0; j < dim + 1; j++) tup[j] += numNodeSkip;
1398: for (c = 0; c < Nk; c++) {
1399: for (j = 0; j < dim + 1; j++) {
1400: ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1401: }
1402: }
1403: PetscDTBaryToIndex(dim + 1, extraSum, tup, &index);
1404: for (j = 0; j < dim; j++) nodeCoords[k * dim + j] = extraNodeCoords[index * dim + j];
1405: }
1406: PetscFree(tup);
1407: PetscFree(extraNodeCoords);
1408: } else {
1409: PetscInt k;
1410: PetscInt *tup;
1412: nodeCoords = extraNodeCoords;
1413: PetscMalloc1(dim + 1, &tup);
1414: for (k = 0; k < nNodes; k++) {
1415: PetscInt j, c;
1417: PetscDTIndexToBary(dim + 1, sum, k, tup);
1418: for (c = 0; c < Nk; c++) {
1419: for (j = 0; j < dim + 1; j++) {
1420: /* barycentric indices can have zeros, but we don't want to push forward zeros because it makes it harder to
1421: * determine which nodes correspond to which under symmetries, so we increase by 1. This is fine
1422: * because the nodeIdx coordinates don't have any meaning other than helping to identify symmetries */
1423: ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1;
1424: }
1425: }
1426: }
1427: PetscFree(tup);
1428: }
1429: PetscQuadratureCreate(PETSC_COMM_SELF, &intNodes);
1430: PetscQuadratureSetData(intNodes, dim, 0, nNodes, nodeCoords, NULL);
1431: MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes * Nk, nNodes * Nk, Nk, NULL, &intMat);
1432: MatSetOption(intMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);
1433: for (j = 0; j < nNodes * Nk; j++) {
1434: PetscInt rem = j % Nk;
1435: PetscInt a, aprev = j - rem;
1436: PetscInt anext = aprev + Nk;
1438: for (a = aprev; a < anext; a++) {
1439: MatSetValue(intMat, j, a, (a == j) ? 1. : 0., INSERT_VALUES);
1440: }
1441: }
1442: MatAssemblyBegin(intMat, MAT_FINAL_ASSEMBLY);
1443: MatAssemblyEnd(intMat, MAT_FINAL_ASSEMBLY);
1444: *iNodes = intNodes;
1445: *iMat = intMat;
1446: *nodeIndices = ni;
1447: return(0);
1448: }
1450: /* once the nodeIndices have been created for the interior of the reference cell, and for all of the boundary cells,
1451: * push forward the boudary dofs and concatenate them into the full node indices for the dual space */
1452: static PetscErrorCode PetscDualSpaceLagrangeCreateAllNodeIdx(PetscDualSpace sp)
1453: {
1454: DM dm;
1455: PetscInt dim, nDofs;
1456: PetscSection section;
1457: PetscInt pStart, pEnd, p;
1458: PetscInt formDegree, Nk;
1459: PetscInt nodeIdxDim, spintdim;
1460: PetscDualSpace_Lag *lag;
1461: PetscLagNodeIndices ni, verti;
1465: lag = (PetscDualSpace_Lag *) sp->data;
1466: verti = lag->vertIndices;
1467: PetscDualSpaceGetDM(sp, &dm);
1468: DMGetDimension(dm, &dim);
1469: PetscDualSpaceGetFormDegree(sp, &formDegree);
1470: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
1471: PetscDualSpaceGetSection(sp, §ion);
1472: PetscSectionGetStorageSize(section, &nDofs);
1473: PetscNew(&ni);
1474: ni->nodeIdxDim = nodeIdxDim = verti->nodeIdxDim;
1475: ni->nodeVecDim = Nk;
1476: ni->nNodes = nDofs;
1477: ni->refct = 1;
1478: PetscMalloc1(nodeIdxDim * nDofs, &(ni->nodeIdx));
1479: PetscMalloc1(Nk * nDofs, &(ni->nodeVec));
1480: DMPlexGetChart(dm, &pStart, &pEnd);
1481: PetscSectionGetDof(section, 0, &spintdim);
1482: if (spintdim) {
1483: PetscArraycpy(ni->nodeIdx, lag->intNodeIndices->nodeIdx, spintdim * nodeIdxDim);
1484: PetscArraycpy(ni->nodeVec, lag->intNodeIndices->nodeVec, spintdim * Nk);
1485: }
1486: for (p = pStart + 1; p < pEnd; p++) {
1487: PetscDualSpace psp = sp->pointSpaces[p];
1488: PetscDualSpace_Lag *plag;
1489: PetscInt dof, off;
1491: PetscSectionGetDof(section, p, &dof);
1492: if (!dof) continue;
1493: plag = (PetscDualSpace_Lag *) psp->data;
1494: PetscSectionGetOffset(section, p, &off);
1495: PetscLagNodeIndicesPushForward(dm, verti, p, plag->vertIndices, plag->intNodeIndices, 0, formDegree, &(ni->nodeIdx[off * nodeIdxDim]), &(ni->nodeVec[off * Nk]));
1496: }
1497: lag->allNodeIndices = ni;
1498: return(0);
1499: }
1501: /* once the (quadrature, Matrix) forms of the dofs have been created for the interior of the
1502: * reference cell and for the boundary cells, jk
1503: * push forward the boundary data and concatenate them into the full (quadrature, matrix) data
1504: * for the dual space */
1505: static PetscErrorCode PetscDualSpaceCreateAllDataFromInteriorData(PetscDualSpace sp)
1506: {
1507: DM dm;
1508: PetscSection section;
1509: PetscInt pStart, pEnd, p, k, Nk, dim, Nc;
1510: PetscInt nNodes;
1511: PetscInt countNodes;
1512: Mat allMat;
1513: PetscQuadrature allNodes;
1514: PetscInt nDofs;
1515: PetscInt maxNzforms, j;
1516: PetscScalar *work;
1517: PetscReal *L, *J, *Jinv, *v0, *pv0;
1518: PetscInt *iwork;
1519: PetscReal *nodes;
1520: PetscErrorCode ierr;
1523: PetscDualSpaceGetDM(sp, &dm);
1524: DMGetDimension(dm, &dim);
1525: PetscDualSpaceGetSection(sp, §ion);
1526: PetscSectionGetStorageSize(section, &nDofs);
1527: DMPlexGetChart(dm, &pStart, &pEnd);
1528: PetscDualSpaceGetFormDegree(sp, &k);
1529: PetscDualSpaceGetNumComponents(sp, &Nc);
1530: PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1531: for (p = pStart, nNodes = 0, maxNzforms = 0; p < pEnd; p++) {
1532: PetscDualSpace psp;
1533: DM pdm;
1534: PetscInt pdim, pNk;
1535: PetscQuadrature intNodes;
1536: Mat intMat;
1538: PetscDualSpaceGetPointSubspace(sp, p, &psp);
1539: if (!psp) continue;
1540: PetscDualSpaceGetDM(psp, &pdm);
1541: DMGetDimension(pdm, &pdim);
1542: if (pdim < PetscAbsInt(k)) continue;
1543: PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
1544: PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
1545: if (intNodes) {
1546: PetscInt nNodesp;
1548: PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, NULL, NULL);
1549: nNodes += nNodesp;
1550: }
1551: if (intMat) {
1552: PetscInt maxNzsp;
1553: PetscInt maxNzformsp;
1555: MatSeqAIJGetMaxRowNonzeros(intMat, &maxNzsp);
1556: if (maxNzsp % pNk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms");
1557: maxNzformsp = maxNzsp / pNk;
1558: maxNzforms = PetscMax(maxNzforms, maxNzformsp);
1559: }
1560: }
1561: MatCreateSeqAIJ(PETSC_COMM_SELF, nDofs, nNodes * Nc, maxNzforms * Nk, NULL, &allMat);
1562: MatSetOption(allMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);
1563: PetscMalloc7(dim, &v0, dim, &pv0, dim * dim, &J, dim * dim, &Jinv, Nk * Nk, &L, maxNzforms * Nk, &work, maxNzforms * Nk, &iwork);
1564: for (j = 0; j < dim; j++) pv0[j] = -1.;
1565: PetscMalloc1(dim * nNodes, &nodes);
1566: for (p = pStart, countNodes = 0; p < pEnd; p++) {
1567: PetscDualSpace psp;
1568: PetscQuadrature intNodes;
1569: DM pdm;
1570: PetscInt pdim, pNk;
1571: PetscInt countNodesIn = countNodes;
1572: PetscReal detJ;
1573: Mat intMat;
1575: PetscDualSpaceGetPointSubspace(sp, p, &psp);
1576: if (!psp) continue;
1577: PetscDualSpaceGetDM(psp, &pdm);
1578: DMGetDimension(pdm, &pdim);
1579: if (pdim < PetscAbsInt(k)) continue;
1580: PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);
1581: if (intNodes == NULL && intMat == NULL) continue;
1582: PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);
1583: if (p) {
1584: DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, Jinv, &detJ);
1585: } else { /* identity */
1586: PetscInt i,j;
1588: for (i = 0; i < dim; i++) for (j = 0; j < dim; j++) J[i * dim + j] = Jinv[i * dim + j] = 0.;
1589: for (i = 0; i < dim; i++) J[i * dim + i] = Jinv[i * dim + i] = 1.;
1590: for (i = 0; i < dim; i++) v0[i] = -1.;
1591: }
1592: if (pdim != dim) { /* compactify Jacobian */
1593: PetscInt i, j;
1595: for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
1596: }
1597: PetscDTAltVPullbackMatrix(pdim, dim, J, k, L);
1598: if (intNodes) { /* push forward quadrature locations by the affine transformation */
1599: PetscInt nNodesp;
1600: const PetscReal *nodesp;
1601: PetscInt j;
1603: PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, &nodesp, NULL);
1604: for (j = 0; j < nNodesp; j++, countNodes++) {
1605: PetscInt d, e;
1607: for (d = 0; d < dim; d++) {
1608: nodes[countNodes * dim + d] = v0[d];
1609: for (e = 0; e < pdim; e++) {
1610: nodes[countNodes * dim + d] += J[d * pdim + e] * (nodesp[j * pdim + e] - pv0[e]);
1611: }
1612: }
1613: }
1614: }
1615: if (intMat) {
1616: PetscInt nrows;
1617: PetscInt off;
1619: PetscSectionGetDof(section, p, &nrows);
1620: PetscSectionGetOffset(section, p, &off);
1621: for (j = 0; j < nrows; j++) {
1622: PetscInt ncols;
1623: const PetscInt *cols;
1624: const PetscScalar *vals;
1625: PetscInt l, d, e;
1626: PetscInt row = j + off;
1628: MatGetRow(intMat, j, &ncols, &cols, &vals);
1629: if (ncols % pNk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms");
1630: for (l = 0; l < ncols / pNk; l++) {
1631: PetscInt blockcol;
1633: for (d = 0; d < pNk; d++) {
1634: if ((cols[l * pNk + d] % pNk) != d) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms");
1635: }
1636: blockcol = cols[l * pNk] / pNk;
1637: for (d = 0; d < Nk; d++) {
1638: iwork[l * Nk + d] = (blockcol + countNodesIn) * Nk + d;
1639: }
1640: for (d = 0; d < Nk; d++) work[l * Nk + d] = 0.;
1641: for (d = 0; d < Nk; d++) {
1642: for (e = 0; e < pNk; e++) {
1643: /* "push forward" dof by pulling back a k-form to be evaluated on the point: multiply on the right by L */
1644: work[l * Nk + d] += vals[l * pNk + e] * L[e * pNk + d];
1645: }
1646: }
1647: }
1648: MatSetValues(allMat, 1, &row, (ncols / pNk) * Nk, iwork, work, INSERT_VALUES);
1649: MatRestoreRow(intMat, j, &ncols, &cols, &vals);
1650: }
1651: }
1652: }
1653: MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);
1654: MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);
1655: PetscQuadratureCreate(PETSC_COMM_SELF, &allNodes);
1656: PetscQuadratureSetData(allNodes, dim, 0, nNodes, nodes, NULL);
1657: PetscFree7(v0, pv0, J, Jinv, L, work, iwork);
1658: MatDestroy(&(sp->allMat));
1659: sp->allMat = allMat;
1660: PetscQuadratureDestroy(&(sp->allNodes));
1661: sp->allNodes = allNodes;
1662: return(0);
1663: }
1665: /* rather than trying to get all data from the functionals, we create
1666: * the functionals from rows of the quadrature -> dof matrix.
1667: *
1668: * Ideally most of the uses of PetscDualSpace in PetscFE will switch
1669: * to using intMat and allMat, so that the individual functionals
1670: * don't need to be constructed at all */
1671: static PetscErrorCode PetscDualSpaceComputeFunctionalsFromAllData(PetscDualSpace sp)
1672: {
1673: PetscQuadrature allNodes;
1674: Mat allMat;
1675: PetscInt nDofs;
1676: PetscInt dim, k, Nk, Nc, f;
1677: DM dm;
1678: PetscInt nNodes, spdim;
1679: const PetscReal *nodes = NULL;
1680: PetscSection section;
1681: PetscErrorCode ierr;
1684: PetscDualSpaceGetDM(sp, &dm);
1685: DMGetDimension(dm, &dim);
1686: PetscDualSpaceGetNumComponents(sp, &Nc);
1687: PetscDualSpaceGetFormDegree(sp, &k);
1688: PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);
1689: PetscDualSpaceGetAllData(sp, &allNodes, &allMat);
1690: nNodes = 0;
1691: if (allNodes) {
1692: PetscQuadratureGetData(allNodes, NULL, NULL, &nNodes, &nodes, NULL);
1693: }
1694: MatGetSize(allMat, &nDofs, NULL);
1695: PetscDualSpaceGetSection(sp, §ion);
1696: PetscSectionGetStorageSize(section, &spdim);
1697: if (spdim != nDofs) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "incompatible all matrix size");
1698: PetscMalloc1(nDofs, &(sp->functional));
1699: for (f = 0; f < nDofs; f++) {
1700: PetscInt ncols, c;
1701: const PetscInt *cols;
1702: const PetscScalar *vals;
1703: PetscReal *nodesf;
1704: PetscReal *weightsf;
1705: PetscInt nNodesf;
1706: PetscInt countNodes;
1708: MatGetRow(allMat, f, &ncols, &cols, &vals);
1709: if (ncols % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "all matrix is not laid out as blocks of k-forms");
1710: for (c = 1, nNodesf = 1; c < ncols; c++) {
1711: if ((cols[c] / Nc) != (cols[c-1] / Nc)) nNodesf++;
1712: }
1713: PetscMalloc1(dim * nNodesf, &nodesf);
1714: PetscMalloc1(Nc * nNodesf, &weightsf);
1715: for (c = 0, countNodes = 0; c < ncols; c++) {
1716: if (!c || ((cols[c] / Nc) != (cols[c-1] / Nc))) {
1717: PetscInt d;
1719: for (d = 0; d < Nc; d++) {
1720: weightsf[countNodes * Nc + d] = 0.;
1721: }
1722: for (d = 0; d < dim; d++) {
1723: nodesf[countNodes * dim + d] = nodes[(cols[c] / Nc) * dim + d];
1724: }
1725: countNodes++;
1726: }
1727: weightsf[(countNodes - 1) * Nc + (cols[c] % Nc)] = PetscRealPart(vals[c]);
1728: }
1729: PetscQuadratureCreate(PETSC_COMM_SELF, &(sp->functional[f]));
1730: PetscQuadratureSetData(sp->functional[f], dim, Nc, nNodesf, nodesf, weightsf);
1731: MatRestoreRow(allMat, f, &ncols, &cols, &vals);
1732: }
1733: return(0);
1734: }
1736: /* take a matrix meant for k-forms and expand it to one for Ncopies */
1737: static PetscErrorCode PetscDualSpaceLagrangeMatrixCreateCopies(Mat A, PetscInt Nk, PetscInt Ncopies, Mat *Abs)
1738: {
1739: PetscInt m, n, i, j, k;
1740: PetscInt maxnnz, *nnz, *iwork;
1741: Mat Ac;
1745: MatGetSize(A, &m, &n);
1746: if (n % Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of columns in A %D is not a multiple of Nk %D", n, Nk);
1747: PetscMalloc1(m * Ncopies, &nnz);
1748: for (i = 0, maxnnz = 0; i < m; i++) {
1749: PetscInt innz;
1750: MatGetRow(A, i, &innz, NULL, NULL);
1751: if (innz % Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "A row %D nnzs is not a multiple of Nk %D", innz, Nk);
1752: for (j = 0; j < Ncopies; j++) nnz[i * Ncopies + j] = innz;
1753: maxnnz = PetscMax(maxnnz, innz);
1754: }
1755: MatCreateSeqAIJ(PETSC_COMM_SELF, m * Ncopies, n * Ncopies, 0, nnz, &Ac);
1756: MatSetOption(Ac, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);
1757: PetscFree(nnz);
1758: PetscMalloc1(maxnnz, &iwork);
1759: for (i = 0; i < m; i++) {
1760: PetscInt innz;
1761: const PetscInt *cols;
1762: const PetscScalar *vals;
1764: MatGetRow(A, i, &innz, &cols, &vals);
1765: for (j = 0; j < innz; j++) iwork[j] = (cols[j] / Nk) * (Nk * Ncopies) + (cols[j] % Nk);
1766: for (j = 0; j < Ncopies; j++) {
1767: PetscInt row = i * Ncopies + j;
1769: MatSetValues(Ac, 1, &row, innz, iwork, vals, INSERT_VALUES);
1770: for (k = 0; k < innz; k++) iwork[k] += Nk;
1771: }
1772: MatRestoreRow(A, i, &innz, &cols, &vals);
1773: }
1774: PetscFree(iwork);
1775: MatAssemblyBegin(Ac, MAT_FINAL_ASSEMBLY);
1776: MatAssemblyEnd(Ac, MAT_FINAL_ASSEMBLY);
1777: *Abs = Ac;
1778: return(0);
1779: }
1781: /* check if a cell is a tensor product of the segment with a facet,
1782: * specifically checking if f and f2 can be the "endpoints" (like the triangles
1783: * at either end of a wedge) */
1784: static PetscErrorCode DMPlexPointIsTensor_Internal_Given(DM dm, PetscInt p, PetscInt f, PetscInt f2, PetscBool *isTensor)
1785: {
1786: PetscInt coneSize, c;
1787: const PetscInt *cone;
1788: const PetscInt *fCone;
1789: const PetscInt *f2Cone;
1790: PetscInt fs[2];
1791: PetscInt meetSize, nmeet;
1792: const PetscInt *meet;
1793: PetscErrorCode ierr;
1796: fs[0] = f;
1797: fs[1] = f2;
1798: DMPlexGetMeet(dm, 2, fs, &meetSize, &meet);
1799: nmeet = meetSize;
1800: DMPlexRestoreMeet(dm, 2, fs, &meetSize, &meet);
1801: /* two points that have a non-empty meet cannot be at opposite ends of a cell */
1802: if (nmeet) {
1803: *isTensor = PETSC_FALSE;
1804: return(0);
1805: }
1806: DMPlexGetConeSize(dm, p, &coneSize);
1807: DMPlexGetCone(dm, p, &cone);
1808: DMPlexGetCone(dm, f, &fCone);
1809: DMPlexGetCone(dm, f2, &f2Cone);
1810: for (c = 0; c < coneSize; c++) {
1811: PetscInt e, ef;
1812: PetscInt d = -1, d2 = -1;
1813: PetscInt dcount, d2count;
1814: PetscInt t = cone[c];
1815: PetscInt tConeSize;
1816: PetscBool tIsTensor;
1817: const PetscInt *tCone;
1819: if (t == f || t == f2) continue;
1820: /* for every other facet in the cone, check that is has
1821: * one ridge in common with each end */
1822: DMPlexGetConeSize(dm, t, &tConeSize);
1823: DMPlexGetCone(dm, t, &tCone);
1825: dcount = 0;
1826: d2count = 0;
1827: for (e = 0; e < tConeSize; e++) {
1828: PetscInt q = tCone[e];
1829: for (ef = 0; ef < coneSize - 2; ef++) {
1830: if (fCone[ef] == q) {
1831: if (dcount) {
1832: *isTensor = PETSC_FALSE;
1833: return(0);
1834: }
1835: d = q;
1836: dcount++;
1837: } else if (f2Cone[ef] == q) {
1838: if (d2count) {
1839: *isTensor = PETSC_FALSE;
1840: return(0);
1841: }
1842: d2 = q;
1843: d2count++;
1844: }
1845: }
1846: }
1847: /* if the whole cell is a tensor with the segment, then this
1848: * facet should be a tensor with the segment */
1849: DMPlexPointIsTensor_Internal_Given(dm, t, d, d2, &tIsTensor);
1850: if (!tIsTensor) {
1851: *isTensor = PETSC_FALSE;
1852: return(0);
1853: }
1854: }
1855: *isTensor = PETSC_TRUE;
1856: return(0);
1857: }
1859: /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair
1860: * that could be the opposite ends */
1861: static PetscErrorCode DMPlexPointIsTensor_Internal(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB)
1862: {
1863: PetscInt coneSize, c, c2;
1864: const PetscInt *cone;
1865: PetscErrorCode ierr;
1868: DMPlexGetConeSize(dm, p, &coneSize);
1869: if (!coneSize) {
1870: if (isTensor) *isTensor = PETSC_FALSE;
1871: if (endA) *endA = -1;
1872: if (endB) *endB = -1;
1873: }
1874: DMPlexGetCone(dm, p, &cone);
1875: for (c = 0; c < coneSize; c++) {
1876: PetscInt f = cone[c];
1877: PetscInt fConeSize;
1879: DMPlexGetConeSize(dm, f, &fConeSize);
1880: if (fConeSize != coneSize - 2) continue;
1882: for (c2 = c + 1; c2 < coneSize; c2++) {
1883: PetscInt f2 = cone[c2];
1884: PetscBool isTensorff2;
1885: PetscInt f2ConeSize;
1887: DMPlexGetConeSize(dm, f2, &f2ConeSize);
1888: if (f2ConeSize != coneSize - 2) continue;
1890: DMPlexPointIsTensor_Internal_Given(dm, p, f, f2, &isTensorff2);
1891: if (isTensorff2) {
1892: if (isTensor) *isTensor = PETSC_TRUE;
1893: if (endA) *endA = f;
1894: if (endB) *endB = f2;
1895: return(0);
1896: }
1897: }
1898: }
1899: if (isTensor) *isTensor = PETSC_FALSE;
1900: if (endA) *endA = -1;
1901: if (endB) *endB = -1;
1902: return(0);
1903: }
1905: /* determine if a cell is a tensor with a segment by looping over pairs of facets to find a pair
1906: * that could be the opposite ends */
1907: static PetscErrorCode DMPlexPointIsTensor(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB)
1908: {
1909: DMPlexInterpolatedFlag interpolated;
1913: DMPlexIsInterpolated(dm, &interpolated);
1914: if (interpolated != DMPLEX_INTERPOLATED_FULL) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Only for interpolated DMPlex's");
1915: DMPlexPointIsTensor_Internal(dm, p, isTensor, endA, endB);
1916: return(0);
1917: }
1919: /* permute a quadrature -> dof matrix so that its rows are in revlex order by nodeIdx */
1920: static PetscErrorCode MatPermuteByNodeIdx(Mat A, PetscLagNodeIndices ni, Mat *Aperm)
1921: {
1922: PetscInt m, n, i, j;
1923: PetscInt nodeIdxDim = ni->nodeIdxDim;
1924: PetscInt nodeVecDim = ni->nodeVecDim;
1925: PetscInt *perm;
1926: IS permIS;
1927: IS id;
1928: PetscInt *nIdxPerm;
1929: PetscReal *nVecPerm;
1933: PetscLagNodeIndicesGetPermutation(ni, &perm);
1934: MatGetSize(A, &m, &n);
1935: PetscMalloc1(nodeIdxDim * m, &nIdxPerm);
1936: PetscMalloc1(nodeVecDim * m, &nVecPerm);
1937: for (i = 0; i < m; i++) for (j = 0; j < nodeIdxDim; j++) nIdxPerm[i * nodeIdxDim + j] = ni->nodeIdx[perm[i] * nodeIdxDim + j];
1938: for (i = 0; i < m; i++) for (j = 0; j < nodeVecDim; j++) nVecPerm[i * nodeVecDim + j] = ni->nodeVec[perm[i] * nodeVecDim + j];
1939: ISCreateGeneral(PETSC_COMM_SELF, m, perm, PETSC_USE_POINTER, &permIS);
1940: ISSetPermutation(permIS);
1941: ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &id);
1942: ISSetPermutation(id);
1943: MatPermute(A, permIS, id, Aperm);
1944: ISDestroy(&permIS);
1945: ISDestroy(&id);
1946: for (i = 0; i < m; i++) perm[i] = i;
1947: PetscFree(ni->nodeIdx);
1948: PetscFree(ni->nodeVec);
1949: ni->nodeIdx = nIdxPerm;
1950: ni->nodeVec = nVecPerm;
1951: return(0);
1952: }
1954: static PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1955: {
1956: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1957: DM dm = sp->dm;
1958: DM dmint = NULL;
1959: PetscInt order;
1960: PetscInt Nc = sp->Nc;
1961: MPI_Comm comm;
1962: PetscBool continuous;
1963: PetscSection section;
1964: PetscInt depth, dim, pStart, pEnd, cStart, cEnd, p, *pStratStart, *pStratEnd, d;
1965: PetscInt formDegree, Nk, Ncopies;
1966: PetscInt tensorf = -1, tensorf2 = -1;
1967: PetscBool tensorCell, tensorSpace;
1968: PetscBool uniform, trimmed;
1969: Petsc1DNodeFamily nodeFamily;
1970: PetscInt numNodeSkip;
1971: DMPlexInterpolatedFlag interpolated;
1972: PetscBool isbdm;
1973: PetscErrorCode ierr;
1976: /* step 1: sanitize input */
1977: PetscObjectGetComm((PetscObject) sp, &comm);
1978: DMGetDimension(dm, &dim);
1979: PetscObjectTypeCompare((PetscObject)sp, "bdm", &isbdm);
1980: if (isbdm) {
1981: sp->k = -(dim-1); /* form degree of H-div */
1982: PetscObjectChangeTypeName((PetscObject)sp, PETSCDUALSPACELAGRANGE);
1983: }
1984: PetscDualSpaceGetFormDegree(sp, &formDegree);
1985: if (PetscAbsInt(formDegree) > dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Form degree must be bounded by dimension");
1986: PetscDTBinomialInt(dim,PetscAbsInt(formDegree),&Nk);
1987: if (sp->Nc <= 0 && lag->numCopies > 0) sp->Nc = Nk * lag->numCopies;
1988: Nc = sp->Nc;
1989: if (Nc % Nk) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "Number of components is not a multiple of form degree size");
1990: if (lag->numCopies <= 0) lag->numCopies = Nc / Nk;
1991: Ncopies = lag->numCopies;
1992: if (Nc / Nk != Ncopies) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "Number of copies * (dim choose k) != Nc");
1993: if (!dim) sp->order = 0;
1994: order = sp->order;
1995: uniform = sp->uniform;
1996: if (!uniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Variable order not supported yet");
1997: if (lag->trimmed && !formDegree) lag->trimmed = PETSC_FALSE; /* trimmed spaces are the same as full spaces for 0-forms */
1998: if (lag->nodeType == PETSCDTNODES_DEFAULT) {
1999: lag->nodeType = PETSCDTNODES_GAUSSJACOBI;
2000: lag->nodeExponent = 0.;
2001: /* trimmed spaces don't include corner vertices, so don't use end nodes by default */
2002: lag->endNodes = lag->trimmed ? PETSC_FALSE : PETSC_TRUE;
2003: }
2004: /* If a trimmed space and the user did choose nodes with endpoints, skip them by default */
2005: if (lag->numNodeSkip < 0) lag->numNodeSkip = (lag->trimmed && lag->endNodes) ? 1 : 0;
2006: numNodeSkip = lag->numNodeSkip;
2007: if (lag->trimmed && !order) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zeroth order trimmed elements");
2008: if (lag->trimmed && PetscAbsInt(formDegree) == dim) { /* convert trimmed n-forms to untrimmed of one polynomial order less */
2009: sp->order--;
2010: order--;
2011: lag->trimmed = PETSC_FALSE;
2012: }
2013: trimmed = lag->trimmed;
2014: if (!order || PetscAbsInt(formDegree) == dim) lag->continuous = PETSC_FALSE;
2015: continuous = lag->continuous;
2016: DMPlexGetDepth(dm, &depth);
2017: DMPlexGetChart(dm, &pStart, &pEnd);
2018: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2019: if (pStart != 0 || cStart != 0) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Expect DM with chart starting at zero and cells first");
2020: if (cEnd != 1) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Use PETSCDUALSPACEREFINED for multi-cell reference meshes");
2021: DMPlexIsInterpolated(dm, &interpolated);
2022: if (interpolated != DMPLEX_INTERPOLATED_FULL) {
2023: DMPlexInterpolate(dm, &dmint);
2024: } else {
2025: PetscObjectReference((PetscObject)dm);
2026: dmint = dm;
2027: }
2028: tensorCell = PETSC_FALSE;
2029: if (dim > 1) {
2030: DMPlexPointIsTensor(dmint, 0, &tensorCell, &tensorf, &tensorf2);
2031: }
2032: lag->tensorCell = tensorCell;
2033: if (dim < 2 || !lag->tensorCell) lag->tensorSpace = PETSC_FALSE;
2034: tensorSpace = lag->tensorSpace;
2035: if (!lag->nodeFamily) {
2036: Petsc1DNodeFamilyCreate(lag->nodeType, lag->nodeExponent, lag->endNodes, &lag->nodeFamily);
2037: }
2038: nodeFamily = lag->nodeFamily;
2039: if (interpolated != DMPLEX_INTERPOLATED_FULL && continuous && (PetscAbsInt(formDegree) > 0 || order > 1)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Reference element won't support all boundary nodes");
2041: /* step 2: construct the boundary spaces */
2042: PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);
2043: PetscCalloc1(pEnd,&(sp->pointSpaces));
2044: for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);}
2045: PetscDualSpaceSectionCreate_Internal(sp, §ion);
2046: sp->pointSection = section;
2047: if (continuous && !(lag->interiorOnly)) {
2048: PetscInt h;
2050: for (p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */
2051: PetscReal v0[3];
2052: DMPolytopeType ptype;
2053: PetscReal J[9], detJ;
2054: PetscInt q;
2056: DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, NULL, &detJ);
2057: DMPlexGetCellType(dm, p, &ptype);
2059: /* compare to previous facets: if computed, reference that dualspace */
2060: for (q = pStratStart[depth - 1]; q < p; q++) {
2061: DMPolytopeType qtype;
2063: DMPlexGetCellType(dm, q, &qtype);
2064: if (qtype == ptype) break;
2065: }
2066: if (q < p) { /* this facet has the same dual space as that one */
2067: PetscObjectReference((PetscObject)sp->pointSpaces[q]);
2068: sp->pointSpaces[p] = sp->pointSpaces[q];
2069: continue;
2070: }
2071: /* if not, recursively compute this dual space */
2072: PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,p,formDegree,Ncopies,PETSC_FALSE,&sp->pointSpaces[p]);
2073: }
2074: for (h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */
2075: PetscInt hd = depth - h;
2076: PetscInt hdim = dim - h;
2078: if (hdim < PetscAbsInt(formDegree)) break;
2079: for (p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
2080: PetscInt suppSize, s;
2081: const PetscInt *supp;
2083: DMPlexGetSupportSize(dm, p, &suppSize);
2084: DMPlexGetSupport(dm, p, &supp);
2085: for (s = 0; s < suppSize; s++) {
2086: DM qdm;
2087: PetscDualSpace qsp, psp;
2088: PetscInt c, coneSize, q;
2089: const PetscInt *cone;
2090: const PetscInt *refCone;
2092: q = supp[0];
2093: qsp = sp->pointSpaces[q];
2094: DMPlexGetConeSize(dm, q, &coneSize);
2095: DMPlexGetCone(dm, q, &cone);
2096: for (c = 0; c < coneSize; c++) if (cone[c] == p) break;
2097: if (c == coneSize) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "cone/suppport mismatch");
2098: PetscDualSpaceGetDM(qsp, &qdm);
2099: DMPlexGetCone(qdm, 0, &refCone);
2100: /* get the equivalent dual space from the support dual space */
2101: PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp);
2102: if (!s) {
2103: PetscObjectReference((PetscObject)psp);
2104: sp->pointSpaces[p] = psp;
2105: }
2106: }
2107: }
2108: }
2109: for (p = 1; p < pEnd; p++) {
2110: PetscInt pspdim;
2111: if (!sp->pointSpaces[p]) continue;
2112: PetscDualSpaceGetInteriorDimension(sp->pointSpaces[p], &pspdim);
2113: PetscSectionSetDof(section, p, pspdim);
2114: }
2115: }
2117: if (Ncopies > 1) {
2118: Mat intMatScalar, allMatScalar;
2119: PetscDualSpace scalarsp;
2120: PetscDualSpace_Lag *scalarlag;
2122: PetscDualSpaceDuplicate(sp, &scalarsp);
2123: /* Setting the number of components to Nk is a space with 1 copy of each k-form */
2124: PetscDualSpaceSetNumComponents(scalarsp, Nk);
2125: PetscDualSpaceSetUp(scalarsp);
2126: PetscDualSpaceGetInteriorData(scalarsp, &(sp->intNodes), &intMatScalar);
2127: PetscObjectReference((PetscObject)(sp->intNodes));
2128: if (intMatScalar) {PetscDualSpaceLagrangeMatrixCreateCopies(intMatScalar, Nk, Ncopies, &(sp->intMat));}
2129: PetscDualSpaceGetAllData(scalarsp, &(sp->allNodes), &allMatScalar);
2130: PetscObjectReference((PetscObject)(sp->allNodes));
2131: PetscDualSpaceLagrangeMatrixCreateCopies(allMatScalar, Nk, Ncopies, &(sp->allMat));
2132: sp->spdim = scalarsp->spdim * Ncopies;
2133: sp->spintdim = scalarsp->spintdim * Ncopies;
2134: scalarlag = (PetscDualSpace_Lag *) scalarsp->data;
2135: PetscLagNodeIndicesReference(scalarlag->vertIndices);
2136: lag->vertIndices = scalarlag->vertIndices;
2137: PetscLagNodeIndicesReference(scalarlag->intNodeIndices);
2138: lag->intNodeIndices = scalarlag->intNodeIndices;
2139: PetscLagNodeIndicesReference(scalarlag->allNodeIndices);
2140: lag->allNodeIndices = scalarlag->allNodeIndices;
2141: PetscDualSpaceDestroy(&scalarsp);
2142: PetscSectionSetDof(section, 0, sp->spintdim);
2143: PetscDualSpaceSectionSetUp_Internal(sp, section);
2144: PetscDualSpaceComputeFunctionalsFromAllData(sp);
2145: PetscFree2(pStratStart, pStratEnd);
2146: DMDestroy(&dmint);
2147: return(0);
2148: }
2150: if (trimmed && !continuous) {
2151: /* the dofs of a trimmed space don't have a nice tensor/lattice structure:
2152: * just construct the continuous dual space and copy all of the data over,
2153: * allocating it all to the cell instead of splitting it up between the boundaries */
2154: PetscDualSpace spcont;
2155: PetscInt spdim, f;
2156: PetscQuadrature allNodes;
2157: PetscDualSpace_Lag *lagc;
2158: Mat allMat;
2160: PetscDualSpaceDuplicate(sp, &spcont);
2161: PetscDualSpaceLagrangeSetContinuity(spcont, PETSC_TRUE);
2162: PetscDualSpaceSetUp(spcont);
2163: PetscDualSpaceGetDimension(spcont, &spdim);
2164: sp->spdim = sp->spintdim = spdim;
2165: PetscSectionSetDof(section, 0, spdim);
2166: PetscDualSpaceSectionSetUp_Internal(sp, section);
2167: PetscMalloc1(spdim, &(sp->functional));
2168: for (f = 0; f < spdim; f++) {
2169: PetscQuadrature fn;
2171: PetscDualSpaceGetFunctional(spcont, f, &fn);
2172: PetscObjectReference((PetscObject)fn);
2173: sp->functional[f] = fn;
2174: }
2175: PetscDualSpaceGetAllData(spcont, &allNodes, &allMat);
2176: PetscObjectReference((PetscObject) allNodes);
2177: PetscObjectReference((PetscObject) allNodes);
2178: sp->allNodes = sp->intNodes = allNodes;
2179: PetscObjectReference((PetscObject) allMat);
2180: PetscObjectReference((PetscObject) allMat);
2181: sp->allMat = sp->intMat = allMat;
2182: lagc = (PetscDualSpace_Lag *) spcont->data;
2183: PetscLagNodeIndicesReference(lagc->vertIndices);
2184: lag->vertIndices = lagc->vertIndices;
2185: PetscLagNodeIndicesReference(lagc->allNodeIndices);
2186: PetscLagNodeIndicesReference(lagc->allNodeIndices);
2187: lag->intNodeIndices = lagc->allNodeIndices;
2188: lag->allNodeIndices = lagc->allNodeIndices;
2189: PetscDualSpaceDestroy(&spcont);
2190: PetscFree2(pStratStart, pStratEnd);
2191: DMDestroy(&dmint);
2192: return(0);
2193: }
2195: /* step 3: construct intNodes, and intMat, and combine it with boundray data to make allNodes and allMat */
2196: if (!tensorSpace) {
2197: if (!tensorCell) {PetscLagNodeIndicesCreateSimplexVertices(dm, &(lag->vertIndices));}
2199: if (trimmed) {
2200: /* there is one dof in the interior of the a trimmed element for each full polynomial of with degree at most
2201: * order + k - dim - 1 */
2202: if (order + PetscAbsInt(formDegree) > dim) {
2203: PetscInt sum = order + PetscAbsInt(formDegree) - dim - 1;
2204: PetscInt nDofs;
2206: PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));
2207: MatGetSize(sp->intMat, &nDofs, NULL);
2208: PetscSectionSetDof(section, 0, nDofs);
2209: }
2210: PetscDualSpaceSectionSetUp_Internal(sp, section);
2211: PetscDualSpaceCreateAllDataFromInteriorData(sp);
2212: PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2213: } else {
2214: if (!continuous) {
2215: /* if discontinuous just construct one node for each set of dofs (a set of dofs is a basis for the k-form
2216: * space) */
2217: PetscInt sum = order;
2218: PetscInt nDofs;
2220: PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));
2221: MatGetSize(sp->intMat, &nDofs, NULL);
2222: PetscSectionSetDof(section, 0, nDofs);
2223: PetscDualSpaceSectionSetUp_Internal(sp, section);
2224: PetscObjectReference((PetscObject)(sp->intNodes));
2225: sp->allNodes = sp->intNodes;
2226: PetscObjectReference((PetscObject)(sp->intMat));
2227: sp->allMat = sp->intMat;
2228: PetscLagNodeIndicesReference(lag->intNodeIndices);
2229: lag->allNodeIndices = lag->intNodeIndices;
2230: } else {
2231: /* there is one dof in the interior of the a full element for each trimmed polynomial of with degree at most
2232: * order + k - dim, but with complementary form degree */
2233: if (order + PetscAbsInt(formDegree) > dim) {
2234: PetscDualSpace trimmedsp;
2235: PetscDualSpace_Lag *trimmedlag;
2236: PetscQuadrature intNodes;
2237: PetscInt trFormDegree = formDegree >= 0 ? formDegree - dim : dim - PetscAbsInt(formDegree);
2238: PetscInt nDofs;
2239: Mat intMat;
2241: PetscDualSpaceDuplicate(sp, &trimmedsp);
2242: PetscDualSpaceLagrangeSetTrimmed(trimmedsp, PETSC_TRUE);
2243: PetscDualSpaceSetOrder(trimmedsp, order + PetscAbsInt(formDegree) - dim);
2244: PetscDualSpaceSetFormDegree(trimmedsp, trFormDegree);
2245: trimmedlag = (PetscDualSpace_Lag *) trimmedsp->data;
2246: trimmedlag->numNodeSkip = numNodeSkip + 1;
2247: PetscDualSpaceSetUp(trimmedsp);
2248: PetscDualSpaceGetAllData(trimmedsp, &intNodes, &intMat);
2249: PetscObjectReference((PetscObject)intNodes);
2250: sp->intNodes = intNodes;
2251: PetscObjectReference((PetscObject)intMat);
2252: sp->intMat = intMat;
2253: MatGetSize(sp->intMat, &nDofs, NULL);
2254: PetscLagNodeIndicesReference(trimmedlag->allNodeIndices);
2255: lag->intNodeIndices = trimmedlag->allNodeIndices;
2256: PetscDualSpaceDestroy(&trimmedsp);
2257: PetscSectionSetDof(section, 0, nDofs);
2258: }
2259: PetscDualSpaceSectionSetUp_Internal(sp, section);
2260: PetscDualSpaceCreateAllDataFromInteriorData(sp);
2261: PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2262: }
2263: }
2264: } else {
2265: PetscQuadrature intNodesTrace = NULL;
2266: PetscQuadrature intNodesFiber = NULL;
2267: PetscQuadrature intNodes = NULL;
2268: PetscLagNodeIndices intNodeIndices = NULL;
2269: Mat intMat = NULL;
2271: if (PetscAbsInt(formDegree) < dim) { /* get the trace k-forms on the first facet, and the 0-forms on the edge,
2272: and wedge them together to create some of the k-form dofs */
2273: PetscDualSpace trace, fiber;
2274: PetscDualSpace_Lag *tracel, *fiberl;
2275: Mat intMatTrace, intMatFiber;
2277: if (sp->pointSpaces[tensorf]) {
2278: PetscObjectReference((PetscObject)(sp->pointSpaces[tensorf]));
2279: trace = sp->pointSpaces[tensorf];
2280: } else {
2281: PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,formDegree,Ncopies,PETSC_TRUE,&trace);
2282: }
2283: PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,0,1,PETSC_TRUE,&fiber);
2284: tracel = (PetscDualSpace_Lag *) trace->data;
2285: fiberl = (PetscDualSpace_Lag *) fiber->data;
2286: PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));
2287: PetscDualSpaceGetInteriorData(trace, &intNodesTrace, &intMatTrace);
2288: PetscDualSpaceGetInteriorData(fiber, &intNodesFiber, &intMatFiber);
2289: if (intNodesTrace && intNodesFiber) {
2290: PetscQuadratureCreateTensor(intNodesTrace, intNodesFiber, &intNodes);
2291: MatTensorAltV(intMatTrace, intMatFiber, dim-1, formDegree, 1, 0, &intMat);
2292: PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, formDegree, fiberl->intNodeIndices, 1, 0, &intNodeIndices);
2293: }
2294: PetscObjectReference((PetscObject) intNodesTrace);
2295: PetscObjectReference((PetscObject) intNodesFiber);
2296: PetscDualSpaceDestroy(&fiber);
2297: PetscDualSpaceDestroy(&trace);
2298: }
2299: if (PetscAbsInt(formDegree) > 0) { /* get the trace (k-1)-forms on the first facet, and the 1-forms on the edge,
2300: and wedge them together to create the remaining k-form dofs */
2301: PetscDualSpace trace, fiber;
2302: PetscDualSpace_Lag *tracel, *fiberl;
2303: PetscQuadrature intNodesTrace2, intNodesFiber2, intNodes2;
2304: PetscLagNodeIndices intNodeIndices2;
2305: Mat intMatTrace, intMatFiber, intMat2;
2306: PetscInt traceDegree = formDegree > 0 ? formDegree - 1 : formDegree + 1;
2307: PetscInt fiberDegree = formDegree > 0 ? 1 : -1;
2309: PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,traceDegree,Ncopies,PETSC_TRUE,&trace);
2310: PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,fiberDegree,1,PETSC_TRUE,&fiber);
2311: tracel = (PetscDualSpace_Lag *) trace->data;
2312: fiberl = (PetscDualSpace_Lag *) fiber->data;
2313: if (!lag->vertIndices) {
2314: PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));
2315: }
2316: PetscDualSpaceGetInteriorData(trace, &intNodesTrace2, &intMatTrace);
2317: PetscDualSpaceGetInteriorData(fiber, &intNodesFiber2, &intMatFiber);
2318: if (intNodesTrace2 && intNodesFiber2) {
2319: PetscQuadratureCreateTensor(intNodesTrace2, intNodesFiber2, &intNodes2);
2320: MatTensorAltV(intMatTrace, intMatFiber, dim-1, traceDegree, 1, fiberDegree, &intMat2);
2321: PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, traceDegree, fiberl->intNodeIndices, 1, fiberDegree, &intNodeIndices2);
2322: if (!intMat) {
2323: intMat = intMat2;
2324: intNodes = intNodes2;
2325: intNodeIndices = intNodeIndices2;
2326: } else {
2327: /* merge the matrices, quadrature points, and nodes */
2328: PetscInt nM;
2329: PetscInt nDof, nDof2;
2330: PetscInt *toMerged = NULL, *toMerged2 = NULL;
2331: PetscQuadrature merged = NULL;
2332: PetscLagNodeIndices intNodeIndicesMerged = NULL;
2333: Mat matMerged = NULL;
2335: MatGetSize(intMat, &nDof, 0);
2336: MatGetSize(intMat2, &nDof2, 0);
2337: PetscQuadraturePointsMerge(intNodes, intNodes2, &merged, &toMerged, &toMerged2);
2338: PetscQuadratureGetData(merged, NULL, NULL, &nM, NULL, NULL);
2339: MatricesMerge(intMat, intMat2, dim, formDegree, nM, toMerged, toMerged2, &matMerged);
2340: PetscLagNodeIndicesMerge(intNodeIndices, intNodeIndices2, &intNodeIndicesMerged);
2341: PetscFree(toMerged);
2342: PetscFree(toMerged2);
2343: MatDestroy(&intMat);
2344: MatDestroy(&intMat2);
2345: PetscQuadratureDestroy(&intNodes);
2346: PetscQuadratureDestroy(&intNodes2);
2347: PetscLagNodeIndicesDestroy(&intNodeIndices);
2348: PetscLagNodeIndicesDestroy(&intNodeIndices2);
2349: intNodes = merged;
2350: intMat = matMerged;
2351: intNodeIndices = intNodeIndicesMerged;
2352: if (!trimmed) {
2353: /* I think users expect that, when a node has a full basis for the k-forms,
2354: * they should be consecutive dofs. That isn't the case for trimmed spaces,
2355: * but is for some of the nodes in untrimmed spaces, so in that case we
2356: * sort them to group them by node */
2357: Mat intMatPerm;
2359: MatPermuteByNodeIdx(intMat, intNodeIndices, &intMatPerm);
2360: MatDestroy(&intMat);
2361: intMat = intMatPerm;
2362: }
2363: }
2364: }
2365: PetscDualSpaceDestroy(&fiber);
2366: PetscDualSpaceDestroy(&trace);
2367: }
2368: PetscQuadratureDestroy(&intNodesTrace);
2369: PetscQuadratureDestroy(&intNodesFiber);
2370: sp->intNodes = intNodes;
2371: sp->intMat = intMat;
2372: lag->intNodeIndices = intNodeIndices;
2373: {
2374: PetscInt nDofs = 0;
2376: if (intMat) {
2377: MatGetSize(intMat, &nDofs, NULL);
2378: }
2379: PetscSectionSetDof(section, 0, nDofs);
2380: }
2381: PetscDualSpaceSectionSetUp_Internal(sp, section);
2382: if (continuous) {
2383: PetscDualSpaceCreateAllDataFromInteriorData(sp);
2384: PetscDualSpaceLagrangeCreateAllNodeIdx(sp);
2385: } else {
2386: PetscObjectReference((PetscObject) intNodes);
2387: sp->allNodes = intNodes;
2388: PetscObjectReference((PetscObject) intMat);
2389: sp->allMat = intMat;
2390: PetscLagNodeIndicesReference(intNodeIndices);
2391: lag->allNodeIndices = intNodeIndices;
2392: }
2393: }
2394: PetscSectionGetStorageSize(section, &sp->spdim);
2395: PetscSectionGetConstrainedStorageSize(section, &sp->spintdim);
2396: PetscDualSpaceComputeFunctionalsFromAllData(sp);
2397: PetscFree2(pStratStart, pStratEnd);
2398: DMDestroy(&dmint);
2399: return(0);
2400: }
2402: /* Create a matrix that represents the transformation that DMPlexVecGetClosure() would need
2403: * to get the representation of the dofs for a mesh point if the mesh point had this orientation
2404: * relative to the cell */
2405: PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat)
2406: {
2407: PetscDualSpace_Lag *lag;
2408: DM dm;
2409: PetscLagNodeIndices vertIndices, intNodeIndices;
2410: PetscLagNodeIndices ni;
2411: PetscInt nodeIdxDim, nodeVecDim, nNodes;
2412: PetscInt formDegree;
2413: PetscInt *perm, *permOrnt;
2414: PetscInt *nnz;
2415: PetscInt n;
2416: PetscInt maxGroupSize;
2417: PetscScalar *V, *W, *work;
2418: Mat A;
2422: if (!sp->spintdim) {
2423: *symMat = NULL;
2424: return(0);
2425: }
2426: lag = (PetscDualSpace_Lag *) sp->data;
2427: vertIndices = lag->vertIndices;
2428: intNodeIndices = lag->intNodeIndices;
2429: PetscDualSpaceGetDM(sp, &dm);
2430: PetscDualSpaceGetFormDegree(sp, &formDegree);
2431: PetscNew(&ni);
2432: ni->refct = 1;
2433: ni->nodeIdxDim = nodeIdxDim = intNodeIndices->nodeIdxDim;
2434: ni->nodeVecDim = nodeVecDim = intNodeIndices->nodeVecDim;
2435: ni->nNodes = nNodes = intNodeIndices->nNodes;
2436: PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));
2437: PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));
2438: /* push forward the dofs by the symmetry of the reference element induced by ornt */
2439: PetscLagNodeIndicesPushForward(dm, vertIndices, 0, vertIndices, intNodeIndices, ornt, formDegree, ni->nodeIdx, ni->nodeVec);
2440: /* get the revlex order for both the original and transformed dofs */
2441: PetscLagNodeIndicesGetPermutation(intNodeIndices, &perm);
2442: PetscLagNodeIndicesGetPermutation(ni, &permOrnt);
2443: PetscMalloc1(nNodes, &nnz);
2444: for (n = 0, maxGroupSize = 0; n < nNodes;) { /* incremented in the loop */
2445: PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2446: PetscInt m, nEnd;
2447: PetscInt groupSize;
2448: /* for each group of dofs that have the same nodeIdx coordinate */
2449: for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2450: PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2451: PetscInt d;
2453: /* compare the oriented permutation indices */
2454: for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break;
2455: if (d < nodeIdxDim) break;
2456: }
2457: /* permOrnt[[n, nEnd)] is a group of dofs that, under the symmetry are at the same location */
2458: #if defined(PETSC_USE_DEBUG)
2459: /* the symmetry had better map the group of dofs with the same permuted nodeIdx
2460: * to a group of dofs with the same size, otherwise we messed up */
2461: {
2462: PetscInt m;
2463: PetscInt *nind = &(intNodeIndices->nodeIdx[perm[n] * nodeIdxDim]);
2465: for (m = n + 1; m < nEnd; m++) {
2466: PetscInt *mind = &(intNodeIndices->nodeIdx[perm[m] * nodeIdxDim]);
2467: PetscInt d;
2469: /* compare the oriented permutation indices */
2470: for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break;
2471: if (d < nodeIdxDim) break;
2472: }
2473: if (m < nEnd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs with same index after symmetry not same block size");
2474: }
2475: #endif
2476: groupSize = nEnd - n;
2477: /* each pushforward dof vector will be expressed in a basis of the unpermuted dofs */
2478: for (m = n; m < nEnd; m++) nnz[permOrnt[m]] = groupSize;
2480: maxGroupSize = PetscMax(maxGroupSize, nEnd - n);
2481: n = nEnd;
2482: }
2483: if (maxGroupSize > nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs are not in blocks that can be solved");
2484: MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes, nNodes, 0, nnz, &A);
2485: PetscFree(nnz);
2486: PetscMalloc3(maxGroupSize * nodeVecDim, &V, maxGroupSize * nodeVecDim, &W, nodeVecDim * 2, &work);
2487: for (n = 0; n < nNodes;) { /* incremented in the loop */
2488: PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2489: PetscInt nEnd;
2490: PetscInt m;
2491: PetscInt groupSize;
2492: for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2493: PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2494: PetscInt d;
2496: /* compare the oriented permutation indices */
2497: for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break;
2498: if (d < nodeIdxDim) break;
2499: }
2500: groupSize = nEnd - n;
2501: /* get all of the vectors from the original and all of the pushforward vectors */
2502: for (m = n; m < nEnd; m++) {
2503: PetscInt d;
2505: for (d = 0; d < nodeVecDim; d++) {
2506: V[(m - n) * nodeVecDim + d] = intNodeIndices->nodeVec[perm[m] * nodeVecDim + d];
2507: W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
2508: }
2509: }
2510: /* now we have to solve for W in terms of V: the systems isn't always square, but the span
2511: * of V and W should always be the same, so the solution of the normal equations works */
2512: {
2513: char transpose = 'N';
2514: PetscBLASInt bm = nodeVecDim;
2515: PetscBLASInt bn = groupSize;
2516: PetscBLASInt bnrhs = groupSize;
2517: PetscBLASInt blda = bm;
2518: PetscBLASInt bldb = bm;
2519: PetscBLASInt blwork = 2 * nodeVecDim;
2520: PetscBLASInt info;
2522: PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&bm,&bn,&bnrhs,V,&blda,W,&bldb,work,&blwork, &info));
2523: if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS");
2524: /* repack */
2525: {
2526: PetscInt i, j;
2528: for (i = 0; i < groupSize; i++) {
2529: for (j = 0; j < groupSize; j++) {
2530: /* notice the different leading dimension */
2531: V[i * groupSize + j] = W[i * nodeVecDim + j];
2532: }
2533: }
2534: }
2535: }
2536: MatSetValues(A, groupSize, &permOrnt[n], groupSize, &perm[n], V, INSERT_VALUES);
2537: n = nEnd;
2538: }
2539: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
2540: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
2541: *symMat = A;
2542: PetscFree3(V,W,work);
2543: PetscLagNodeIndicesDestroy(&ni);
2544: return(0);
2545: }
2547: #define BaryIndex(perEdge,a,b,c) (((b)*(2*perEdge+1-(b)))/2)+(c)
2549: #define CartIndex(perEdge,a,b) (perEdge*(a)+b)
2551: /* the existing interface for symmetries is insufficient for all cases:
2552: * - it should be sufficient for form degrees that are scalar (0 and n)
2553: * - it should be sufficient for hypercube dofs
2554: * - it isn't sufficient for simplex cells with non-scalar form degrees if
2555: * there are any dofs in the interior
2556: *
2557: * We compute the general transformation matrices, and if they fit, we return them,
2558: * otherwise we error (but we should probably change the interface to allow for
2559: * these symmetries)
2560: */
2561: static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
2562: {
2563: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2564: PetscInt dim, order, Nc;
2565: PetscErrorCode ierr;
2568: PetscDualSpaceGetOrder(sp,&order);
2569: PetscDualSpaceGetNumComponents(sp,&Nc);
2570: DMGetDimension(sp->dm,&dim);
2571: if (!lag->symComputed) { /* store symmetries */
2572: PetscInt pStart, pEnd, p;
2573: PetscInt numPoints;
2574: PetscInt numFaces;
2575: PetscInt spintdim;
2576: PetscInt ***symperms;
2577: PetscScalar ***symflips;
2579: DMPlexGetChart(sp->dm, &pStart, &pEnd);
2580: numPoints = pEnd - pStart;
2581: DMPlexGetConeSize(sp->dm, 0, &numFaces);
2582: PetscCalloc1(numPoints,&symperms);
2583: PetscCalloc1(numPoints,&symflips);
2584: spintdim = sp->spintdim;
2585: /* The nodal symmetry behavior is not present when tensorSpace != tensorCell: someone might want this for the "S"
2586: * family of FEEC spaces. Most used in particular are discontinuous polynomial L2 spaces in tensor cells, where
2587: * the symmetries are not necessary for FE assembly. So for now we assume this is the case and don't return
2588: * symmetries if tensorSpace != tensorCell */
2589: if (spintdim && 0 < dim && dim < 3 && (lag->tensorSpace == lag->tensorCell)) { /* compute self symmetries */
2590: PetscInt **cellSymperms;
2591: PetscScalar **cellSymflips;
2592: PetscInt ornt;
2593: PetscInt nCopies = Nc / lag->intNodeIndices->nodeVecDim;
2594: PetscInt nNodes = lag->intNodeIndices->nNodes;
2596: lag->numSelfSym = 2 * numFaces;
2597: lag->selfSymOff = numFaces;
2598: PetscCalloc1(2*numFaces,&cellSymperms);
2599: PetscCalloc1(2*numFaces,&cellSymflips);
2600: /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
2601: symperms[0] = &cellSymperms[numFaces];
2602: symflips[0] = &cellSymflips[numFaces];
2603: if (lag->intNodeIndices->nodeVecDim * nCopies != Nc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs");
2604: if (nNodes * nCopies != spintdim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs");
2605: 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 */
2606: Mat symMat;
2607: PetscInt *perm;
2608: PetscScalar *flips;
2609: PetscInt i;
2611: if (!ornt) continue;
2612: PetscMalloc1(spintdim, &perm);
2613: PetscCalloc1(spintdim, &flips);
2614: for (i = 0; i < spintdim; i++) perm[i] = -1;
2615: PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, ornt, &symMat);
2616: for (i = 0; i < nNodes; i++) {
2617: PetscInt ncols;
2618: PetscInt j, k;
2619: const PetscInt *cols;
2620: const PetscScalar *vals;
2621: PetscBool nz_seen = PETSC_FALSE;
2623: MatGetRow(symMat, i, &ncols, &cols, &vals);
2624: for (j = 0; j < ncols; j++) {
2625: if (PetscAbsScalar(vals[j]) > PETSC_SMALL) {
2626: if (nz_seen) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2627: nz_seen = PETSC_TRUE;
2628: if (PetscAbsReal(PetscAbsScalar(vals[j]) - PetscRealConstant(1.)) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2629: if (PetscAbsReal(PetscImaginaryPart(vals[j])) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2630: if (perm[cols[j] * nCopies] >= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips");
2631: for (k = 0; k < nCopies; k++) {
2632: perm[cols[j] * nCopies + k] = i * nCopies + k;
2633: }
2634: if (PetscRealPart(vals[j]) < 0.) {
2635: for (k = 0; k < nCopies; k++) {
2636: flips[i * nCopies + k] = -1.;
2637: }
2638: } else {
2639: for (k = 0; k < nCopies; k++) {
2640: flips[i * nCopies + k] = 1.;
2641: }
2642: }
2643: }
2644: }
2645: MatRestoreRow(symMat, i, &ncols, &cols, &vals);
2646: }
2647: MatDestroy(&symMat);
2648: /* if there were no sign flips, keep NULL */
2649: for (i = 0; i < spintdim; i++) if (flips[i] != 1.) break;
2650: if (i == spintdim) {
2651: PetscFree(flips);
2652: flips = NULL;
2653: }
2654: /* if the permutation is identity, keep NULL */
2655: for (i = 0; i < spintdim; i++) if (perm[i] != i) break;
2656: if (i == spintdim) {
2657: PetscFree(perm);
2658: perm = NULL;
2659: }
2660: symperms[0][ornt] = perm;
2661: symflips[0][ornt] = flips;
2662: }
2663: /* if no orientations produced non-identity permutations, keep NULL */
2664: for (ornt = -numFaces; ornt < numFaces; ornt++) if (symperms[0][ornt]) break;
2665: if (ornt == numFaces) {
2666: PetscFree(cellSymperms);
2667: symperms[0] = NULL;
2668: }
2669: /* if no orientations produced sign flips, keep NULL */
2670: for (ornt = -numFaces; ornt < numFaces; ornt++) if (symflips[0][ornt]) break;
2671: if (ornt == numFaces) {
2672: PetscFree(cellSymflips);
2673: symflips[0] = NULL;
2674: }
2675: }
2676: { /* get the symmetries of closure points */
2677: PetscInt closureSize = 0;
2678: PetscInt *closure = NULL;
2679: PetscInt r;
2681: DMPlexGetTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);
2682: for (r = 0; r < closureSize; r++) {
2683: PetscDualSpace psp;
2684: PetscInt point = closure[2 * r];
2685: PetscInt pspintdim;
2686: const PetscInt ***psymperms = NULL;
2687: const PetscScalar ***psymflips = NULL;
2689: if (!point) continue;
2690: PetscDualSpaceGetPointSubspace(sp, point, &psp);
2691: if (!psp) continue;
2692: PetscDualSpaceGetInteriorDimension(psp, &pspintdim);
2693: if (!pspintdim) continue;
2694: PetscDualSpaceGetSymmetries(psp,&psymperms,&psymflips);
2695: symperms[r] = (PetscInt **) (psymperms ? psymperms[0] : NULL);
2696: symflips[r] = (PetscScalar **) (psymflips ? psymflips[0] : NULL);
2697: }
2698: DMPlexRestoreTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);
2699: }
2700: for (p = 0; p < pEnd; p++) if (symperms[p]) break;
2701: if (p == pEnd) {
2702: PetscFree(symperms);
2703: symperms = NULL;
2704: }
2705: for (p = 0; p < pEnd; p++) if (symflips[p]) break;
2706: if (p == pEnd) {
2707: PetscFree(symflips);
2708: symflips = NULL;
2709: }
2710: lag->symperms = symperms;
2711: lag->symflips = symflips;
2712: lag->symComputed = PETSC_TRUE;
2713: }
2714: if (perms) *perms = (const PetscInt ***) lag->symperms;
2715: if (flips) *flips = (const PetscScalar ***) lag->symflips;
2716: return(0);
2717: }
2719: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2720: {
2721: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2726: *continuous = lag->continuous;
2727: return(0);
2728: }
2730: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2731: {
2732: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2736: lag->continuous = continuous;
2737: return(0);
2738: }
2740: /*@
2741: PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity
2743: Not Collective
2745: Input Parameter:
2746: . sp - the PetscDualSpace
2748: Output Parameter:
2749: . continuous - flag for element continuity
2751: Level: intermediate
2753: .seealso: PetscDualSpaceLagrangeSetContinuity()
2754: @*/
2755: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2756: {
2762: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2763: return(0);
2764: }
2766: /*@
2767: PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous
2769: Logically Collective on sp
2771: Input Parameters:
2772: + sp - the PetscDualSpace
2773: - continuous - flag for element continuity
2775: Options Database:
2776: . -petscdualspace_lagrange_continuity <bool>
2778: Level: intermediate
2780: .seealso: PetscDualSpaceLagrangeGetContinuity()
2781: @*/
2782: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2783: {
2789: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2790: return(0);
2791: }
2793: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
2794: {
2795: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2798: *tensor = lag->tensorSpace;
2799: return(0);
2800: }
2802: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
2803: {
2804: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2807: lag->tensorSpace = tensor;
2808: return(0);
2809: }
2811: static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Lagrange(PetscDualSpace sp, PetscBool *trimmed)
2812: {
2813: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2816: *trimmed = lag->trimmed;
2817: return(0);
2818: }
2820: static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp, PetscBool trimmed)
2821: {
2822: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2825: lag->trimmed = trimmed;
2826: return(0);
2827: }
2829: static PetscErrorCode PetscDualSpaceLagrangeGetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
2830: {
2831: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2834: if (nodeType) *nodeType = lag->nodeType;
2835: if (boundary) *boundary = lag->endNodes;
2836: if (exponent) *exponent = lag->nodeExponent;
2837: return(0);
2838: }
2840: static PetscErrorCode PetscDualSpaceLagrangeSetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
2841: {
2842: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2845: if (nodeType == PETSCDTNODES_GAUSSJACOBI && exponent <= -1.) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Exponent must be > -1");
2846: lag->nodeType = nodeType;
2847: lag->endNodes = boundary;
2848: lag->nodeExponent = exponent;
2849: return(0);
2850: }
2852: /*@
2853: PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space
2855: Not collective
2857: Input Parameter:
2858: . sp - The PetscDualSpace
2860: Output Parameter:
2861: . tensor - Whether the dual space has tensor layout (vs. simplicial)
2863: Level: intermediate
2865: .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
2866: @*/
2867: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
2868: {
2874: PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));
2875: return(0);
2876: }
2878: /*@
2879: PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space
2881: Not collective
2883: Input Parameters:
2884: + sp - The PetscDualSpace
2885: - tensor - Whether the dual space has tensor layout (vs. simplicial)
2887: Level: intermediate
2889: .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
2890: @*/
2891: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
2892: {
2897: PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));
2898: return(0);
2899: }
2901: /*@
2902: PetscDualSpaceLagrangeGetTrimmed - Get the trimmed nature of the dual space
2904: Not collective
2906: Input Parameter:
2907: . sp - The PetscDualSpace
2909: Output Parameter:
2910: . 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)
2912: Level: intermediate
2914: .seealso: PetscDualSpaceLagrangeSetTrimmed(), PetscDualSpaceCreate()
2915: @*/
2916: PetscErrorCode PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp, PetscBool *trimmed)
2917: {
2923: PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTrimmed_C",(PetscDualSpace,PetscBool *),(sp,trimmed));
2924: return(0);
2925: }
2927: /*@
2928: PetscDualSpaceLagrangeSetTrimmed - Set the trimmed nature of the dual space
2930: Not collective
2932: Input Parameters:
2933: + sp - The PetscDualSpace
2934: - 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)
2936: Level: intermediate
2938: .seealso: PetscDualSpaceLagrangeGetTrimmed(), PetscDualSpaceCreate()
2939: @*/
2940: PetscErrorCode PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp, PetscBool trimmed)
2941: {
2946: PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTrimmed_C",(PetscDualSpace,PetscBool),(sp,trimmed));
2947: return(0);
2948: }
2950: /*@
2951: PetscDualSpaceLagrangeGetNodeType - Get a description of how nodes are laid out for Lagrange polynomials in this
2952: dual space
2954: Not collective
2956: Input Parameter:
2957: . sp - The PetscDualSpace
2959: Output Parameters:
2960: + nodeType - The type of nodes
2961: . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that
2962: include the boundary are Gauss-Lobatto-Jacobi nodes)
2963: - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function
2964: '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type
2966: Level: advanced
2968: .seealso: PetscDTNodeType, PetscDualSpaceLagrangeSetNodeType()
2969: @*/
2970: PetscErrorCode PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
2971: {
2979: PetscTryMethod(sp,"PetscDualSpaceLagrangeGetNodeType_C",(PetscDualSpace,PetscDTNodeType *,PetscBool *,PetscReal *),(sp,nodeType,boundary,exponent));
2980: return(0);
2981: }
2983: /*@
2984: PetscDualSpaceLagrangeSetNodeType - Set a description of how nodes are laid out for Lagrange polynomials in this
2985: dual space
2987: Logically collective
2989: Input Parameters:
2990: + sp - The PetscDualSpace
2991: . nodeType - The type of nodes
2992: . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that
2993: include the boundary are Gauss-Lobatto-Jacobi nodes)
2994: - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function
2995: '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type
2997: Level: advanced
2999: .seealso: PetscDTNodeType, PetscDualSpaceLagrangeGetNodeType()
3000: @*/
3001: PetscErrorCode PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
3002: {
3007: PetscTryMethod(sp,"PetscDualSpaceLagrangeSetNodeType_C",(PetscDualSpace,PetscDTNodeType,PetscBool,PetscReal),(sp,nodeType,boundary,exponent));
3008: return(0);
3009: }
3012: static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
3013: {
3015: sp->ops->destroy = PetscDualSpaceDestroy_Lagrange;
3016: sp->ops->view = PetscDualSpaceView_Lagrange;
3017: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
3018: sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange;
3019: sp->ops->setup = PetscDualSpaceSetUp_Lagrange;
3020: sp->ops->createheightsubspace = NULL;
3021: sp->ops->createpointsubspace = NULL;
3022: sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Lagrange;
3023: sp->ops->apply = PetscDualSpaceApplyDefault;
3024: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
3025: sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
3026: sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
3027: sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
3028: return(0);
3029: }
3031: /*MC
3032: PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals
3034: Level: intermediate
3036: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
3037: M*/
3038: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
3039: {
3040: PetscDualSpace_Lag *lag;
3041: PetscErrorCode ierr;
3045: PetscNewLog(sp,&lag);
3046: sp->data = lag;
3048: lag->tensorCell = PETSC_FALSE;
3049: lag->tensorSpace = PETSC_FALSE;
3050: lag->continuous = PETSC_TRUE;
3051: lag->numCopies = PETSC_DEFAULT;
3052: lag->numNodeSkip = PETSC_DEFAULT;
3053: lag->nodeType = PETSCDTNODES_DEFAULT;
3055: PetscDualSpaceInitialize_Lagrange(sp);
3056: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
3057: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
3058: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
3059: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
3060: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Lagrange);
3061: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Lagrange);
3062: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Lagrange);
3063: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Lagrange);
3064: return(0);
3065: }