Actual source code: dspacelagrange.c

petsc-3.10.2 2018-10-09
Report Typos and Errors
  1:  #include <petsc/private/petscfeimpl.h>
  2:  #include <petscdmplex.h>

  4: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor)
  5: {
  6:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

  9:   *tensor = lag->tensorSpace;
 10:   return(0);
 11: }

 13: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor)
 14: {
 15:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;

 18:   lag->tensorSpace = tensor;
 19:   return(0);
 20: }

 22: /*
 23:   LatticePointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to at most 'max'.
 24:                                        Ordering is lexicographic with lowest index as least significant in ordering.
 25:                                        e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,0}.

 27:   Input Parameters:
 28: + len - The length of the tuple
 29: . max - The maximum sum
 30: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition

 32:   Output Parameter:
 33: . tup - A tuple of len integers whos sum is at most 'max'
 34: */
 35: static PetscErrorCode LatticePointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
 36: {
 38:   while (len--) {
 39:     max -= tup[len];
 40:     if (!max) {
 41:       tup[len] = 0;
 42:       break;
 43:     }
 44:   }
 45:   tup[++len]++;
 46:   return(0);
 47: }

 49: /*
 50:   TensorPointLexicographic_Internal - Returns all tuples of size 'len' with nonnegative integers that are all less than or equal to 'max'.
 51:                                       Ordering is lexicographic with lowest index as least significant in ordering.
 52:                                       e.g. for len == 2 and max == 2, this will return, in order, {0,0}, {1,0}, {2,0}, {0,1}, {1,1}, {2,1}, {0,2}, {1,2}, {2,2}.

 54:   Input Parameters:
 55: + len - The length of the tuple
 56: . max - The maximum value
 57: - tup - A tuple of length len+1: tup[len] > 0 indicates a stopping condition

 59:   Output Parameter:
 60: . tup - A tuple of len integers whos sum is at most 'max'
 61: */
 62: static PetscErrorCode TensorPointLexicographic_Internal(PetscInt len, PetscInt max, PetscInt tup[])
 63: {
 64:   PetscInt       i;

 67:   for (i = 0; i < len; i++) {
 68:     if (tup[i] < max) {
 69:       break;
 70:     } else {
 71:       tup[i] = 0;
 72:     }
 73:   }
 74:   tup[i]++;
 75:   return(0);
 76: }


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

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

 83: static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
 84: {

 86:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
 87:   PetscInt           dim, order, p, Nc;
 88:   PetscErrorCode     ierr;

 91:   PetscDualSpaceGetOrder(sp,&order);
 92:   PetscDualSpaceGetNumComponents(sp,&Nc);
 93:   DMGetDimension(sp->dm,&dim);
 94:   if (!dim || !lag->continuous || order < 3) return(0);
 95:   if (dim > 3) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Lagrange symmetries not implemented for dim = %D > 3",dim);
 96:   if (!lag->symmetries) { /* store symmetries */
 97:     PetscDualSpace hsp;
 98:     DM             K;
 99:     PetscInt       numPoints = 1, d;
100:     PetscInt       numFaces;
101:     PetscInt       ***symmetries;
102:     const PetscInt ***hsymmetries;

104:     if (lag->simplexCell) {
105:       numFaces = 1 + dim;
106:       for (d = 0; d < dim; d++) numPoints = numPoints * 2 + 1;
107:     }
108:     else {
109:       numPoints = PetscPowInt(3,dim);
110:       numFaces  = 2 * dim;
111:     }
112:     PetscCalloc1(numPoints,&symmetries);
113:     if (0 < dim && dim < 3) { /* compute self symmetries */
114:       PetscInt **cellSymmetries;

116:       lag->numSelfSym = 2 * numFaces;
117:       lag->selfSymOff = numFaces;
118:       PetscCalloc1(2*numFaces,&cellSymmetries);
119:       /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */
120:       symmetries[0] = &cellSymmetries[numFaces];
121:       if (dim == 1) {
122:         PetscInt dofPerEdge = order - 1;

124:         if (dofPerEdge > 1) {
125:           PetscInt i, j, *reverse;

127:           PetscMalloc1(dofPerEdge*Nc,&reverse);
128:           for (i = 0; i < dofPerEdge; i++) {
129:             for (j = 0; j < Nc; j++) {
130:               reverse[i*Nc + j] = Nc * (dofPerEdge - 1 - i) + j;
131:             }
132:           }
133:           symmetries[0][-2] = reverse;

135:           /* yes, this is redundant, but it makes it easier to cleanup if I don't have to worry about what not to free */
136:           PetscMalloc1(dofPerEdge*Nc,&reverse);
137:           for (i = 0; i < dofPerEdge; i++) {
138:             for (j = 0; j < Nc; j++) {
139:               reverse[i*Nc + j] = Nc * (dofPerEdge - 1 - i) + j;
140:             }
141:           }
142:           symmetries[0][1] = reverse;
143:         }
144:       } else {
145:         PetscInt dofPerEdge = lag->simplexCell ? (order - 2) : (order - 1), s;
146:         PetscInt dofPerFace;

148:         if (dofPerEdge > 1) {
149:           for (s = -numFaces; s < numFaces; s++) {
150:             PetscInt *sym, i, j, k, l;

152:             if (!s) continue;
153:             if (lag->simplexCell) {
154:               dofPerFace = (dofPerEdge * (dofPerEdge + 1))/2;
155:               PetscMalloc1(Nc*dofPerFace,&sym);
156:               for (j = 0, l = 0; j < dofPerEdge; j++) {
157:                 for (k = 0; k < dofPerEdge - j; k++, l++) {
158:                   i = dofPerEdge - 1 - j - k;
159:                   switch (s) {
160:                   case -3:
161:                     sym[Nc*l] = BaryIndex(dofPerEdge,i,k,j);
162:                     break;
163:                   case -2:
164:                     sym[Nc*l] = BaryIndex(dofPerEdge,j,i,k);
165:                     break;
166:                   case -1:
167:                     sym[Nc*l] = BaryIndex(dofPerEdge,k,j,i);
168:                     break;
169:                   case 1:
170:                     sym[Nc*l] = BaryIndex(dofPerEdge,k,i,j);
171:                     break;
172:                   case 2:
173:                     sym[Nc*l] = BaryIndex(dofPerEdge,j,k,i);
174:                     break;
175:                   }
176:                 }
177:               }
178:             } else {
179:               dofPerFace = dofPerEdge * dofPerEdge;
180:               PetscMalloc1(Nc*dofPerFace,&sym);
181:               for (j = 0, l = 0; j < dofPerEdge; j++) {
182:                 for (k = 0; k < dofPerEdge; k++, l++) {
183:                   switch (s) {
184:                   case -4:
185:                     sym[Nc*l] = CartIndex(dofPerEdge,k,j);
186:                     break;
187:                   case -3:
188:                     sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),k);
189:                     break;
190:                   case -2:
191:                     sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),(dofPerEdge - 1 - j));
192:                     break;
193:                   case -1:
194:                     sym[Nc*l] = CartIndex(dofPerEdge,j,(dofPerEdge - 1 - k));
195:                     break;
196:                   case 1:
197:                     sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - k),j);
198:                     break;
199:                   case 2:
200:                     sym[Nc*l] = CartIndex(dofPerEdge,(dofPerEdge - 1 - j),(dofPerEdge - 1 - k));
201:                     break;
202:                   case 3:
203:                     sym[Nc*l] = CartIndex(dofPerEdge,k,(dofPerEdge - 1 - j));
204:                     break;
205:                   }
206:                 }
207:               }
208:             }
209:             for (i = 0; i < dofPerFace; i++) {
210:               sym[Nc*i] *= Nc;
211:               for (j = 1; j < Nc; j++) {
212:                 sym[Nc*i+j] = sym[Nc*i] + j;
213:               }
214:             }
215:             symmetries[0][s] = sym;
216:           }
217:         }
218:       }
219:     }
220:     PetscDualSpaceGetHeightSubspace(sp,1,&hsp);
221:     PetscDualSpaceGetSymmetries(hsp,&hsymmetries,NULL);
222:     if (hsymmetries) {
223:       PetscBool      *seen;
224:       const PetscInt *cone;
225:       PetscInt       KclosureSize, *Kclosure = NULL;

227:       PetscDualSpaceGetDM(sp,&K);
228:       PetscCalloc1(numPoints,&seen);
229:       DMPlexGetCone(K,0,&cone);
230:       DMPlexGetTransitiveClosure(K,0,PETSC_TRUE,&KclosureSize,&Kclosure);
231:       for (p = 0; p < numFaces; p++) {
232:         PetscInt closureSize, *closure = NULL, q;

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

238:           if(!seen[point]) {
239:             for (r = 0; r < KclosureSize; r++) {
240:               if (Kclosure[2 * r] == point) break;
241:             }
242:             seen[point] = PETSC_TRUE;
243:             symmetries[r] = (PetscInt **) hsymmetries[q];
244:           }
245:         }
246:         DMPlexRestoreTransitiveClosure(K,cone[p],PETSC_TRUE,&closureSize,&closure);
247:       }
248:       DMPlexRestoreTransitiveClosure(K,0,PETSC_TRUE,&KclosureSize,&Kclosure);
249:       PetscFree(seen);
250:     }
251:     lag->symmetries = symmetries;
252:   }
253:   if (perms) *perms = (const PetscInt ***) lag->symmetries;
254:   return(0);
255: }

