Actual source code: dspacebdm.c

petsc-master 2020-01-21
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:         PetscTabulation  T;
352:         PetscReal       *B, n[3];
353:         const PetscReal *fqpoints, *fqweights;
354:         PetscInt         faceDim = PetscMax(dim-1, 1), Nq, q, fdim, fb;

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

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

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

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

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

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

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

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

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

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

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

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

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

498:   *numDof = bdm->numDof;
499:   return(0);
500: }

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

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

516:     PetscDualSpaceGetDM(sp, &dm);
517:     DMGetDimension(dm, &dim);
518:     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);
519:     if (height <= bdm->height) {*bdsp = bdm->subspaces[height-1];}
520:     else                       {*bdsp = NULL;}
521:   }
522:   return(0);
523: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

748:   Level: intermediate

750: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
751: M*/

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

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

764:   bdm->numDof      = NULL;
765:   bdm->simplexCell = PETSC_TRUE;

767:   PetscDualSpaceInitialize_BDM(sp);
768:   return(0);
769: }