Actual source code: dspacebdm.c

petsc-master 2019-06-17
Report Typos and Errors
  1:  #include <petsc/private/petscfeimpl.h>
  2:  #include <petscdmplex.h>

  4: /*
  5: Let's work out BDM_1:

  7: The model basis is
  8:   \phi(x, y) = / a + b x + c y \
  9:                \ d + e x + f y /
 10: which is a 6 dimensional space. There are also six dual basis functions,
 11:   \psi_0(v) = \int^1_{-1} dx v(x, -1) \cdot <0, -1> (1-x)/2
 12:   \psi_1(v) = \int^1_{-1} dx v(x, -1) \cdot <0, -1> (1+x)/2
 13:   \psi_2(v) = 1/2 \int^1_{-1} ds v(-s, s) \cdot <1, 1> (1-s)/2 TODO I think the 1/2 is wrong here
 14:   \psi_3(v) = 1/2 \int^1_{-1} ds v(-s, s) \cdot <1, 1> (1+s)/2
 15:   \psi_4(v) = -\int^1_{-1} dy v(-1, y) \cdot <-1, 0> (1+y)/2
 16:   \psi_5(v) = -\int^1_{-1} dy v(-1, y) \cdot <-1, 0> (1-y)/2
 17: So we do the integrals
 18:   \psi_0(\phi) = \int^1_{-1} dx (f - d - e x) (1-x)/2 = (f - d) + e/3
 19:   \psi_1(\phi) = \int^1_{-1} dx (f - d - e x) (1+x)/2 = (f - d) - e/3
 20:   \psi_2(\phi) = \int^1_{-1} ds (a - b s + c s + d - e s + f s) (1-s)/2 = (a + d)/2 - (c + f - b - e)/6
 21:   \psi_3(\phi) = \int^1_{-1} ds (a - b s + c s + d - e s + f s) (1+s)/2 = (a + d)/2 + (c + f - b - e)/6
 22:   \psi_4(\phi) = \int^1_{-1} dy (b - a - c y) (1+y)/2 = (a - b) + c/3
 23:   \psi_5(\phi) = \int^1_{-1} dy (b - a - c y) (1-y)/2 = (a - b) - c/3
 24: so the nodal basis is
 25:   \phi_0 = / -(1+x)/2        \
 26:            \ 1/2 + 3/2 x + y /
 27:   \phi_1 = / 1+x                \
 28:            \ -1 - 3/2 x - 1/2 y /
 29:   \phi_2 = / 1+x      \
 30:            \ -(1+y)/2 /
 31:   \phi_3 = / -(1+x)/2 \
 32:            \ (1+y)    /
 33:   \phi_4 = / -1 - 1/2 x - 3/2 y \
 34:            \ (1+y)             /
 35:   \phi_5 = / 1/2 + x + 3/2 y \
 36:            \ -(1+y)/2        /
 37: */

 39: static PetscErrorCode PetscDualSpaceDestroy_BDM(PetscDualSpace sp)
 40: {
 41:   PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
 42:   PetscInt            i;
 43:   PetscErrorCode      ierr;

 46:   if (bdm->symmetries) {
 47:     PetscInt    **selfSyms = bdm->symmetries[0];
 48:     PetscScalar **selfFlps = bdm->flips[0];

 50:     if (selfSyms) {
 51:       PetscInt    **fSyms = &selfSyms[-bdm->selfSymOff], i;
 52:       PetscScalar **fFlps = &selfFlps[-bdm->selfSymOff];

 54:       for (i = 0; i < bdm->numSelfSym; i++) {
 55:         PetscFree(fSyms[i]);
 56:         PetscFree(fFlps[i]);
 57:       }
 58:       PetscFree(fSyms);
 59:       PetscFree(fFlps);
 60:     }
 61:     PetscFree(bdm->symmetries);
 62:     PetscFree(bdm->flips);
 63:   }
 64:   for (i = 0; i < bdm->height; i++) {
 65:     PetscDualSpaceDestroy(&bdm->subspaces[i]);
 66:   }
 67:   PetscFree(bdm->subspaces);
 68:   PetscFree(bdm->numDof);
 69:   PetscFree(bdm);
 70:   return(0);
 71: }

 73: static PetscErrorCode PetscDualSpaceBDMView_Ascii(PetscDualSpace sp, PetscViewer viewer)
 74: {

 78:   PetscViewerASCIIPrintf(viewer, "BDM(%D) dual space\n", sp->order);
 79:   return(0);
 80: }

 82: static PetscErrorCode PetscDualSpaceView_BDM(PetscDualSpace sp, PetscViewer viewer)
 83: {
 84:   PetscBool      iascii;

 90:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 91:   if (iascii) {PetscDualSpaceBDMView_Ascii(sp, viewer);}
 92:   return(0);
 93: }

 95: static PetscErrorCode PetscDualSpaceSetFromOptions_BDM(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
 96: {

100:   PetscOptionsHead(PetscOptionsObject,"PetscDualSpace BDM Options");
101:   PetscOptionsTail();
102:   return(0);
103: }

105: static PetscErrorCode PetscDualSpaceDuplicate_BDM(PetscDualSpace sp, PetscDualSpace *spNew)
106: {
107:   PetscInt       order, Nc;
108:   const char    *name;

112:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
113:   PetscObjectGetName((PetscObject) sp,     &name);
114:   PetscObjectSetName((PetscObject) *spNew,  name);
115:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACEBDM);
116:   PetscDualSpaceGetOrder(sp, &order);
117:   PetscDualSpaceSetOrder(*spNew, order);
118:   PetscDualSpaceGetNumComponents(sp, &Nc);
119:   PetscDualSpaceSetNumComponents(*spNew, Nc);
120:   return(0);
121: }