257: /*@C
258:   PetscDualSpaceGetSymmetries - Returns a description of the symmetries of this basis

260:   Not collective

262:   Input Parameter:
263: . sp - the PetscDualSpace object

265:   Output Parameters:
266: + perms - Permutations of the local degrees of freedom, parameterized by the point orientation
267: - flips - Sign reversal of the local degrees of freedom, parameterized by the point orientation

269:   Note: The permutation and flip arrays are organized in the following way
270: $ perms[p][ornt][dof # on point] = new local dof #
271: $ flips[p][ornt][dof # on point] = reversal or not

273:   Level: developer

275: .seealso: PetscDualSpaceSetSymmetries()
276: @*/
277: PetscErrorCode PetscDualSpaceGetSymmetries(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
278: {

283:   if (perms) {
285:     *perms = NULL;
286:   }
287:   if (flips) {
289:     *flips = NULL;
290:   }
291:   if (sp->ops->getsymmetries) {
292:     (sp->ops->getsymmetries)(sp,perms,flips);
293:   }
294:   return(0);
295: }

297: static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer)
298: {
299:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
300:   PetscErrorCode      ierr;

303:   PetscViewerASCIIPrintf(viewer, "%s %sLagrange dual space of order %D", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "Tensor " : "", sp->order, sp->Nc);
304:   if (sp->Nc > 1) {PetscViewerASCIIPrintf(viewer, " with %D components", sp->Nc);}
305:   PetscViewerASCIIPrintf(viewer, "\n");
306:   return(0);
307: }

