Actual source code: dspacelagrange.c
petsc-3.11.4 2019-09-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\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "");
304: return(0);
305: }
307: PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
308: {
309: PetscBool iascii;
315: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
316: if (iascii) {PetscDualSpaceLagrangeView_Ascii(sp, viewer);}
317: return(0);
318: }
320: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
321: {
322: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
323: PetscReal D = 1.0;
324: PetscInt n, i;
325: PetscErrorCode ierr;
328: *dim = -1; /* Ensure that the compiler knows *dim is set. */
329: DMGetDimension(sp->dm, &n);
330: if (!lag->tensorSpace) {
331: for (i = 1; i <= n; ++i) {
332: D *= ((PetscReal) (order+i))/i;
333: }
334: *dim = (PetscInt) (D + 0.5);
335: } else {
336: *dim = 1;
337: for (i = 0; i < n; ++i) *dim *= (order+1);
338: }
339: *dim *= sp->Nc;
340: return(0);
341: }
343: static PetscErrorCode PetscDualSpaceCreateHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
344: {
345: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
346: PetscBool continuous, tensor;
347: PetscInt order;
348: PetscErrorCode ierr;
353: PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
354: PetscDualSpaceGetOrder(sp,&order);
355: if (height == 0) {
356: PetscObjectReference((PetscObject)sp);
357: *bdsp = sp;
358: } else if (continuous == PETSC_FALSE || !order) {
359: *bdsp = NULL;
360: } else {
361: DM dm, K;
362: PetscInt dim;
364: PetscDualSpaceGetDM(sp,&dm);
365: DMGetDimension(dm,&dim);
366: 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);}
367: PetscDualSpaceDuplicate(sp,bdsp);
368: PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, lag->simplexCell, &K);
369: PetscDualSpaceSetDM(*bdsp, K);
370: DMDestroy(&K);
371: PetscDualSpaceLagrangeGetTensor(sp,&tensor);
372: PetscDualSpaceLagrangeSetTensor(*bdsp,tensor);
373: PetscDualSpaceSetUp(*bdsp);
374: }
375: return(0);
376: }
378: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
379: {
380: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
381: DM dm = sp->dm;
382: PetscInt order = sp->order;
383: PetscInt Nc = sp->Nc;
384: PetscBool continuous;
385: PetscSection csection;
386: Vec coordinates;
387: PetscReal *qpoints, *qweights;
388: PetscInt depth, dim, pdimMax, pStart, pEnd, p, *pStratStart, *pStratEnd, coneSize, d, f = 0, c;
389: PetscBool simplex, tensorSpace;
390: PetscErrorCode ierr;
393: /* Classify element type */
394: if (!order) lag->continuous = PETSC_FALSE;
395: continuous = lag->continuous;
396: DMGetDimension(dm, &dim);
397: DMPlexGetDepth(dm, &depth);
398: DMPlexGetChart(dm, &pStart, &pEnd);
399: PetscCalloc1(dim+1, &lag->numDof);
400: PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);
401: for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);}
402: DMPlexGetConeSize(dm, pStratStart[depth], &coneSize);
403: DMGetCoordinateSection(dm, &csection);
404: DMGetCoordinatesLocal(dm, &coordinates);
405: if (depth == 1) {
406: if (coneSize == dim+1) simplex = PETSC_TRUE;
407: else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
408: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
409: } else if (depth == dim) {
410: if (coneSize == dim+1) simplex = PETSC_TRUE;
411: else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
412: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
413: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support cell-vertex meshes or interpolated meshes");
414: lag->simplexCell = simplex;
415: 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");
416: tensorSpace = lag->tensorSpace;
417: lag->height = 0;
418: lag->subspaces = NULL;
419: if (continuous && sp->order > 0 && dim > 0) {
420: PetscInt i;
422: lag->height = dim;
423: PetscMalloc1(dim,&lag->subspaces);
424: PetscDualSpaceCreateHeightSubspace_Lagrange(sp,1,&lag->subspaces[0]);
425: PetscDualSpaceSetUp(lag->subspaces[0]);
426: for (i = 1; i < dim; i++) {
427: PetscDualSpaceGetHeightSubspace(lag->subspaces[i-1],1,&lag->subspaces[i]);
428: PetscObjectReference((PetscObject)(lag->subspaces[i]));
429: }
430: }
431: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
432: pdimMax *= (pStratEnd[depth] - pStratStart[depth]);
433: PetscMalloc1(pdimMax, &sp->functional);
434: if (!dim) {
435: for (c = 0; c < Nc; ++c) {
436: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
437: PetscCalloc1(Nc, &qweights);
438: PetscQuadratureSetOrder(sp->functional[f], 0);
439: PetscQuadratureSetData(sp->functional[f], 0, Nc, 1, NULL, qweights);
440: qweights[c] = 1.0;
441: ++f;
442: lag->numDof[0]++;
443: }
444: } else {
445: PetscInt *tup;
446: PetscReal *v0, *hv0, *J, *invJ, detJ, hdetJ;
447: PetscSection section;
449: PetscSectionCreate(PETSC_COMM_SELF,§ion);
450: PetscSectionSetChart(section,pStart,pEnd);
451: PetscCalloc5(dim+1,&tup,dim,&v0,dim,&hv0,dim*dim,&J,dim*dim,&invJ);
452: for (p = pStart; p < pEnd; p++) {
453: PetscInt pointDim, d, nFunc = 0;
454: PetscDualSpace hsp;
456: DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);
457: for (d = 0; d < depth; d++) {if (p >= pStratStart[d] && p < pStratEnd[d]) break;}
458: pointDim = (depth == 1 && d == 1) ? dim : d;
459: hsp = ((pointDim < dim) && lag->subspaces) ? lag->subspaces[dim - pointDim - 1] : NULL;
460: if (hsp) {
461: PetscDualSpace_Lag *hlag = (PetscDualSpace_Lag *) hsp->data;
462: DM hdm;
464: PetscDualSpaceGetDM(hsp,&hdm);
465: DMPlexComputeCellGeometryFEM(hdm, 0, NULL, hv0, NULL, NULL, &hdetJ);
466: nFunc = lag->numDof[pointDim] = hlag->numDof[pointDim];
467: }
468: if (pointDim == dim) {
469: /* Cells, create for self */
470: PetscInt orderEff = continuous ? (!tensorSpace ? order-1-dim : order-2) : order;
471: PetscReal denom = continuous ? order : (!tensorSpace ? order+1+dim : order+2);
472: PetscReal numer = (!simplex || !tensorSpace) ? 2. : (2./dim);
473: PetscReal dx = numer/denom;
474: PetscInt cdim, d, d2;
476: if (orderEff < 0) continue;
477: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, orderEff, &cdim);
478: PetscMemzero(tup,(dim+1)*sizeof(PetscInt));
479: if (!tensorSpace) {
480: while (!tup[dim]) {
481: for (c = 0; c < Nc; ++c) {
482: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
483: PetscMalloc1(dim, &qpoints);
484: PetscCalloc1(Nc, &qweights);
485: PetscQuadratureSetOrder(sp->functional[f], 0);
486: PetscQuadratureSetData(sp->functional[f], dim, Nc, 1, qpoints, qweights);
487: for (d = 0; d < dim; ++d) {
488: qpoints[d] = v0[d];
489: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
490: }
491: qweights[c] = 1.0;
492: ++f;
493: }
494: LatticePointLexicographic_Internal(dim, orderEff, tup);
495: }
496: } else {
497: while (!tup[dim]) {
498: for (c = 0; c < Nc; ++c) {
499: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
500: PetscMalloc1(dim, &qpoints);
501: PetscCalloc1(Nc, &qweights);
502: PetscQuadratureSetOrder(sp->functional[f], 0);
503: PetscQuadratureSetData(sp->functional[f], dim, Nc, 1, qpoints, qweights);
504: for (d = 0; d < dim; ++d) {
505: qpoints[d] = v0[d];
506: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
507: }
508: qweights[c] = 1.0;
509: ++f;
510: }
511: TensorPointLexicographic_Internal(dim, orderEff, tup);
512: }
513: }
514: lag->numDof[dim] = cdim;
515: } else { /* transform functionals from subspaces */
516: PetscInt q;
518: for (q = 0; q < nFunc; q++, f++) {
519: PetscQuadrature fn;
520: PetscInt fdim, Nc, c, nPoints, i;
521: const PetscReal *points;
522: const PetscReal *weights;
523: PetscReal *qpoints;
524: PetscReal *qweights;
526: PetscDualSpaceGetFunctional(hsp, q, &fn);
527: PetscQuadratureGetData(fn,&fdim,&Nc,&nPoints,&points,&weights);
528: if (fdim != pointDim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected height dual space dim %D, got %D",pointDim,fdim);
529: PetscMalloc1(nPoints * dim, &qpoints);
530: PetscCalloc1(nPoints * Nc, &qweights);
531: for (i = 0; i < nPoints; i++) {
532: PetscInt j, k;
533: PetscReal *qp = &qpoints[i * dim];
535: for (c = 0; c < Nc; ++c) qweights[i*Nc+c] = weights[i*Nc+c];
536: for (j = 0; j < dim; ++j) qp[j] = v0[j];
537: for (j = 0; j < dim; ++j) {
538: for (k = 0; k < pointDim; k++) qp[j] += J[dim * j + k] * (points[pointDim * i + k] - hv0[k]);
539: }
540: }
541: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
542: PetscQuadratureSetOrder(sp->functional[f],0);
543: PetscQuadratureSetData(sp->functional[f],dim,Nc,nPoints,qpoints,qweights);
544: }
545: }
546: PetscSectionSetDof(section,p,lag->numDof[pointDim]);
547: }
548: PetscFree5(tup,v0,hv0,J,invJ);
549: PetscSectionSetUp(section);
550: { /* reorder to closure order */
551: PetscInt *key, count;
552: PetscQuadrature *reorder = NULL;
554: PetscCalloc1(f,&key);
555: PetscMalloc1(f*sp->Nc,&reorder);
557: for (p = pStratStart[depth], count = 0; p < pStratEnd[depth]; p++) {
558: PetscInt *closure = NULL, closureSize, c;
560: DMPlexGetTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
561: for (c = 0; c < closureSize; c++) {
562: PetscInt point = closure[2 * c], dof, off, i;
564: PetscSectionGetDof(section,point,&dof);
565: PetscSectionGetOffset(section,point,&off);
566: for (i = 0; i < dof; i++) {
567: PetscInt fi = i + off;
568: if (!key[fi]) {
569: key[fi] = 1;
570: reorder[count++] = sp->functional[fi];
571: }
572: }
573: }
574: DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
575: }
576: PetscFree(sp->functional);
577: sp->functional = reorder;
578: PetscFree(key);
579: }
580: PetscSectionDestroy(§ion);
581: }
582: 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);
583: PetscFree2(pStratStart, pStratEnd);
584: if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
585: return(0);
586: }
588: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
589: {
590: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
591: PetscInt i;
592: PetscErrorCode ierr;
595: if (lag->symmetries) {
596: PetscInt **selfSyms = lag->symmetries[0];
598: if (selfSyms) {
599: PetscInt i, **allocated = &selfSyms[-lag->selfSymOff];
601: for (i = 0; i < lag->numSelfSym; i++) {
602: PetscFree(allocated[i]);
603: }
604: PetscFree(allocated);
605: }
606: PetscFree(lag->symmetries);
607: }
608: for (i = 0; i < lag->height; i++) {
609: PetscDualSpaceDestroy(&lag->subspaces[i]);
610: }
611: PetscFree(lag->subspaces);
612: PetscFree(lag->numDof);
613: PetscFree(lag);
614: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
615: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
616: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);
617: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);
618: return(0);
619: }
621: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
622: {
623: PetscInt order, Nc;
624: PetscBool cont, tensor;
625: const char *name;
629: PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
630: PetscObjectGetName((PetscObject) sp, &name);
631: PetscObjectSetName((PetscObject) *spNew, name);
632: PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
633: PetscDualSpaceGetOrder(sp, &order);
634: PetscDualSpaceSetOrder(*spNew, order);
635: PetscDualSpaceGetNumComponents(sp, &Nc);
636: PetscDualSpaceSetNumComponents(*spNew, Nc);
637: PetscDualSpaceLagrangeGetContinuity(sp, &cont);
638: PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
639: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
640: PetscDualSpaceLagrangeSetTensor(*spNew, tensor);
641: return(0);
642: }
644: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
645: {
646: PetscBool continuous, tensor, flg;
650: PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
651: PetscDualSpaceLagrangeGetTensor(sp, &tensor);
652: PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
653: PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
654: if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
655: PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetContinuity", tensor, &tensor, &flg);
656: if (flg) {PetscDualSpaceLagrangeSetTensor(sp, tensor);}
657: PetscOptionsTail();
658: return(0);
659: }
661: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
662: {
663: DM K;
664: const PetscInt *numDof;
665: PetscInt spatialDim, Nc, size = 0, d;
666: PetscErrorCode ierr;
669: PetscDualSpaceGetDM(sp, &K);
670: PetscDualSpaceGetNumDof(sp, &numDof);
671: DMGetDimension(K, &spatialDim);
672: DMPlexGetHeightStratum(K, 0, NULL, &Nc);
673: if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
674: for (d = 0; d <= spatialDim; ++d) {
675: PetscInt pStart, pEnd;
677: DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
678: size += (pEnd-pStart)*numDof[d];
679: }
680: *dim = size;
681: return(0);
682: }
684: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
685: {
686: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
689: *numDof = lag->numDof;
690: return(0);
691: }
693: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
694: {
695: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
700: *continuous = lag->continuous;
701: return(0);
702: }
704: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
705: {
706: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
710: lag->continuous = continuous;
711: return(0);
712: }
714: /*@
715: PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity
717: Not Collective
719: Input Parameter:
720: . sp - the PetscDualSpace
722: Output Parameter:
723: . continuous - flag for element continuity
725: Level: intermediate
727: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
728: .seealso: PetscDualSpaceLagrangeSetContinuity()
729: @*/
730: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
731: {
737: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
738: return(0);
739: }
741: /*@
742: PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous
744: Logically Collective on PetscDualSpace
746: Input Parameters:
747: + sp - the PetscDualSpace
748: - continuous - flag for element continuity
750: Options Database:
751: . -petscdualspace_lagrange_continuity <bool>
753: Level: intermediate
755: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
756: .seealso: PetscDualSpaceLagrangeGetContinuity()
757: @*/
758: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
759: {
765: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
766: return(0);
767: }
769: PetscErrorCode PetscDualSpaceGetHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
770: {
771: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
772: PetscErrorCode ierr;
777: if (height == 0) {
778: *bdsp = sp;
779: }
780: else {
781: DM dm;
782: PetscInt dim;
784: PetscDualSpaceGetDM(sp,&dm);
785: DMGetDimension(dm,&dim);
786: 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);}
787: if (height <= lag->height) {
788: *bdsp = lag->subspaces[height-1];
789: }
790: else {
791: *bdsp = NULL;
792: }
793: }
794: return(0);
795: }
797: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
798: {
800: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
801: sp->ops->setup = PetscDualSpaceSetUp_Lagrange;
802: sp->ops->view = PetscDualSpaceView_Lagrange;
803: sp->ops->destroy = PetscDualSpaceDestroy_Lagrange;
804: sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange;
805: sp->ops->getdimension = PetscDualSpaceGetDimension_Lagrange;
806: sp->ops->getnumdof = PetscDualSpaceGetNumDof_Lagrange;
807: sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_Lagrange;
808: sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Lagrange;
809: sp->ops->apply = PetscDualSpaceApplyDefault;
810: sp->ops->applyall = PetscDualSpaceApplyAllDefault;
811: sp->ops->createallpoints = PetscDualSpaceCreateAllPointsDefault;
812: return(0);
813: }
815: /*MC
816: PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals
818: Level: intermediate
820: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
821: M*/
823: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
824: {
825: PetscDualSpace_Lag *lag;
826: PetscErrorCode ierr;
830: PetscNewLog(sp,&lag);
831: sp->data = lag;
833: lag->numDof = NULL;
834: lag->simplexCell = PETSC_TRUE;
835: lag->tensorSpace = PETSC_FALSE;
836: lag->continuous = PETSC_TRUE;
838: PetscDualSpaceInitialize_Lagrange(sp);
839: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
840: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
841: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
842: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
843: return(0);
844: }