Actual source code: dspacelagrange.c

petsc-master 2019-10-20
Report Typos and Errors
  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,&section);
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(&section);
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: }