Actual source code: dspacebdm.c
petsc-3.12.5 2020-03-29
1: #include <petsc/private/petscfeimpl.h>
2: #include <petscdmplex.h>
4: /*
5: Let's work out BDM_1:
7: The model basis is
8: \phi(x, y) = / a + b x + c y \
9: \ d + e x + f y /
10: which is a 6 dimensional space. There are also six dual basis functions,
11: \psi_0(v) = \int^1_{-1} dx v(x, -1) \cdot <0, -1> (1-x)/2
12: \psi_1(v) = \int^1_{-1} dx v(x, -1) \cdot <0, -1> (1+x)/2
13: \psi_2(v) = 1/2 \int^1_{-1} ds v(-s, s) \cdot <1, 1> (1-s)/2 TODO I think the 1/2 is wrong here
14: \psi_3(v) = 1/2 \int^1_{-1} ds v(-s, s) \cdot <1, 1> (1+s)/2
15: \psi_4(v) = -\int^1_{-1} dy v(-1, y) \cdot <-1, 0> (1+y)/2
16: \psi_5(v) = -\int^1_{-1} dy v(-1, y) \cdot <-1, 0> (1-y)/2
17: So we do the integrals
18: \psi_0(\phi) = \int^1_{-1} dx (f - d - e x) (1-x)/2 = (f - d) + e/3
19: \psi_1(\phi) = \int^1_{-1} dx (f - d - e x) (1+x)/2 = (f - d) - e/3
20: \psi_2(\phi) = \int^1_{-1} ds (a - b s + c s + d - e s + f s) (1-s)/2 = (a + d)/2 - (c + f - b - e)/6
21: \psi_3(\phi) = \int^1_{-1} ds (a - b s + c s + d - e s + f s) (1+s)/2 = (a + d)/2 + (c + f - b - e)/6
22: \psi_4(\phi) = \int^1_{-1} dy (b - a - c y) (1+y)/2 = (a - b) + c/3
23: \psi_5(\phi) = \int^1_{-1} dy (b - a - c y) (1-y)/2 = (a - b) - c/3
24: so the nodal basis is
25: \phi_0 = / -(1+x)/2 \
26: \ 1/2 + 3/2 x + y /
27: \phi_1 = / 1+x \
28: \ -1 - 3/2 x - 1/2 y /
29: \phi_2 = / 1+x \
30: \ -(1+y)/2 /
31: \phi_3 = / -(1+x)/2 \
32: \ (1+y) /
33: \phi_4 = / -1 - 1/2 x - 3/2 y \
34: \ (1+y) /
35: \phi_5 = / 1/2 + x + 3/2 y \
36: \ -(1+y)/2 /
37: */
39: static PetscErrorCode PetscDualSpaceDestroy_BDM(PetscDualSpace sp)
40: {
41: PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
42: PetscInt i;
43: PetscErrorCode ierr;
46: if (bdm->symmetries) {
47: PetscInt **selfSyms = bdm->symmetries[0];
48: PetscScalar **selfFlps = bdm->flips[0];
50: if (selfSyms) {
51: PetscInt **fSyms = &selfSyms[-bdm->selfSymOff], i;
52: PetscScalar **fFlps = &selfFlps[-bdm->selfSymOff];
54: for (i = 0; i < bdm->numSelfSym; i++) {
55: PetscFree(fSyms[i]);
56: PetscFree(fFlps[i]);
57: }
58: PetscFree(fSyms);
59: PetscFree(fFlps);
60: }
61: PetscFree(bdm->symmetries);
62: PetscFree(bdm->flips);
63: }
64: for (i = 0; i < bdm->height; i++) {
65: PetscDualSpaceDestroy(&bdm->subspaces[i]);
66: }
67: PetscFree(bdm->subspaces);
68: PetscFree(bdm->numDof);
69: PetscFree(bdm);
70: return(0);
71: }
73: static PetscErrorCode PetscDualSpaceBDMView_Ascii(PetscDualSpace sp, PetscViewer viewer)
74: {
78: PetscViewerASCIIPrintf(viewer, "BDM(%D) dual space\n", sp->order);
79: return(0);
80: }
82: static PetscErrorCode PetscDualSpaceView_BDM(PetscDualSpace sp, PetscViewer viewer)
83: {
84: PetscBool iascii;
90: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
91: if (iascii) {PetscDualSpaceBDMView_Ascii(sp, viewer);}
92: return(0);
93: }
95: static PetscErrorCode PetscDualSpaceSetFromOptions_BDM(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
96: {
100: PetscOptionsHead(PetscOptionsObject,"PetscDualSpace BDM Options");
101: PetscOptionsTail();
102: return(0);
103: }
105: static PetscErrorCode PetscDualSpaceDuplicate_BDM(PetscDualSpace sp, PetscDualSpace *spNew)
106: {
107: PetscInt order, Nc;
108: const char *name;
112: PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
113: PetscObjectGetName((PetscObject) sp, &name);
114: PetscObjectSetName((PetscObject) *spNew, name);
115: PetscDualSpaceSetType(*spNew, PETSCDUALSPACEBDM);
116: PetscDualSpaceGetOrder(sp, &order);
117: PetscDualSpaceSetOrder(*spNew, order);
118: PetscDualSpaceGetNumComponents(sp, &Nc);
119: PetscDualSpaceSetNumComponents(*spNew, Nc);
120: return(0);
121: }
123: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_BDM(PetscDualSpace sp, PetscInt order, PetscInt *dim)
124: {
125: PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
126: PetscReal D = 1.0;
127: PetscInt n, d;
128: PetscErrorCode ierr;
131: *dim = -1;
132: DMGetDimension(sp->dm, &n);
133: if (!n) {*dim = 0; return(0);}
134: if (bdm->simplexCell) {
135: for (d = 1; d <= n; ++d) {
136: D *= ((PetscReal) (order+d))/d;
137: }
138: *dim = (PetscInt) (D + 0.5);
139: } else {
140: *dim = 1;
141: for (d = 0; d < n; ++d) *dim *= (order+1);
142: }
143: if (!bdm->faceSpace) {
144: *dim *= sp->Nc;
145: }
146: return(0);
147: }
149: static PetscErrorCode PetscDualSpaceCreateHeightSubspace_BDM(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
150: {
151: PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
152: PetscInt order;
153: PetscErrorCode ierr;
156: PetscDualSpaceGetOrder(sp, &order);
157: if (!height) {
158: PetscObjectReference((PetscObject) sp);
159: *bdsp = sp;
160: } else if (!order) {
161: *bdsp = NULL;
162: } else if (height == 1) {
163: DM dm, K;
164: PetscInt dim;
166: PetscDualSpaceGetDM(sp, &dm);
167: DMGetDimension(dm, &dim);
168: if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for dual space at height %d for dimension %d reference element\n", height, dim);
169: PetscDualSpaceDuplicate(sp, bdsp);
170: PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, bdm->simplexCell, &K);
171: PetscDualSpaceSetDM(*bdsp, K);
172: DMDestroy(&K);
173: ((PetscDualSpace_BDM *) (*bdsp)->data)->faceSpace = PETSC_TRUE;
174: PetscDualSpaceSetUp(*bdsp);
175: } else {
176: *bdsp = NULL;
177: }
178: return(0);
179: }
181: static PetscErrorCode PetscDualSpaceBDMCreateFaceFE(PetscDualSpace sp, PetscBool tensor, PetscInt faceDim, PetscInt order, PetscFE *fe)
182: {
183: DM K;
184: PetscSpace P;
185: PetscDualSpace Q;
186: PetscQuadrature q;
187: const PetscInt Nc = 1;
188: PetscInt quadPointsPerEdge;
189: PetscErrorCode ierr;
191: /* Create space */
192: PetscSpaceCreate(PETSC_COMM_SELF, &P);
193: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
194: PetscSpacePolynomialSetTensor(P, tensor);
195: PetscSpaceSetNumComponents(P, Nc);
196: PetscSpaceSetNumVariables(P, faceDim);
197: PetscSpaceSetDegree(P, order, order);
198: PetscSpaceSetUp(P);
199: /* Create dual space */
200: PetscDualSpaceCreate(PETSC_COMM_SELF, &Q);
201: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
202: PetscDualSpaceCreateReferenceCell(Q, faceDim, tensor ? PETSC_FALSE : PETSC_TRUE, &K);
203: PetscDualSpaceSetDM(Q, K);
204: DMDestroy(&K);
205: PetscDualSpaceSetNumComponents(Q, Nc);
206: PetscDualSpaceSetOrder(Q, order);
207: PetscDualSpaceLagrangeSetTensor(Q, tensor);
208: PetscDualSpaceSetUp(Q);
209: /* Create element */
210: PetscFECreate(PETSC_COMM_SELF, fe);
211: PetscFESetType(*fe, PETSCFEBASIC);
212: PetscFESetBasisSpace(*fe, P);
213: PetscFESetDualSpace(*fe, Q);
214: PetscFESetNumComponents(*fe, Nc);
215: PetscFESetUp(*fe);
216: PetscSpaceDestroy(&P);
217: PetscDualSpaceDestroy(&Q);
218: /* Create quadrature (with specified order if given) */
219: quadPointsPerEdge = PetscMax(order + 1, 1);
220: if (tensor) {PetscDTGaussTensorQuadrature(faceDim, 1, quadPointsPerEdge, -1.0, 1.0, &q);}
221: else {PetscDTGaussJacobiQuadrature(faceDim, 1, quadPointsPerEdge, -1.0, 1.0, &q);}
222: PetscFESetQuadrature(*fe, q);
223: PetscQuadratureDestroy(&q);
224: return(0);
225: }
227: static PetscErrorCode PetscDualSpaceBDMCreateCellFE(PetscDualSpace sp, PetscBool tensor, PetscInt dim, PetscInt Nc, PetscInt order, PetscFE *fe)
228: {
229: DM K;
230: PetscSpace P;
231: PetscDualSpace Q;
232: PetscQuadrature q;
233: PetscInt quadPointsPerEdge;
234: PetscErrorCode ierr;
236: /* Create space */
237: PetscSpaceCreate(PETSC_COMM_SELF, &P);
238: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
239: PetscSpacePolynomialSetTensor(P, tensor);
240: PetscSpaceSetNumComponents(P, Nc);
241: PetscSpaceSetNumVariables(P, dim);
242: PetscSpaceSetDegree(P, order, order);
243: PetscSpaceSetUp(P);
244: /* Create dual space */
245: /* TODO Needs NED1 dual space */
246: PetscDualSpaceCreate(PETSC_COMM_SELF, &Q);
247: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
248: PetscDualSpaceCreateReferenceCell(Q, dim, tensor ? PETSC_FALSE : PETSC_TRUE, &K);
249: PetscDualSpaceSetDM(Q, K);
250: DMDestroy(&K);
251: PetscDualSpaceSetNumComponents(Q, Nc);
252: PetscDualSpaceSetOrder(Q, order);
253: PetscDualSpaceLagrangeSetTensor(Q, tensor);
254: PetscDualSpaceSetUp(Q);
255: /* Create element */
256: PetscFECreate(PETSC_COMM_SELF, fe);
257: PetscFESetBasisSpace(*fe, P);
258: PetscFESetDualSpace(*fe, Q);
259: PetscFESetNumComponents(*fe, Nc);
260: PetscFESetUp(*fe);
261: PetscSpaceDestroy(&P);
262: PetscDualSpaceDestroy(&Q);
263: /* Create quadrature (with specified order if given) */
264: quadPointsPerEdge = PetscMax(order + 1, 1);
265: if (tensor) {PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);}
266: else {PetscDTGaussJacobiQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);}
267: PetscFESetQuadrature(*fe, q);
268: PetscQuadratureDestroy(&q);
269: return(0);
270: }
272: static PetscErrorCode PetscDualSpaceSetUp_BDM(PetscDualSpace sp)
273: {
274: PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
275: DM dm = sp->dm;
276: PetscInt order = sp->order;
277: PetscInt Nc = sp->Nc;
278: PetscBool faceSp = bdm->faceSpace;
279: MPI_Comm comm;
280: PetscSection csection;
281: Vec coordinates;
282: PetscInt depth, dim, cdim, pdimMax, pStart, pEnd, p, *pStratStart, *pStratEnd, coneSize, d, f = 0, c;
283: PetscBool simplex;
284: PetscErrorCode ierr;
287: PetscObjectGetComm((PetscObject) sp, &comm);
288: DMGetDimension(dm, &dim);
289: DMGetCoordinateDim(dm, &cdim);
290: DMPlexGetDepth(dm, &depth);
291: DMPlexGetChart(dm, &pStart, &pEnd);
292: if (depth != dim) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "BDM element requires interpolated meshes, but depth %D != topological dimension %D", depth, dim);
293: if (!order) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "BDM elements not defined for order 0");
294: if (!faceSp && Nc != cdim) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "BDM element has %D components != %D space dimension", Nc, cdim);
295: PetscCalloc1(dim+1, &bdm->numDof);
296: PetscMalloc2(depth+1, &pStratStart, depth+1, &pStratEnd);
297: for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);}
298: DMPlexGetConeSize(dm, pStratStart[depth], &coneSize);
299: DMGetCoordinateSection(dm, &csection);
300: DMGetCoordinatesLocal(dm, &coordinates);
301: if (coneSize == dim+1) simplex = PETSC_TRUE;
302: else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
303: else SETERRQ(comm, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
304: bdm->simplexCell = simplex;
305: bdm->height = 0;
306: bdm->subspaces = NULL;
307: /* Create the height 1 subspace for every dimension */
308: if (order > 0 && dim > 0) {
309: PetscInt i;
311: bdm->height = dim;
312: PetscMalloc1(dim, &bdm->subspaces);
313: PetscDualSpaceCreateHeightSubspace_BDM(sp, 1, &bdm->subspaces[0]);
314: PetscDualSpaceSetUp(bdm->subspaces[0]);
315: for (i = 1; i < dim; i++) {
316: PetscDualSpaceGetHeightSubspace(bdm->subspaces[i-1], 1, &bdm->subspaces[i]);
317: PetscObjectReference((PetscObject) bdm->subspaces[i]);
318: }
319: }
320: PetscDualSpaceGetDimension_SingleCell_BDM(sp, sp->order, &pdimMax);
321: pdimMax *= (pStratEnd[depth] - pStratStart[depth]);
322: PetscMalloc1(pdimMax, &sp->functional);
323: if (!dim) {
324: bdm->numDof[0] = 0;
325: } else {
326: PetscFE faceFE, cellFE;
327: PetscSection section;
328: CellRefiner cellRefiner;
329: PetscInt faceDim = PetscMax(dim-1, 1), faceNum = 0;
330: PetscReal *v0 = NULL, *J = NULL, *detJ = NULL;
332: PetscSectionCreate(PETSC_COMM_SELF, §ion);
333: PetscSectionSetChart(section, pStart, pEnd);
334: if (!faceSp) {
335: DMPlexGetCellRefiner_Internal(dm, &cellRefiner);
336: CellRefinerGetAffineFaceTransforms_Internal(cellRefiner, NULL, &v0, &J, NULL, &detJ);
337: }
338: /* Create P_q(f) */
339: PetscDualSpaceBDMCreateFaceFE(sp, simplex ? PETSC_FALSE : PETSC_TRUE, faceDim, order, &faceFE);
340: /* Create NED^1_{q-1}(T) = P^d_{q-2} + S_{q-1}(T) */
341: PetscDualSpaceBDMCreateCellFE(sp, simplex ? PETSC_FALSE : PETSC_TRUE, faceDim, Nc, order-1, &cellFE);
342: for (p = pStart; p < pEnd; p++) {
343: PetscBool isFace, isCell;
344: PetscInt d;
346: for (d = 0; d < depth; d++) {if (p >= pStratStart[d] && p < pStratEnd[d]) break;}
347: isFace = ((d == dim) && faceSp) || ((d == dim-1) && !faceSp) ? PETSC_TRUE : PETSC_FALSE;
348: isCell = ((d == dim) && !faceSp) ? PETSC_TRUE : PETSC_FALSE;
349: if (isFace) {
350: PetscQuadrature fq;
351: PetscReal *B, n[3];
352: const PetscReal *fqpoints, *fqweights;
353: PetscInt faceDim = PetscMax(dim-1, 1), Nq, q, fdim, fb;
355: if (cdim == 1) {n[0] = 0.; n[1] = 1.;}
356: else {DMPlexComputeCellGeometryFVM(dm, p, NULL, NULL, n);}
357: PetscFEGetDefaultTabulation(faceFE, &B, NULL, NULL);
358: PetscFEGetQuadrature(faceFE, &fq);
359: PetscQuadratureGetData(fq, NULL, NULL, &Nq, &fqpoints, &fqweights);
360: /* Create a dual basis vector for each basis function */
361: PetscFEGetDimension(faceFE, &fdim);
362: for (fb = 0; fb < fdim; ++fb, ++f) {
363: PetscReal *qpoints, *qweights;
365: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
366: PetscMalloc1(Nq*dim, &qpoints);
367: PetscCalloc1(Nq*Nc, &qweights);
368: PetscQuadratureSetOrder(sp->functional[f], order);
369: PetscQuadratureSetData(sp->functional[f], dim, Nc, Nq, qpoints, qweights);
370: for (q = 0; q < Nq; ++q) {
371: PetscInt g, h;
373: if (faceDim < dim) {
374: /* Transform quadrature points from face coordinates to cell coordinates */
375: for (g = 0; g < dim; ++g) {
376: qpoints[q*dim+g] = v0[faceNum*dim+g];
377: for (h = 0; h < faceDim; ++h) qpoints[q*dim+g] += J[faceNum*dim*faceDim+g*faceDim+h] * fqpoints[q*faceDim+h];
378: }
379: } else {
380: for (g = 0; g < dim; ++g) qpoints[q*dim+g] = fqpoints[q*faceDim+g];
381: }
382: /* Make Radon measure for integral \hat n p ds */
383: for (c = 0; c < Nc; ++c) qweights[q*Nc+c] = B[q*fdim+fb]*n[c]*fqweights[q]*(detJ ? detJ[faceNum] : 1.0);
384: }
385: }
386: bdm->numDof[d] = fdim;
387: PetscSectionSetDof(section, p, bdm->numDof[d]);
388: ++faceNum;
389: }
390: if (order < 2) continue;
391: if (isCell) {
392: PetscSpace csp;
393: PetscQuadrature cq;
394: PetscReal *B;
395: const PetscReal *cqpoints, *cqweights;
396: PetscInt Nq, q, cdim, cb;
398: PetscFEGetBasisSpace(cellFE, &csp);
399: PetscFEGetQuadrature(cellFE, &cq);
400: PetscQuadratureGetData(cq, NULL, NULL, &Nq, &cqpoints, &cqweights);
401: /* Create a dual basis vector for each basis function */
402: PetscSpaceGetDimension(csp, &cdim);
403: PetscMalloc1(Nq*cdim*Nc, &B);
404: PetscSpaceEvaluate(csp, Nq, cqpoints, B, NULL, NULL);
405: for (cb = 0; cb < cdim; ++cb, ++f) {
406: PetscReal *qpoints, *qweights;
408: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
409: PetscMalloc1(Nq*dim, &qpoints);
410: PetscCalloc1(Nq*Nc, &qweights);
411: PetscQuadratureSetOrder(sp->functional[f], order);
412: PetscQuadratureSetData(sp->functional[f], dim, Nc, Nq, qpoints, qweights);
413: PetscArraycpy(qpoints, cqpoints, Nq*dim);
414: for (q = 0; q < Nq; ++q) {
415: /* Make Radon measure for integral p dx */
416: for (c = 0; c < Nc; ++c) qweights[q*Nc+c] = B[(q*cdim+cb)*Nc+c]*cqweights[q*Nc+c];
417: }
418: }
419: PetscFree(B);
420: bdm->numDof[d] = cdim;
421: PetscSectionSetDof(section, p, bdm->numDof[d]);
422: }
423: }
424: PetscFEDestroy(&faceFE);
425: PetscFEDestroy(&cellFE);
426: PetscFree(v0);
427: PetscFree(J);
428: PetscFree(detJ);
429: PetscSectionSetUp(section);
430: { /* reorder to closure order */
431: PetscQuadrature *reorder = NULL;
432: PetscInt *key, count;
434: PetscCalloc1(f, &key);
435: PetscMalloc1(f, &reorder);
436: for (p = pStratStart[depth], count = 0; p < pStratEnd[depth]; p++) {
437: PetscInt *closure = NULL, closureSize, c;
439: DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
440: for (c = 0; c < closureSize*2; c += 2) {
441: PetscInt point = closure[c], dof, off, i;
443: PetscSectionGetDof(section, point, &dof);
444: PetscSectionGetOffset(section, point, &off);
445: for (i = 0; i < dof; ++i) {
446: PetscInt fi = off + i;
448: if (!key[fi]) {
449: key[fi] = 1;
450: reorder[count++] = sp->functional[fi];
451: }
452: }
453: }
454: DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
455: }
456: PetscFree(key);
457: PetscFree(sp->functional);
458: sp->functional = reorder;
459: }
460: PetscSectionDestroy(§ion);
461: }
462: if (pStratEnd[depth] == 1 && f != pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %D not equal to dimension %D", f, pdimMax);
463: if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %D is greater than max size %D", f, pdimMax);
464: PetscFree2(pStratStart, pStratEnd);
465: return(0);
466: }
468: static PetscErrorCode PetscDualSpaceGetDimension_BDM(PetscDualSpace sp, PetscInt *dim)
469: {
470: DM K;
471: const PetscInt *numDof;
472: PetscInt spatialDim, cEnd, size = 0, d;
473: PetscErrorCode ierr;
476: PetscDualSpaceGetDM(sp, &K);
477: PetscDualSpaceGetNumDof(sp, &numDof);
478: DMGetDimension(K, &spatialDim);
479: DMPlexGetHeightStratum(K, 0, NULL, &cEnd);
480: if (cEnd == 1) {PetscDualSpaceGetDimension_SingleCell_BDM(sp, sp->order, dim); return(0);}
481: for (d = 0; d <= spatialDim; ++d) {
482: PetscInt pStart, pEnd;
484: DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
485: size += (pEnd-pStart)*numDof[d];
486: }
487: *dim = size;
488: return(0);
489: }
491: static PetscErrorCode PetscDualSpaceGetNumDof_BDM(PetscDualSpace sp, const PetscInt **numDof)
492: {
493: PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
496: *numDof = bdm->numDof;
497: return(0);
498: }
500: static PetscErrorCode PetscDualSpaceGetHeightSubspace_BDM(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
501: {
502: PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
503: PetscErrorCode ierr;
508: if (height == 0) {
509: *bdsp = sp;
510: } else {
511: DM dm;
512: PetscInt dim;
514: PetscDualSpaceGetDM(sp, &dm);
515: DMGetDimension(dm, &dim);
516: if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for dual space at height %D for dimension %D reference element\n", height, dim);
517: if (height <= bdm->height) {*bdsp = bdm->subspaces[height-1];}
518: else {*bdsp = NULL;}
519: }
520: return(0);
521: }
523: #define BaryIndex(perEdge,a,b,c) (((b)*(2*perEdge+1-(b)))/2)+(c)
525: #define CartIndex(perEdge,a,b) (perEdge*(a)+b)
527: static PetscErrorCode PetscDualSpaceGetSymmetries_BDM(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****rots)
528: {
530: PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
531: PetscInt dim, order, p, Nc;
532: PetscErrorCode ierr;
535: PetscDualSpaceGetOrder(sp, &order);
536: PetscDualSpaceGetNumComponents(sp, &Nc);
537: DMGetDimension(sp->dm, &dim);
538: if (dim < 1) return(0);
539: if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "BDM symmetries not implemented for dim = %D > 3", dim);
540: if (!bdm->symmetries) { /* store symmetries */
541: PetscInt ***symmetries, **cellSymmetries;
542: PetscScalar ***flips, **cellFlips;
543: PetscInt numPoints, numFaces, d;
545: if (bdm->simplexCell) {
546: numPoints = 1;
547: for (d = 0; d < dim; d++) numPoints = numPoints * 2 + 1;
548: numFaces = 1 + dim;
549: } else {
550: numPoints = PetscPowInt(3,dim);
551: numFaces = 2 * dim;
552: }
553: PetscCalloc1(numPoints, &symmetries);
554: PetscCalloc1(numPoints, &flips);
555: /* compute self symmetries */
556: bdm->numSelfSym = 2 * numFaces;
557: bdm->selfSymOff = numFaces;
558: PetscCalloc1(2*numFaces, &cellSymmetries);
559: PetscCalloc1(2*numFaces, &cellFlips);
560: /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
561: symmetries[0] = &cellSymmetries[bdm->selfSymOff];
562: flips[0] = &cellFlips[bdm->selfSymOff];
563: switch (dim) {
564: case 1: /* Edge symmetries */
565: {
566: PetscScalar *invert;
567: PetscInt *reverse;
568: PetscInt dofPerEdge = order+1, eNc = 1 /* ??? */, i, j;
570: PetscMalloc1(dofPerEdge*eNc, &reverse);
571: PetscMalloc1(dofPerEdge*eNc, &invert);
572: for (i = 0; i < dofPerEdge; ++i) {
573: for (j = 0; j < eNc; ++j) {
574: reverse[i*eNc + j] = eNc * (dofPerEdge - 1 - i) + j;
575: invert[i*eNc + j] = -1.0;
576: }
577: }
578: symmetries[0][-2] = reverse;
579: flips[0][-2] = invert;
581: /* yes, this is redundant, but it makes it easier to cleanup if I don't have to worry about what not to free */
582: PetscMalloc1(dofPerEdge*eNc, &reverse);
583: PetscMalloc1(dofPerEdge*eNc, &invert);
584: for (i = 0; i < dofPerEdge; i++) {
585: for (j = 0; j < eNc; j++) {
586: reverse[i*eNc + j] = eNc * (dofPerEdge - 1 - i) + j;
587: invert[i*eNc + j] = -1.0;
588: }
589: }
590: symmetries[0][1] = reverse;
591: flips[0][1] = invert;
592: break;
593: }
594: case 2: /* Face symmetries */
595: {
596: PetscInt dofPerEdge = bdm->simplexCell ? (order - 2) : (order - 1), s;
597: PetscInt dofPerFace;
599: for (s = -numFaces; s < numFaces; s++) {
600: PetscScalar *flp;
601: PetscInt *sym, fNc = 1, i, j, k, l;
603: if (!s) continue;
604: if (bdm->simplexCell) {
605: dofPerFace = (dofPerEdge * (dofPerEdge + 1))/2;
606: PetscMalloc1(fNc*dofPerFace,&sym);
607: PetscMalloc1(fNc*dofPerFace,&flp);
608: for (j = 0, l = 0; j < dofPerEdge; j++) {
609: for (k = 0; k < dofPerEdge - j; k++, l++) {
610: i = dofPerEdge - 1 - j - k;
611: switch (s) {
612: case -3:
613: sym[fNc*l] = BaryIndex(dofPerEdge,i,k,j);
614: break;
615: case -2:
616: sym[fNc*l] = BaryIndex(dofPerEdge,j,i,k);
617: break;
618: case -1:
619: sym[fNc*l] = BaryIndex(dofPerEdge,k,j,i);
620: break;
621: case 1:
622: sym[fNc*l] = BaryIndex(dofPerEdge,k,i,j);
623: break;
624: case 2:
625: sym[fNc*l] = BaryIndex(dofPerEdge,j,k,i);
626: break;
627: }
628: flp[fNc*l] = s < 0 ? -1.0 : 1.0;
629: }
630: }
631: } else {
632: dofPerFace = dofPerEdge * dofPerEdge;
633: PetscMalloc1(fNc*dofPerFace,&sym);
634: PetscMalloc1(fNc*dofPerFace,&flp);
635: for (j = 0, l = 0; j < dofPerEdge; j++) {
636: for (k = 0; k < dofPerEdge; k++, l++) {
637: switch (s) {
638: case -4:
639: sym[fNc*l] = CartIndex(dofPerEdge,k,j);
640: break;
641: case -3:
642: sym[fNc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),k);
643: break;
644: case -2:
645: sym[fNc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),(dofPerEdge - 1 - j));
646: break;
647: case -1:
648: sym[fNc*l] = CartIndex(dofPerEdge,j,(dofPerEdge - 1 - k));
649: break;
650: case 1:
651: sym[fNc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),j);
652: break;
653: case 2:
654: sym[fNc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),(dofPerEdge - 1 - k));
655: break;
656: case 3:
657: sym[fNc*l] = CartIndex(dofPerEdge,k,(dofPerEdge - 1 - j));
658: break;
659: }
660: flp[fNc*l] = s < 0 ? -1.0 : 1.0;
661: }
662: }
663: }
664: for (i = 0; i < dofPerFace; i++) {
665: sym[fNc*i] *= fNc;
666: for (j = 1; j < fNc; j++) {
667: sym[fNc*i+j] = sym[fNc*i] + j;
668: flp[fNc*i+j] = flp[fNc*i];
669: }
670: }
671: symmetries[0][s] = sym;
672: flips[0][s] = flp;
673: }
674: break;
675: }
676: default: SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "No symmetries for point of dimension %D", dim);
677: }
678: /* Copy subspace symmetries */
679: {
680: PetscDualSpace hsp;
681: DM K;
682: const PetscInt ***hsymmetries;
683: const PetscScalar ***hflips;
685: PetscDualSpaceGetHeightSubspace(sp, 1, &hsp);
686: PetscDualSpaceGetSymmetries(hsp, &hsymmetries, &hflips);
687: if (hsymmetries || hflips) {
688: PetscBool *seen;
689: const PetscInt *cone;
690: PetscInt KclosureSize, *Kclosure = NULL;
692: PetscDualSpaceGetDM(sp, &K);
693: PetscCalloc1(numPoints, &seen);
694: DMPlexGetCone(K, 0, &cone);
695: DMPlexGetTransitiveClosure(K, 0, PETSC_TRUE, &KclosureSize, &Kclosure);
696: for (p = 0; p < numFaces; ++p) {
697: PetscInt closureSize, *closure = NULL, q;
699: DMPlexGetTransitiveClosure(K, cone[p], PETSC_TRUE, &closureSize, &closure);
700: for (q = 0; q < closureSize; ++q) {
701: PetscInt point = closure[q*2], r;
703: if (!seen[point]) {
704: for (r = 0; r < KclosureSize; ++r) {
705: if (Kclosure[r*2] == point) break;
706: }
707: seen[point] = PETSC_TRUE;
708: symmetries[r] = (PetscInt **) hsymmetries[q];
709: flips[r] = (PetscScalar **) hflips[q];
710: }
711: }
712: DMPlexRestoreTransitiveClosure(K, cone[p], PETSC_TRUE, &closureSize, &closure);
713: }
714: DMPlexRestoreTransitiveClosure(K, 0, PETSC_TRUE, &KclosureSize, &Kclosure);
715: PetscFree(seen);
716: }
717: }
718: bdm->symmetries = symmetries;
719: bdm->flips = flips;
720: }
721: if (perms) *perms = (const PetscInt ***) bdm->symmetries;
722: if (rots) *rots = (const PetscScalar ***) bdm->flips;
723: return(0);
724: }
726: static PetscErrorCode PetscDualSpaceInitialize_BDM(PetscDualSpace sp)
727: {
729: sp->ops->destroy = PetscDualSpaceDestroy_BDM;
730: sp->ops->view = PetscDualSpaceView_BDM;
731: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_BDM;
732: sp->ops->duplicate = PetscDualSpaceDuplicate_BDM;
733: sp->ops->setup = PetscDualSpaceSetUp_BDM;
734: sp->ops->getdimension = PetscDualSpaceGetDimension_BDM;
735: sp->ops->getnumdof = PetscDualSpaceGetNumDof_BDM;
736: sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_BDM;
737: sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_BDM;
738: sp->ops->apply = PetscDualSpaceApplyDefault;
739: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
740: sp->ops->createallpoints = PetscDualSpaceCreateAllPointsDefault;
741: return(0);
742: }
743: /*MC
744: PETSCDUALSPACEBDM = "bdm" - A PetscDualSpace object that encapsulates a dual space for Brezzi-Douglas-Marini elements
746: Level: intermediate
748: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
749: M*/
751: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_BDM(PetscDualSpace sp)
752: {
753: PetscDualSpace_BDM *bdm;
754: PetscErrorCode ierr;
758: PetscNewLog(sp, &bdm);
759: sp->data = bdm;
760: sp->k = 3;
762: bdm->numDof = NULL;
763: bdm->simplexCell = PETSC_TRUE;
765: PetscDualSpaceInitialize_BDM(sp);
766: return(0);
767: }