123: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_BDM(PetscDualSpace sp, PetscInt order, PetscInt *dim)
124: {
125:   PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
126:   PetscReal           D   = 1.0;
127:   PetscInt            n, d;
128:   PetscErrorCode      ierr;

131:   *dim = -1;
132:   DMGetDimension(sp->dm, &n);
133:   if (!n) {*dim = 0; return(0);}
134:   if (bdm->simplexCell) {
135:     for (d = 1; d <= n; ++d) {
136:       D *= ((PetscReal) (order+d))/d;
137:     }
138:     *dim = (PetscInt) (D + 0.5);
139:   } else {
140:     *dim = 1;
141:     for (d = 0; d < n; ++d) *dim *= (order+1);
142:   }
143:   if (!bdm->faceSpace) {
144:     *dim *= sp->Nc;
145:   }
146:   return(0);
147: }

149: static PetscErrorCode PetscDualSpaceCreateHeightSubspace_BDM(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
150: {
151:   PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
152:   PetscInt            order;
153:   PetscErrorCode      ierr;

156:   PetscDualSpaceGetOrder(sp, &order);
157:   if (!height) {
158:     PetscObjectReference((PetscObject) sp);
159:     *bdsp = sp;
160:   } else if (!order) {
161:     *bdsp = NULL;
162:   } else if (height == 1) {
163:     DM       dm, K;
164:     PetscInt dim;

166:     PetscDualSpaceGetDM(sp, &dm);
167:     DMGetDimension(dm, &dim);
168:     if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for dual space at height %d for dimension %d reference element\n", height, dim);
169:     PetscDualSpaceDuplicate(sp, bdsp);
170:     PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, bdm->simplexCell, &K);
171:     PetscDualSpaceSetDM(*bdsp, K);
172:     DMDestroy(&K);
173:     ((PetscDualSpace_BDM *) (*bdsp)->data)->faceSpace = PETSC_TRUE;
174:     PetscDualSpaceSetUp(*bdsp);
175:   } else {
176:     *bdsp = NULL;
177:   }
178:   return(0);
179: }

181: static PetscErrorCode PetscDualSpaceBDMCreateFaceFE(PetscDualSpace sp, PetscBool tensor, PetscInt faceDim, PetscInt order, PetscFE *fe)
182: {
183:   DM              K;
184:   PetscSpace      P;
185:   PetscDualSpace  Q;
186:   PetscQuadrature q;
187:   const PetscInt  Nc = 1;
188:   PetscInt        quadPointsPerEdge;
189:   PetscErrorCode  ierr;

191:   /* Create space */
192:   PetscSpaceCreate(PETSC_COMM_SELF, &P);
193:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
194:   PetscSpacePolynomialSetTensor(P, tensor);
195:   PetscSpaceSetNumComponents(P, Nc);
196:   PetscSpaceSetNumVariables(P, faceDim);
197:   PetscSpaceSetDegree(P, order, order);
198:   PetscSpaceSetUp(P);
199:   /* Create dual space */
200:   PetscDualSpaceCreate(PETSC_COMM_SELF, &Q);
201:   PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
202:   PetscDualSpaceCreateReferenceCell(Q, faceDim, tensor ? PETSC_FALSE : PETSC_TRUE, &K);
203:   PetscDualSpaceSetDM(Q, K);
204:   DMDestroy(&K);
205:   PetscDualSpaceSetNumComponents(Q, Nc);
206:   PetscDualSpaceSetOrder(Q, order);
207:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
208:   PetscDualSpaceSetUp(Q);
209:   /* Create element */
210:   PetscFECreate(PETSC_COMM_SELF, fe);
211:   PetscFESetType(*fe, PETSCFEBASIC);
212:   PetscFESetBasisSpace(*fe, P);
213:   PetscFESetDualSpace(*fe, Q);
214:   PetscFESetNumComponents(*fe, Nc);
215:   PetscFESetUp(*fe);
216:   PetscSpaceDestroy(&P);
217:   PetscDualSpaceDestroy(&Q);
218:   /* Create quadrature (with specified order if given) */
219:   quadPointsPerEdge = PetscMax(order + 1, 1);
220:   if (tensor) {PetscDTGaussTensorQuadrature(faceDim, 1, quadPointsPerEdge, -1.0, 1.0, &q);}
221:   else        {PetscDTGaussJacobiQuadrature(faceDim, 1, quadPointsPerEdge, -1.0, 1.0, &q);}
222:   PetscFESetQuadrature(*fe, q);
223:   PetscQuadratureDestroy(&q);
224:   return(0);
225: }

