Actual source code: dspacelagrange.c

petsc-3.11.3 2019-06-26
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\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "");
304:   return(0);
305: }

307: PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer)
308: {
309:   PetscBool      iascii;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

621: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
622: {
623:   PetscInt       order, Nc;
624:   PetscBool      cont, tensor;
625:   const char    *name;

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

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

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

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

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

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

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

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

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

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

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

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

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

717:   Not Collective

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

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

725:   Level: intermediate

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

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

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

744:   Logically Collective on PetscDualSpace

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

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

753:   Level: intermediate

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

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

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

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

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

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

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

818:   Level: intermediate

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

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

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

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

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