309: PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
310: {
311:   PetscBool      iascii;

317:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
318:   if (iascii) {PetscDualSpaceLagrangeView_Ascii(sp, viewer);}
319:   return(0);
320: }

322: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
323: {
324:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
325:   PetscReal           D   = 1.0;
326:   PetscInt            n, i;
327:   PetscErrorCode      ierr;

330:   *dim = -1;                    /* Ensure that the compiler knows *dim is set. */
331:   DMGetDimension(sp->dm, &n);
332:   if (!lag->tensorSpace) {
333:     for (i = 1; i <= n; ++i) {
334:       D *= ((PetscReal) (order+i))/i;
335:     }
336:     *dim = (PetscInt) (D + 0.5);
337:   } else {
338:     *dim = 1;
339:     for (i = 0; i < n; ++i) *dim *= (order+1);
340:   }
341:   *dim *= sp->Nc;
342:   return(0);
343: }

345: static PetscErrorCode PetscDualSpaceCreateHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
346: {
347:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
348:   PetscBool          continuous, tensor;
349:   PetscInt           order;
350:   PetscErrorCode     ierr;

355:   PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
356:   PetscDualSpaceGetOrder(sp,&order);
357:   if (height == 0) {
358:     PetscObjectReference((PetscObject)sp);
359:     *bdsp = sp;
360:   } else if (continuous == PETSC_FALSE || !order) {
361:     *bdsp = NULL;
362:   } else {
363:     DM dm, K;
364:     PetscInt dim;

366:     PetscDualSpaceGetDM(sp,&dm);
367:     DMGetDimension(dm,&dim);
368:     if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Asked for dual space at height %d for dimension %d reference element\n",height,dim);}
369:     PetscDualSpaceDuplicate(sp,bdsp);
370:     PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, lag->simplexCell, &K);
371:     PetscDualSpaceSetDM(*bdsp, K);
372:     DMDestroy(&K);
373:     PetscDualSpaceLagrangeGetTensor(sp,&tensor);
374:     PetscDualSpaceLagrangeSetTensor(*bdsp,tensor);
375:     PetscDualSpaceSetUp(*bdsp);
376:   }
377:   return(0);
378: }