227: static PetscErrorCode PetscDualSpaceBDMCreateCellFE(PetscDualSpace sp, PetscBool tensor, PetscInt dim, PetscInt Nc, PetscInt order, PetscFE *fe)
228: {
229:   DM              K;
230:   PetscSpace      P;
231:   PetscDualSpace  Q;
232:   PetscQuadrature q;
233:   PetscInt        quadPointsPerEdge;
234:   PetscErrorCode  ierr;

236:   /* Create space */
237:   PetscSpaceCreate(PETSC_COMM_SELF, &P);
238:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
239:   PetscSpacePolynomialSetTensor(P, tensor);
240:   PetscSpaceSetNumComponents(P, Nc);
241:   PetscSpaceSetNumVariables(P, dim);
242:   PetscSpaceSetDegree(P, order, order);
243:   PetscSpaceSetUp(P);
244:   /* Create dual space */
245:   /* TODO Needs NED1 dual space */
246:   PetscDualSpaceCreate(PETSC_COMM_SELF, &Q);
247:   PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
248:   PetscDualSpaceCreateReferenceCell(Q, dim, tensor ? PETSC_FALSE : PETSC_TRUE, &K);
249:   PetscDualSpaceSetDM(Q, K);
250:   DMDestroy(&K);
251:   PetscDualSpaceSetNumComponents(Q, Nc);
252:   PetscDualSpaceSetOrder(Q, order);
253:   PetscDualSpaceLagrangeSetTensor(Q, tensor);
254:   PetscDualSpaceSetUp(Q);
255:   /* Create element */
256:   PetscFECreate(PETSC_COMM_SELF, fe);
257:   PetscFESetBasisSpace(*fe, P);
258:   PetscFESetDualSpace(*fe, Q);
259:   PetscFESetNumComponents(*fe, Nc);
260:   PetscFESetUp(*fe);
261:   PetscSpaceDestroy(&P);
262:   PetscDualSpaceDestroy(&Q);
263:   /* Create quadrature (with specified order if given) */
264:   quadPointsPerEdge = PetscMax(order + 1, 1);
265:   if (tensor) {PetscDTGaussTensorQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);}
266:   else        {PetscDTGaussJacobiQuadrature(dim, 1, quadPointsPerEdge, -1.0, 1.0, &q);}
267:   PetscFESetQuadrature(*fe, q);
268:   PetscQuadratureDestroy(&q);
269:   return(0);
270: }

