Actual source code: dspacelagrange.c
petsc-3.10.5 2019-03-28
1: #include <petsc/private/petscfeimpl.h>
2: #include <petscdmplex.h>
4: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
5: {
6: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
9: *tensor = lag->tensorSpace;
10: return(0);
11: }
13: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
14: {
15: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
18: lag->tensorSpace = tensor;
19: return(0);
20: }
22: /*
23: LatticePointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to at most 'max'.
24: Ordering is lexicographic with lowest index as least significant in ordering.
25: e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,0}.
27: Input Parameters:
28: + len - The length of the tuple
29: . max - The maximum sum
30: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition
32: Output Parameter:
33: . tup - A tuple of len integers whos sum is at most 'max'
34: */
35: static PetscErrorCode LatticePointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
36: {
38: while (len--) {
39: max -= tup[len];
40: if (!max) {
41: tup[len] = 0;
42: break;
43: }
44: }
45: tup[++len]++;
46: return(0);
47: }
49: /*
50: TensorPointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that are all less than or equal to 'max'.
51: Ordering is lexicographic with lowest index as least significant in ordering.
52: e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,1}, {0,2}, {1,2}, {2,2}.
54: Input Parameters:
55: + len - The length of the tuple
56: . max - The maximum value
57: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition
59: Output Parameter:
60: . tup - A tuple of len integers whos sum is at most 'max'
61: */
62: static PetscErrorCode TensorPointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
63: {
64: PetscInt i;
67: for (i = 0; i < len; i++) {
68: if (tup[i] < max) {
69: break;
70: } else {
71: tup[i] = 0;
72: }
73: }
74: tup[i]++;
75: return(0);
76: }
79: #define BaryIndex(perEdge,a,b,c) (((b)*(2*perEdge+1-(b)))/2)+(c)
81: #define CartIndex(perEdge,a,b) (perEdge*(a)+b)
83: static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
84: {
86: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
87: PetscInt dim, order, p, Nc;
88: PetscErrorCode ierr;
91: PetscDualSpaceGetOrder(sp,&order);
92: PetscDualSpaceGetNumComponents(sp,&Nc);
93: DMGetDimension(sp->dm,&dim);
94: if (!dim || !lag->continuous || order < 3) return(0);
95: if (dim > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Lagrange symmetries not implemented for dim = %D > 3",dim);
96: if (!lag->symmetries) { /* store symmetries */
97: PetscDualSpace hsp;
98: DM K;
99: PetscInt numPoints = 1, d;
100: PetscInt numFaces;
101: PetscInt ***symmetries;
102: const PetscInt ***hsymmetries;
104: if (lag->simplexCell) {
105: numFaces = 1 + dim;
106: for (d = 0; d < dim; d++) numPoints = numPoints * 2 + 1;
107: }
108: else {
109: numPoints = PetscPowInt(3,dim);
110: numFaces = 2 * dim;
111: }
112: PetscCalloc1(numPoints,&symmetries);
113: if (0 < dim && dim < 3) { /* compute self symmetries */
114: PetscInt **cellSymmetries;
116: lag->numSelfSym = 2 * numFaces;
117: lag->selfSymOff = numFaces;
118: PetscCalloc1(2*numFaces,&cellSymmetries);
119: /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
120: symmetries[0] = &cellSymmetries[numFaces];
121: if (dim == 1) {
122: PetscInt dofPerEdge = order - 1;
124: if (dofPerEdge > 1) {
125: PetscInt i, j, *reverse;
127: PetscMalloc1(dofPerEdge*Nc,&reverse);
128: for (i = 0; i < dofPerEdge; i++) {
129: for (j = 0; j < Nc; j++) {
130: reverse[i*Nc + j] = Nc * (dofPerEdge - 1 - i) + j;
131: }
132: }
133: symmetries[0][-2] = reverse;
135: /* yes, this is redundant, but it makes it easier to cleanup if I don't have to worry about what not to free */
136: PetscMalloc1(dofPerEdge*Nc,&reverse);
137: for (i = 0; i < dofPerEdge; i++) {
138: for (j = 0; j < Nc; j++) {
139: reverse[i*Nc + j] = Nc * (dofPerEdge - 1 - i) + j;
140: }
141: }
142: symmetries[0][1] = reverse;
143: }
144: } else {
145: PetscInt dofPerEdge = lag->simplexCell ? (order - 2) : (order - 1), s;
146: PetscInt dofPerFace;
148: if (dofPerEdge > 1) {
149: for (s = -numFaces; s < numFaces; s++) {
150: PetscInt *sym, i, j, k, l;
152: if (!s) continue;
153: if (lag->simplexCell) {
154: dofPerFace = (dofPerEdge * (dofPerEdge + 1))/2;
155: PetscMalloc1(Nc*dofPerFace,&sym);
156: for (j = 0, l = 0; j < dofPerEdge; j++) {
157: for (k = 0; k < dofPerEdge - j; k++, l++) {
158: i = dofPerEdge - 1 - j - k;
159: switch (s) {
160: case -3:
161: sym[Nc*l] = BaryIndex(dofPerEdge,i,k,j);
162: break;
163: case -2:
164: sym[Nc*l] = BaryIndex(dofPerEdge,j,i,k);
165: break;
166: case -1:
167: sym[Nc*l] = BaryIndex(dofPerEdge,k,j,i);
168: break;
169: case 1:
170: sym[Nc*l] = BaryIndex(dofPerEdge,k,i,j);
171: break;
172: case 2:
173: sym[Nc*l] = BaryIndex(dofPerEdge,j,k,i);
174: break;
175: }
176: }
177: }
178: } else {
179: dofPerFace = dofPerEdge * dofPerEdge;
180: PetscMalloc1(Nc*dofPerFace,&sym);
181: for (j = 0, l = 0; j < dofPerEdge; j++) {
182: for (k = 0; k < dofPerEdge; k++, l++) {
183: switch (s) {
184: case -4:
185: sym[Nc*l] = CartIndex(dofPerEdge,k,j);
186: break;
187: case -3:
188: sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),k);
189: break;
190: case -2:
191: sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),(dofPerEdge - 1 - j));
192: break;
193: case -1:
194: sym[Nc*l] = CartIndex(dofPerEdge,j,(dofPerEdge - 1 - k));
195: break;
196: case 1:
197: sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),j);
198: break;
199: case 2:
200: sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),(dofPerEdge - 1 - k));
201: break;
202: case 3:
203: sym[Nc*l] = CartIndex(dofPerEdge,k,(dofPerEdge - 1 - j));
204: break;
205: }
206: }
207: }
208: }
209: for (i = 0; i < dofPerFace; i++) {
210: sym[Nc*i] *= Nc;
211: for (j = 1; j < Nc; j++) {
212: sym[Nc*i+j] = sym[Nc*i] + j;
213: }
214: }
215: symmetries[0][s] = sym;
216: }
217: }
218: }
219: }
220: PetscDualSpaceGetHeightSubspace(sp,1,&hsp);
221: PetscDualSpaceGetSymmetries(hsp,&hsymmetries,NULL);
222: if (hsymmetries) {
223: PetscBool *seen;
224: const PetscInt *cone;
225: PetscInt KclosureSize, *Kclosure = NULL;
227: PetscDualSpaceGetDM(sp,&K);
228: PetscCalloc1(numPoints,&seen);
229: DMPlexGetCone(K,0,&cone);
230: DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&KclosureSize,&Kclosure);
231: for (p = 0; p < numFaces; p++) {
232: PetscInt closureSize, *closure = NULL, q;
234: DMPlexGetTransitiveClosure(K,cone[p],PETSC_TRUE,&closureSize,&closure);
235: for (q = 0; q < closureSize; q++) {
236: PetscInt point = closure[2*q], r;
238: if(!seen[point]) {
239: for (r = 0; r < KclosureSize; r++) {
240: if (Kclosure[2 * r] == point) break;
241: }
242: seen[point] = PETSC_TRUE;
243: symmetries[r] = (PetscInt **) hsymmetries[q];
244: }
245: }
246: DMPlexRestoreTransitiveClosure(K,cone[p],PETSC_TRUE,&closureSize,&closure);
247: }
248: DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&KclosureSize,&Kclosure);
249: PetscFree(seen);
250: }
251: lag->symmetries = symmetries;
252: }
253: if (perms) *perms = (const PetscInt ***) lag->symmetries;
254: return(0);
255: }
257: /*@C
258: PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis
260: Not collective
262: Input Parameter:
263: . sp - the PetscDualSpace object
265: Output Parameters:
266: + perms - Permutations of the local degrees of freedom, parameterized by the point orientation
267: - flips - Sign reversal of the local degrees of freedom, parameterized by the point orientation
269: Note: The permutation and flip arrays are organized in the following way
270: $ perms[p][ornt][dof # on point] = new local dof #
271: $ flips[p][ornt][dof # on point] = reversal or not
273: Level: developer
275: .seealso: PetscDualSpaceSetSymmetries()
276: @*/
277: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
278: {
283: if (perms) {
285: *perms = NULL;
286: }
287: if (flips) {
289: *flips = NULL;
290: }
291: if (sp->ops->getsymmetries) {
292: (sp->ops->getsymmetries)(sp,perms,flips);
293: }
294: return(0);
295: }
297: static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
298: {
299: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
300: PetscErrorCode ierr;
303: PetscViewerASCIIPrintf(viewer, "%s %sLagrange dual space of order %D", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "Tensor " : "", sp->order, sp->Nc);
304: if (sp->Nc > 1) {PetscViewerASCIIPrintf(viewer, " with %D components", sp->Nc);}
305: PetscViewerASCIIPrintf(viewer, "\n");
306: return(0);
307: }
309: PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
310: {
311: PetscBool iascii;
317: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
318: if (iascii) {PetscDualSpaceLagrangeView_Ascii(sp, viewer);}
319: return(0);
320: }
322: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
323: {
324: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
325: PetscReal D = 1.0;
326: PetscInt n, i;
327: PetscErrorCode ierr;
330: *dim = -1; /* Ensure that the compiler knows *dim is set. */
331: DMGetDimension(sp->dm, &n);
332: if (!lag->tensorSpace) {
333: for (i = 1; i <= n; ++i) {
334: D *= ((PetscReal) (order+i))/i;
335: }
336: *dim = (PetscInt) (D + 0.5);
337: } else {
338: *dim = 1;
339: for (i = 0; i < n; ++i) *dim *= (order+1);
340: }
341: *dim *= sp->Nc;
342: return(0);
343: }
345: static PetscErrorCode PetscDualSpaceCreateHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
346: {
347: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
348: PetscBool continuous, tensor;
349: PetscInt order;
350: PetscErrorCode ierr;
355: PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
356: PetscDualSpaceGetOrder(sp,&order);
357: if (height == 0) {
358: PetscObjectReference((PetscObject)sp);
359: *bdsp = sp;
360: } else if (continuous == PETSC_FALSE || !order) {
361: *bdsp = NULL;
362: } else {
363: DM dm, K;
364: PetscInt dim;
366: PetscDualSpaceGetDM(sp,&dm);
367: DMGetDimension(dm,&dim);
368: 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);}
369: PetscDualSpaceDuplicate(sp,bdsp);
370: PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, lag->simplexCell, &K);
371: PetscDualSpaceSetDM(*bdsp, K);
372: DMDestroy(&K);
373: PetscDualSpaceLagrangeGetTensor(sp,&tensor);
374: PetscDualSpaceLagrangeSetTensor(*bdsp,tensor);
375: PetscDualSpaceSetUp(*bdsp);
376: }
377: return(0);
378: }
380: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
381: {
382: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
383: DM dm = sp->dm;
384: PetscInt order = sp->order;
385: PetscInt Nc = sp->Nc;
386: PetscBool continuous;
387: PetscSection csection;
388: Vec coordinates;
389: PetscReal *qpoints, *qweights;
390: PetscInt depth, dim, pdimMax, pStart, pEnd, p, *pStratStart, *pStratEnd, coneSize, d, f = 0, c;
391: PetscBool simplex, tensorSpace;
392: PetscErrorCode ierr;
395: /* Classify element type */
396: if (!order) lag->continuous = PETSC_FALSE;
397: continuous = lag->continuous;
398: DMGetDimension(dm, &dim);
399: DMPlexGetDepth(dm, &depth);
400: DMPlexGetChart(dm, &pStart, &pEnd);
401: PetscCalloc1(dim+1, &lag->numDof);
402: PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);
403: for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);}
404: DMPlexGetConeSize(dm, pStratStart[depth], &coneSize);
405: DMGetCoordinateSection(dm, &csection);
406: DMGetCoordinatesLocal(dm, &coordinates);
407: if (depth == 1) {
408: if (coneSize == dim+1) simplex = PETSC_TRUE;
409: else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
410: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
411: } else if (depth == dim) {
412: if (coneSize == dim+1) simplex = PETSC_TRUE;
413: else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
414: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
415: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support cell-vertex meshes or interpolated meshes");
416: lag->simplexCell = simplex;
417: 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");
418: tensorSpace = lag->tensorSpace;
419: lag->height = 0;
420: lag->subspaces = NULL;
421: if (continuous && sp->order > 0 && dim > 0) {
422: PetscInt i;
424: lag->height = dim;
425: PetscMalloc1(dim,&lag->subspaces);
426: PetscDualSpaceCreateHeightSubspace_Lagrange(sp,1,&lag->subspaces[0]);
427: PetscDualSpaceSetUp(lag->subspaces[0]);
428: for (i = 1; i < dim; i++) {
429: PetscDualSpaceGetHeightSubspace(lag->subspaces[i-1],1,&lag->subspaces[i]);
430: PetscObjectReference((PetscObject)(lag->subspaces[i]));
431: }
432: }
433: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
434: pdimMax *= (pStratEnd[depth] - pStratStart[depth]);
435: PetscMalloc1(pdimMax, &sp->functional);
436: if (!dim) {
437: for (c = 0; c < Nc; ++c) {
438: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
439: PetscCalloc1(Nc, &qweights);
440: PetscQuadratureSetOrder(sp->functional[f], 0);
441: PetscQuadratureSetData(sp->functional[f], 0, Nc, 1, NULL, qweights);
442: qweights[c] = 1.0;
443: ++f;
444: lag->numDof[0]++;
445: }
446: } else {
447: PetscInt *tup;
448: PetscReal *v0, *hv0, *J, *invJ, detJ, hdetJ;
449: PetscSection section;
451: PetscSectionCreate(PETSC_COMM_SELF,§ion);
452: PetscSectionSetChart(section,pStart,pEnd);
453: PetscCalloc5(dim+1,&tup,dim,&v0,dim,&hv0,dim*dim,&J,dim*dim,&invJ);
454: for (p = pStart; p < pEnd; p++) {
455: PetscInt pointDim, d, nFunc = 0;
456: PetscDualSpace hsp;
458: DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);
459: for (d = 0; d < depth; d++) {if (p >= pStratStart[d] && p < pStratEnd[d]) break;}
460: pointDim = (depth == 1 && d == 1) ? dim : d;
461: hsp = ((pointDim < dim) && lag->subspaces) ? lag->subspaces[dim - pointDim - 1] : NULL;
462: if (hsp) {
463: PetscDualSpace_Lag *hlag = (PetscDualSpace_Lag *) hsp->data;
464: DM hdm;
466: PetscDualSpaceGetDM(hsp,&hdm);
467: DMPlexComputeCellGeometryFEM(hdm, 0, NULL, hv0, NULL, NULL, &hdetJ);
468: nFunc = lag->numDof[pointDim] = hlag->numDof[pointDim];
469: }
470: if (pointDim == dim) {
471: /* Cells, create for self */
472: PetscInt orderEff = continuous ? (!tensorSpace ? order-1-dim : order-2) : order;
473: PetscReal denom = continuous ? order : (!tensorSpace ? order+1+dim : order+2);
474: PetscReal numer = (!simplex || !tensorSpace) ? 2. : (2./dim);
475: PetscReal dx = numer/denom;
476: PetscInt cdim, d, d2;
478: if (orderEff < 0) continue;
479: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, orderEff, &cdim);
480: PetscMemzero(tup,(dim+1)*sizeof(PetscInt));
481: if (!tensorSpace) {
482: while (!tup[dim]) {
483: for (c = 0; c < Nc; ++c) {
484: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
485: PetscMalloc1(dim, &qpoints);
486: PetscCalloc1(Nc, &qweights);
487: PetscQuadratureSetOrder(sp->functional[f], 0);
488: PetscQuadratureSetData(sp->functional[f], dim, Nc, 1, qpoints, qweights);
489: for (d = 0; d < dim; ++d) {
490: qpoints[d] = v0[d];
491: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
492: }
493: qweights[c] = 1.0;
494: ++f;
495: }
496: LatticePointLexicographic_Internal(dim, orderEff, tup);
497: }
498: } else {
499: while (!tup[dim]) {
500: for (c = 0; c < Nc; ++c) {
501: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
502: PetscMalloc1(dim, &qpoints);
503: PetscCalloc1(Nc, &qweights);
504: PetscQuadratureSetOrder(sp->functional[f], 0);
505: PetscQuadratureSetData(sp->functional[f], dim, Nc, 1, qpoints, qweights);
506: for (d = 0; d < dim; ++d) {
507: qpoints[d] = v0[d];
508: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
509: }
510: qweights[c] = 1.0;
511: ++f;
512: }
513: TensorPointLexicographic_Internal(dim, orderEff, tup);
514: }
515: }
516: lag->numDof[dim] = cdim;
517: } else { /* transform functionals from subspaces */
518: PetscInt q;
520: for (q = 0; q < nFunc; q++, f++) {
521: PetscQuadrature fn;
522: PetscInt fdim, Nc, c, nPoints, i;
523: const PetscReal *points;
524: const PetscReal *weights;
525: PetscReal *qpoints;
526: PetscReal *qweights;
528: PetscDualSpaceGetFunctional(hsp, q, &fn);
529: PetscQuadratureGetData(fn,&fdim,&Nc,&nPoints,&points,&weights);
530: if (fdim != pointDim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected height dual space dim %D, got %D",pointDim,fdim);
531: PetscMalloc1(nPoints * dim, &qpoints);
532: PetscCalloc1(nPoints * Nc, &qweights);
533: for (i = 0; i < nPoints; i++) {
534: PetscInt j, k;
535: PetscReal *qp = &qpoints[i * dim];
537: for (c = 0; c < Nc; ++c) qweights[i*Nc+c] = weights[i*Nc+c];
538: for (j = 0; j < dim; ++j) qp[j] = v0[j];
539: for (j = 0; j < dim; ++j) {
540: for (k = 0; k < pointDim; k++) qp[j] += J[dim * j + k] * (points[pointDim * i + k] - hv0[k]);
541: }
542: }
543: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
544: PetscQuadratureSetOrder(sp->functional[f],0);
545: PetscQuadratureSetData(sp->functional[f],dim,Nc,nPoints,qpoints,qweights);
546: }
547: }
548: PetscSectionSetDof(section,p,lag->numDof[pointDim]);
549: }
550: PetscFree5(tup,v0,hv0,J,invJ);
551: PetscSectionSetUp(section);
552: { /* reorder to closure order */
553: PetscInt *key, count;
554: PetscQuadrature *reorder = NULL;
556: PetscCalloc1(f,&key);
557: PetscMalloc1(f*sp->Nc,&reorder);
559: for (p = pStratStart[depth], count = 0; p < pStratEnd[depth]; p++) {
560: PetscInt *closure = NULL, closureSize, c;
562: DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
563: for (c = 0; c < closureSize; c++) {
564: PetscInt point = closure[2 * c], dof, off, i;
566: PetscSectionGetDof(section,point,&dof);
567: PetscSectionGetOffset(section,point,&off);
568: for (i = 0; i < dof; i++) {
569: PetscInt fi = i + off;
570: if (!key[fi]) {
571: key[fi] = 1;
572: reorder[count++] = sp->functional[fi];
573: }
574: }
575: }
576: DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
577: }
578: PetscFree(sp->functional);
579: sp->functional = reorder;
580: PetscFree(key);
581: }
582: PetscSectionDestroy(§ion);
583: }
584: 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);
585: PetscFree2(pStratStart, pStratEnd);
586: if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
587: return(0);
588: }
590: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
591: {
592: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
593: PetscInt i;
594: PetscErrorCode ierr;
597: if (lag->symmetries) {
598: PetscInt **selfSyms = lag->symmetries[0];
600: if (selfSyms) {
601: PetscInt i, **allocated = &selfSyms[-lag->selfSymOff];
603: for (i = 0; i < lag->numSelfSym; i++) {
604: PetscFree(allocated[i]);
605: }
606: PetscFree(allocated);
607: }
608: PetscFree(lag->symmetries);
609: }
610: for (i = 0; i < lag->height; i++) {
611: PetscDualSpaceDestroy(&lag->subspaces[i]);
612: }
613: PetscFree(lag->subspaces);
614: PetscFree(lag->numDof);
615: PetscFree(lag);
616: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
617: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
618: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);
619: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);
620: return(0);
621: }
623: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
624: {
625: PetscInt order, Nc;
626: PetscBool cont, tensor;
630: PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
631: PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
632: PetscDualSpaceGetOrder(sp, &order);
633: PetscDualSpaceSetOrder(*spNew, order);
634: PetscDualSpaceGetNumComponents(sp, &Nc);
635: PetscDualSpaceSetNumComponents(*spNew, Nc);
636: PetscDualSpaceLagrangeGetContinuity(sp, &cont);
637: PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
638: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
639: PetscDualSpaceLagrangeSetTensor(*spNew, tensor);
640: return(0);
641: }
643: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
644: {
645: PetscBool continuous, tensor, flg;
649: PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
650: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
651: PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
652: PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
653: if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
654: PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetContinuity", tensor, &tensor, &flg);
655: if (flg) {PetscDualSpaceLagrangeSetTensor(sp, tensor);}
656: PetscOptionsTail();
657: return(0);
658: }
660: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
661: {
662: DM K;
663: const PetscInt *numDof;
664: PetscInt spatialDim, Nc, size = 0, d;
665: PetscErrorCode ierr;
668: PetscDualSpaceGetDM(sp, &K);
669: PetscDualSpaceGetNumDof(sp, &numDof);
670: DMGetDimension(K, &spatialDim);
671: DMPlexGetHeightStratum(K, 0, NULL, &Nc);
672: if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
673: for (d = 0; d <= spatialDim; ++d) {
674: PetscInt pStart, pEnd;
676: DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
677: size += (pEnd-pStart)*numDof[d];
678: }
679: *dim = size;
680: return(0);
681: }
683: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
684: {
685: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
688: *numDof = lag->numDof;
689: return(0);
690: }
692: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
693: {
694: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
699: *continuous = lag->continuous;
700: return(0);
701: }
703: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
704: {
705: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
709: lag->continuous = continuous;
710: return(0);
711: }
713: /*@
714: PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity
716: Not Collective
718: Input Parameter:
719: . sp - the PetscDualSpace
721: Output Parameter:
722: . continuous - flag for element continuity
724: Level: intermediate
726: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
727: .seealso: PetscDualSpaceLagrangeSetContinuity()
728: @*/
729: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
730: {
736: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
737: return(0);
738: }
740: /*@
741: PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous
743: Logically Collective on PetscDualSpace
745: Input Parameters:
746: + sp - the PetscDualSpace
747: - continuous - flag for element continuity
749: Options Database:
750: . -petscdualspace_lagrange_continuity <bool>
752: Level: intermediate
754: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
755: .seealso: PetscDualSpaceLagrangeGetContinuity()
756: @*/
757: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
758: {
764: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
765: return(0);
766: }
768: PetscErrorCode PetscDualSpaceGetHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
769: {
770: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
771: PetscErrorCode ierr;
776: if (height == 0) {
777: *bdsp = sp;
778: }
779: else {
780: DM dm;
781: PetscInt dim;
783: PetscDualSpaceGetDM(sp,&dm);
784: DMGetDimension(dm,&dim);
785: 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);}
786: if (height <= lag->height) {
787: *bdsp = lag->subspaces[height-1];
788: }
789: else {
790: *bdsp = NULL;
791: }
792: }
793: return(0);
794: }
796: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
797: {
799: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
800: sp->ops->setup = PetscDualSpaceSetUp_Lagrange;
801: sp->ops->view = PetscDualSpaceView_Lagrange;
802: sp->ops->destroy = PetscDualSpaceDestroy_Lagrange;
803: sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange;
804: sp->ops->getdimension = PetscDualSpaceGetDimension_Lagrange;
805: sp->ops->getnumdof = PetscDualSpaceGetNumDof_Lagrange;
806: sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_Lagrange;
807: sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Lagrange;
808: sp->ops->apply = PetscDualSpaceApplyDefault;
809: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
810: sp->ops->createallpoints = PetscDualSpaceCreateAllPointsDefault;
811: return(0);
812: }
814: /*MC
815: PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals
817: Level: intermediate
819: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
820: M*/
822: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
823: {
824: PetscDualSpace_Lag *lag;
825: PetscErrorCode ierr;
829: PetscNewLog(sp,&lag);
830: sp->data = lag;
832: lag->numDof = NULL;
833: lag->simplexCell = PETSC_TRUE;
834: lag->tensorSpace = PETSC_FALSE;
835: lag->continuous = PETSC_TRUE;
837: PetscDualSpaceInitialize_Lagrange(sp);
838: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
839: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
840: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
841: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
842: return(0);
843: }