380: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
381: {
382:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
383:   DM                  dm    = sp->dm;
384:   PetscInt            order = sp->order;
385:   PetscInt            Nc    = sp->Nc;
386:   PetscBool           continuous;
387:   PetscSection        csection;
388:   Vec                 coordinates;
389:   PetscReal          *qpoints, *qweights;
390:   PetscInt            depth, dim, pdimMax, pStart, pEnd, p, *pStratStart, *pStratEnd, coneSize, d, f = 0, c;
391:   PetscBool           simplex, tensorSpace;
392:   PetscErrorCode      ierr;

395:   /* Classify element type */
396:   if (!order) lag->continuous = PETSC_FALSE;
397:   continuous = lag->continuous;
398:   DMGetDimension(dm, &dim);
399:   DMPlexGetDepth(dm, &depth);
400:   DMPlexGetChart(dm, &pStart, &pEnd);
401:   PetscCalloc1(dim+1, &lag->numDof);
402:   PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);
403:   for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);}
404:   DMPlexGetConeSize(dm, pStratStart[depth], &coneSize);
405:   DMGetCoordinateSection(dm, &csection);
406:   DMGetCoordinatesLocal(dm, &coordinates);
407:   if (depth == 1) {
408:     if      (coneSize == dim+1)    simplex = PETSC_TRUE;
409:     else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
410:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
411:   } else if (depth == dim) {
412:     if      (coneSize == dim+1)   simplex = PETSC_TRUE;
413:     else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
414:     else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
415:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support cell-vertex meshes or interpolated meshes");
416:   lag->simplexCell = simplex;
417:   if (dim > 1 && continuous && lag->simplexCell == lag->tensorSpace) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Mismatching simplex/tensor cells and spaces only allowed for discontinuous elements");
418:   tensorSpace    = lag->tensorSpace;
419:   lag->height    = 0;
420:   lag->subspaces = NULL;
421:   if (continuous && sp->order > 0 && dim > 0) {
422:     PetscInt i;

424:     lag->height = dim;
425:     PetscMalloc1(dim,&lag->subspaces);
426:     PetscDualSpaceCreateHeightSubspace_Lagrange(sp,1,&lag->subspaces[0]);
427:     PetscDualSpaceSetUp(lag->subspaces[0]);
428:     for (i = 1; i < dim; i++) {
429:       PetscDualSpaceGetHeightSubspace(lag->subspaces[i-1],1,&lag->subspaces[i]);
430:       PetscObjectReference((PetscObject)(lag->subspaces[i]));
431:     }
432:   }
433:   PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
434:   pdimMax *= (pStratEnd[depth] - pStratStart[depth]);
435:   PetscMalloc1(pdimMax, &sp->functional);
436:   if (!dim) {
437:     for (c = 0; c < Nc; ++c) {
438:       PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
439:       PetscCalloc1(Nc, &qweights);
440:       PetscQuadratureSetOrder(sp->functional[f], 0);
441:       PetscQuadratureSetData(sp->functional[f], 0, Nc, 1, NULL, qweights);
442:       qweights[c] = 1.0;
443:       ++f;
444:       lag->numDof[0]++;
445:     }
446:   } else {
447:     PetscInt     *tup;
448:     PetscReal    *v0, *hv0, *J, *invJ, detJ, hdetJ;
449:     PetscSection section;

451:     PetscSectionCreate(PETSC_COMM_SELF,&section);
452:     PetscSectionSetChart(section,pStart,pEnd);
453:     PetscCalloc5(dim+1,&tup,dim,&v0,dim,&hv0,dim*dim,&J,dim*dim,&invJ);
454:     for (p = pStart; p < pEnd; p++) {
455:       PetscInt       pointDim, d, nFunc = 0;
456:       PetscDualSpace hsp;

458:       DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);
459:       for (d = 0; d < depth; d++) {if (p >= pStratStart[d] && p < pStratEnd[d]) break;}
460:       pointDim = (depth == 1 && d == 1) ? dim : d;
461:       hsp = ((pointDim < dim) && lag->subspaces) ? lag->subspaces[dim - pointDim - 1] : NULL;
462:       if (hsp) {
463:         PetscDualSpace_Lag *hlag = (PetscDualSpace_Lag *) hsp->data;
464:         DM                 hdm;

466:         PetscDualSpaceGetDM(hsp,&hdm);
467:         DMPlexComputeCellGeometryFEM(hdm, 0, NULL, hv0, NULL, NULL, &hdetJ);
468:         nFunc = lag->numDof[pointDim] = hlag->numDof[pointDim];
469:       }
470:       if (pointDim == dim) {
471:         /* Cells, create for self */
472:         PetscInt     orderEff = continuous ? (!tensorSpace ? order-1-dim : order-2) : order;
473:         PetscReal    denom    = continuous ? order : (!tensorSpace ? order+1+dim : order+2);
474:         PetscReal    numer    = (!simplex || !tensorSpace) ? 2. : (2./dim);
475:         PetscReal    dx = numer/denom;
476:         PetscInt     cdim, d, d2;

478:         if (orderEff < 0) continue;
479:         PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, orderEff, &cdim);
480:         PetscMemzero(tup,(dim+1)*sizeof(PetscInt));
481:         if (!tensorSpace) {
482:           while (!tup[dim]) {
483:             for (c = 0; c < Nc; ++c) {
484:               PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
485:               PetscMalloc1(dim, &qpoints);
486:               PetscCalloc1(Nc,  &qweights);
487:               PetscQuadratureSetOrder(sp->functional[f], 0);
488:               PetscQuadratureSetData(sp->functional[f], dim, Nc, 1, qpoints, qweights);
489:               for (d = 0; d < dim; ++d) {
490:                 qpoints[d] = v0[d];
491:                 for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
492:               }
493:               qweights[c] = 1.0;
494:               ++f;
495:             }
496:             LatticePointLexicographic_Internal(dim, orderEff, tup);
497:           }
498:         } else {
499:           while (!tup[dim]) {
500:             for (c = 0; c < Nc; ++c) {
501:               PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
502:               PetscMalloc1(dim, &qpoints);
503:               PetscCalloc1(Nc,  &qweights);
504:               PetscQuadratureSetOrder(sp->functional[f], 0);
505:               PetscQuadratureSetData(sp->functional[f], dim, Nc, 1, qpoints, qweights);
506:               for (d = 0; d < dim; ++d) {
507:                 qpoints[d] = v0[d];
508:                 for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
509:               }
510:               qweights[c] = 1.0;
511:               ++f;
512:             }
513:             TensorPointLexicographic_Internal(dim, orderEff, tup);
514:           }
515:         }
516:         lag->numDof[dim] = cdim;
517:       } else { /* transform functionals from subspaces */
518:         PetscInt q;

520:         for (q = 0; q < nFunc; q++, f++) {
521:           PetscQuadrature fn;
522:           PetscInt        fdim, Nc, c, nPoints, i;
523:           const PetscReal *points;
524:           const PetscReal *weights;
525:           PetscReal       *qpoints;
526:           PetscReal       *qweights;

528:           PetscDualSpaceGetFunctional(hsp, q, &fn);
529:           PetscQuadratureGetData(fn,&fdim,&Nc,&nPoints,&points,&weights);
530:           if (fdim != pointDim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expected height dual space dim %D, got %D",pointDim,fdim);
531:           PetscMalloc1(nPoints * dim, &qpoints);
532:           PetscCalloc1(nPoints * Nc,  &qweights);
533:           for (i = 0; i < nPoints; i++) {
534:             PetscInt  j, k;
535:             PetscReal *qp = &qpoints[i * dim];

537:             for (c = 0; c < Nc; ++c) qweights[i*Nc+c] = weights[i*Nc+c];
538:             for (j = 0; j < dim; ++j) qp[j] = v0[j];
539:             for (j = 0; j < dim; ++j) {
540:               for (k = 0; k < pointDim; k++) qp[j] += J[dim * j + k] * (points[pointDim * i + k] - hv0[k]);
541:             }
542:           }
543:           PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
544:           PetscQuadratureSetOrder(sp->functional[f],0);
545:           PetscQuadratureSetData(sp->functional[f],dim,Nc,nPoints,qpoints,qweights);
546:         }
547:       }
548:       PetscSectionSetDof(section,p,lag->numDof[pointDim]);
549:     }
550:     PetscFree5(tup,v0,hv0,J,invJ);
551:     PetscSectionSetUp(section);
552:     { /* reorder to closure order */
553:       PetscInt *key, count;
554:       PetscQuadrature *reorder = NULL;

556:       PetscCalloc1(f,&key);
557:       PetscMalloc1(f*sp->Nc,&reorder);

559:       for (p = pStratStart[depth], count = 0; p < pStratEnd[depth]; p++) {
560:         PetscInt *closure = NULL, closureSize, c;

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

566:           PetscSectionGetDof(section,point,&dof);
567:           PetscSectionGetOffset(section,point,&off);
568:           for (i = 0; i < dof; i++) {
569:             PetscInt fi = i + off;
570:             if (!key[fi]) {
571:               key[fi] = 1;
572:               reorder[count++] = sp->functional[fi];
573:             }
574:           }
575:         }
576:         DMPlexRestoreTransitiveClosure(dm,p,PETSC_TRUE,&closureSize,&closure);
577:       }
578:       PetscFree(sp->functional);
579:       sp->functional = reorder;
580:       PetscFree(key);
581:     }
582:     PetscSectionDestroy(&section);
583:   }
584:   if (pStratEnd[depth] == 1 && f != pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d not equal to dimension %d", f, pdimMax);
585:   PetscFree2(pStratStart, pStratEnd);
586:   if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
587:   return(0);
588: }