272: static PetscErrorCode PetscDualSpaceSetUp_BDM(PetscDualSpace sp)
273: {
274:   PetscDualSpace_BDM *bdm     = (PetscDualSpace_BDM *) sp->data;
275:   DM                  dm      = sp->dm;
276:   PetscInt            order   = sp->order;
277:   PetscInt            Nc      = sp->Nc;
278:   PetscBool           faceSp  = bdm->faceSpace;
279:   MPI_Comm            comm;
280:   PetscSection        csection;
281:   Vec                 coordinates;
282:   PetscInt            depth, dim, cdim, pdimMax, pStart, pEnd, p, *pStratStart, *pStratEnd, coneSize, d, f = 0, c;
283:   PetscBool           simplex;
284:   PetscErrorCode      ierr;

287:   PetscObjectGetComm((PetscObject) sp, &comm);
288:   DMGetDimension(dm, &dim);
289:   DMGetCoordinateDim(dm, &cdim);
290:   DMPlexGetDepth(dm, &depth);
291:   DMPlexGetChart(dm, &pStart, &pEnd);
292:   if (depth != dim) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "BDM element requires interpolated meshes, but depth %D != topological dimension %D", depth, dim);
293:   if (!order)       SETERRQ(comm, PETSC_ERR_ARG_WRONG, "BDM elements not defined for order 0");
294:   if (!faceSp && Nc != cdim) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "BDM element has %D components != %D space dimension", Nc, cdim);
295:   PetscCalloc1(dim+1, &bdm->numDof);
296:   PetscMalloc2(depth+1, &pStratStart, depth+1, &pStratEnd);
297:   for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);}
298:   DMPlexGetConeSize(dm, pStratStart[depth], &coneSize);
299:   DMGetCoordinateSection(dm, &csection);
300:   DMGetCoordinatesLocal(dm, &coordinates);
301:   if      (coneSize == dim+1)   simplex = PETSC_TRUE;
302:   else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
303:   else SETERRQ(comm, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
304:   bdm->simplexCell = simplex;
305:   bdm->height      = 0;
306:   bdm->subspaces   = NULL;
307:   /* Create the height 1 subspace for every dimension */
308:   if (order > 0 && dim > 0) {
309:     PetscInt i;

311:     bdm->height = dim;
312:     PetscMalloc1(dim, &bdm->subspaces);
313:     PetscDualSpaceCreateHeightSubspace_BDM(sp, 1, &bdm->subspaces[0]);
314:     PetscDualSpaceSetUp(bdm->subspaces[0]);
315:     for (i = 1; i < dim; i++) {
316:       PetscDualSpaceGetHeightSubspace(bdm->subspaces[i-1], 1, &bdm->subspaces[i]);
317:       PetscObjectReference((PetscObject) bdm->subspaces[i]);
318:     }
319:   }
320:   PetscDualSpaceGetDimension_SingleCell_BDM(sp, sp->order, &pdimMax);
321:   pdimMax *= (pStratEnd[depth] - pStratStart[depth]);
322:   PetscMalloc1(pdimMax, &sp->functional);
323:   if (!dim) {
324:     bdm->numDof[0] = 0;
325:   } else {
326:     PetscFE      faceFE, cellFE;
327:     PetscSection section;
328:     CellRefiner  cellRefiner;
329:     PetscInt     faceDim = PetscMax(dim-1, 1), faceNum = 0;
330:     PetscReal   *v0 = NULL, *J = NULL, *detJ = NULL;

332:     PetscSectionCreate(PETSC_COMM_SELF, &section);
333:     PetscSectionSetChart(section, pStart, pEnd);
334:     if (!faceSp) {
335:       DMPlexGetCellRefiner_Internal(dm, &cellRefiner);
336:       CellRefinerGetAffineFaceTransforms_Internal(cellRefiner, NULL, &v0, &J, NULL, &detJ);
337:     }
338:     /* Create P_q(f) */
339:     PetscDualSpaceBDMCreateFaceFE(sp, simplex ? PETSC_FALSE : PETSC_TRUE, faceDim, order, &faceFE);
340:     /* Create NED^1_{q-1}(T) = P^d_{q-2} + S_{q-1}(T) */
341:     PetscDualSpaceBDMCreateCellFE(sp, simplex ? PETSC_FALSE : PETSC_TRUE, faceDim, Nc, order-1, &cellFE);
342:     for (p = pStart; p < pEnd; p++) {
343:       PetscBool isFace, isCell;
344:       PetscInt  d;

346:       for (d = 0; d < depth; d++) {if (p >= pStratStart[d] && p < pStratEnd[d]) break;}
347:       isFace = ((d == dim) &&  faceSp) || ((d == dim-1) && !faceSp) ? PETSC_TRUE : PETSC_FALSE;
348:       isCell = ((d == dim) && !faceSp) ? PETSC_TRUE : PETSC_FALSE;
349:       if (isFace) {
350:         PetscQuadrature  fq;
351:         PetscReal       *B, n[3];
352:         const PetscReal *fqpoints, *fqweights;
353:         PetscInt         faceDim = PetscMax(dim-1, 1), Nq, q, fdim, fb;

355:         if (cdim == 1) {n[0] = 0.; n[1] = 1.;}
356:         else           {DMPlexComputeCellGeometryFVM(dm, p, NULL, NULL, n);}
357:         PetscFEGetDefaultTabulation(faceFE, &B, NULL, NULL);
358:         PetscFEGetQuadrature(faceFE, &fq);
359:         PetscQuadratureGetData(fq, NULL, NULL, &Nq, &fqpoints, &fqweights);
360:         /* Create a dual basis vector for each basis function */
361:         PetscFEGetDimension(faceFE, &fdim);
362:         for (fb = 0; fb < fdim; ++fb, ++f) {
363:           PetscReal *qpoints, *qweights;

365:           PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
366:           PetscMalloc1(Nq*dim, &qpoints);
367:           PetscCalloc1(Nq*Nc,  &qweights);
368:           PetscQuadratureSetOrder(sp->functional[f], order);
369:           PetscQuadratureSetData(sp->functional[f], dim, Nc, Nq, qpoints, qweights);
370:           for (q = 0; q < Nq; ++q) {
371:             PetscInt g, h;

373:             if (faceDim < dim) {
374:               /* Transform quadrature points from face coordinates to cell coordinates */
375:               for (g = 0; g < dim; ++g) {
376:                 qpoints[q*dim+g] = v0[faceNum*dim+g];
377:                 for (h = 0; h < faceDim; ++h) qpoints[q*dim+g] += J[faceNum*dim*faceDim+g*faceDim+h] * fqpoints[q*faceDim+h];
378:               }
379:             } else {
380:               for (g = 0; g < dim; ++g) qpoints[q*dim+g] = fqpoints[q*faceDim+g];
381:             }
382:             /* Make Radon measure for integral \hat n p ds */
383:             for (c = 0; c < Nc; ++c) qweights[q*Nc+c] = B[q*fdim+fb]*n[c]*fqweights[q]*(detJ ? detJ[faceNum] : 1.0);
384:           }
385:         }
386:         bdm->numDof[d] = fdim;
387:         PetscSectionSetDof(section, p, bdm->numDof[d]);
388:         ++faceNum;
389:       }
390:       if (order < 2) continue;
391:       if (isCell) {
392:         PetscSpace       csp;
393:         PetscQuadrature  cq;
394:         PetscReal       *B;
395:         const PetscReal *cqpoints, *cqweights;
396:         PetscInt         Nq, q, cdim, cb;

398:         PetscFEGetBasisSpace(cellFE, &csp);
399:         PetscFEGetQuadrature(cellFE, &cq);
400:         PetscQuadratureGetData(cq, NULL, NULL, &Nq, &cqpoints, &cqweights);
401:         /* Create a dual basis vector for each basis function */
402:         PetscSpaceGetDimension(csp, &cdim);
403:         PetscMalloc1(Nq*cdim*Nc, &B);
404:         PetscSpaceEvaluate(csp, Nq, cqpoints, B, NULL, NULL);
405:         for (cb = 0; cb < cdim; ++cb, ++f) {
406:           PetscReal *qpoints, *qweights;

408:           PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
409:           PetscMalloc1(Nq*dim, &qpoints);
410:           PetscCalloc1(Nq*Nc,  &qweights);
411:           PetscQuadratureSetOrder(sp->functional[f], order);
412:           PetscQuadratureSetData(sp->functional[f], dim, Nc, Nq, qpoints, qweights);
413:           PetscMemcpy(qpoints, cqpoints, Nq*dim * sizeof(PetscReal));
414:           for (q = 0; q < Nq; ++q) {
415:             /* Make Radon measure for integral p dx */
416:             for (c = 0; c < Nc; ++c) qweights[q*Nc+c] = B[(q*cdim+cb)*Nc+c]*cqweights[q*Nc+c];
417:           }
418:         }
419:         PetscFree(B);
420:         bdm->numDof[d] = cdim;
421:         PetscSectionSetDof(section, p, bdm->numDof[d]);
422:       }
423:     }
424:     PetscFEDestroy(&faceFE);
425:     PetscFEDestroy(&cellFE);
426:     PetscFree(v0);
427:     PetscFree(J);
428:     PetscFree(detJ);
429:     PetscSectionSetUp(section);
430:     { /* reorder to closure order */
431:       PetscQuadrature *reorder = NULL;
432:       PetscInt        *key, count;

434:       PetscCalloc1(f, &key);
435:       PetscMalloc1(f, &reorder);
436:       for (p = pStratStart[depth], count = 0; p < pStratEnd[depth]; p++) {
437:         PetscInt *closure = NULL, closureSize, c;

439:         DMPlexGetTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
440:         for (c = 0; c < closureSize*2; c += 2) {
441:           PetscInt point = closure[c], dof, off, i;

443:           PetscSectionGetDof(section, point, &dof);
444:           PetscSectionGetOffset(section, point, &off);
445:           for (i = 0; i < dof; ++i) {
446:             PetscInt fi = off + i;

448:             if (!key[fi]) {
449:               key[fi] = 1;
450:               reorder[count++] = sp->functional[fi];
451:             }
452:           }
453:         }
454:         DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);
455:       }
456:       PetscFree(key);
457:       PetscFree(sp->functional);
458:       sp->functional = reorder;
459:     }
460:     PetscSectionDestroy(&section);
461:   }
462:   if (pStratEnd[depth] == 1 && f != pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %D not equal to dimension %D", f, pdimMax);
463:   if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %D is greater than max size %D", f, pdimMax);
464:   PetscFree2(pStratStart, pStratEnd);
465:   return(0);
466: }

