Actual source code: dspacelagrange.c
petsc-3.12.5 2020-03-29
1: #include <petsc/private/petscfeimpl.h>
2: #include <petscdmplex.h>
4: static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
5: {
6: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
7: PetscInt i;
8: PetscErrorCode ierr;
11: if (lag->symmetries) {
12: PetscInt **selfSyms = lag->symmetries[0];
14: if (selfSyms) {
15: PetscInt i, **allocated = &selfSyms[-lag->selfSymOff];
17: for (i = 0; i < lag->numSelfSym; i++) {
18: PetscFree(allocated[i]);
19: }
20: PetscFree(allocated);
21: }
22: PetscFree(lag->symmetries);
23: }
24: for (i = 0; i < lag->height; i++) {
25: PetscDualSpaceDestroy(&lag->subspaces[i]);
26: }
27: PetscFree(lag->subspaces);
28: PetscFree(lag->numDof);
29: PetscFree(lag);
30: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
31: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
32: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);
33: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);
34: return(0);
35: }
37: static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
38: {
39: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
40: PetscErrorCode ierr;
43: PetscViewerASCIIPrintf(viewer, "%s %sLagrange dual space\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "");
44: return(0);
45: }
47: static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
48: {
49: PetscBool iascii;
55: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
56: if (iascii) {PetscDualSpaceLagrangeView_Ascii(sp, viewer);}
57: return(0);
58: }
60: static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
61: {
62: PetscBool continuous, tensor, flg;
66: PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
67: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
68: PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
69: PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
70: if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
71: PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetContinuity", tensor, &tensor, &flg);
72: if (flg) {PetscDualSpaceLagrangeSetTensor(sp, tensor);}
73: PetscOptionsTail();
74: return(0);
75: }
77: static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
78: {
79: PetscInt order, Nc;
80: PetscBool cont, tensor;
81: const char *name;
85: PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
86: PetscObjectGetName((PetscObject) sp, &name);
87: PetscObjectSetName((PetscObject) *spNew, name);
88: PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
89: PetscDualSpaceGetOrder(sp, &order);
90: PetscDualSpaceSetOrder(*spNew, order);
91: PetscDualSpaceGetNumComponents(sp, &Nc);
92: PetscDualSpaceSetNumComponents(*spNew, Nc);
93: PetscDualSpaceLagrangeGetContinuity(sp, &cont);
94: PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
95: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
96: PetscDualSpaceLagrangeSetTensor(*spNew, tensor);
97: return(0);
98: }
100: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
101: {
102: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
103: PetscReal D = 1.0;
104: PetscInt n, d;
105: PetscErrorCode ierr;
108: *dim = -1;
109: DMGetDimension(sp->dm, &n);
110: if (!lag->tensorSpace) {
111: for (d = 1; d <= n; ++d) {
112: D *= ((PetscReal) (order+d))/d;
113: }
114: *dim = (PetscInt) (D + 0.5);
115: } else {
116: *dim = 1;
117: for (d = 0; d < n; ++d) *dim *= (order+1);
118: }
119: *dim *= sp->Nc;
120: return(0);
121: }
123: static PetscErrorCode PetscDualSpaceCreateHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
124: {
125: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
126: PetscBool continuous, tensor;
127: PetscInt order;
128: PetscErrorCode ierr;
131: PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
132: PetscDualSpaceGetOrder(sp,&order);
133: if (height == 0) {
134: PetscObjectReference((PetscObject)sp);
135: *bdsp = sp;
136: } else if (continuous == PETSC_FALSE || !order) {
137: *bdsp = NULL;
138: } else {
139: DM dm, K;
140: PetscInt dim;
142: PetscDualSpaceGetDM(sp,&dm);
143: DMGetDimension(dm,&dim);
144: 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);
145: PetscDualSpaceDuplicate(sp,bdsp);
146: PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, lag->simplexCell, &K);
147: PetscDualSpaceSetDM(*bdsp, K);
148: DMDestroy(&K);
149: PetscDualSpaceLagrangeGetTensor(sp,&tensor);
150: PetscDualSpaceLagrangeSetTensor(*bdsp,tensor);
151: PetscDualSpaceSetUp(*bdsp);
152: }
153: return(0);
154: }
156: static PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
157: {
158: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
159: DM dm = sp->dm;
160: PetscInt order = sp->order;
161: PetscInt Nc = sp->Nc;
162: MPI_Comm comm;
163: PetscBool continuous;
164: PetscSection csection;
165: Vec coordinates;
166: PetscReal *qpoints, *qweights;
167: PetscInt depth, dim, pdimMax, pStart, pEnd, p, *pStratStart, *pStratEnd, coneSize, d, f = 0, c;
168: PetscBool simplex, tensorSpace;
169: PetscErrorCode ierr;
172: PetscObjectGetComm((PetscObject) sp, &comm);
173: if (!order) lag->continuous = PETSC_FALSE;
174: continuous = lag->continuous;
175: DMGetDimension(dm, &dim);
176: DMPlexGetDepth(dm, &depth);
177: DMPlexGetChart(dm, &pStart, &pEnd);
178: PetscCalloc1(dim+1, &lag->numDof);
179: PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);
180: for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);}
181: DMPlexGetConeSize(dm, pStratStart[depth], &coneSize);
182: DMGetCoordinateSection(dm, &csection);
183: DMGetCoordinatesLocal(dm, &coordinates);
184: if (depth == 1) {
185: if (coneSize == dim+1) simplex = PETSC_TRUE;
186: else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
187: else SETERRQ(comm, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
188: } else if (depth == dim) {
189: if (coneSize == dim+1) simplex = PETSC_TRUE;
190: else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
191: else SETERRQ(comm, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
192: } else SETERRQ(comm, PETSC_ERR_SUP, "Only support cell-vertex meshes or fully interpolated meshes");
193: lag->simplexCell = simplex;
194: if (dim > 1 && continuous && lag->simplexCell == lag->tensorSpace) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Mismatching simplex/tensor cells and spaces only allowed for discontinuous elements");
195: tensorSpace = lag->tensorSpace;
196: lag->height = 0;
197: lag->subspaces = NULL;
198: if (continuous && order > 0 && dim > 0) {
199: PetscInt i;
201: lag->height = dim;
202: PetscMalloc1(dim,&lag->subspaces);
203: PetscDualSpaceCreateHeightSubspace_Lagrange(sp,1,&lag->subspaces[0]);
204: PetscDualSpaceSetUp(lag->subspaces[0]);
205: for (i = 1; i < dim; i++) {
206: PetscDualSpaceGetHeightSubspace(lag->subspaces[i-1],1,&lag->subspaces[i]);
207: PetscObjectReference((PetscObject)(lag->subspaces[i]));
208: }
209: }
210: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
211: pdimMax *= (pStratEnd[depth] - pStratStart[depth]);
212: PetscMalloc1(pdimMax, &sp->functional);
213: if (!dim) {
214: for (c = 0; c < Nc; ++c) {
215: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
216: PetscCalloc1(Nc, &qweights);
217: PetscQuadratureSetOrder(sp->functional[f], 0);
218: PetscQuadratureSetData(sp->functional[f], 0, Nc, 1, NULL, qweights);
219: qweights[c] = 1.0;
220: ++f;
221: lag->numDof[0]++;
222: }
223: } else {
224: PetscSection section;
225: PetscReal *v0, *hv0, *J, *invJ, detJ, hdetJ;
226: PetscInt *tup;
228: PetscSectionCreate(PETSC_COMM_SELF,§ion);
229: PetscSectionSetChart(section,pStart,pEnd);
230: PetscCalloc5(dim+1,&tup,dim,&v0,dim,&hv0,dim*dim,&J,dim*dim,&invJ);
231: for (p = pStart; p < pEnd; p++) {
232: PetscInt pointDim, d, nFunc = 0;
233: PetscDualSpace hsp;
235: DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);
236: for (d = 0; d < depth; d++) {if (p >= pStratStart[d] && p < pStratEnd[d]) break;}
237: pointDim = (depth == 1 && d == 1) ? dim : d;
238: hsp = ((pointDim < dim) && lag->subspaces) ? lag->subspaces[dim - pointDim - 1] : NULL;
239: if (hsp) {
240: PetscDualSpace_Lag *hlag = (PetscDualSpace_Lag *) hsp->data;
241: DM hdm;
243: PetscDualSpaceGetDM(hsp,&hdm);
244: DMPlexComputeCellGeometryFEM(hdm, 0, NULL, hv0, NULL, NULL, &hdetJ);
245: nFunc = lag->numDof[pointDim] = hlag->numDof[pointDim];
246: }
247: if (pointDim == dim) {
248: /* Cells, create for self */
249: PetscInt orderEff = continuous ? (!tensorSpace ? order-1-dim : order-2) : order;
250: PetscReal denom = continuous ? order : (!tensorSpace ? order+1+dim : order+2);
251: PetscReal numer = (!simplex || !tensorSpace) ? 2. : (2./dim);
252: PetscReal dx = numer/denom;
253: PetscInt cdim, d, d2;
255: if (orderEff < 0) continue;
256: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, orderEff, &cdim);
257: PetscArrayzero(tup,dim+1);
258: if (!tensorSpace) {
259: while (!tup[dim]) {
260: for (c = 0; c < Nc; ++c) {
261: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
262: PetscMalloc1(dim, &qpoints);
263: PetscCalloc1(Nc, &qweights);
264: PetscQuadratureSetOrder(sp->functional[f], 0);
265: PetscQuadratureSetData(sp->functional[f], dim, Nc, 1, qpoints, qweights);
266: for (d = 0; d < dim; ++d) {
267: qpoints[d] = v0[d];
268: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
269: }
270: qweights[c] = 1.0;
271: ++f;
272: }
273: PetscDualSpaceLatticePointLexicographic_Internal(dim, orderEff, tup);
274: }
275: } else {
276: while (!tup[dim]) {
277: for (c = 0; c < Nc; ++c) {
278: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
279: PetscMalloc1(dim, &qpoints);
280: PetscCalloc1(Nc, &qweights);
281: PetscQuadratureSetOrder(sp->functional[f], 0);
282: PetscQuadratureSetData(sp->functional[f], dim, Nc, 1, qpoints, qweights);
283: for (d = 0; d < dim; ++d) {
284: qpoints[d] = v0[d];
285: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
286: }
287: qweights[c] = 1.0;
288: ++f;
289: }
290: PetscDualSpaceTensorPointLexicographic_Internal(dim, orderEff, tup);
291: }
292: }
293: lag->numDof[dim] = cdim;
294: } else { /* transform functionals from subspaces */
295: PetscInt q;
297: for (q = 0; q < nFunc; q++, f++) {
298: PetscQuadrature fn;
299: PetscInt fdim, Nc, c, nPoints, i;
300: const PetscReal *points;
301: const PetscReal *weights;
302: PetscReal *qpoints;
303: PetscReal *qweights;
305: PetscDualSpaceGetFunctional(hsp, q, &fn);
306: PetscQuadratureGetData(fn,&fdim,&Nc,&nPoints,&points,&weights);
307: if (fdim != pointDim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected height dual space dim %D, got %D",pointDim,fdim);
308: PetscMalloc1(nPoints * dim, &qpoints);
309: PetscCalloc1(nPoints * Nc, &qweights);
310: for (i = 0; i < nPoints; i++) {
311: PetscInt j, k;
312: PetscReal *qp = &qpoints[i * dim];
314: for (c = 0; c < Nc; ++c) qweights[i*Nc+c] = weights[i*Nc+c];
315: for (j = 0; j < dim; ++j) qp[j] = v0[j];
316: for (j = 0; j < dim; ++j) {
317: for (k = 0; k < pointDim; k++) qp[j] += J[dim * j + k] * (points[pointDim * i + k] - hv0[k]);
318: }
319: }
320: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
321: PetscQuadratureSetOrder(sp->functional[f],0);
322: PetscQuadratureSetData(sp->functional[f],dim,Nc,nPoints,qpoints,qweights);
323: }
324: }
325: PetscSectionSetDof(section,p,lag->numDof[pointDim]);
326: }
327: PetscFree5(tup,v0,hv0,J,invJ);
328: PetscSectionSetUp(section);
329: { /* reorder to closure order */
330: PetscInt *key, count;
331: PetscQuadrature *reorder = NULL;
333: PetscCalloc1(f,&key);
334: PetscMalloc1(f*sp->Nc,&reorder);
336: for (p = pStratStart[depth], count = 0; p < pStratEnd[depth]; p++) {
337: PetscInt *closure = NULL, closureSize, c;
339: DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
340: for (c = 0; c < closureSize; c++) {
341: PetscInt point = closure[2 * c], dof, off, i;
343: PetscSectionGetDof(section,point,&dof);
344: PetscSectionGetOffset(section,point,&off);
345: for (i = 0; i < dof; i++) {
346: PetscInt fi = i + off;
347: if (!key[fi]) {
348: key[fi] = 1;
349: reorder[count++] = sp->functional[fi];
350: }
351: }
352: }
353: DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
354: }
355: PetscFree(sp->functional);
356: sp->functional = reorder;
357: PetscFree(key);
358: }
359: PetscSectionDestroy(§ion);
360: }
361: 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);
362: if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %D is greater than max size %D", f, pdimMax);
363: PetscFree2(pStratStart, pStratEnd);
364: return(0);
365: }
367: static PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
368: {
369: DM K;
370: const PetscInt *numDof;
371: PetscInt spatialDim, Nc, size = 0, d;
372: PetscErrorCode ierr;
375: PetscDualSpaceGetDM(sp, &K);
376: PetscDualSpaceGetNumDof(sp, &numDof);
377: DMGetDimension(K, &spatialDim);
378: DMPlexGetHeightStratum(K, 0, NULL, &Nc);
379: if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
380: for (d = 0; d <= spatialDim; ++d) {
381: PetscInt pStart, pEnd;
383: DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
384: size += (pEnd-pStart)*numDof[d];
385: }
386: *dim = size;
387: return(0);
388: }
390: static PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
391: {
392: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
395: *numDof = lag->numDof;
396: return(0);
397: }
399: static PetscErrorCode PetscDualSpaceGetHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
400: {
401: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
402: PetscErrorCode ierr;
405: if (height == 0) {
406: *bdsp = sp;
407: } else {
408: DM dm;
409: PetscInt dim;
411: PetscDualSpaceGetDM(sp,&dm);
412: DMGetDimension(dm,&dim);
413: 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);
414: if (height <= lag->height) {*bdsp = lag->subspaces[height-1];}
415: else {*bdsp = NULL;}
416: }
417: return(0);
418: }
420: #define BaryIndex(perEdge,a,b,c) (((b)*(2*perEdge+1-(b)))/2)+(c)
422: #define CartIndex(perEdge,a,b) (perEdge*(a)+b)
424: static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
425: {
427: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
428: PetscInt dim, order, p, Nc;
429: PetscErrorCode ierr;
432: PetscDualSpaceGetOrder(sp,&order);
433: PetscDualSpaceGetNumComponents(sp,&Nc);
434: DMGetDimension(sp->dm,&dim);
435: if (!dim || !lag->continuous || order < 3) return(0);
436: if (dim > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Lagrange symmetries not implemented for dim = %D > 3",dim);
437: if (!lag->symmetries) { /* store symmetries */
438: PetscDualSpace hsp;
439: DM K;
440: PetscInt numPoints = 1, d;
441: PetscInt numFaces;
442: PetscInt ***symmetries;
443: const PetscInt ***hsymmetries;
445: if (lag->simplexCell) {
446: numFaces = 1 + dim;
447: for (d = 0; d < dim; d++) numPoints = numPoints * 2 + 1;
448: } else {
449: numPoints = PetscPowInt(3,dim);
450: numFaces = 2 * dim;
451: }
452: PetscCalloc1(numPoints,&symmetries);
453: if (0 < dim && dim < 3) { /* compute self symmetries */
454: PetscInt **cellSymmetries;
456: lag->numSelfSym = 2 * numFaces;
457: lag->selfSymOff = numFaces;
458: PetscCalloc1(2*numFaces,&cellSymmetries);
459: /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
460: symmetries[0] = &cellSymmetries[numFaces];
461: if (dim == 1) {
462: PetscInt dofPerEdge = order - 1;
464: if (dofPerEdge > 1) {
465: PetscInt i, j, *reverse;
467: PetscMalloc1(dofPerEdge*Nc,&reverse);
468: for (i = 0; i < dofPerEdge; i++) {
469: for (j = 0; j < Nc; j++) {
470: reverse[i*Nc + j] = Nc * (dofPerEdge - 1 - i) + j;
471: }
472: }
473: symmetries[0][-2] = reverse;
475: /* yes, this is redundant, but it makes it easier to cleanup if I don't have to worry about what not to free */
476: PetscMalloc1(dofPerEdge*Nc,&reverse);
477: for (i = 0; i < dofPerEdge; i++) {
478: for (j = 0; j < Nc; j++) {
479: reverse[i*Nc + j] = Nc * (dofPerEdge - 1 - i) + j;
480: }
481: }
482: symmetries[0][1] = reverse;
483: }
484: } else {
485: PetscInt dofPerEdge = lag->simplexCell ? (order - 2) : (order - 1), s;
486: PetscInt dofPerFace;
488: if (dofPerEdge > 1) {
489: for (s = -numFaces; s < numFaces; s++) {
490: PetscInt *sym, i, j, k, l;
492: if (!s) continue;
493: if (lag->simplexCell) {
494: dofPerFace = (dofPerEdge * (dofPerEdge + 1))/2;
495: PetscMalloc1(Nc*dofPerFace,&sym);
496: for (j = 0, l = 0; j < dofPerEdge; j++) {
497: for (k = 0; k < dofPerEdge - j; k++, l++) {
498: i = dofPerEdge - 1 - j - k;
499: switch (s) {
500: case -3:
501: sym[Nc*l] = BaryIndex(dofPerEdge,i,k,j);
502: break;
503: case -2:
504: sym[Nc*l] = BaryIndex(dofPerEdge,j,i,k);
505: break;
506: case -1:
507: sym[Nc*l] = BaryIndex(dofPerEdge,k,j,i);
508: break;
509: case 1:
510: sym[Nc*l] = BaryIndex(dofPerEdge,k,i,j);
511: break;
512: case 2:
513: sym[Nc*l] = BaryIndex(dofPerEdge,j,k,i);
514: break;
515: }
516: }
517: }
518: } else {
519: dofPerFace = dofPerEdge * dofPerEdge;
520: PetscMalloc1(Nc*dofPerFace,&sym);
521: for (j = 0, l = 0; j < dofPerEdge; j++) {
522: for (k = 0; k < dofPerEdge; k++, l++) {
523: switch (s) {
524: case -4:
525: sym[Nc*l] = CartIndex(dofPerEdge,k,j);
526: break;
527: case -3:
528: sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),k);
529: break;
530: case -2:
531: sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),(dofPerEdge - 1 - j));
532: break;
533: case -1:
534: sym[Nc*l] = CartIndex(dofPerEdge,j,(dofPerEdge - 1 - k));
535: break;
536: case 1:
537: sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),j);
538: break;
539: case 2:
540: sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),(dofPerEdge - 1 - k));
541: break;
542: case 3:
543: sym[Nc*l] = CartIndex(dofPerEdge,k,(dofPerEdge - 1 - j));
544: break;
545: }
546: }
547: }
548: }
549: for (i = 0; i < dofPerFace; i++) {
550: sym[Nc*i] *= Nc;
551: for (j = 1; j < Nc; j++) {
552: sym[Nc*i+j] = sym[Nc*i] + j;
553: }
554: }
555: symmetries[0][s] = sym;
556: }
557: }
558: }
559: }
560: PetscDualSpaceGetHeightSubspace(sp,1,&hsp);
561: PetscDualSpaceGetSymmetries(hsp,&hsymmetries,NULL);
562: if (hsymmetries) {
563: PetscBool *seen;
564: const PetscInt *cone;
565: PetscInt KclosureSize, *Kclosure = NULL;
567: PetscDualSpaceGetDM(sp,&K);
568: PetscCalloc1(numPoints,&seen);
569: DMPlexGetCone(K,0,&cone);
570: DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&KclosureSize,&Kclosure);
571: for (p = 0; p < numFaces; p++) {
572: PetscInt closureSize, *closure = NULL, q;
574: DMPlexGetTransitiveClosure(K,cone[p],PETSC_TRUE,&closureSize,&closure);
575: for (q = 0; q < closureSize; q++) {
576: PetscInt point = closure[2*q], r;
578: if(!seen[point]) {
579: for (r = 0; r < KclosureSize; r++) {
580: if (Kclosure[2 * r] == point) break;
581: }
582: seen[point] = PETSC_TRUE;
583: symmetries[r] = (PetscInt **) hsymmetries[q];
584: }
585: }
586: DMPlexRestoreTransitiveClosure(K,cone[p],PETSC_TRUE,&closureSize,&closure);
587: }
588: DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&KclosureSize,&Kclosure);
589: PetscFree(seen);
590: }
591: lag->symmetries = symmetries;
592: }
593: if (perms) *perms = (const PetscInt ***) lag->symmetries;
594: return(0);
595: }
597: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
598: {
599: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
604: *continuous = lag->continuous;
605: return(0);
606: }
608: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
609: {
610: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
614: lag->continuous = continuous;
615: return(0);
616: }
618: /*@
619: PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity
621: Not Collective
623: Input Parameter:
624: . sp - the PetscDualSpace
626: Output Parameter:
627: . continuous - flag for element continuity
629: Level: intermediate
631: .seealso: PetscDualSpaceLagrangeSetContinuity()
632: @*/
633: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
634: {
640: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
641: return(0);
642: }
644: /*@
645: PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous
647: Logically Collective on sp
649: Input Parameters:
650: + sp - the PetscDualSpace
651: - continuous - flag for element continuity
653: Options Database:
654: . -petscdualspace_lagrange_continuity <bool>
656: Level: intermediate
658: .seealso: PetscDualSpaceLagrangeGetContinuity()
659: @*/
660: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
661: {
667: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
668: return(0);
669: }
671: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
672: {
673: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
676: *tensor = lag->tensorSpace;
677: return(0);
678: }
680: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
681: {
682: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
685: lag->tensorSpace = tensor;
686: return(0);
687: }
689: /*@
690: PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space
692: Not collective
694: Input Parameter:
695: . sp - The PetscDualSpace
697: Output Parameter:
698: . tensor - Whether the dual space has tensor layout (vs. simplicial)
700: Level: intermediate
702: .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate()
703: @*/
704: PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor)
705: {
711: PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));
712: return(0);
713: }
715: /*@
716: PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space
718: Not collective
720: Input Parameters:
721: + sp - The PetscDualSpace
722: - tensor - Whether the dual space has tensor layout (vs. simplicial)
724: Level: intermediate
726: .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate()
727: @*/
728: PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor)
729: {
734: PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));
735: return(0);
736: }
738: static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
739: {
741: sp->ops->destroy = PetscDualSpaceDestroy_Lagrange;
742: sp->ops->view = PetscDualSpaceView_Lagrange;
743: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
744: sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange;
745: sp->ops->setup = PetscDualSpaceSetUp_Lagrange;
746: sp->ops->getdimension = PetscDualSpaceGetDimension_Lagrange;
747: sp->ops->getnumdof = PetscDualSpaceGetNumDof_Lagrange;
748: sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_Lagrange;
749: sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Lagrange;
750: sp->ops->apply = PetscDualSpaceApplyDefault;
751: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
752: sp->ops->createallpoints = PetscDualSpaceCreateAllPointsDefault;
753: return(0);
754: }
756: /*MC
757: PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals
759: Level: intermediate
761: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
762: M*/
764: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
765: {
766: PetscDualSpace_Lag *lag;
767: PetscErrorCode ierr;
771: PetscNewLog(sp,&lag);
772: sp->data = lag;
774: lag->numDof = NULL;
775: lag->simplexCell = PETSC_TRUE;
776: lag->tensorSpace = PETSC_FALSE;
777: lag->continuous = PETSC_TRUE;
779: PetscDualSpaceInitialize_Lagrange(sp);
780: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
781: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
782: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
783: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
784: return(0);
785: }