590: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
591: {
592:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
593:   PetscInt            i;
594:   PetscErrorCode      ierr;

597:   if (lag->symmetries) {
598:     PetscInt **selfSyms = lag->symmetries[0];

600:     if (selfSyms) {
601:       PetscInt i, **allocated = &selfSyms[-lag->selfSymOff];

603:       for (i = 0; i < lag->numSelfSym; i++) {
604:         PetscFree(allocated[i]);
605:       }
606:       PetscFree(allocated);
607:     }
608:     PetscFree(lag->symmetries);
609:   }
610:   for (i = 0; i < lag->height; i++) {
611:     PetscDualSpaceDestroy(&lag->subspaces[i]);
612:   }
613:   PetscFree(lag->subspaces);
614:   PetscFree(lag->numDof);
615:   PetscFree(lag);
616:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
617:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
618:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);
619:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);
620:   return(0);
621: }

623: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
624: {
625:   PetscInt       order, Nc;
626:   PetscBool      cont, tensor;

630:   PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
631:   PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
632:   PetscDualSpaceGetOrder(sp, &order);
633:   PetscDualSpaceSetOrder(*spNew, order);
634:   PetscDualSpaceGetNumComponents(sp, &Nc);
635:   PetscDualSpaceSetNumComponents(*spNew, Nc);
636:   PetscDualSpaceLagrangeGetContinuity(sp, &cont);
637:   PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
638:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
639:   PetscDualSpaceLagrangeSetTensor(*spNew, tensor);
640:   return(0);
641: }