468: static PetscErrorCode PetscDualSpaceGetDimension_BDM(PetscDualSpace sp, PetscInt *dim)
469: {
470:   DM              K;
471:   const PetscInt *numDof;
472:   PetscInt        spatialDim, cEnd, size = 0, d;
473:   PetscErrorCode  ierr;

476:   PetscDualSpaceGetDM(sp, &K);
477:   PetscDualSpaceGetNumDof(sp, &numDof);
478:   DMGetDimension(K, &spatialDim);
479:   DMPlexGetHeightStratum(K, 0, NULL, &cEnd);
480:   if (cEnd == 1) {PetscDualSpaceGetDimension_SingleCell_BDM(sp, sp->order, dim); return(0);}
481:   for (d = 0; d <= spatialDim; ++d) {
482:     PetscInt pStart, pEnd;

484:     DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
485:     size += (pEnd-pStart)*numDof[d];
486:   }
487:   *dim = size;
488:   return(0);
489: }

491: static PetscErrorCode PetscDualSpaceGetNumDof_BDM(PetscDualSpace sp, const PetscInt **numDof)
492: {
493:   PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;

496:   *numDof = bdm->numDof;
497:   return(0);
498: }

500: static PetscErrorCode PetscDualSpaceGetHeightSubspace_BDM(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
501: {
502:   PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
503:   PetscErrorCode      ierr;

508:   if (height == 0) {
509:     *bdsp = sp;
510:   } else {
511:     DM       dm;
512:     PetscInt dim;

514:     PetscDualSpaceGetDM(sp, &dm);
515:     DMGetDimension(dm, &dim);
516:     if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for dual space at height %D for dimension %D reference element\n", height, dim);
517:     if (height <= bdm->height) {*bdsp = bdm->subspaces[height-1];}
518:     else                       {*bdsp = NULL;}
519:   }
520:   return(0);
521: }

523: #define BaryIndex(perEdge,a,b,c) (((b)*(2*perEdge+1-(b)))/2)+(c)

525: #define CartIndex(perEdge,a,b) (perEdge*(a)+b)

527: static PetscErrorCode PetscDualSpaceGetSymmetries_BDM(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****rots)
528: {

530:   PetscDualSpace_BDM *bdm = (PetscDualSpace_BDM *) sp->data;
531:   PetscInt            dim, order, p, Nc;
532:   PetscErrorCode      ierr;

535:   PetscDualSpaceGetOrder(sp, &order);
536:   PetscDualSpaceGetNumComponents(sp, &Nc);
537:   DMGetDimension(sp->dm, &dim);
538:   if (dim < 1) return(0);
539:   if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "BDM symmetries not implemented for dim = %D > 3", dim);
540:   if (!bdm->symmetries) { /* store symmetries */
541:     PetscInt    ***symmetries, **cellSymmetries;
542:     PetscScalar ***flips,      **cellFlips;
543:     PetscInt       numPoints, numFaces, d;

545:     if (bdm->simplexCell) {
546:       numPoints = 1;
547:       for (d = 0; d < dim; d++) numPoints = numPoints * 2 + 1;
548:       numFaces  = 1 + dim;
549:     } else {
550:       numPoints = PetscPowInt(3,dim);
551:       numFaces  = 2 * dim;
552:     }
553:     PetscCalloc1(numPoints, &symmetries);
554:     PetscCalloc1(numPoints, &flips);
555:     /* compute self symmetries */
556:     bdm->numSelfSym = 2 * numFaces;
557:     bdm->selfSymOff = numFaces;
558:     PetscCalloc1(2*numFaces, &cellSymmetries);
559:     PetscCalloc1(2*numFaces, &cellFlips);
560:     /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
561:     symmetries[0] = &cellSymmetries[bdm->selfSymOff];
562:     flips[0]      = &cellFlips[bdm->selfSymOff];
563:     switch (dim) {
564:     case 1: /* Edge symmetries */
565:     {
566:       PetscScalar *invert;
567:       PetscInt    *reverse;
568:       PetscInt     dofPerEdge = order+1, eNc = 1 /* ??? */, i, j;

570:       PetscMalloc1(dofPerEdge*eNc, &reverse);
571:       PetscMalloc1(dofPerEdge*eNc, &invert);
572:       for (i = 0; i < dofPerEdge; ++i) {
573:         for (j = 0; j < eNc; ++j) {
574:           reverse[i*eNc + j] = eNc * (dofPerEdge - 1 - i) + j;
575:           invert[i*eNc + j]  = -1.0;
576:         }
577:       }
578:       symmetries[0][-2] = reverse;
579:       flips[0][-2]      = invert;

581:       /* yes, this is redundant, but it makes it easier to cleanup if I don't have to worry about what not to free */
582:       PetscMalloc1(dofPerEdge*eNc, &reverse);
583:       PetscMalloc1(dofPerEdge*eNc, &invert);
584:       for (i = 0; i < dofPerEdge; i++) {
585:         for (j = 0; j < eNc; j++) {
586:           reverse[i*eNc + j] = eNc * (dofPerEdge - 1 - i) + j;
587:           invert[i*eNc + j]  = -1.0;
588:         }
589:       }
590:       symmetries[0][1] = reverse;
591:       flips[0][1]      = invert;
592:       break;
593:     }
594:     case 2: /* Face symmetries  */
595:     {
596:       PetscInt dofPerEdge = bdm->simplexCell ? (order - 2) : (order - 1), s;
597:       PetscInt dofPerFace;

599:       for (s = -numFaces; s < numFaces; s++) {
600:         PetscScalar *flp;
601:         PetscInt    *sym, fNc = 1, i, j, k, l;

603:         if (!s) continue;
604:         if (bdm->simplexCell) {
605:           dofPerFace = (dofPerEdge * (dofPerEdge + 1))/2;
606:           PetscMalloc1(fNc*dofPerFace,&sym);
607:           PetscMalloc1(fNc*dofPerFace,&flp);
608:           for (j = 0, l = 0; j < dofPerEdge; j++) {
609:             for (k = 0; k < dofPerEdge - j; k++, l++) {
610:               i = dofPerEdge - 1 - j - k;
611:               switch (s) {
612:               case -3:
613:                 sym[fNc*l] = BaryIndex(dofPerEdge,i,k,j);
614:                 break;
615:               case -2:
616:                 sym[fNc*l] = BaryIndex(dofPerEdge,j,i,k);
617:                 break;
618:               case -1:
619:                 sym[fNc*l] = BaryIndex(dofPerEdge,k,j,i);
620:                 break;
621:               case 1:
622:                 sym[fNc*l] = BaryIndex(dofPerEdge,k,i,j);
623:                 break;
624:               case 2:
625:                 sym[fNc*l] = BaryIndex(dofPerEdge,j,k,i);
626:                 break;
627:               }
628:               flp[fNc*l] = s < 0 ? -1.0 : 1.0;
629:             }
630:           }
631:         } else {
632:           dofPerFace = dofPerEdge * dofPerEdge;
633:           PetscMalloc1(fNc*dofPerFace,&sym);
634:           PetscMalloc1(fNc*dofPerFace,&flp);
635:           for (j = 0, l = 0; j < dofPerEdge; j++) {
636:             for (k = 0; k < dofPerEdge; k++, l++) {
637:               switch (s) {
638:               case -4:
639:                 sym[fNc*l] = CartIndex(dofPerEdge,k,j);
640:                 break;
641:               case -3:
642:                 sym[fNc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),k);
643:                 break;
644:               case -2:
645:                 sym[fNc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),(dofPerEdge - 1 - j));
646:                 break;
647:               case -1:
648:                 sym[fNc*l] = CartIndex(dofPerEdge,j,(dofPerEdge - 1 - k));
649:                 break;
650:               case 1:
651:                 sym[fNc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),j);
652:                 break;
653:               case 2:
654:                 sym[fNc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),(dofPerEdge - 1 - k));
655:                 break;
656:               case 3:
657:                 sym[fNc*l] = CartIndex(dofPerEdge,k,(dofPerEdge - 1 - j));
658:                 break;
659:               }
660:               flp[fNc*l] = s < 0 ? -1.0 : 1.0;
661:             }
662:           }
663:         }
664:         for (i = 0; i < dofPerFace; i++) {
665:           sym[fNc*i] *= fNc;
666:           for (j = 1; j < fNc; j++) {
667:             sym[fNc*i+j] = sym[fNc*i] + j;
668:             flp[fNc*i+j] = flp[fNc*i];
669:           }
670:         }
671:         symmetries[0][s] = sym;
672:         flips[0][s]      = flp;
673:       }
674:       break;
675:     }
676:     default: SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "No symmetries for point of dimension %D", dim);
677:     }
678:     /* Copy subspace symmetries */
679:     {
680:       PetscDualSpace       hsp;
681:       DM                   K;
682:       const PetscInt    ***hsymmetries;
683:       const PetscScalar ***hflips;

685:       PetscDualSpaceGetHeightSubspace(sp, 1, &hsp);
686:       PetscDualSpaceGetSymmetries(hsp, &hsymmetries, &hflips);
687:       if (hsymmetries || hflips) {
688:         PetscBool      *seen;
689:         const PetscInt *cone;
690:         PetscInt        KclosureSize, *Kclosure = NULL;

692:         PetscDualSpaceGetDM(sp, &K);
693:         PetscCalloc1(numPoints, &seen);
694:         DMPlexGetCone(K, 0, &cone);
695:         DMPlexGetTransitiveClosure(K, 0, PETSC_TRUE, &KclosureSize, &Kclosure);
696:         for (p = 0; p < numFaces; ++p) {
697:           PetscInt closureSize, *closure = NULL, q;

699:           DMPlexGetTransitiveClosure(K, cone[p], PETSC_TRUE, &closureSize, &closure);
700:           for (q = 0; q < closureSize; ++q) {
701:             PetscInt point = closure[q*2], r;

703:             if (!seen[point]) {
704:               for (r = 0; r < KclosureSize; ++r) {
705:                 if (Kclosure[r*2] == point) break;
706:               }
707:               seen[point] = PETSC_TRUE;
708:               symmetries[r] = (PetscInt **)    hsymmetries[q];
709:               flips[r]      = (PetscScalar **) hflips[q];
710:             }
711:           }
712:           DMPlexRestoreTransitiveClosure(K, cone[p], PETSC_TRUE, &closureSize, &closure);
713:         }
714:         DMPlexRestoreTransitiveClosure(K, 0, PETSC_TRUE, &KclosureSize, &Kclosure);
715:         PetscFree(seen);
716:       }
717:     }
718:     bdm->symmetries = symmetries;
719:     bdm->flips      = flips;
720:   }
721:   if (perms) *perms = (const PetscInt ***)    bdm->symmetries;
722:   if (rots)  *rots  = (const PetscScalar ***) bdm->flips;
723:   return(0);
724: }

