1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2 #include <petscdmplex.h>
3 #include <petscblaslapack.h>
4
5 PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]);
6
7 struct _n_Petsc1DNodeFamily
8 {
9 PetscInt refct;
10 PetscDTNodeType nodeFamily;
11 PetscReal gaussJacobiExp;
12 PetscInt nComputed;
13 PetscReal **nodesets;
14 PetscBool endpoints;
15 };
16
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 */
Petsc1DNodeFamilyCreate(PetscDTNodeType family,PetscReal gaussJacobiExp,PetscBool endpoints,Petsc1DNodeFamily * nf)19 static PetscErrorCode Petsc1DNodeFamilyCreate(PetscDTNodeType family, PetscReal gaussJacobiExp, PetscBool endpoints, Petsc1DNodeFamily *nf)
20 {
21 Petsc1DNodeFamily f;
22 PetscErrorCode ierr;
23
24 PetscFunctionBegin;
25 ierr = PetscNew(&f);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
42 }
43
Petsc1DNodeFamilyReference(Petsc1DNodeFamily nf)44 static PetscErrorCode Petsc1DNodeFamilyReference(Petsc1DNodeFamily nf)
45 {
46 PetscFunctionBegin;
47 if (nf) nf->refct++;
48 PetscFunctionReturn(0);
49 }
50
Petsc1DNodeFamilyDestroy(Petsc1DNodeFamily * nf)51 static PetscErrorCode Petsc1DNodeFamilyDestroy(Petsc1DNodeFamily *nf) {
52 PetscInt i, nc;
53 PetscErrorCode ierr;
54
55 PetscFunctionBegin;
56 if (!(*nf)) PetscFunctionReturn(0);
57 if (--(*nf)->refct > 0) {
58 *nf = NULL;
59 PetscFunctionReturn(0);
60 }
61 nc = (*nf)->nComputed;
62 for (i = 0; i < nc; i++) {
63 ierr = PetscFree((*nf)->nodesets[i]);CHKERRQ(ierr);
64 }
65 ierr = PetscFree((*nf)->nodesets);CHKERRQ(ierr);
66 ierr = PetscFree(*nf);CHKERRQ(ierr);
67 *nf = NULL;
68 PetscFunctionReturn(0);
69 }
70
Petsc1DNodeFamilyGetNodeSets(Petsc1DNodeFamily f,PetscInt degree,PetscReal *** nodesets)71 static PetscErrorCode Petsc1DNodeFamilyGetNodeSets(Petsc1DNodeFamily f, PetscInt degree, PetscReal ***nodesets)
72 {
73 PetscInt nc;
74 PetscErrorCode ierr;
75
76 PetscFunctionBegin;
77 nc = f->nComputed;
78 if (degree >= nc) {
79 PetscInt i, j;
80 PetscReal **new_nodesets;
81 PetscReal *w;
82
83 ierr = PetscMalloc1(degree + 1, &new_nodesets);CHKERRQ(ierr);
84 ierr = PetscArraycpy(new_nodesets, f->nodesets, nc);CHKERRQ(ierr);
85 ierr = PetscFree(f->nodesets);CHKERRQ(ierr);
86 f->nodesets = new_nodesets;
87 ierr = PetscMalloc1(degree + 1, &w);CHKERRQ(ierr);
88 for (i = nc; i < degree + 1; i++) {
89 ierr = PetscMalloc1(i + 1, &(f->nodesets[i]));CHKERRQ(ierr);
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 ierr = PetscDTGaussLobattoJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w);CHKERRQ(ierr);
106 } else {
107 ierr = PetscDTGaussJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w);CHKERRQ(ierr);
108 }
109 break;
110 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family");
111 }
112 }
113 }
114 ierr = PetscFree(w);CHKERRQ(ierr);
115 f->nComputed = degree + 1;
116 }
117 *nodesets = f->nodesets;
118 PetscFunctionReturn(0);
119 }
120
121 /* http://arxiv.org/abs/2002.09421 for details */
PetscNodeRecursive_Internal(PetscInt dim,PetscInt degree,PetscReal ** nodesets,PetscInt tup[],PetscReal node[])122 static PetscErrorCode PetscNodeRecursive_Internal(PetscInt dim, PetscInt degree, PetscReal **nodesets, PetscInt tup[], PetscReal node[])
123 {
124 PetscReal w;
125 PetscInt i, j;
126 PetscErrorCode ierr;
127
128 PetscFunctionBeginHot;
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]];
137
138 for (j = 0; j < dim+1; j++) tup[dim+1+j] = tup[j+(j>=i)];
139 ierr = PetscNodeRecursive_Internal(dim-1,degree-tup[i],nodesets,&tup[dim+1],&node[dim+1]);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
146 }
147
148 /* compute simplex nodes for the biunit simplex from the 1D node family */
Petsc1DNodeFamilyComputeSimplexNodes(Petsc1DNodeFamily f,PetscInt dim,PetscInt degree,PetscReal points[])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;
158 PetscErrorCode ierr;
159
160 PetscFunctionBegin;
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) PetscFunctionReturn(0);
164 ierr = PetscCalloc1(dim+2, &tup);CHKERRQ(ierr);
165 k = 0;
166 ierr = PetscDTBinomialInt(degree + dim, dim, &npoints);CHKERRQ(ierr);
167 ierr = Petsc1DNodeFamilyGetNodeSets(f, degree, &nodesets);CHKERRQ(ierr);
168 worksize = ((dim + 2) * (dim + 3)) / 2;
169 ierr = PetscMalloc2(worksize, &nodework, worksize, &tupwork);CHKERRQ(ierr);
170 /* loop over the tuples of length dim with sum at most degree */
171 for (k = 0; k < npoints; k++) {
172 PetscInt i;
173
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 ierr = PetscNodeRecursive_Internal(dim, degree, nodesets, tupwork, nodework);CHKERRQ(ierr);
199 for (i = 0; i < dim; i++) points[dim*k + i] = nodework[i + 1];
200 break;
201 }
202 ierr = PetscDualSpaceLatticePointLexicographic_Internal(dim, degree, &tup[1]);CHKERRQ(ierr);
203 }
204 /* map from unit simplex to biunit simplex */
205 for (k = 0; k < npoints * dim; k++) points[k] = points[k] * 2. - 1.;
206 ierr = PetscFree2(nodework, tupwork);CHKERRQ(ierr);
207 ierr = PetscFree(tup);
208 PetscFunctionReturn(0);
209 }
210
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 };
329
330 /* this is just here so I can access the values in tests/ex1.c outside the library */
PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices ni,PetscInt * nodeIdxDim,PetscInt * nodeVecDim,PetscInt * nNodes,const PetscInt * nodeIdx[],const PetscReal * nodeVec[])331 PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices ni, PetscInt *nodeIdxDim, PetscInt *nodeVecDim, PetscInt *nNodes, const PetscInt *nodeIdx[], const PetscReal *nodeVec[])
332 {
333 PetscFunctionBegin;
334 *nodeIdxDim = ni->nodeIdxDim;
335 *nodeVecDim = ni->nodeVecDim;
336 *nNodes = ni->nNodes;
337 *nodeIdx = ni->nodeIdx;
338 *nodeVec = ni->nodeVec;
339 PetscFunctionReturn(0);
340 }
341
PetscLagNodeIndicesReference(PetscLagNodeIndices ni)342 static PetscErrorCode PetscLagNodeIndicesReference(PetscLagNodeIndices ni)
343 {
344 PetscFunctionBegin;
345 if (ni) ni->refct++;
346 PetscFunctionReturn(0);
347 }
348
PetscLagNodeIndicesDestroy(PetscLagNodeIndices * ni)349 static PetscErrorCode PetscLagNodeIndicesDestroy(PetscLagNodeIndices *ni) {
350 PetscErrorCode ierr;
351
352 PetscFunctionBegin;
353 if (!(*ni)) PetscFunctionReturn(0);
354 if (--(*ni)->refct > 0) {
355 *ni = NULL;
356 PetscFunctionReturn(0);
357 }
358 ierr = PetscFree((*ni)->nodeIdx);CHKERRQ(ierr);
359 ierr = PetscFree((*ni)->nodeVec);CHKERRQ(ierr);
360 ierr = PetscFree((*ni)->perm);CHKERRQ(ierr);
361 ierr = PetscFree(*ni);CHKERRQ(ierr);
362 *ni = NULL;
363 PetscFunctionReturn(0);
364 }
365
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 */
PetscLagNodeIndicesComputeVertexOrder(DM dm,PetscLagNodeIndices ni,PetscBool sortIdx)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;
390
391 PetscFunctionBegin;
392 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
393 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
394 nVerts = vEnd - vStart;
395 ierr = PetscMalloc1(nVerts, &closureOrder);CHKERRQ(ierr);
396 ierr = PetscMalloc1(nVerts, &invClosureOrder);CHKERRQ(ierr);
397 ierr = PetscMalloc1(nVerts, &revlexOrder);CHKERRQ(ierr);
398 if (sortIdx) { /* bubble sort nodeIdx into revlex order */
399 PetscInt nodeIdxDim = ni->nodeIdxDim;
400 PetscInt *idxOrder;
401
402 ierr = PetscMalloc1(nVerts * nodeIdxDim, &newNodeIdx);CHKERRQ(ierr);
403 ierr = PetscMalloc1(nVerts, &idxOrder);CHKERRQ(ierr);
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;
410
411 for (d = nodeIdxDim - 1; d >= 0; d--) if ((diff = (iv[d] - iw[d]))) break;
412 if (diff > 0) {
413 PetscInt swap = idxOrder[v];
414
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 ierr = PetscFree(ni->nodeIdx);CHKERRQ(ierr);
426 ni->nodeIdx = newNodeIdx;
427 newNodeIdx = NULL;
428 ierr = PetscFree(idxOrder);CHKERRQ(ierr);
429 }
430 ierr = DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
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 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
435 ierr = DMGetCoordinatesLocal(dm, &coordVec);CHKERRQ(ierr);
436 ierr = VecGetArrayRead(coordVec, &coords);CHKERRQ(ierr);
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;
444
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];
448
449 revlexOrder[v] = revlexOrder[w];
450 revlexOrder[w] = swap;
451 }
452 }
453 }
454 ierr = VecRestoreArrayRead(coordVec, &coords);CHKERRQ(ierr);
455 ierr = PetscMalloc1(ni->nodeIdxDim * nVerts, &newNodeIdx);CHKERRQ(ierr);
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 ierr = PetscFree(ni->nodeIdx);CHKERRQ(ierr);
463 ni->nodeIdx = newNodeIdx;
464 ni->perm = invClosureOrder;
465 ierr = PetscFree(revlexOrder);CHKERRQ(ierr);
466 ierr = PetscFree(closureOrder);CHKERRQ(ierr);
467 PetscFunctionReturn(0);
468 }
469
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 */
PetscLagNodeIndicesCreateSimplexVertices(DM dm,PetscLagNodeIndices * nodeIndices)472 static PetscErrorCode PetscLagNodeIndicesCreateSimplexVertices(DM dm, PetscLagNodeIndices *nodeIndices)
473 {
474 PetscLagNodeIndices ni;
475 PetscInt dim, d;
476
477 PetscErrorCode ierr;
478
479 PetscFunctionBegin;
480 ierr = PetscNew(&ni);CHKERRQ(ierr);
481 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
482 ni->nodeIdxDim = dim + 1;
483 ni->nodeVecDim = 0;
484 ni->nNodes = dim + 1;
485 ni->refct = 1;
486 ierr = PetscCalloc1((dim + 1)*(dim + 1), &(ni->nodeIdx));CHKERRQ(ierr);
487 for (d = 0; d < dim + 1; d++) ni->nodeIdx[d*(dim + 2)] = 1;
488 ierr = PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_FALSE);CHKERRQ(ierr);
489 *nodeIndices = ni;
490 PetscFunctionReturn(0);
491 }
492
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 */
PetscLagNodeIndicesCreateTensorVertices(DM dm,PetscLagNodeIndices facetni,PetscLagNodeIndices * nodeIndices)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;
504
505 PetscErrorCode ierr;
506
507 PetscFunctionBegin;
508 ierr = PetscNew(&ni);CHKERRQ(ierr);
509 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
510 ni->nodeIdxDim = nodeIdxDim = subNodeIdxDim + 2;
511 ni->nodeVecDim = 0;
512 ni->nNodes = nVerts = 2 * nSubVerts;
513 ni->refct = 1;
514 ierr = PetscCalloc1(nodeIdxDim * nVerts, &(ni->nodeIdx));CHKERRQ(ierr);
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 ierr = PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_TRUE);CHKERRQ(ierr);
525 *nodeIndices = ni;
526 PetscFunctionReturn(0);
527 }
528
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 */
PetscLagNodeIndicesPushForward(DM dm,PetscLagNodeIndices vert,PetscInt p,PetscLagNodeIndices vertp,PetscLagNodeIndices nodep,PetscInt ornt,PetscInt formDegree,PetscInt pfNodeIdx[],PetscReal pfNodeVec[])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;
568
569 PetscFunctionBegin;
570 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
571 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
572 ierr = DMGetCoordinatesLocal(dm, &coordVec);CHKERRQ(ierr);
573 ierr = DMPlexGetPointDepth(dm, p, &pdepth);CHKERRQ(ierr);
574 pdim = pdepth != depth ? pdepth != 0 ? pdepth : 0 : dim;
575 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
576 ierr = DMGetWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);CHKERRQ(ierr);
577 ierr = DMPlexGetTransitiveClosure_Internal(dm, p, ornt, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
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 ierr = DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
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 ierr = DMGetWorkArray(dm, dim * dim, MPIU_REAL, &J);CHKERRQ(ierr);
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;
614
615 ierr = DMPlexGetTransitiveClosure_Internal(dm, p, 0, PETSC_TRUE, &closureSize2, &closure2);CHKERRQ(ierr);
616 ierr = PetscMalloc1(dim * nSubVert, &newCoords);CHKERRQ(ierr);
617 ierr = VecGetArrayRead(coordVec, &oldCoords);CHKERRQ(ierr);
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 ierr = VecRestoreArrayRead(coordVec, &oldCoords);CHKERRQ(ierr);
625 ierr = DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize2, &closure2);CHKERRQ(ierr);
626 ierr = VecPlaceArray(coordVec, newCoords);CHKERRQ(ierr);
627 }
628 ierr = DMPlexComputeCellGeometryAffineFEM(dm, p, NULL, J, NULL, &detJ);CHKERRQ(ierr);
629 if (ornt != 0) {
630 ierr = VecResetArray(coordVec);CHKERRQ(ierr);
631 ierr = PetscFree(newCoords);CHKERRQ(ierr);
632 }
633 ierr = DMRestoreWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);CHKERRQ(ierr);
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 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);CHKERRQ(ierr);
640 ierr = PetscDTBinomialInt(pdim, PetscAbsInt(formDegree), &pNk);CHKERRQ(ierr);
641 ierr = DMGetWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);CHKERRQ(ierr);
642 ierr = PetscDTAltVPullbackMatrix(pdim, dim, J, formDegree, Jstar);CHKERRQ(ierr);
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 ierr = DMRestoreWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);CHKERRQ(ierr);
651 ierr = DMRestoreWorkArray(dm, dim * dim, MPIU_REAL, &J);CHKERRQ(ierr);
652 PetscFunctionReturn(0);
653 }
654
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 */
PetscLagNodeIndicesTensor(PetscLagNodeIndices tracei,PetscInt dimT,PetscInt kT,PetscLagNodeIndices fiberi,PetscInt dimF,PetscInt kF,PetscLagNodeIndices * nodeIndices)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;
671 PetscErrorCode ierr;
672
673 PetscFunctionBegin;
674 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);CHKERRQ(ierr);
675 ierr = PetscDTBinomialInt(dimT, PetscAbsInt(kT), &NkT);CHKERRQ(ierr);
676 ierr = PetscDTBinomialInt(dimF, PetscAbsInt(kF), &NkF);CHKERRQ(ierr);
677 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kT), &MkT);CHKERRQ(ierr);
678 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kF), &MkF);CHKERRQ(ierr);
679 ierr = PetscNew(&ni);CHKERRQ(ierr);
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 ierr = PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));CHKERRQ(ierr);
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;
689
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 }
694
695 /* now wedge together the push-forward vectors */
696 ierr = PetscMalloc1(nNodes * Nk, &(ni->nodeVec));CHKERRQ(ierr);
697 ierr = PetscCalloc2(dimT*dim, &projT, dimF*dim, &projF);CHKERRQ(ierr);
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 ierr = PetscMalloc2(MkT*NkT, &projTstar, MkF*NkF, &projFstar);CHKERRQ(ierr);
701 ierr = PetscDTAltVPullbackMatrix(dim, dimT, projT, kT, projTstar);CHKERRQ(ierr);
702 ierr = PetscDTAltVPullbackMatrix(dim, dimF, projF, kF, projFstar);CHKERRQ(ierr);
703 ierr = PetscMalloc6(MkT, &workT, MkT, &workT2, MkF, &workF, MkF, &workF2, Nk, &work, Nk, &work2);CHKERRQ(ierr);
704 ierr = PetscMalloc1(Nk * MkT, &wedgeMat);CHKERRQ(ierr);
705 sign = (PetscAbsInt(kT * kF) & 1) ? -1. : 1.;
706 for (l = 0, j = 0; j < fiberi->nNodes; j++) {
707 PetscInt d, e;
708
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 ierr = PetscDTAltVStar(dim, PetscAbsInt(kF), 1, workF2, workF);CHKERRQ(ierr);
719 }
720 /* Compute the matrix that wedges this form with one of the trace k-form */
721 ierr = PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kF), PetscAbsInt(kT), workF, wedgeMat);CHKERRQ(ierr);
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 ierr = PetscDTAltVStar(dim, PetscAbsInt(kT), 1, workT2, workT);CHKERRQ(ierr);
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 ierr = PetscDTAltVStar(dim, PetscAbsInt(formDegree), -1, work2, work);CHKERRQ(ierr);
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 ierr = PetscFree(wedgeMat);CHKERRQ(ierr);
750 ierr = PetscFree6(workT, workT2, workF, workF2, work, work2);CHKERRQ(ierr);
751 ierr = PetscFree2(projTstar, projFstar);CHKERRQ(ierr);
752 ierr = PetscFree2(projT, projF);CHKERRQ(ierr);
753 *nodeIndices = ni;
754 PetscFunctionReturn(0);
755 }
756
757 /* simple union of two sets of nodes */
PetscLagNodeIndicesMerge(PetscLagNodeIndices niA,PetscLagNodeIndices niB,PetscLagNodeIndices * nodeIndices)758 static PetscErrorCode PetscLagNodeIndicesMerge(PetscLagNodeIndices niA, PetscLagNodeIndices niB, PetscLagNodeIndices *nodeIndices)
759 {
760 PetscLagNodeIndices ni;
761 PetscInt nodeIdxDim, nodeVecDim, nNodes;
762 PetscErrorCode ierr;
763
764 PetscFunctionBegin;
765 ierr = PetscNew(&ni);CHKERRQ(ierr);
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 ierr = PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));CHKERRQ(ierr);
773 ierr = PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));CHKERRQ(ierr);
774 ierr = PetscArraycpy(ni->nodeIdx, niA->nodeIdx, niA->nNodes * nodeIdxDim);CHKERRQ(ierr);
775 ierr = PetscArraycpy(ni->nodeVec, niA->nodeVec, niA->nNodes * nodeVecDim);CHKERRQ(ierr);
776 ierr = PetscArraycpy(&(ni->nodeIdx[niA->nNodes * nodeIdxDim]), niB->nodeIdx, niB->nNodes * nodeIdxDim);CHKERRQ(ierr);
777 ierr = PetscArraycpy(&(ni->nodeVec[niA->nNodes * nodeVecDim]), niB->nodeVec, niB->nNodes * nodeVecDim);CHKERRQ(ierr);
778 *nodeIndices = ni;
779 PetscFunctionReturn(0);
780 }
781
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 }
795
796 PETSCTUPINTCOMPREVLEX(3)
797 PETSCTUPINTCOMPREVLEX(4)
798 PETSCTUPINTCOMPREVLEX(5)
799 PETSCTUPINTCOMPREVLEX(6)
800 PETSCTUPINTCOMPREVLEX(7)
801
PetscTupIntCompRevlex_N(const void * a,const void * b)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 }
815
816 /* The nodes are not necessarily in revlex order wrt nodeIdx: get the permutation
817 * that puts them in that order */
PetscLagNodeIndicesGetPermutation(PetscLagNodeIndices ni,PetscInt * perm[])818 static PetscErrorCode PetscLagNodeIndicesGetPermutation(PetscLagNodeIndices ni, PetscInt *perm[])
819 {
820 PetscErrorCode ierr;
821
822 PetscFunctionBegin;
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 *);
830
831 ierr = PetscMalloc1((nodeIdxDim + 2) * m, &sorter);CHKERRQ(ierr);
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 ierr = PetscMalloc1(m, &prm);CHKERRQ(ierr);
859 for (i = 0; i < m; i++) prm[i] = sorter[(nodeIdxDim + 2) * i + 1];
860 ni->perm = prm;
861 ierr = PetscFree(sorter);
862 }
863 *perm = ni->perm;
864 PetscFunctionReturn(0);
865 }
866
PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)867 static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
868 {
869 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
870 PetscErrorCode ierr;
871
872 PetscFunctionBegin;
873 if (lag->symperms) {
874 PetscInt **selfSyms = lag->symperms[0];
875
876 if (selfSyms) {
877 PetscInt i, **allocated = &selfSyms[-lag->selfSymOff];
878
879 for (i = 0; i < lag->numSelfSym; i++) {
880 ierr = PetscFree(allocated[i]);CHKERRQ(ierr);
881 }
882 ierr = PetscFree(allocated);CHKERRQ(ierr);
883 }
884 ierr = PetscFree(lag->symperms);CHKERRQ(ierr);
885 }
886 if (lag->symflips) {
887 PetscScalar **selfSyms = lag->symflips[0];
888
889 if (selfSyms) {
890 PetscInt i;
891 PetscScalar **allocated = &selfSyms[-lag->selfSymOff];
892
893 for (i = 0; i < lag->numSelfSym; i++) {
894 ierr = PetscFree(allocated[i]);CHKERRQ(ierr);
895 }
896 ierr = PetscFree(allocated);CHKERRQ(ierr);
897 }
898 ierr = PetscFree(lag->symflips);CHKERRQ(ierr);
899 }
900 ierr = Petsc1DNodeFamilyDestroy(&(lag->nodeFamily));CHKERRQ(ierr);
901 ierr = PetscLagNodeIndicesDestroy(&(lag->vertIndices));CHKERRQ(ierr);
902 ierr = PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));CHKERRQ(ierr);
903 ierr = PetscLagNodeIndicesDestroy(&(lag->allNodeIndices));CHKERRQ(ierr);
904 ierr = PetscFree(lag);CHKERRQ(ierr);
905 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);CHKERRQ(ierr);
906 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);CHKERRQ(ierr);
907 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);CHKERRQ(ierr);
908 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);CHKERRQ(ierr);
909 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL);CHKERRQ(ierr);
910 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL);CHKERRQ(ierr);
911 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL);CHKERRQ(ierr);
912 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL);CHKERRQ(ierr);
913 PetscFunctionReturn(0);
914 }
915
PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp,PetscViewer viewer)916 static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
917 {
918 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
919 PetscErrorCode ierr;
920
921 PetscFunctionBegin;
922 ierr = PetscViewerASCIIPrintf(viewer, "%s %s%sLagrange dual space\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "", lag->trimmed ? "trimmed " : "");CHKERRQ(ierr);
923 PetscFunctionReturn(0);
924 }
925
PetscDualSpaceView_Lagrange(PetscDualSpace sp,PetscViewer viewer)926 static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
927 {
928 PetscBool iascii;
929 PetscErrorCode ierr;
930
931 PetscFunctionBegin;
932 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
933 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
934 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
935 if (iascii) {ierr = PetscDualSpaceLagrangeView_Ascii(sp, viewer);CHKERRQ(ierr);}
936 PetscFunctionReturn(0);
937 }
938
PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems * PetscOptionsObject,PetscDualSpace sp)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;
945 PetscErrorCode ierr;
946
947 PetscFunctionBegin;
948 ierr = PetscDualSpaceLagrangeGetContinuity(sp, &continuous);CHKERRQ(ierr);
949 ierr = PetscDualSpaceLagrangeGetTensor(sp, &tensor);CHKERRQ(ierr);
950 ierr = PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);CHKERRQ(ierr);
951 ierr = PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &nodeEndpoints, &nodeExponent);CHKERRQ(ierr);
952 if (nodeType == PETSCDTNODES_DEFAULT) nodeType = PETSCDTNODES_GAUSSJACOBI;
953 ierr = PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");CHKERRQ(ierr);
954 ierr = PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);CHKERRQ(ierr);
955 if (flg) {ierr = PetscDualSpaceLagrangeSetContinuity(sp, continuous);CHKERRQ(ierr);}
956 ierr = PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetTensor", tensor, &tensor, &flg);CHKERRQ(ierr);
957 if (flg) {ierr = PetscDualSpaceLagrangeSetTensor(sp, tensor);CHKERRQ(ierr);}
958 ierr = PetscOptionsBool("-petscdualspace_lagrange_trimmed", "Flag for trimmed dual space", "PetscDualSpaceLagrangeSetTrimmed", trimmed, &trimmed, &flg);CHKERRQ(ierr);
959 if (flg) {ierr = PetscDualSpaceLagrangeSetTrimmed(sp, trimmed);CHKERRQ(ierr);}
960 ierr = PetscOptionsEnum("-petscdualspace_lagrange_node_type", "Lagrange node location type", "PetscDualSpaceLagrangeSetNodeType", PetscDTNodeTypes, (PetscEnum)nodeType, (PetscEnum *)&nodeType, &flg);CHKERRQ(ierr);
961 ierr = PetscOptionsBool("-petscdualspace_lagrange_node_endpoints", "Flag for nodes that include endpoints", "PetscDualSpaceLagrangeSetNodeType", nodeEndpoints, &nodeEndpoints, &flg2);CHKERRQ(ierr);
962 flg3 = PETSC_FALSE;
963 if (nodeType == PETSCDTNODES_GAUSSJACOBI) {
964 ierr = PetscOptionsReal("-petscdualspace_lagrange_node_exponent", "Gauss-Jacobi weight function exponent", "PetscDualSpaceLagrangeSetNodeType", nodeExponent, &nodeExponent, &flg3);CHKERRQ(ierr);
965 }
966 if (flg || flg2 || flg3) {ierr = PetscDualSpaceLagrangeSetNodeType(sp, nodeType, nodeEndpoints, nodeExponent);CHKERRQ(ierr);}
967 ierr = PetscOptionsTail();CHKERRQ(ierr);
968 PetscFunctionReturn(0);
969 }
970
PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp,PetscDualSpace spNew)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;
978
979 PetscFunctionBegin;
980 ierr = PetscDualSpaceLagrangeGetContinuity(sp, &cont);CHKERRQ(ierr);
981 ierr = PetscDualSpaceLagrangeSetContinuity(spNew, cont);CHKERRQ(ierr);
982 ierr = PetscDualSpaceLagrangeGetTensor(sp, &tensor);CHKERRQ(ierr);
983 ierr = PetscDualSpaceLagrangeSetTensor(spNew, tensor);CHKERRQ(ierr);
984 ierr = PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);CHKERRQ(ierr);
985 ierr = PetscDualSpaceLagrangeSetTrimmed(spNew, trimmed);CHKERRQ(ierr);
986 ierr = PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &boundary, &exponent);CHKERRQ(ierr);
987 ierr = PetscDualSpaceLagrangeSetNodeType(spNew, nodeType, boundary, exponent);CHKERRQ(ierr);
988 if (lag->nodeFamily) {
989 PetscDualSpace_Lag *lagnew = (PetscDualSpace_Lag *) spNew->data;
990
991 ierr = Petsc1DNodeFamilyReference(lag->nodeFamily);CHKERRQ(ierr);
992 lagnew->nodeFamily = lag->nodeFamily;
993 }
994 PetscFunctionReturn(0);
995 }
996
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 */
PetscDualSpaceCreateEdgeSubspace_Lagrange(PetscDualSpace sp,PetscInt order,PetscInt k,PetscInt Nc,PetscBool interiorOnly,PetscDualSpace * bdsp)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;
1004
1005 PetscFunctionBegin;
1006 ierr = PetscDualSpaceDuplicate(sp,bdsp);CHKERRQ(ierr);
1007 ierr = PetscDualSpaceSetFormDegree(*bdsp, k);CHKERRQ(ierr);
1008 ierr = PetscDualSpaceCreateReferenceCell(*bdsp, 1, PETSC_TRUE, &K);CHKERRQ(ierr);
1009 ierr = PetscDualSpaceSetDM(*bdsp, K);CHKERRQ(ierr);
1010 ierr = DMDestroy(&K);CHKERRQ(ierr);
1011 ierr = PetscDualSpaceSetOrder(*bdsp, order);CHKERRQ(ierr);
1012 ierr = PetscDualSpaceSetNumComponents(*bdsp, Nc);CHKERRQ(ierr);
1013 newlag = (PetscDualSpace_Lag *) (*bdsp)->data;
1014 newlag->interiorOnly = interiorOnly;
1015 ierr = PetscDualSpaceSetUp(*bdsp);CHKERRQ(ierr);
1016 PetscFunctionReturn(0);
1017 }
1018
1019 /* just the points, weights aren't handled */
PetscQuadratureCreateTensor(PetscQuadrature trace,PetscQuadrature fiber,PetscQuadrature * product)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;
1030
1031 PetscFunctionBegin;
1032 ierr = PetscQuadratureGetData(trace, &dimTrace, NULL, &numPointsTrace, &pointsTrace, NULL);CHKERRQ(ierr);
1033 ierr = PetscQuadratureGetData(fiber, &dimFiber, NULL, &numPointsFiber, &pointsFiber, NULL);CHKERRQ(ierr);
1034 dim = dimTrace + dimFiber;
1035 numPoints = numPointsFiber * numPointsTrace;
1036 ierr = PetscMalloc1(numPoints * dim, &points);CHKERRQ(ierr);
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 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, product);CHKERRQ(ierr);
1044 ierr = PetscQuadratureSetData(*product, dim, 0, numPoints, points, NULL);CHKERRQ(ierr);
1045 PetscFunctionReturn(0);
1046 }
1047
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 */
MatTensorAltV(Mat trace,Mat fiber,PetscInt dimTrace,PetscInt kTrace,PetscInt dimFiber,PetscInt kFiber,Mat * product)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 */
1065 PetscErrorCode ierr;
1066
1067 PetscFunctionBegin;
1068 ierr = MatGetSize(trace, &mTrace, &nTrace);CHKERRQ(ierr);
1069 ierr = PetscDTBinomialInt(dimTrace, PetscAbsInt(kTrace), &NkTrace);CHKERRQ(ierr);
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 ierr = MatGetSize(fiber, &mFiber, &nFiber);CHKERRQ(ierr);
1072 ierr = PetscDTBinomialInt(dimFiber, PetscAbsInt(kFiber), &NkFiber);CHKERRQ(ierr);
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 ierr = PetscMalloc2(mTrace, &nnzTrace, mFiber, &nnzFiber);CHKERRQ(ierr);
1075 for (i = 0; i < mTrace; i++) {
1076 ierr = MatGetRow(trace, i, &(nnzTrace[i]), NULL, NULL);CHKERRQ(ierr);
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 ierr = MatGetRow(fiber, i, &(nnzFiber[i]), NULL, NULL);CHKERRQ(ierr);
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 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr);
1086 m = mTrace * mFiber;
1087 ierr = PetscMalloc1(m, &nnz);CHKERRQ(ierr);
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 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &prod);CHKERRQ(ierr);
1091 ierr = PetscFree(nnz);CHKERRQ(ierr);
1092 ierr = PetscFree2(nnzTrace,nnzFiber);CHKERRQ(ierr);
1093 /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1094 ierr = MatSetOption(prod, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr);
1095 /* compute pullbacks */
1096 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kTrace), &dT);CHKERRQ(ierr);
1097 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kFiber), &dF);CHKERRQ(ierr);
1098 ierr = PetscMalloc4(dimTrace * dim, &projT, dimFiber * dim, &projF, dT * NkTrace, &projTstar, dF * NkFiber, &projFstar);CHKERRQ(ierr);
1099 ierr = PetscArrayzero(projT, dimTrace * dim);CHKERRQ(ierr);
1100 for (i = 0; i < dimTrace; i++) projT[i * (dim + 1)] = 1.;
1101 ierr = PetscArrayzero(projF, dimFiber * dim);CHKERRQ(ierr);
1102 for (i = 0; i < dimFiber; i++) projF[i * (dim + 1) + dimTrace] = 1.;
1103 ierr = PetscDTAltVPullbackMatrix(dim, dimTrace, projT, kTrace, projTstar);CHKERRQ(ierr);
1104 ierr = PetscDTAltVPullbackMatrix(dim, dimFiber, projF, kFiber, projFstar);CHKERRQ(ierr);
1105 ierr = PetscMalloc5(dT, &workT, dF, &workF, Nk, &work, Nk, &workstar, Nk, &workS);CHKERRQ(ierr);
1106 ierr = PetscMalloc2(dT, &workT2, dF, &workF2);CHKERRQ(ierr);
1107 ierr = PetscMalloc1(Nk * dT, &wedgeMat);CHKERRQ(ierr);
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;
1113
1114 ierr = MatGetRow(fiber, iF, &ncolsF, &colsF, &valsF);CHKERRQ(ierr);
1115 nformsF = ncolsF / NkFiber;
1116 for (iT = 0; iT < mTrace; iT++, i++) {
1117 PetscInt ncolsT, nformsT;
1118 const PetscInt *colsT;
1119 const PetscScalar *valsT;
1120
1121 ierr = MatGetRow(trace, iT, &ncolsT, &colsT, &valsT);CHKERRQ(ierr);
1122 nformsT = ncolsT / NkTrace;
1123 for (j = 0, jF = 0; jF < nformsF; jF++) {
1124 PetscInt colF = colsF[jF * NkFiber] / NkFiber;
1125
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 ierr = PetscDTAltVStar(dim, PetscAbsInt(kFiber), 1, workF2, workF);CHKERRQ(ierr);
1134 }
1135 ierr = PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kFiber), PetscAbsInt(kTrace), workF, wedgeMat);CHKERRQ(ierr);
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;
1140
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 ierr = PetscDTAltVStar(dim, PetscAbsInt(kTrace), 1, workT2, workT);CHKERRQ(ierr);
1149 }
1150
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 ierr = PetscDTAltVStar(dim, PetscAbsInt(k), -1, work, workstar);CHKERRQ(ierr);
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 ierr = MatSetValue(prod, i, col * Nk + l, vals[l], INSERT_VALUES);CHKERRQ(ierr);
1174 } /* Nk */
1175 } /* jT */
1176 } /* jF */
1177 ierr = MatRestoreRow(trace, iT, &ncolsT, &colsT, &valsT);CHKERRQ(ierr);
1178 } /* iT */
1179 ierr = MatRestoreRow(fiber, iF, &ncolsF, &colsF, &valsF);CHKERRQ(ierr);
1180 } /* iF */
1181 ierr = PetscFree(wedgeMat);CHKERRQ(ierr);
1182 ierr = PetscFree4(projT, projF, projTstar, projFstar);CHKERRQ(ierr);
1183 ierr = PetscFree2(workT2, workF2);CHKERRQ(ierr);
1184 ierr = PetscFree5(workT, workF, work, workstar, workS);CHKERRQ(ierr);
1185 ierr = MatAssemblyBegin(prod, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1186 ierr = MatAssemblyEnd(prod, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1187 *product = prod;
1188 PetscFunctionReturn(0);
1189 }
1190
1191 /* Union of quadrature points, with an attempt to identify commont points in the two sets */
PetscQuadraturePointsMerge(PetscQuadrature quadA,PetscQuadrature quadB,PetscQuadrature * quadJoint,PetscInt * aToJoint[],PetscInt * bToJoint[])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;
1202
1203 PetscFunctionBegin;
1204 ierr = PetscQuadratureGetData(quadA, &dimA, NULL, &nA, &pointsA, NULL);CHKERRQ(ierr);
1205 ierr = PetscQuadratureGetData(quadB, &dimB, NULL, &nB, &pointsB, NULL);CHKERRQ(ierr);
1206 if (dimA != dimB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Quadrature points must be in the same dimension");
1207 nJoint = nA;
1208 ierr = PetscMalloc1(nA, &aToJ);CHKERRQ(ierr);
1209 for (i = 0; i < nA; i++) aToJ[i] = i;
1210 ierr = PetscMalloc1(nB, &bToJ);CHKERRQ(ierr);
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 ierr = PetscMalloc1(nJoint * dimA, &pointsJoint);CHKERRQ(ierr);
1227 ierr = PetscArraycpy(pointsJoint, pointsA, nA * dimA);CHKERRQ(ierr);
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 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &qJ);CHKERRQ(ierr);
1234 ierr = PetscQuadratureSetData(qJ, dimA, 0, nJoint, pointsJoint, NULL);CHKERRQ(ierr);
1235 *quadJoint = qJ;
1236 PetscFunctionReturn(0);
1237 }
1238
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 */
MatricesMerge(Mat matA,Mat matB,PetscInt dim,PetscInt k,PetscInt numMerged,const PetscInt aToMerged[],const PetscInt bToMerged[],Mat * matMerged)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;
1248 PetscErrorCode ierr;
1249
1250 PetscFunctionBegin;
1251 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr);
1252 ierr = MatGetSize(matA, &mA, &nA);CHKERRQ(ierr);
1253 if (nA % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matA column space not a multiple of k-form size");
1254 ierr = MatGetSize(matB, &mB, &nB);CHKERRQ(ierr);
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 ierr = PetscMalloc1(m, &nnz);CHKERRQ(ierr);
1259 maxnnz = 0;
1260 for (i = 0; i < mA; i++) {
1261 ierr = MatGetRow(matA, i, &(nnz[i]), NULL, NULL);CHKERRQ(ierr);
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 ierr = MatGetRow(matB, i, &(nnz[i+mA]), NULL, NULL);CHKERRQ(ierr);
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 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &M);CHKERRQ(ierr);
1271 ierr = PetscFree(nnz);CHKERRQ(ierr);
1272 /* reasoning about which points each dof needs depends on having zeros computed at points preserved */
1273 ierr = MatSetOption(M, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr);
1274 ierr = PetscMalloc1(maxnnz, &work);CHKERRQ(ierr);
1275 for (i = 0; i < mA; i++) {
1276 const PetscInt *cols;
1277 const PetscScalar *vals;
1278 PetscInt nCols;
1279 ierr = MatGetRow(matA, i, &nCols, &cols, &vals);CHKERRQ(ierr);
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 ierr = MatSetValuesBlocked(M, 1, &i, nCols, work, vals, INSERT_VALUES);CHKERRQ(ierr);
1285 ierr = MatRestoreRow(matA, i, &nCols, &cols, &vals);CHKERRQ(ierr);
1286 }
1287 for (i = 0; i < mB; i++) {
1288 const PetscInt *cols;
1289 const PetscScalar *vals;
1290
1291 PetscInt row = i + mA;
1292 PetscInt nCols;
1293 ierr = MatGetRow(matB, i, &nCols, &cols, &vals);CHKERRQ(ierr);
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 ierr = MatSetValuesBlocked(M, 1, &row, nCols, work, vals, INSERT_VALUES);CHKERRQ(ierr);
1299 ierr = MatRestoreRow(matB, i, &nCols, &cols, &vals);CHKERRQ(ierr);
1300 }
1301 ierr = PetscFree(work);CHKERRQ(ierr);
1302 ierr = MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1303 ierr = MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1304 *matMerged = M;
1305 PetscFunctionReturn(0);
1306 }
1307
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 */
PetscDualSpaceCreateFacetSubspace_Lagrange(PetscDualSpace sp,DM K,PetscInt f,PetscInt k,PetscInt Ncopies,PetscBool interiorOnly,PetscDualSpace * bdsp)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;
1318
1319 PetscFunctionBegin;
1320 ierr = PetscDualSpaceGetDM(sp,&dm);CHKERRQ(ierr);
1321 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1322 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr);
1323 ierr = PetscDualSpaceDuplicate(sp,bdsp);CHKERRQ(ierr);
1324 ierr = PetscDualSpaceSetFormDegree(*bdsp,k);CHKERRQ(ierr);
1325 if (!K) {
1326 PetscBool isSimplex;
1327
1328
1329 if (depth == dim) {
1330 PetscInt coneSize;
1331
1332 pointDim = dim - 1;
1333 ierr = DMPlexGetConeSize(dm,f,&coneSize);CHKERRQ(ierr);
1334 isSimplex = (PetscBool) (coneSize == dim);
1335 ierr = PetscDualSpaceCreateReferenceCell(*bdsp, dim-1, isSimplex, &K);CHKERRQ(ierr);
1336 } else if (depth == 1) {
1337 pointDim = 0;
1338 ierr = PetscDualSpaceCreateReferenceCell(*bdsp, 0, PETSC_TRUE, &K);CHKERRQ(ierr);
1339 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
1340 } else {
1341 ierr = PetscObjectReference((PetscObject)K);CHKERRQ(ierr);
1342 ierr = DMGetDimension(K, &pointDim);CHKERRQ(ierr);
1343 }
1344 ierr = PetscDualSpaceSetDM(*bdsp, K);CHKERRQ(ierr);
1345 ierr = DMDestroy(&K);CHKERRQ(ierr);
1346 ierr = PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew);CHKERRQ(ierr);
1347 Ncnew = Nknew * Ncopies;
1348 ierr = PetscDualSpaceSetNumComponents(*bdsp, Ncnew);CHKERRQ(ierr);
1349 newlag = (PetscDualSpace_Lag *) (*bdsp)->data;
1350 newlag->interiorOnly = interiorOnly;
1351 ierr = PetscDualSpaceSetUp(*bdsp);CHKERRQ(ierr);
1352 PetscFunctionReturn(0);
1353 }
1354
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 * */
PetscDualSpaceLagrangeCreateSimplexNodeMat(Petsc1DNodeFamily nodeFamily,PetscInt dim,PetscInt sum,PetscInt Nk,PetscInt numNodeSkip,PetscQuadrature * iNodes,Mat * iMat,PetscLagNodeIndices * nodeIndices)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;
1370 PetscErrorCode ierr;
1371
1372 PetscFunctionBegin;
1373 ierr = PetscDTBinomialInt(dim + sum, dim, &nNodes);CHKERRQ(ierr);
1374 ierr = PetscDTBinomialInt(dim + extraSum, dim, &nExtraNodes);CHKERRQ(ierr);
1375
1376 ierr = PetscMalloc1(dim * nExtraNodes, &extraNodeCoords);CHKERRQ(ierr);
1377 ierr = PetscNew(&ni);CHKERRQ(ierr);
1378 ni->nodeIdxDim = dim + 1;
1379 ni->nodeVecDim = Nk;
1380 ni->nNodes = nNodes * Nk;
1381 ni->refct = 1;
1382 ierr = PetscMalloc1(nNodes * Nk * (dim + 1), &(ni->nodeIdx));CHKERRQ(ierr);
1383 ierr = PetscMalloc1(nNodes * Nk * Nk, &(ni->nodeVec));CHKERRQ(ierr);
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 ierr = Petsc1DNodeFamilyComputeSimplexNodes(nodeFamily, dim, extraSum, extraNodeCoords);CHKERRQ(ierr);
1386 if (numNodeSkip) {
1387 PetscInt k;
1388 PetscInt *tup;
1389
1390 ierr = PetscMalloc1(dim * nNodes, &nodeCoords);CHKERRQ(ierr);
1391 ierr = PetscMalloc1(dim + 1, &tup);CHKERRQ(ierr);
1392 for (k = 0; k < nNodes; k++) {
1393 PetscInt j, c;
1394 PetscInt index;
1395
1396 ierr = PetscDTIndexToBary(dim + 1, sum, k, tup);CHKERRQ(ierr);
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 ierr = PetscDTBaryToIndex(dim + 1, extraSum, tup, &index);CHKERRQ(ierr);
1404 for (j = 0; j < dim; j++) nodeCoords[k * dim + j] = extraNodeCoords[index * dim + j];
1405 }
1406 ierr = PetscFree(tup);CHKERRQ(ierr);
1407 ierr = PetscFree(extraNodeCoords);CHKERRQ(ierr);
1408 } else {
1409 PetscInt k;
1410 PetscInt *tup;
1411
1412 nodeCoords = extraNodeCoords;
1413 ierr = PetscMalloc1(dim + 1, &tup);CHKERRQ(ierr);
1414 for (k = 0; k < nNodes; k++) {
1415 PetscInt j, c;
1416
1417 ierr = PetscDTIndexToBary(dim + 1, sum, k, tup);CHKERRQ(ierr);
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 ierr = PetscFree(tup);CHKERRQ(ierr);
1428 }
1429 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &intNodes);CHKERRQ(ierr);
1430 ierr = PetscQuadratureSetData(intNodes, dim, 0, nNodes, nodeCoords, NULL);CHKERRQ(ierr);
1431 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes * Nk, nNodes * Nk, Nk, NULL, &intMat);CHKERRQ(ierr);
1432 ierr = MatSetOption(intMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);CHKERRQ(ierr);
1433 for (j = 0; j < nNodes * Nk; j++) {
1434 PetscInt rem = j % Nk;
1435 PetscInt a, aprev = j - rem;
1436 PetscInt anext = aprev + Nk;
1437
1438 for (a = aprev; a < anext; a++) {
1439 ierr = MatSetValue(intMat, j, a, (a == j) ? 1. : 0., INSERT_VALUES);CHKERRQ(ierr);
1440 }
1441 }
1442 ierr = MatAssemblyBegin(intMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1443 ierr = MatAssemblyEnd(intMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1444 *iNodes = intNodes;
1445 *iMat = intMat;
1446 *nodeIndices = ni;
1447 PetscFunctionReturn(0);
1448 }
1449
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 */
PetscDualSpaceLagrangeCreateAllNodeIdx(PetscDualSpace sp)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;
1462 PetscErrorCode ierr;
1463
1464 PetscFunctionBegin;
1465 lag = (PetscDualSpace_Lag *) sp->data;
1466 verti = lag->vertIndices;
1467 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr);
1468 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1469 ierr = PetscDualSpaceGetFormDegree(sp, &formDegree);CHKERRQ(ierr);
1470 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);CHKERRQ(ierr);
1471 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr);
1472 ierr = PetscSectionGetStorageSize(section, &nDofs);CHKERRQ(ierr);
1473 ierr = PetscNew(&ni);CHKERRQ(ierr);
1474 ni->nodeIdxDim = nodeIdxDim = verti->nodeIdxDim;
1475 ni->nodeVecDim = Nk;
1476 ni->nNodes = nDofs;
1477 ni->refct = 1;
1478 ierr = PetscMalloc1(nodeIdxDim * nDofs, &(ni->nodeIdx));CHKERRQ(ierr);
1479 ierr = PetscMalloc1(Nk * nDofs, &(ni->nodeVec));CHKERRQ(ierr);
1480 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
1481 ierr = PetscSectionGetDof(section, 0, &spintdim);CHKERRQ(ierr);
1482 if (spintdim) {
1483 ierr = PetscArraycpy(ni->nodeIdx, lag->intNodeIndices->nodeIdx, spintdim * nodeIdxDim);CHKERRQ(ierr);
1484 ierr = PetscArraycpy(ni->nodeVec, lag->intNodeIndices->nodeVec, spintdim * Nk);CHKERRQ(ierr);
1485 }
1486 for (p = pStart + 1; p < pEnd; p++) {
1487 PetscDualSpace psp = sp->pointSpaces[p];
1488 PetscDualSpace_Lag *plag;
1489 PetscInt dof, off;
1490
1491 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
1492 if (!dof) continue;
1493 plag = (PetscDualSpace_Lag *) psp->data;
1494 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
1495 ierr = PetscLagNodeIndicesPushForward(dm, verti, p, plag->vertIndices, plag->intNodeIndices, 0, formDegree, &(ni->nodeIdx[off * nodeIdxDim]), &(ni->nodeVec[off * Nk]));CHKERRQ(ierr);
1496 }
1497 lag->allNodeIndices = ni;
1498 PetscFunctionReturn(0);
1499 }
1500
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 */
PetscDualSpaceCreateAllDataFromInteriorData(PetscDualSpace sp)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;
1521
1522 PetscFunctionBegin;
1523 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr);
1524 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1525 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr);
1526 ierr = PetscSectionGetStorageSize(section, &nDofs);CHKERRQ(ierr);
1527 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
1528 ierr = PetscDualSpaceGetFormDegree(sp, &k);CHKERRQ(ierr);
1529 ierr = PetscDualSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr);
1530 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr);
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;
1537
1538 ierr = PetscDualSpaceGetPointSubspace(sp, p, &psp);CHKERRQ(ierr);
1539 if (!psp) continue;
1540 ierr = PetscDualSpaceGetDM(psp, &pdm);CHKERRQ(ierr);
1541 ierr = DMGetDimension(pdm, &pdim);CHKERRQ(ierr);
1542 if (pdim < PetscAbsInt(k)) continue;
1543 ierr = PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);CHKERRQ(ierr);
1544 ierr = PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);CHKERRQ(ierr);
1545 if (intNodes) {
1546 PetscInt nNodesp;
1547
1548 ierr = PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, NULL, NULL);CHKERRQ(ierr);
1549 nNodes += nNodesp;
1550 }
1551 if (intMat) {
1552 PetscInt maxNzsp;
1553 PetscInt maxNzformsp;
1554
1555 ierr = MatSeqAIJGetMaxRowNonzeros(intMat, &maxNzsp);CHKERRQ(ierr);
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 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nDofs, nNodes * Nc, maxNzforms * Nk, NULL, &allMat);CHKERRQ(ierr);
1562 ierr = MatSetOption(allMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);CHKERRQ(ierr);
1563 ierr = PetscMalloc7(dim, &v0, dim, &pv0, dim * dim, &J, dim * dim, &Jinv, Nk * Nk, &L, maxNzforms * Nk, &work, maxNzforms * Nk, &iwork);CHKERRQ(ierr);
1564 for (j = 0; j < dim; j++) pv0[j] = -1.;
1565 ierr = PetscMalloc1(dim * nNodes, &nodes);CHKERRQ(ierr);
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;
1574
1575 ierr = PetscDualSpaceGetPointSubspace(sp, p, &psp);CHKERRQ(ierr);
1576 if (!psp) continue;
1577 ierr = PetscDualSpaceGetDM(psp, &pdm);CHKERRQ(ierr);
1578 ierr = DMGetDimension(pdm, &pdim);CHKERRQ(ierr);
1579 if (pdim < PetscAbsInt(k)) continue;
1580 ierr = PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);CHKERRQ(ierr);
1581 if (intNodes == NULL && intMat == NULL) continue;
1582 ierr = PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);CHKERRQ(ierr);
1583 if (p) {
1584 ierr = DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, Jinv, &detJ);CHKERRQ(ierr);
1585 } else { /* identity */
1586 PetscInt i,j;
1587
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;
1594
1595 for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j];
1596 }
1597 ierr = PetscDTAltVPullbackMatrix(pdim, dim, J, k, L);CHKERRQ(ierr);
1598 if (intNodes) { /* push forward quadrature locations by the affine transformation */
1599 PetscInt nNodesp;
1600 const PetscReal *nodesp;
1601 PetscInt j;
1602
1603 ierr = PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, &nodesp, NULL);CHKERRQ(ierr);
1604 for (j = 0; j < nNodesp; j++, countNodes++) {
1605 PetscInt d, e;
1606
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;
1618
1619 ierr = PetscSectionGetDof(section, p, &nrows);CHKERRQ(ierr);
1620 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr);
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;
1627
1628 ierr = MatGetRow(intMat, j, &ncols, &cols, &vals);CHKERRQ(ierr);
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;
1632
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 ierr = MatSetValues(allMat, 1, &row, (ncols / pNk) * Nk, iwork, work, INSERT_VALUES);CHKERRQ(ierr);
1649 ierr = MatRestoreRow(intMat, j, &ncols, &cols, &vals);CHKERRQ(ierr);
1650 }
1651 }
1652 }
1653 ierr = MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1654 ierr = MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1655 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &allNodes);CHKERRQ(ierr);
1656 ierr = PetscQuadratureSetData(allNodes, dim, 0, nNodes, nodes, NULL);CHKERRQ(ierr);
1657 ierr = PetscFree7(v0, pv0, J, Jinv, L, work, iwork);CHKERRQ(ierr);
1658 ierr = MatDestroy(&(sp->allMat));CHKERRQ(ierr);
1659 sp->allMat = allMat;
1660 ierr = PetscQuadratureDestroy(&(sp->allNodes));CHKERRQ(ierr);
1661 sp->allNodes = allNodes;
1662 PetscFunctionReturn(0);
1663 }
1664
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 */
PetscDualSpaceComputeFunctionalsFromAllData(PetscDualSpace sp)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;
1682
1683 PetscFunctionBegin;
1684 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr);
1685 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1686 ierr = PetscDualSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr);
1687 ierr = PetscDualSpaceGetFormDegree(sp, &k);CHKERRQ(ierr);
1688 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr);
1689 ierr = PetscDualSpaceGetAllData(sp, &allNodes, &allMat);CHKERRQ(ierr);
1690 nNodes = 0;
1691 if (allNodes) {
1692 ierr = PetscQuadratureGetData(allNodes, NULL, NULL, &nNodes, &nodes, NULL);CHKERRQ(ierr);
1693 }
1694 ierr = MatGetSize(allMat, &nDofs, NULL);CHKERRQ(ierr);
1695 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr);
1696 ierr = PetscSectionGetStorageSize(section, &spdim);CHKERRQ(ierr);
1697 if (spdim != nDofs) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "incompatible all matrix size");
1698 ierr = PetscMalloc1(nDofs, &(sp->functional));CHKERRQ(ierr);
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;
1707
1708 ierr = MatGetRow(allMat, f, &ncols, &cols, &vals);CHKERRQ(ierr);
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 ierr = PetscMalloc1(dim * nNodesf, &nodesf);CHKERRQ(ierr);
1714 ierr = PetscMalloc1(Nc * nNodesf, &weightsf);CHKERRQ(ierr);
1715 for (c = 0, countNodes = 0; c < ncols; c++) {
1716 if (!c || ((cols[c] / Nc) != (cols[c-1] / Nc))) {
1717 PetscInt d;
1718
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 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &(sp->functional[f]));CHKERRQ(ierr);
1730 ierr = PetscQuadratureSetData(sp->functional[f], dim, Nc, nNodesf, nodesf, weightsf);CHKERRQ(ierr);
1731 ierr = MatRestoreRow(allMat, f, &ncols, &cols, &vals);CHKERRQ(ierr);
1732 }
1733 PetscFunctionReturn(0);
1734 }
1735
1736 /* take a matrix meant for k-forms and expand it to one for Ncopies */
PetscDualSpaceLagrangeMatrixCreateCopies(Mat A,PetscInt Nk,PetscInt Ncopies,Mat * Abs)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;
1742 PetscErrorCode ierr;
1743
1744 PetscFunctionBegin;
1745 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr);
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 ierr = PetscMalloc1(m * Ncopies, &nnz);CHKERRQ(ierr);
1748 for (i = 0, maxnnz = 0; i < m; i++) {
1749 PetscInt innz;
1750 ierr = MatGetRow(A, i, &innz, NULL, NULL);CHKERRQ(ierr);
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 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, m * Ncopies, n * Ncopies, 0, nnz, &Ac);CHKERRQ(ierr);
1756 ierr = MatSetOption(Ac, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr);
1757 ierr = PetscFree(nnz);CHKERRQ(ierr);
1758 ierr = PetscMalloc1(maxnnz, &iwork);CHKERRQ(ierr);
1759 for (i = 0; i < m; i++) {
1760 PetscInt innz;
1761 const PetscInt *cols;
1762 const PetscScalar *vals;
1763
1764 ierr = MatGetRow(A, i, &innz, &cols, &vals);CHKERRQ(ierr);
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;
1768
1769 ierr = MatSetValues(Ac, 1, &row, innz, iwork, vals, INSERT_VALUES);CHKERRQ(ierr);
1770 for (k = 0; k < innz; k++) iwork[k] += Nk;
1771 }
1772 ierr = MatRestoreRow(A, i, &innz, &cols, &vals);CHKERRQ(ierr);
1773 }
1774 ierr = PetscFree(iwork);CHKERRQ(ierr);
1775 ierr = MatAssemblyBegin(Ac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1776 ierr = MatAssemblyEnd(Ac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1777 *Abs = Ac;
1778 PetscFunctionReturn(0);
1779 }
1780
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) */
DMPlexPointIsTensor_Internal_Given(DM dm,PetscInt p,PetscInt f,PetscInt f2,PetscBool * isTensor)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;
1794
1795 PetscFunctionBegin;
1796 fs[0] = f;
1797 fs[1] = f2;
1798 ierr = DMPlexGetMeet(dm, 2, fs, &meetSize, &meet);CHKERRQ(ierr);
1799 nmeet = meetSize;
1800 ierr = DMPlexRestoreMeet(dm, 2, fs, &meetSize, &meet);CHKERRQ(ierr);
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 PetscFunctionReturn(0);
1805 }
1806 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
1807 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
1808 ierr = DMPlexGetCone(dm, f, &fCone);CHKERRQ(ierr);
1809 ierr = DMPlexGetCone(dm, f2, &f2Cone);CHKERRQ(ierr);
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;
1818
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 ierr = DMPlexGetConeSize(dm, t, &tConeSize);CHKERRQ(ierr);
1823 ierr = DMPlexGetCone(dm, t, &tCone);CHKERRQ(ierr);
1824
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 PetscFunctionReturn(0);
1834 }
1835 d = q;
1836 dcount++;
1837 } else if (f2Cone[ef] == q) {
1838 if (d2count) {
1839 *isTensor = PETSC_FALSE;
1840 PetscFunctionReturn(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 ierr = DMPlexPointIsTensor_Internal_Given(dm, t, d, d2, &tIsTensor);CHKERRQ(ierr);
1850 if (!tIsTensor) {
1851 *isTensor = PETSC_FALSE;
1852 PetscFunctionReturn(0);
1853 }
1854 }
1855 *isTensor = PETSC_TRUE;
1856 PetscFunctionReturn(0);
1857 }
1858
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 */
DMPlexPointIsTensor_Internal(DM dm,PetscInt p,PetscBool * isTensor,PetscInt * endA,PetscInt * endB)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;
1866
1867 PetscFunctionBegin;
1868 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr);
1869 if (!coneSize) {
1870 if (isTensor) *isTensor = PETSC_FALSE;
1871 if (endA) *endA = -1;
1872 if (endB) *endB = -1;
1873 }
1874 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr);
1875 for (c = 0; c < coneSize; c++) {
1876 PetscInt f = cone[c];
1877 PetscInt fConeSize;
1878
1879 ierr = DMPlexGetConeSize(dm, f, &fConeSize);CHKERRQ(ierr);
1880 if (fConeSize != coneSize - 2) continue;
1881
1882 for (c2 = c + 1; c2 < coneSize; c2++) {
1883 PetscInt f2 = cone[c2];
1884 PetscBool isTensorff2;
1885 PetscInt f2ConeSize;
1886
1887 ierr = DMPlexGetConeSize(dm, f2, &f2ConeSize);CHKERRQ(ierr);
1888 if (f2ConeSize != coneSize - 2) continue;
1889
1890 ierr = DMPlexPointIsTensor_Internal_Given(dm, p, f, f2, &isTensorff2);CHKERRQ(ierr);
1891 if (isTensorff2) {
1892 if (isTensor) *isTensor = PETSC_TRUE;
1893 if (endA) *endA = f;
1894 if (endB) *endB = f2;
1895 PetscFunctionReturn(0);
1896 }
1897 }
1898 }
1899 if (isTensor) *isTensor = PETSC_FALSE;
1900 if (endA) *endA = -1;
1901 if (endB) *endB = -1;
1902 PetscFunctionReturn(0);
1903 }
1904
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 */
DMPlexPointIsTensor(DM dm,PetscInt p,PetscBool * isTensor,PetscInt * endA,PetscInt * endB)1907 static PetscErrorCode DMPlexPointIsTensor(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB)
1908 {
1909 DMPlexInterpolatedFlag interpolated;
1910 PetscErrorCode ierr;
1911
1912 PetscFunctionBegin;
1913 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr);
1914 if (interpolated != DMPLEX_INTERPOLATED_FULL) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Only for interpolated DMPlex's");
1915 ierr = DMPlexPointIsTensor_Internal(dm, p, isTensor, endA, endB);CHKERRQ(ierr);
1916 PetscFunctionReturn(0);
1917 }
1918
1919 /* permute a quadrature -> dof matrix so that its rows are in revlex order by nodeIdx */
MatPermuteByNodeIdx(Mat A,PetscLagNodeIndices ni,Mat * Aperm)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;
1930 PetscErrorCode ierr;
1931
1932 PetscFunctionBegin;
1933 ierr = PetscLagNodeIndicesGetPermutation(ni, &perm);CHKERRQ(ierr);
1934 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr);
1935 ierr = PetscMalloc1(nodeIdxDim * m, &nIdxPerm);CHKERRQ(ierr);
1936 ierr = PetscMalloc1(nodeVecDim * m, &nVecPerm);CHKERRQ(ierr);
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 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, perm, PETSC_USE_POINTER, &permIS);CHKERRQ(ierr);
1940 ierr = ISSetPermutation(permIS);CHKERRQ(ierr);
1941 ierr = ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &id);CHKERRQ(ierr);
1942 ierr = ISSetPermutation(id);CHKERRQ(ierr);
1943 ierr = MatPermute(A, permIS, id, Aperm);CHKERRQ(ierr);
1944 ierr = ISDestroy(&permIS);CHKERRQ(ierr);
1945 ierr = ISDestroy(&id);CHKERRQ(ierr);
1946 for (i = 0; i < m; i++) perm[i] = i;
1947 ierr = PetscFree(ni->nodeIdx);CHKERRQ(ierr);
1948 ierr = PetscFree(ni->nodeVec);CHKERRQ(ierr);
1949 ni->nodeIdx = nIdxPerm;
1950 ni->nodeVec = nVecPerm;
1951 PetscFunctionReturn(0);
1952 }
1953
PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)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;
1974
1975 PetscFunctionBegin;
1976 /* step 1: sanitize input */
1977 ierr = PetscObjectGetComm((PetscObject) sp, &comm);CHKERRQ(ierr);
1978 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
1979 ierr = PetscObjectTypeCompare((PetscObject)sp, PETSCDUALSPACEBDM, &isbdm);CHKERRQ(ierr);
1980 if (isbdm) {
1981 sp->k = -(dim-1); /* form degree of H-div */
1982 ierr = PetscObjectChangeTypeName((PetscObject)sp, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr);
1983 }
1984 ierr = PetscDualSpaceGetFormDegree(sp, &formDegree);CHKERRQ(ierr);
1985 if (PetscAbsInt(formDegree) > dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Form degree must be bounded by dimension");
1986 ierr = PetscDTBinomialInt(dim,PetscAbsInt(formDegree),&Nk);CHKERRQ(ierr);
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 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
2017 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
2018 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
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 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr);
2022 if (interpolated != DMPLEX_INTERPOLATED_FULL) {
2023 ierr = DMPlexInterpolate(dm, &dmint);CHKERRQ(ierr);
2024 } else {
2025 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);
2026 dmint = dm;
2027 }
2028 tensorCell = PETSC_FALSE;
2029 if (dim > 1) {
2030 ierr = DMPlexPointIsTensor(dmint, 0, &tensorCell, &tensorf, &tensorf2);CHKERRQ(ierr);
2031 }
2032 lag->tensorCell = tensorCell;
2033 if (dim < 2 || !lag->tensorCell) lag->tensorSpace = PETSC_FALSE;
2034 tensorSpace = lag->tensorSpace;
2035 if (!lag->nodeFamily) {
2036 ierr = Petsc1DNodeFamilyCreate(lag->nodeType, lag->nodeExponent, lag->endNodes, &lag->nodeFamily);CHKERRQ(ierr);
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");
2040
2041 /* step 2: construct the boundary spaces */
2042 ierr = PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);CHKERRQ(ierr);
2043 ierr = PetscCalloc1(pEnd,&(sp->pointSpaces));CHKERRQ(ierr);
2044 for (d = 0; d <= depth; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);CHKERRQ(ierr);}
2045 ierr = PetscDualSpaceSectionCreate_Internal(sp, §ion);CHKERRQ(ierr);
2046 sp->pointSection = section;
2047 if (continuous && !(lag->interiorOnly)) {
2048 PetscInt h;
2049
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;
2055
2056 ierr = DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, NULL, &detJ);CHKERRQ(ierr);
2057 ierr = DMPlexGetCellType(dm, p, &ptype);CHKERRQ(ierr);
2058
2059 /* compare to previous facets: if computed, reference that dualspace */
2060 for (q = pStratStart[depth - 1]; q < p; q++) {
2061 DMPolytopeType qtype;
2062
2063 ierr = DMPlexGetCellType(dm, q, &qtype);CHKERRQ(ierr);
2064 if (qtype == ptype) break;
2065 }
2066 if (q < p) { /* this facet has the same dual space as that one */
2067 ierr = PetscObjectReference((PetscObject)sp->pointSpaces[q]);CHKERRQ(ierr);
2068 sp->pointSpaces[p] = sp->pointSpaces[q];
2069 continue;
2070 }
2071 /* if not, recursively compute this dual space */
2072 ierr = PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,p,formDegree,Ncopies,PETSC_FALSE,&sp->pointSpaces[p]);CHKERRQ(ierr);
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;
2077
2078 if (hdim < PetscAbsInt(formDegree)) break;
2079 for (p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
2080 PetscInt suppSize, s;
2081 const PetscInt *supp;
2082
2083 ierr = DMPlexGetSupportSize(dm, p, &suppSize);CHKERRQ(ierr);
2084 ierr = DMPlexGetSupport(dm, p, &supp);CHKERRQ(ierr);
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;
2091
2092 q = supp[0];
2093 qsp = sp->pointSpaces[q];
2094 ierr = DMPlexGetConeSize(dm, q, &coneSize);CHKERRQ(ierr);
2095 ierr = DMPlexGetCone(dm, q, &cone);CHKERRQ(ierr);
2096 for (c = 0; c < coneSize; c++) if (cone[c] == p) break;
2097 if (c == coneSize) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "cone/support mismatch");
2098 ierr = PetscDualSpaceGetDM(qsp, &qdm);CHKERRQ(ierr);
2099 ierr = DMPlexGetCone(qdm, 0, &refCone);CHKERRQ(ierr);
2100 /* get the equivalent dual space from the support dual space */
2101 ierr = PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp);CHKERRQ(ierr);
2102 if (!s) {
2103 ierr = PetscObjectReference((PetscObject)psp);CHKERRQ(ierr);
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 ierr = PetscDualSpaceGetInteriorDimension(sp->pointSpaces[p], &pspdim);CHKERRQ(ierr);
2113 ierr = PetscSectionSetDof(section, p, pspdim);CHKERRQ(ierr);
2114 }
2115 }
2116
2117 if (Ncopies > 1) {
2118 Mat intMatScalar, allMatScalar;
2119 PetscDualSpace scalarsp;
2120 PetscDualSpace_Lag *scalarlag;
2121
2122 ierr = PetscDualSpaceDuplicate(sp, &scalarsp);CHKERRQ(ierr);
2123 /* Setting the number of components to Nk is a space with 1 copy of each k-form */
2124 ierr = PetscDualSpaceSetNumComponents(scalarsp, Nk);CHKERRQ(ierr);
2125 ierr = PetscDualSpaceSetUp(scalarsp);CHKERRQ(ierr);
2126 ierr = PetscDualSpaceGetInteriorData(scalarsp, &(sp->intNodes), &intMatScalar);CHKERRQ(ierr);
2127 ierr = PetscObjectReference((PetscObject)(sp->intNodes));CHKERRQ(ierr);
2128 if (intMatScalar) {ierr = PetscDualSpaceLagrangeMatrixCreateCopies(intMatScalar, Nk, Ncopies, &(sp->intMat));CHKERRQ(ierr);}
2129 ierr = PetscDualSpaceGetAllData(scalarsp, &(sp->allNodes), &allMatScalar);CHKERRQ(ierr);
2130 ierr = PetscObjectReference((PetscObject)(sp->allNodes));CHKERRQ(ierr);
2131 ierr = PetscDualSpaceLagrangeMatrixCreateCopies(allMatScalar, Nk, Ncopies, &(sp->allMat));CHKERRQ(ierr);
2132 sp->spdim = scalarsp->spdim * Ncopies;
2133 sp->spintdim = scalarsp->spintdim * Ncopies;
2134 scalarlag = (PetscDualSpace_Lag *) scalarsp->data;
2135 ierr = PetscLagNodeIndicesReference(scalarlag->vertIndices);CHKERRQ(ierr);
2136 lag->vertIndices = scalarlag->vertIndices;
2137 ierr = PetscLagNodeIndicesReference(scalarlag->intNodeIndices);CHKERRQ(ierr);
2138 lag->intNodeIndices = scalarlag->intNodeIndices;
2139 ierr = PetscLagNodeIndicesReference(scalarlag->allNodeIndices);CHKERRQ(ierr);
2140 lag->allNodeIndices = scalarlag->allNodeIndices;
2141 ierr = PetscDualSpaceDestroy(&scalarsp);CHKERRQ(ierr);
2142 ierr = PetscSectionSetDof(section, 0, sp->spintdim);CHKERRQ(ierr);
2143 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr);
2144 ierr = PetscDualSpaceComputeFunctionalsFromAllData(sp);CHKERRQ(ierr);
2145 ierr = PetscFree2(pStratStart, pStratEnd);CHKERRQ(ierr);
2146 ierr = DMDestroy(&dmint);CHKERRQ(ierr);
2147 PetscFunctionReturn(0);
2148 }
2149
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;
2159
2160 ierr = PetscDualSpaceDuplicate(sp, &spcont);CHKERRQ(ierr);
2161 ierr = PetscDualSpaceLagrangeSetContinuity(spcont, PETSC_TRUE);CHKERRQ(ierr);
2162 ierr = PetscDualSpaceSetUp(spcont);CHKERRQ(ierr);
2163 ierr = PetscDualSpaceGetDimension(spcont, &spdim);CHKERRQ(ierr);
2164 sp->spdim = sp->spintdim = spdim;
2165 ierr = PetscSectionSetDof(section, 0, spdim);CHKERRQ(ierr);
2166 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr);
2167 ierr = PetscMalloc1(spdim, &(sp->functional));CHKERRQ(ierr);
2168 for (f = 0; f < spdim; f++) {
2169 PetscQuadrature fn;
2170
2171 ierr = PetscDualSpaceGetFunctional(spcont, f, &fn);CHKERRQ(ierr);
2172 ierr = PetscObjectReference((PetscObject)fn);CHKERRQ(ierr);
2173 sp->functional[f] = fn;
2174 }
2175 ierr = PetscDualSpaceGetAllData(spcont, &allNodes, &allMat);CHKERRQ(ierr);
2176 ierr = PetscObjectReference((PetscObject) allNodes);CHKERRQ(ierr);
2177 ierr = PetscObjectReference((PetscObject) allNodes);CHKERRQ(ierr);
2178 sp->allNodes = sp->intNodes = allNodes;
2179 ierr = PetscObjectReference((PetscObject) allMat);CHKERRQ(ierr);
2180 ierr = PetscObjectReference((PetscObject) allMat);CHKERRQ(ierr);
2181 sp->allMat = sp->intMat = allMat;
2182 lagc = (PetscDualSpace_Lag *) spcont->data;
2183 ierr = PetscLagNodeIndicesReference(lagc->vertIndices);CHKERRQ(ierr);
2184 lag->vertIndices = lagc->vertIndices;
2185 ierr = PetscLagNodeIndicesReference(lagc->allNodeIndices);CHKERRQ(ierr);
2186 ierr = PetscLagNodeIndicesReference(lagc->allNodeIndices);CHKERRQ(ierr);
2187 lag->intNodeIndices = lagc->allNodeIndices;
2188 lag->allNodeIndices = lagc->allNodeIndices;
2189 ierr = PetscDualSpaceDestroy(&spcont);CHKERRQ(ierr);
2190 ierr = PetscFree2(pStratStart, pStratEnd);CHKERRQ(ierr);
2191 ierr = DMDestroy(&dmint);CHKERRQ(ierr);
2192 PetscFunctionReturn(0);
2193 }
2194
2195 /* step 3: construct intNodes, and intMat, and combine it with boundray data to make allNodes and allMat */
2196 if (!tensorSpace) {
2197 if (!tensorCell) {ierr = PetscLagNodeIndicesCreateSimplexVertices(dm, &(lag->vertIndices));CHKERRQ(ierr);}
2198
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;
2205
2206 ierr = PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));CHKERRQ(ierr);
2207 ierr = MatGetSize(sp->intMat, &nDofs, NULL);CHKERRQ(ierr);
2208 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr);
2209 }
2210 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr);
2211 ierr = PetscDualSpaceCreateAllDataFromInteriorData(sp);CHKERRQ(ierr);
2212 ierr = PetscDualSpaceLagrangeCreateAllNodeIdx(sp);CHKERRQ(ierr);
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;
2219
2220 ierr = PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));CHKERRQ(ierr);
2221 ierr = MatGetSize(sp->intMat, &nDofs, NULL);CHKERRQ(ierr);
2222 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr);
2223 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr);
2224 ierr = PetscObjectReference((PetscObject)(sp->intNodes));CHKERRQ(ierr);
2225 sp->allNodes = sp->intNodes;
2226 ierr = PetscObjectReference((PetscObject)(sp->intMat));CHKERRQ(ierr);
2227 sp->allMat = sp->intMat;
2228 ierr = PetscLagNodeIndicesReference(lag->intNodeIndices);CHKERRQ(ierr);
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;
2240
2241 ierr = PetscDualSpaceDuplicate(sp, &trimmedsp);CHKERRQ(ierr);
2242 ierr = PetscDualSpaceLagrangeSetTrimmed(trimmedsp, PETSC_TRUE);CHKERRQ(ierr);
2243 ierr = PetscDualSpaceSetOrder(trimmedsp, order + PetscAbsInt(formDegree) - dim);CHKERRQ(ierr);
2244 ierr = PetscDualSpaceSetFormDegree(trimmedsp, trFormDegree);CHKERRQ(ierr);
2245 trimmedlag = (PetscDualSpace_Lag *) trimmedsp->data;
2246 trimmedlag->numNodeSkip = numNodeSkip + 1;
2247 ierr = PetscDualSpaceSetUp(trimmedsp);CHKERRQ(ierr);
2248 ierr = PetscDualSpaceGetAllData(trimmedsp, &intNodes, &intMat);CHKERRQ(ierr);
2249 ierr = PetscObjectReference((PetscObject)intNodes);CHKERRQ(ierr);
2250 sp->intNodes = intNodes;
2251 ierr = PetscObjectReference((PetscObject)intMat);CHKERRQ(ierr);
2252 sp->intMat = intMat;
2253 ierr = MatGetSize(sp->intMat, &nDofs, NULL);CHKERRQ(ierr);
2254 ierr = PetscLagNodeIndicesReference(trimmedlag->allNodeIndices);CHKERRQ(ierr);
2255 lag->intNodeIndices = trimmedlag->allNodeIndices;
2256 ierr = PetscDualSpaceDestroy(&trimmedsp);CHKERRQ(ierr);
2257 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr);
2258 }
2259 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr);
2260 ierr = PetscDualSpaceCreateAllDataFromInteriorData(sp);CHKERRQ(ierr);
2261 ierr = PetscDualSpaceLagrangeCreateAllNodeIdx(sp);CHKERRQ(ierr);
2262 }
2263 }
2264 } else {
2265 PetscQuadrature intNodesTrace = NULL;
2266 PetscQuadrature intNodesFiber = NULL;
2267 PetscQuadrature intNodes = NULL;
2268 PetscLagNodeIndices intNodeIndices = NULL;
2269 Mat intMat = NULL;
2270
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;
2276
2277 if (sp->pointSpaces[tensorf]) {
2278 ierr = PetscObjectReference((PetscObject)(sp->pointSpaces[tensorf]));CHKERRQ(ierr);
2279 trace = sp->pointSpaces[tensorf];
2280 } else {
2281 ierr = PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,formDegree,Ncopies,PETSC_TRUE,&trace);CHKERRQ(ierr);
2282 }
2283 ierr = PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,0,1,PETSC_TRUE,&fiber);CHKERRQ(ierr);
2284 tracel = (PetscDualSpace_Lag *) trace->data;
2285 fiberl = (PetscDualSpace_Lag *) fiber->data;
2286 ierr = PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));CHKERRQ(ierr);
2287 ierr = PetscDualSpaceGetInteriorData(trace, &intNodesTrace, &intMatTrace);CHKERRQ(ierr);
2288 ierr = PetscDualSpaceGetInteriorData(fiber, &intNodesFiber, &intMatFiber);CHKERRQ(ierr);
2289 if (intNodesTrace && intNodesFiber) {
2290 ierr = PetscQuadratureCreateTensor(intNodesTrace, intNodesFiber, &intNodes);CHKERRQ(ierr);
2291 ierr = MatTensorAltV(intMatTrace, intMatFiber, dim-1, formDegree, 1, 0, &intMat);CHKERRQ(ierr);
2292 ierr = PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, formDegree, fiberl->intNodeIndices, 1, 0, &intNodeIndices);CHKERRQ(ierr);
2293 }
2294 ierr = PetscObjectReference((PetscObject) intNodesTrace);CHKERRQ(ierr);
2295 ierr = PetscObjectReference((PetscObject) intNodesFiber);CHKERRQ(ierr);
2296 ierr = PetscDualSpaceDestroy(&fiber);CHKERRQ(ierr);
2297 ierr = PetscDualSpaceDestroy(&trace);CHKERRQ(ierr);
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;
2308
2309 ierr = PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,traceDegree,Ncopies,PETSC_TRUE,&trace);CHKERRQ(ierr);
2310 ierr = PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,fiberDegree,1,PETSC_TRUE,&fiber);CHKERRQ(ierr);
2311 tracel = (PetscDualSpace_Lag *) trace->data;
2312 fiberl = (PetscDualSpace_Lag *) fiber->data;
2313 if (!lag->vertIndices) {
2314 ierr = PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));CHKERRQ(ierr);
2315 }
2316 ierr = PetscDualSpaceGetInteriorData(trace, &intNodesTrace2, &intMatTrace);CHKERRQ(ierr);
2317 ierr = PetscDualSpaceGetInteriorData(fiber, &intNodesFiber2, &intMatFiber);CHKERRQ(ierr);
2318 if (intNodesTrace2 && intNodesFiber2) {
2319 ierr = PetscQuadratureCreateTensor(intNodesTrace2, intNodesFiber2, &intNodes2);CHKERRQ(ierr);
2320 ierr = MatTensorAltV(intMatTrace, intMatFiber, dim-1, traceDegree, 1, fiberDegree, &intMat2);CHKERRQ(ierr);
2321 ierr = PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, traceDegree, fiberl->intNodeIndices, 1, fiberDegree, &intNodeIndices2);CHKERRQ(ierr);
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;
2334
2335 ierr = MatGetSize(intMat, &nDof, NULL);CHKERRQ(ierr);
2336 ierr = MatGetSize(intMat2, &nDof2, NULL);CHKERRQ(ierr);
2337 ierr = PetscQuadraturePointsMerge(intNodes, intNodes2, &merged, &toMerged, &toMerged2);CHKERRQ(ierr);
2338 ierr = PetscQuadratureGetData(merged, NULL, NULL, &nM, NULL, NULL);CHKERRQ(ierr);
2339 ierr = MatricesMerge(intMat, intMat2, dim, formDegree, nM, toMerged, toMerged2, &matMerged);CHKERRQ(ierr);
2340 ierr = PetscLagNodeIndicesMerge(intNodeIndices, intNodeIndices2, &intNodeIndicesMerged);CHKERRQ(ierr);
2341 ierr = PetscFree(toMerged);CHKERRQ(ierr);
2342 ierr = PetscFree(toMerged2);CHKERRQ(ierr);
2343 ierr = MatDestroy(&intMat);CHKERRQ(ierr);
2344 ierr = MatDestroy(&intMat2);CHKERRQ(ierr);
2345 ierr = PetscQuadratureDestroy(&intNodes);CHKERRQ(ierr);
2346 ierr = PetscQuadratureDestroy(&intNodes2);CHKERRQ(ierr);
2347 ierr = PetscLagNodeIndicesDestroy(&intNodeIndices);CHKERRQ(ierr);
2348 ierr = PetscLagNodeIndicesDestroy(&intNodeIndices2);CHKERRQ(ierr);
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;
2358
2359 ierr = MatPermuteByNodeIdx(intMat, intNodeIndices, &intMatPerm);CHKERRQ(ierr);
2360 ierr = MatDestroy(&intMat);CHKERRQ(ierr);
2361 intMat = intMatPerm;
2362 }
2363 }
2364 }
2365 ierr = PetscDualSpaceDestroy(&fiber);CHKERRQ(ierr);
2366 ierr = PetscDualSpaceDestroy(&trace);CHKERRQ(ierr);
2367 }
2368 ierr = PetscQuadratureDestroy(&intNodesTrace);CHKERRQ(ierr);
2369 ierr = PetscQuadratureDestroy(&intNodesFiber);CHKERRQ(ierr);
2370 sp->intNodes = intNodes;
2371 sp->intMat = intMat;
2372 lag->intNodeIndices = intNodeIndices;
2373 {
2374 PetscInt nDofs = 0;
2375
2376 if (intMat) {
2377 ierr = MatGetSize(intMat, &nDofs, NULL);CHKERRQ(ierr);
2378 }
2379 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr);
2380 }
2381 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr);
2382 if (continuous) {
2383 ierr = PetscDualSpaceCreateAllDataFromInteriorData(sp);CHKERRQ(ierr);
2384 ierr = PetscDualSpaceLagrangeCreateAllNodeIdx(sp);CHKERRQ(ierr);
2385 } else {
2386 ierr = PetscObjectReference((PetscObject) intNodes);CHKERRQ(ierr);
2387 sp->allNodes = intNodes;
2388 ierr = PetscObjectReference((PetscObject) intMat);CHKERRQ(ierr);
2389 sp->allMat = intMat;
2390 ierr = PetscLagNodeIndicesReference(intNodeIndices);CHKERRQ(ierr);
2391 lag->allNodeIndices = intNodeIndices;
2392 }
2393 }
2394 ierr = PetscSectionGetStorageSize(section, &sp->spdim);CHKERRQ(ierr);
2395 ierr = PetscSectionGetConstrainedStorageSize(section, &sp->spintdim);CHKERRQ(ierr);
2396 ierr = PetscDualSpaceComputeFunctionalsFromAllData(sp);CHKERRQ(ierr);
2397 ierr = PetscFree2(pStratStart, pStratEnd);CHKERRQ(ierr);
2398 ierr = DMDestroy(&dmint);CHKERRQ(ierr);
2399 PetscFunctionReturn(0);
2400 }
2401
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 */
PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp,PetscInt ornt,Mat * symMat)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;
2419 PetscErrorCode ierr;
2420
2421 PetscFunctionBegin;
2422 if (!sp->spintdim) {
2423 *symMat = NULL;
2424 PetscFunctionReturn(0);
2425 }
2426 lag = (PetscDualSpace_Lag *) sp->data;
2427 vertIndices = lag->vertIndices;
2428 intNodeIndices = lag->intNodeIndices;
2429 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr);
2430 ierr = PetscDualSpaceGetFormDegree(sp, &formDegree);CHKERRQ(ierr);
2431 ierr = PetscNew(&ni);CHKERRQ(ierr);
2432 ni->refct = 1;
2433 ni->nodeIdxDim = nodeIdxDim = intNodeIndices->nodeIdxDim;
2434 ni->nodeVecDim = nodeVecDim = intNodeIndices->nodeVecDim;
2435 ni->nNodes = nNodes = intNodeIndices->nNodes;
2436 ierr = PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));CHKERRQ(ierr);
2437 ierr = PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));CHKERRQ(ierr);
2438 /* push forward the dofs by the symmetry of the reference element induced by ornt */
2439 ierr = PetscLagNodeIndicesPushForward(dm, vertIndices, 0, vertIndices, intNodeIndices, ornt, formDegree, ni->nodeIdx, ni->nodeVec);CHKERRQ(ierr);
2440 /* get the revlex order for both the original and transformed dofs */
2441 ierr = PetscLagNodeIndicesGetPermutation(intNodeIndices, &perm);CHKERRQ(ierr);
2442 ierr = PetscLagNodeIndicesGetPermutation(ni, &permOrnt);CHKERRQ(ierr);
2443 ierr = PetscMalloc1(nNodes, &nnz);CHKERRQ(ierr);
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;
2452
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
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 if (PetscDefined(USE_DEBUG)) {
2462 PetscInt m;
2463 PetscInt *nind = &(intNodeIndices->nodeIdx[perm[n] * nodeIdxDim]);
2464
2465 for (m = n + 1; m < nEnd; m++) {
2466 PetscInt *mind = &(intNodeIndices->nodeIdx[perm[m] * nodeIdxDim]);
2467 PetscInt d;
2468
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 groupSize = nEnd - n;
2476 /* each pushforward dof vector will be expressed in a basis of the unpermuted dofs */
2477 for (m = n; m < nEnd; m++) nnz[permOrnt[m]] = groupSize;
2478
2479 maxGroupSize = PetscMax(maxGroupSize, nEnd - n);
2480 n = nEnd;
2481 }
2482 if (maxGroupSize > nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs are not in blocks that can be solved");
2483 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes, nNodes, 0, nnz, &A);CHKERRQ(ierr);
2484 ierr = PetscFree(nnz);CHKERRQ(ierr);
2485 ierr = PetscMalloc3(maxGroupSize * nodeVecDim, &V, maxGroupSize * nodeVecDim, &W, nodeVecDim * 2, &work);CHKERRQ(ierr);
2486 for (n = 0; n < nNodes;) { /* incremented in the loop */
2487 PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]);
2488 PetscInt nEnd;
2489 PetscInt m;
2490 PetscInt groupSize;
2491 for (nEnd = n + 1; nEnd < nNodes; nEnd++) {
2492 PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]);
2493 PetscInt d;
2494
2495 /* compare the oriented permutation indices */
2496 for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break;
2497 if (d < nodeIdxDim) break;
2498 }
2499 groupSize = nEnd - n;
2500 /* get all of the vectors from the original and all of the pushforward vectors */
2501 for (m = n; m < nEnd; m++) {
2502 PetscInt d;
2503
2504 for (d = 0; d < nodeVecDim; d++) {
2505 V[(m - n) * nodeVecDim + d] = intNodeIndices->nodeVec[perm[m] * nodeVecDim + d];
2506 W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d];
2507 }
2508 }
2509 /* now we have to solve for W in terms of V: the systems isn't always square, but the span
2510 * of V and W should always be the same, so the solution of the normal equations works */
2511 {
2512 char transpose = 'N';
2513 PetscBLASInt bm = nodeVecDim;
2514 PetscBLASInt bn = groupSize;
2515 PetscBLASInt bnrhs = groupSize;
2516 PetscBLASInt blda = bm;
2517 PetscBLASInt bldb = bm;
2518 PetscBLASInt blwork = 2 * nodeVecDim;
2519 PetscBLASInt info;
2520
2521 PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&bm,&bn,&bnrhs,V,&blda,W,&bldb,work,&blwork, &info));
2522 if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS");
2523 /* repack */
2524 {
2525 PetscInt i, j;
2526
2527 for (i = 0; i < groupSize; i++) {
2528 for (j = 0; j < groupSize; j++) {
2529 /* notice the different leading dimension */
2530 V[i * groupSize + j] = W[i * nodeVecDim + j];
2531 }
2532 }
2533 }
2534 }
2535 ierr = MatSetValues(A, groupSize, &permOrnt[n], groupSize, &perm[n], V, INSERT_VALUES);CHKERRQ(ierr);
2536 n = nEnd;
2537 }
2538 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2539 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
2540 *symMat = A;
2541 ierr = PetscFree3(V,W,work);CHKERRQ(ierr);
2542 ierr = PetscLagNodeIndicesDestroy(&ni);CHKERRQ(ierr);
2543 PetscFunctionReturn(0);
2544 }
2545
2546 #define BaryIndex(perEdge,a,b,c) (((b)*(2*perEdge+1-(b)))/2)+(c)
2547
2548 #define CartIndex(perEdge,a,b) (perEdge*(a)+b)
2549
2550 /* the existing interface for symmetries is insufficient for all cases:
2551 * - it should be sufficient for form degrees that are scalar (0 and n)
2552 * - it should be sufficient for hypercube dofs
2553 * - it isn't sufficient for simplex cells with non-scalar form degrees if
2554 * there are any dofs in the interior
2555 *
2556 * We compute the general transformation matrices, and if they fit, we return them,
2557 * otherwise we error (but we should probably change the interface to allow for
2558 * these symmetries)
2559 */
PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp,const PetscInt **** perms,const PetscScalar **** flips)2560 static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
2561 {
2562 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2563 PetscInt dim, order, Nc;
2564 PetscErrorCode ierr;
2565
2566 PetscFunctionBegin;
2567 ierr = PetscDualSpaceGetOrder(sp,&order);CHKERRQ(ierr);
2568 ierr = PetscDualSpaceGetNumComponents(sp,&Nc);CHKERRQ(ierr);
2569 ierr = DMGetDimension(sp->dm,&dim);CHKERRQ(ierr);
2570 if (!lag->symComputed) { /* store symmetries */
2571 PetscInt pStart, pEnd, p;
2572 PetscInt numPoints;
2573 PetscInt numFaces;
2574 PetscInt spintdim;
2575 PetscInt ***symperms;
2576 PetscScalar ***symflips;
2577
2578 ierr = DMPlexGetChart(sp->dm, &pStart, &pEnd);CHKERRQ(ierr);
2579 numPoints = pEnd - pStart;
2580 ierr = DMPlexGetConeSize(sp->dm, 0, &numFaces);CHKERRQ(ierr);
2581 ierr = PetscCalloc1(numPoints,&symperms);CHKERRQ(ierr);
2582 ierr = PetscCalloc1(numPoints,&symflips);CHKERRQ(ierr);
2583 spintdim = sp->spintdim;
2584 /* The nodal symmetry behavior is not present when tensorSpace != tensorCell: someone might want this for the "S"
2585 * family of FEEC spaces. Most used in particular are discontinuous polynomial L2 spaces in tensor cells, where
2586 * the symmetries are not necessary for FE assembly. So for now we assume this is the case and don't return
2587 * symmetries if tensorSpace != tensorCell */
2588 if (spintdim && 0 < dim && dim < 3 && (lag->tensorSpace == lag->tensorCell)) { /* compute self symmetries */
2589 PetscInt **cellSymperms;
2590 PetscScalar **cellSymflips;
2591 PetscInt ornt;
2592 PetscInt nCopies = Nc / lag->intNodeIndices->nodeVecDim;
2593 PetscInt nNodes = lag->intNodeIndices->nNodes;
2594
2595 lag->numSelfSym = 2 * numFaces;
2596 lag->selfSymOff = numFaces;
2597 ierr = PetscCalloc1(2*numFaces,&cellSymperms);CHKERRQ(ierr);
2598 ierr = PetscCalloc1(2*numFaces,&cellSymflips);CHKERRQ(ierr);
2599 /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
2600 symperms[0] = &cellSymperms[numFaces];
2601 symflips[0] = &cellSymflips[numFaces];
2602 if (lag->intNodeIndices->nodeVecDim * nCopies != Nc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs");
2603 if (nNodes * nCopies != spintdim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs");
2604 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 */
2605 Mat symMat;
2606 PetscInt *perm;
2607 PetscScalar *flips;
2608 PetscInt i;
2609
2610 if (!ornt) continue;
2611 ierr = PetscMalloc1(spintdim, &perm);CHKERRQ(ierr);
2612 ierr = PetscCalloc1(spintdim, &flips);CHKERRQ(ierr);
2613 for (i = 0; i < spintdim; i++) perm[i] = -1;
2614 ierr = PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, ornt, &symMat);CHKERRQ(ierr);
2615 for (i = 0; i < nNodes; i++) {
2616 PetscInt ncols;
2617 PetscInt j, k;
2618 const PetscInt *cols;
2619 const PetscScalar *vals;
2620 PetscBool nz_seen = PETSC_FALSE;
2621
2622 ierr = MatGetRow(symMat, i, &ncols, &cols, &vals);CHKERRQ(ierr);
2623 for (j = 0; j < ncols; j++) {
2624 if (PetscAbsScalar(vals[j]) > PETSC_SMALL) {
2625 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");
2626 nz_seen = PETSC_TRUE;
2627 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");
2628 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");
2629 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");
2630 for (k = 0; k < nCopies; k++) {
2631 perm[cols[j] * nCopies + k] = i * nCopies + k;
2632 }
2633 if (PetscRealPart(vals[j]) < 0.) {
2634 for (k = 0; k < nCopies; k++) {
2635 flips[i * nCopies + k] = -1.;
2636 }
2637 } else {
2638 for (k = 0; k < nCopies; k++) {
2639 flips[i * nCopies + k] = 1.;
2640 }
2641 }
2642 }
2643 }
2644 ierr = MatRestoreRow(symMat, i, &ncols, &cols, &vals);CHKERRQ(ierr);
2645 }
2646 ierr = MatDestroy(&symMat);CHKERRQ(ierr);
2647 /* if there were no sign flips, keep NULL */
2648 for (i = 0; i < spintdim; i++) if (flips[i] != 1.) break;
2649 if (i == spintdim) {
2650 ierr = PetscFree(flips);CHKERRQ(ierr);
2651 flips = NULL;
2652 }
2653 /* if the permutation is identity, keep NULL */
2654 for (i = 0; i < spintdim; i++) if (perm[i] != i) break;
2655 if (i == spintdim) {
2656 ierr = PetscFree(perm);CHKERRQ(ierr);
2657 perm = NULL;
2658 }
2659 symperms[0][ornt] = perm;
2660 symflips[0][ornt] = flips;
2661 }
2662 /* if no orientations produced non-identity permutations, keep NULL */
2663 for (ornt = -numFaces; ornt < numFaces; ornt++) if (symperms[0][ornt]) break;
2664 if (ornt == numFaces) {
2665 ierr = PetscFree(cellSymperms);CHKERRQ(ierr);
2666 symperms[0] = NULL;
2667 }
2668 /* if no orientations produced sign flips, keep NULL */
2669 for (ornt = -numFaces; ornt < numFaces; ornt++) if (symflips[0][ornt]) break;
2670 if (ornt == numFaces) {
2671 ierr = PetscFree(cellSymflips);CHKERRQ(ierr);
2672 symflips[0] = NULL;
2673 }
2674 }
2675 { /* get the symmetries of closure points */
2676 PetscInt closureSize = 0;
2677 PetscInt *closure = NULL;
2678 PetscInt r;
2679
2680 ierr = DMPlexGetTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
2681 for (r = 0; r < closureSize; r++) {
2682 PetscDualSpace psp;
2683 PetscInt point = closure[2 * r];
2684 PetscInt pspintdim;
2685 const PetscInt ***psymperms = NULL;
2686 const PetscScalar ***psymflips = NULL;
2687
2688 if (!point) continue;
2689 ierr = PetscDualSpaceGetPointSubspace(sp, point, &psp);CHKERRQ(ierr);
2690 if (!psp) continue;
2691 ierr = PetscDualSpaceGetInteriorDimension(psp, &pspintdim);CHKERRQ(ierr);
2692 if (!pspintdim) continue;
2693 ierr = PetscDualSpaceGetSymmetries(psp,&psymperms,&psymflips);CHKERRQ(ierr);
2694 symperms[r] = (PetscInt **) (psymperms ? psymperms[0] : NULL);
2695 symflips[r] = (PetscScalar **) (psymflips ? psymflips[0] : NULL);
2696 }
2697 ierr = DMPlexRestoreTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
2698 }
2699 for (p = 0; p < pEnd; p++) if (symperms[p]) break;
2700 if (p == pEnd) {
2701 ierr = PetscFree(symperms);CHKERRQ(ierr);
2702 symperms = NULL;
2703 }
2704 for (p = 0; p < pEnd; p++) if (symflips[p]) break;
2705 if (p == pEnd) {
2706 ierr = PetscFree(symflips);CHKERRQ(ierr);
2707 symflips = NULL;
2708 }
2709 lag->symperms = symperms;
2710 lag->symflips = symflips;
2711 lag->symComputed = PETSC_TRUE;
2712 }
2713 if (perms) *perms = (const PetscInt ***) lag->symperms;
2714 if (flips) *flips = (const PetscScalar ***) lag->symflips;
2715 PetscFunctionReturn(0);
2716 }
2717
PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp,PetscBool * continuous)2718 static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2719 {
2720 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2721
2722 PetscFunctionBegin;
2723 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
2724 PetscValidPointer(continuous, 2);
2725 *continuous = lag->continuous;
2726 PetscFunctionReturn(0);
2727 }
2728
PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp,PetscBool continuous)2729 static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2730 {
2731 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2732
2733 PetscFunctionBegin;
2734 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
2735 lag->continuous = continuous;
2736 PetscFunctionReturn(0);
2737 }
2738
2739 /*@
2740 PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity
2741
2742 Not Collective
2743
2744 Input Parameter:
2745 . sp - the PetscDualSpace
2746
2747 Output Parameter:
2748 . continuous - flag for element continuity
2749
2750 Level: intermediate
2751
2752 .seealso: PetscDualSpaceLagrangeSetContinuity()
2753 @*/
PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp,PetscBool * continuous)2754 PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2755 {
2756 PetscErrorCode ierr;
2757
2758 PetscFunctionBegin;
2759 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
2760 PetscValidPointer(continuous, 2);
2761 ierr = PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));CHKERRQ(ierr);
2762 PetscFunctionReturn(0);
2763 }
2764
2765 /*@
2766 PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous
2767
2768 Logically Collective on sp
2769
2770 Input Parameters:
2771 + sp - the PetscDualSpace
2772 - continuous - flag for element continuity
2773
2774 Options Database:
2775 . -petscdualspace_lagrange_continuity <bool>
2776
2777 Level: intermediate
2778
2779 .seealso: PetscDualSpaceLagrangeGetContinuity()
2780 @*/
PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp,PetscBool continuous)2781 PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2782 {
2783 PetscErrorCode ierr;
2784
2785 PetscFunctionBegin;
2786 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
2787 PetscValidLogicalCollectiveBool(sp, continuous, 2);
2788 ierr = PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));CHKERRQ(ierr);
2789 PetscFunctionReturn(0);
2790 }
2791
PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp,PetscBool * tensor)2792 static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
2793 {
2794 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2795
2796 PetscFunctionBegin;
2797 *tensor = lag->tensorSpace;
2798 PetscFunctionReturn(0);
2799 }
2800
PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp,PetscBool tensor)2801 static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
2802 {
2803 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2804
2805 PetscFunctionBegin;
2806 lag->tensorSpace = tensor;
2807 PetscFunctionReturn(0);
2808 }
2809
PetscDualSpaceLagrangeGetTrimmed_Lagrange(PetscDualSpace sp,PetscBool * trimmed)2810 static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Lagrange(PetscDualSpace sp, PetscBool *trimmed)
2811 {
2812 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2813
2814 PetscFunctionBegin;
2815 *trimmed = lag->trimmed;
2816 PetscFunctionReturn(0);
2817 }
2818
PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp,PetscBool trimmed)2819 static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp, PetscBool trimmed)
2820 {
2821 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2822
2823 PetscFunctionBegin;
2824 lag->trimmed = trimmed;
2825 PetscFunctionReturn(0);
2826 }
2827
PetscDualSpaceLagrangeGetNodeType_Lagrange(PetscDualSpace sp,PetscDTNodeType * nodeType,PetscBool * boundary,PetscReal * exponent)2828 static PetscErrorCode PetscDualSpaceLagrangeGetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
2829 {
2830 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2831
2832 PetscFunctionBegin;
2833 if (nodeType) *nodeType = lag->nodeType;
2834 if (boundary) *boundary = lag->endNodes;
2835 if (exponent) *exponent = lag->nodeExponent;
2836 PetscFunctionReturn(0);
2837 }
2838
PetscDualSpaceLagrangeSetNodeType_Lagrange(PetscDualSpace sp,PetscDTNodeType nodeType,PetscBool boundary,PetscReal exponent)2839 static PetscErrorCode PetscDualSpaceLagrangeSetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
2840 {
2841 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
2842
2843 PetscFunctionBegin;
2844 if (nodeType == PETSCDTNODES_GAUSSJACOBI && exponent <= -1.) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Exponent must be > -1");
2845 lag->nodeType = nodeType;
2846 lag->endNodes = boundary;
2847 lag->nodeExponent = exponent;
2848 PetscFunctionReturn(0);
2849 }
2850
2851 /*@
2852 PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space
2853
2854 Not collective
2855
2856 Input Parameter:
2857 . sp - The PetscDualSpace
2858
2859 Output Parameter:
2860 . tensor - Whether the dual space has tensor layout (vs. simplicial)
2861
2862 Level: intermediate
2863
2864 .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
2865 @*/
PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp,PetscBool * tensor)2866 PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
2867 {
2868 PetscErrorCode ierr;
2869
2870 PetscFunctionBegin;
2871 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
2872 PetscValidPointer(tensor, 2);
2873 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));CHKERRQ(ierr);
2874 PetscFunctionReturn(0);
2875 }
2876
2877 /*@
2878 PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space
2879
2880 Not collective
2881
2882 Input Parameters:
2883 + sp - The PetscDualSpace
2884 - tensor - Whether the dual space has tensor layout (vs. simplicial)
2885
2886 Level: intermediate
2887
2888 .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
2889 @*/
PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp,PetscBool tensor)2890 PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
2891 {
2892 PetscErrorCode ierr;
2893
2894 PetscFunctionBegin;
2895 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
2896 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));CHKERRQ(ierr);
2897 PetscFunctionReturn(0);
2898 }
2899
2900 /*@
2901 PetscDualSpaceLagrangeGetTrimmed - Get the trimmed nature of the dual space
2902
2903 Not collective
2904
2905 Input Parameter:
2906 . sp - The PetscDualSpace
2907
2908 Output Parameter:
2909 . 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)
2910
2911 Level: intermediate
2912
2913 .seealso: PetscDualSpaceLagrangeSetTrimmed(), PetscDualSpaceCreate()
2914 @*/
PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp,PetscBool * trimmed)2915 PetscErrorCode PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp, PetscBool *trimmed)
2916 {
2917 PetscErrorCode ierr;
2918
2919 PetscFunctionBegin;
2920 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
2921 PetscValidPointer(trimmed, 2);
2922 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTrimmed_C",(PetscDualSpace,PetscBool *),(sp,trimmed));CHKERRQ(ierr);
2923 PetscFunctionReturn(0);
2924 }
2925
2926 /*@
2927 PetscDualSpaceLagrangeSetTrimmed - Set the trimmed nature of the dual space
2928
2929 Not collective
2930
2931 Input Parameters:
2932 + sp - The PetscDualSpace
2933 - 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)
2934
2935 Level: intermediate
2936
2937 .seealso: PetscDualSpaceLagrangeGetTrimmed(), PetscDualSpaceCreate()
2938 @*/
PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp,PetscBool trimmed)2939 PetscErrorCode PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp, PetscBool trimmed)
2940 {
2941 PetscErrorCode ierr;
2942
2943 PetscFunctionBegin;
2944 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
2945 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTrimmed_C",(PetscDualSpace,PetscBool),(sp,trimmed));CHKERRQ(ierr);
2946 PetscFunctionReturn(0);
2947 }
2948
2949 /*@
2950 PetscDualSpaceLagrangeGetNodeType - Get a description of how nodes are laid out for Lagrange polynomials in this
2951 dual space
2952
2953 Not collective
2954
2955 Input Parameter:
2956 . sp - The PetscDualSpace
2957
2958 Output Parameters:
2959 + nodeType - The type of nodes
2960 . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that
2961 include the boundary are Gauss-Lobatto-Jacobi nodes)
2962 - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function
2963 '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type
2964
2965 Level: advanced
2966
2967 .seealso: PetscDTNodeType, PetscDualSpaceLagrangeSetNodeType()
2968 @*/
PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp,PetscDTNodeType * nodeType,PetscBool * boundary,PetscReal * exponent)2969 PetscErrorCode PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent)
2970 {
2971 PetscErrorCode ierr;
2972
2973 PetscFunctionBegin;
2974 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
2975 if (nodeType) PetscValidPointer(nodeType, 2);
2976 if (boundary) PetscValidPointer(boundary, 3);
2977 if (exponent) PetscValidPointer(exponent, 4);
2978 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeGetNodeType_C",(PetscDualSpace,PetscDTNodeType *,PetscBool *,PetscReal *),(sp,nodeType,boundary,exponent));CHKERRQ(ierr);
2979 PetscFunctionReturn(0);
2980 }
2981
2982 /*@
2983 PetscDualSpaceLagrangeSetNodeType - Set a description of how nodes are laid out for Lagrange polynomials in this
2984 dual space
2985
2986 Logically collective
2987
2988 Input Parameters:
2989 + sp - The PetscDualSpace
2990 . nodeType - The type of nodes
2991 . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that
2992 include the boundary are Gauss-Lobatto-Jacobi nodes)
2993 - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function
2994 '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type
2995
2996 Level: advanced
2997
2998 .seealso: PetscDTNodeType, PetscDualSpaceLagrangeGetNodeType()
2999 @*/
PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp,PetscDTNodeType nodeType,PetscBool boundary,PetscReal exponent)3000 PetscErrorCode PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent)
3001 {
3002 PetscErrorCode ierr;
3003
3004 PetscFunctionBegin;
3005 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
3006 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetNodeType_C",(PetscDualSpace,PetscDTNodeType,PetscBool,PetscReal),(sp,nodeType,boundary,exponent));CHKERRQ(ierr);
3007 PetscFunctionReturn(0);
3008 }
3009
3010
PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)3011 static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
3012 {
3013 PetscFunctionBegin;
3014 sp->ops->destroy = PetscDualSpaceDestroy_Lagrange;
3015 sp->ops->view = PetscDualSpaceView_Lagrange;
3016 sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
3017 sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange;
3018 sp->ops->setup = PetscDualSpaceSetUp_Lagrange;
3019 sp->ops->createheightsubspace = NULL;
3020 sp->ops->createpointsubspace = NULL;
3021 sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Lagrange;
3022 sp->ops->apply = PetscDualSpaceApplyDefault;
3023 sp->ops->applyall = PetscDualSpaceApplyAllDefault;
3024 sp->ops->applyint = PetscDualSpaceApplyInteriorDefault;
3025 sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault;
3026 sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault;
3027 PetscFunctionReturn(0);
3028 }
3029
3030 /*MC
3031 PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals
3032
3033 Level: intermediate
3034
3035 .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
3036 M*/
PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)3037 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
3038 {
3039 PetscDualSpace_Lag *lag;
3040 PetscErrorCode ierr;
3041
3042 PetscFunctionBegin;
3043 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1);
3044 ierr = PetscNewLog(sp,&lag);CHKERRQ(ierr);
3045 sp->data = lag;
3046
3047 lag->tensorCell = PETSC_FALSE;
3048 lag->tensorSpace = PETSC_FALSE;
3049 lag->continuous = PETSC_TRUE;
3050 lag->numCopies = PETSC_DEFAULT;
3051 lag->numNodeSkip = PETSC_DEFAULT;
3052 lag->nodeType = PETSCDTNODES_DEFAULT;
3053
3054 ierr = PetscDualSpaceInitialize_Lagrange(sp);CHKERRQ(ierr);
3055 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);CHKERRQ(ierr);
3056 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);CHKERRQ(ierr);
3057 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);CHKERRQ(ierr);
3058 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);CHKERRQ(ierr);
3059 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Lagrange);CHKERRQ(ierr);
3060 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Lagrange);CHKERRQ(ierr);
3061 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Lagrange);CHKERRQ(ierr);
3062 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Lagrange);CHKERRQ(ierr);
3063 PetscFunctionReturn(0);
3064 }
3065