643: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
644: {
645:   PetscBool      continuous, tensor, flg;

649:   PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
650:   PetscDualSpaceLagrangeGetTensor(sp, &tensor);
651:   PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
652:   PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
653:   if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
654:   PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetContinuity", tensor, &tensor, &flg);
655:   if (flg) {PetscDualSpaceLagrangeSetTensor(sp, tensor);}
656:   PetscOptionsTail();
657:   return(0);
658: }

660: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
661: {
662:   DM              K;
663:   const PetscInt *numDof;
664:   PetscInt        spatialDim, Nc, size = 0, d;
665:   PetscErrorCode  ierr;

668:   PetscDualSpaceGetDM(sp, &K);
669:   PetscDualSpaceGetNumDof(sp, &numDof);
670:   DMGetDimension(K, &spatialDim);
671:   DMPlexGetHeightStratum(K, 0, NULL, &Nc);
672:   if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
673:   for (d = 0; d <= spatialDim; ++d) {
674:     PetscInt pStart, pEnd;

676:     DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
677:     size += (pEnd-pStart)*numDof[d];
678:   }
679:   *dim = size;
680:   return(0);
681: }

683: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
684: {
685:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

688:   *numDof = lag->numDof;
689:   return(0);
690: }

692: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
693: {
694:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

699:   *continuous = lag->continuous;
700:   return(0);
701: }