726: static PetscErrorCode PetscDualSpaceInitialize_BDM(PetscDualSpace sp)
727: {
729:   sp->ops->destroy           = PetscDualSpaceDestroy_BDM;
730:   sp->ops->view              = PetscDualSpaceView_BDM;
731:   sp->ops->setfromoptions    = PetscDualSpaceSetFromOptions_BDM;
732:   sp->ops->duplicate         = PetscDualSpaceDuplicate_BDM;
733:   sp->ops->setup             = PetscDualSpaceSetUp_BDM;
734:   sp->ops->getdimension      = PetscDualSpaceGetDimension_BDM;
735:   sp->ops->getnumdof         = PetscDualSpaceGetNumDof_BDM;
736:   sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_BDM;
737:   sp->ops->getsymmetries     = PetscDualSpaceGetSymmetries_BDM;
738:   sp->ops->apply             = PetscDualSpaceApplyDefault;
739:   sp->ops->applyall          = PetscDualSpaceApplyAllDefault;
740:   sp->ops->createallpoints   = PetscDualSpaceCreateAllPointsDefault;
741:   return(0);
742: }
743: /*MC
744:   PETSCDUALSPACEBDM = "bdm" - A PetscDualSpace object that encapsulates a dual space for Brezzi-Douglas-Marini elements

746:   Level: intermediate

748: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
749: M*/

751: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_BDM(PetscDualSpace sp)
752: {
753:   PetscDualSpace_BDM *bdm;
754:   PetscErrorCode      ierr;

758:   PetscNewLog(sp, &bdm);
759:   sp->data = bdm;
760:   sp->k    = 3;

762:   bdm->numDof      = NULL;
763:   bdm->simplexCell = PETSC_TRUE;

765:   PetscDualSpaceInitialize_BDM(sp);
766:   return(0);
767: }