703: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
704: {
705:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;

709:   lag->continuous = continuous;
710:   return(0);
711: }

713: /*@
714:   PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity

716:   Not Collective

718:   Input Parameter:
719: . sp         - the PetscDualSpace

721:   Output Parameter:
722: . continuous - flag for element continuity

724:   Level: intermediate

726: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
727: .seealso: PetscDualSpaceLagrangeSetContinuity()
728: @*/
729: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
730: {

736:   PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
737:   return(0);
738: }

740: /*@
741:   PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous

743:   Logically Collective on PetscDualSpace

745:   Input Parameters:
746: + sp         - the PetscDualSpace
747: - continuous - flag for element continuity

749:   Options Database:
750: . -petscdualspace_lagrange_continuity <bool>

752:   Level: intermediate

754: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
755: .seealso: PetscDualSpaceLagrangeGetContinuity()
756: @*/
757: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
758: {

764:   PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
765:   return(0);
766: }

768: PetscErrorCode PetscDualSpaceGetHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
769: {
770:   PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
771:   PetscErrorCode     ierr;

776:   if (height == 0) {
777:     *bdsp = sp;
778:   }
779:   else {
780:     DM dm;
781:     PetscInt dim;

783:     PetscDualSpaceGetDM(sp,&dm);
784:     DMGetDimension(dm,&dim);
785:     if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Asked for dual space at height %d for dimension %d reference element\n",height,dim);}
786:     if (height <= lag->height) {
787:       *bdsp = lag->subspaces[height-1];
788:     }
789:     else {
790:       *bdsp = NULL;
791:     }
792:   }
793:   return(0);
794: }

796: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
797: {
799:   sp->ops->setfromoptions    = PetscDualSpaceSetFromOptions_Lagrange;
800:   sp->ops->setup             = PetscDualSpaceSetUp_Lagrange;
801:   sp->ops->view              = PetscDualSpaceView_Lagrange;
802:   sp->ops->destroy           = PetscDualSpaceDestroy_Lagrange;
803:   sp->ops->duplicate         = PetscDualSpaceDuplicate_Lagrange;
804:   sp->ops->getdimension      = PetscDualSpaceGetDimension_Lagrange;
805:   sp->ops->getnumdof         = PetscDualSpaceGetNumDof_Lagrange;
806:   sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_Lagrange;
807:   sp->ops->getsymmetries     = PetscDualSpaceGetSymmetries_Lagrange;
808:   sp->ops->apply             = PetscDualSpaceApplyDefault;
809:   sp->ops->applyall          = PetscDualSpaceApplyAllDefault;
810:   sp->ops->createallpoints   = PetscDualSpaceCreateAllPointsDefault;
811:   return(0);
812: }

814: /*MC
815:   PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals

817:   Level: intermediate

819: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
820: M*/

822: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
823: {
824:   PetscDualSpace_Lag *lag;
825:   PetscErrorCode      ierr;

829:   PetscNewLog(sp,&lag);
830:   sp->data = lag;

832:   lag->numDof      = NULL;
833:   lag->simplexCell = PETSC_TRUE;
834:   lag->tensorSpace = PETSC_FALSE;
835:   lag->continuous  = PETSC_TRUE;

837:   PetscDualSpaceInitialize_Lagrange(sp);
838:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
839:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
840:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);
841:   PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);
842:   return(0);
843: }