Actual source code: dmpleximpl.h

petsc-3.7.6 2017-04-24
Report Typos and Errors
  1: #if !defined(_PLEXIMPL_H)
  2: #define _PLEXIMPL_H

  4: #include <petscmat.h>       /*I      "petscmat.h"          I*/
  5: #include <petscdmplex.h> /*I      "petscdmplex.h"    I*/
  6: #include <petscbt.h>
  7: #include <petscsf.h>
  8: #include <petsc/private/dmimpl.h>
  9: #include <petsc/private/isimpl.h>     /* for inline access to atlasOff */
 10: #include <../src/sys/utils/hash.h>

 12: PETSC_EXTERN PetscLogEvent DMPLEX_Interpolate, PETSCPARTITIONER_Partition, DMPLEX_Distribute, DMPLEX_DistributeCones, DMPLEX_DistributeLabels, DMPLEX_DistributeSF, DMPLEX_DistributeOverlap, DMPLEX_DistributeField, DMPLEX_DistributeData, DMPLEX_Migrate, DMPLEX_GlobalToNaturalBegin, DMPLEX_GlobalToNaturalEnd, DMPLEX_NaturalToGlobalBegin, DMPLEX_NaturalToGlobalEnd, DMPLEX_Stratify, DMPLEX_Preallocate, DMPLEX_ResidualFEM, DMPLEX_JacobianFEM, DMPLEX_InterpolatorFEM, DMPLEX_InjectorFEM, DMPLEX_IntegralFEM, DMPLEX_CreateGmsh;

 14: PETSC_EXTERN PetscBool      PetscPartitionerRegisterAllCalled;
 15: PETSC_EXTERN PetscErrorCode PetscPartitionerRegisterAll(void);
 16: PETSC_EXTERN PetscErrorCode PetscPartitionerSetTypeFromOptions_Internal(PetscPartitioner);

 18: typedef enum {REFINER_NOOP = 0,
 19:               REFINER_SIMPLEX_1D,
 20:               REFINER_SIMPLEX_2D,
 21:               REFINER_HYBRID_SIMPLEX_2D,
 22:               REFINER_HEX_2D,
 23:               REFINER_HYBRID_HEX_2D,
 24:               REFINER_SIMPLEX_3D,
 25:               REFINER_HYBRID_SIMPLEX_3D,
 26:               REFINER_HEX_3D,
 27:               REFINER_HYBRID_HEX_3D} CellRefiner;

 29: typedef struct _PetscPartitionerOps *PetscPartitionerOps;
 30: struct _PetscPartitionerOps {
 31:   PetscErrorCode (*setfromoptions)(PetscPartitioner);
 32:   PetscErrorCode (*setup)(PetscPartitioner);
 33:   PetscErrorCode (*view)(PetscPartitioner,PetscViewer);
 34:   PetscErrorCode (*destroy)(PetscPartitioner);
 35:   PetscErrorCode (*partition)(PetscPartitioner, DM, PetscInt, PetscInt, PetscInt[], PetscInt[], PetscSection, IS *);
 36: };

 38: struct _p_PetscPartitioner {
 39:   PETSCHEADER(struct _PetscPartitionerOps);
 40:   void           *data;             /* Implementation object */
 41:   PetscInt        height;           /* Height of points to partition into non-overlapping subsets */
 42: };

 44: typedef struct {
 45:   PetscInt dummy;
 46: } PetscPartitioner_Chaco;

 48: typedef struct {
 49:   PetscInt dummy;
 50: } PetscPartitioner_ParMetis;

 52: typedef struct {
 53:   PetscSection section;   /* Sizes for each partition */
 54:   IS           partition; /* Points in each partition */
 55: } PetscPartitioner_Shell;

 57: typedef struct {
 58:   PetscInt dummy;
 59: } PetscPartitioner_Simple;

 61: typedef struct {
 62:   PetscInt dummy;
 63: } PetscPartitioner_Gather;

 65: /* Utility struct to store the contents of a Gmsh file in memory */
 66: typedef struct {
 67:   PetscInt dim;      /* Entity dimension */
 68:   PetscInt id;       /* Element number */
 69:   PetscInt numNodes; /* Size of node array */
 70:   int nodes[8];      /* Node array */
 71:   PetscInt numTags;  /* Size of tag array */
 72:   int tags[4];       /* Tag array */
 73: } GmshElement;

 75: /* Utility struct to store the contents of a Fluent file in memory */
 76: typedef struct {
 77:   int   index;    /* Type of section */
 78:   int   zoneID;
 79:   int   first;
 80:   int   last;
 81:   int   type;
 82:   int   nd;       /* Either ND or element-type */
 83:   void *data;
 84: } FluentSection;

 86: struct _PetscGridHash {
 87:   PetscInt     dim;
 88:   PetscReal    lower[3];    /* The lower-left corner */
 89:   PetscReal    upper[3];    /* The upper-right corner */
 90:   PetscReal    extent[3];   /* The box size */
 91:   PetscReal    h[3];        /* The subbox size */
 92:   PetscInt     n[3];        /* The number of subboxes */
 93:   PetscSection cellSection; /* Offsets for cells in each subbox*/
 94:   IS           cells;       /* List of cells in each subbox */
 95:   DMLabel      cellsSparse; /* Sparse storage for cell map */
 96: };

 98: typedef struct {
 99:   PetscInt             refct;

101:   /* Sieve */
102:   PetscSection         coneSection;       /* Layout of cones (inedges for DAG) */
103:   PetscInt             maxConeSize;       /* Cached for fast lookup */
104:   PetscInt            *cones;             /* Cone for each point */
105:   PetscInt            *coneOrientations;  /* Orientation of each cone point, means cone traveral should start on point 'o', and if negative start on -(o+1) and go in reverse */
106:   PetscSection         supportSection;    /* Layout of cones (inedges for DAG) */
107:   PetscInt             maxSupportSize;    /* Cached for fast lookup */
108:   PetscInt            *supports;          /* Cone for each point */
109:   PetscBool            refinementUniform; /* Flag for uniform cell refinement */
110:   PetscReal            refinementLimit;   /* Maximum volume for refined cell */
111:   PetscErrorCode     (*refinementFunc)(const PetscReal [], PetscReal *); /* Function giving the maximum volume for refined cell */
112:   PetscInt             hybridPointMax[8]; /* Allow segregation of some points, each dimension has a divider (used in VTK output and refinement) */

114:   PetscInt            *facesTmp;          /* Work space for faces operation */

116:   /* Hierarchy */
117:   PetscBool            regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */

119:   /* Generation */
120:   char                *tetgenOpts;
121:   char                *triangleOpts;
122:   PetscPartitioner     partitioner;

124:   /* Submesh */
125:   DMLabel              subpointMap;       /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */

127:   /* Labels and numbering */
128:   PetscObjectState     depthState;        /* State of depth label, so that we can determine if a user changes it */
129:   IS                   globalVertexNumbers;
130:   IS                   globalCellNumbers;

132:   /* Constraints */
133:   PetscSection         anchorSection;      /* maps constrained points to anchor points */
134:   IS                   anchorIS;           /* anchors indexed by the above section */
135:   PetscErrorCode     (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
136:   PetscErrorCode     (*computeanchormatrix)(DM,PetscSection,PetscSection,Mat);

138:   /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
139:   PetscSection         parentSection;     /* dof == 1 if point has parent */
140:   PetscInt            *parents;           /* point to parent */
141:   PetscInt            *childIDs;          /* point to child ID */
142:   PetscSection         childSection;      /* inverse of parent section */
143:   PetscInt            *children;          /* point to children */
144:   DM                   referenceTree;     /* reference tree to which child ID's refer */
145:   PetscErrorCode      (*getchildsymmetry)(DM,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*);

147:   /* Adjacency */
148:   PetscBool            useCone;           /* Use cone() first when defining adjacency */
149:   PetscBool            useClosure;        /* Use the transitive closure when defining adjacency */
150:   PetscBool            useAnchors;        /* Replace constrained points with their anchors in adjacency lists */

152:   /* Projection */
153:   PetscInt             maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */

155:   /* Output */
156:   PetscInt             vtkCellHeight;            /* The height of cells for output, default is 0 */
157:   PetscReal            scale[NUM_PETSC_UNITS];   /* The scale for each SI unit */

159:   /* Geometry */
160:   PetscReal            minradius;         /* Minimum distance from cell centroid to face */
161:   PetscBool            useHashLocation;   /* Use grid hashing for point location */
162:   PetscGridHash        lbox;              /* Local box for searching */

164:   /* Debugging */
165:   PetscBool            printSetValues;
166:   PetscInt             printFEM;
167:   PetscReal            printTol;
168: } DM_Plex;

170: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
171: PETSC_EXTERN PetscErrorCode DMPlexVTKGetCellType(DM,PetscInt,PetscInt,PetscInt*);
172: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
173: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec,PetscViewer);
174: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
175: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
176: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec,PetscViewer);
177: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
178: PETSC_EXTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
179: #if defined(PETSC_HAVE_HDF5)
180: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
181: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
182: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
183: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
184: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
185: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
186: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
187: #endif

189: PETSC_EXTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
190: PETSC_EXTERN PetscErrorCode DMPlexGetFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
191: PETSC_EXTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,PetscInt,PetscInt,const PetscInt[], PetscInt*,PetscInt*,const PetscInt*[]);
192: PETSC_EXTERN PetscErrorCode DMPlexRestoreFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
193: PETSC_EXTERN PetscErrorCode DMPlexRefineUniform_Internal(DM,CellRefiner,DM*);
194: PETSC_EXTERN PetscErrorCode DMPlexGetCellRefiner_Internal(DM,CellRefiner*);
195: PETSC_EXTERN PetscErrorCode CellRefinerGetAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
196: PETSC_EXTERN PetscErrorCode CellRefinerRestoreAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
197: PETSC_EXTERN PetscErrorCode CellRefinerInCellTest_Internal(CellRefiner, const PetscReal[], PetscBool *);
198: PETSC_EXTERN PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer, PetscInt, PetscBool, PetscBool, GmshElement **);
199: PETSC_EXTERN PetscErrorCode DMPlexInvertCell_Internal(PetscInt, PetscInt, PetscInt[]);
200: PETSC_EXTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, const PetscScalar[], InsertMode);
201: PETSC_EXTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
202: PETSC_EXTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems*,DM);
203: PETSC_EXTERN PetscErrorCode DMPlexLabelComplete_Internal(DM,DMLabel,PetscBool);
204: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM,DM,const char *,DM*);
205: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM,DM,PetscSF,PetscInt *,Mat);
206: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM,DM,PetscSF,PetscInt *,Mat);
207: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM,PetscSection,PetscInt,PetscInt,const PetscInt*,const PetscScalar*,PetscInt*,PetscInt*,PetscInt**,PetscScalar**,PetscInt*,PetscBool);
208: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt []);
209: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt [],PetscBool,PetscInt,PetscInt []);
210: PETSC_EXTERN PetscErrorCode DMPlexLocatePoint_Internal(DM,PetscInt,const PetscScalar [],PetscInt,PetscInt *);

212: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
213: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);

217: /* invert dihedral symmetry: return a^-1,
218:  * using the representation described in
219:  * DMPlexGetConeOrientation() */
220: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
221: {
222:   return (a <= 0) ? a : (N - a);
223: }

227: /* invert dihedral symmetry: return b * a,
228:  * using the representation described in
229:  * DMPlexGetConeOrientation() */
230: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
231: {
232:   if (!N) return 0;
233:   return  (a >= 0) ?
234:          ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
235:          ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
236: }

240: /* swap dihedral symmetries: return b * a^-1,
241:  * using the representation described in
242:  * DMPlexGetConeOrientation() */
243: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
244: {
245:   return DihedralCompose(N,DihedralInvert(N,a),b);
246: }

248: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, PetscInt, PetscInt, PetscReal, Vec, Vec, Vec, void *);
249: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, PetscInt, PetscInt, PetscReal, PetscReal, Vec, Vec, Mat, Mat,void *);

253: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
254: {
255:   const PetscReal invDet = 1.0/detJ;

257:   invJ[0] =  invDet*J[3];
258:   invJ[1] = -invDet*J[1];
259:   invJ[2] = -invDet*J[2];
260:   invJ[3] =  invDet*J[0];
261:   (void)PetscLogFlops(5.0);
262: }

266: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
267: {
268:   const PetscReal invDet = 1.0/detJ;

270:   invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
271:   invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
272:   invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
273:   invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
274:   invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
275:   invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
276:   invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
277:   invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
278:   invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
279:   (void)PetscLogFlops(37.0);
280: }

284: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, PetscReal J[])
285: {
286:   *detJ = J[0]*J[3] - J[1]*J[2];
287:   (void)PetscLogFlops(3.0);
288: }

292: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, PetscReal J[])
293: {
294:   *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
295:            J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
296:            J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
297:   (void)PetscLogFlops(12.0);
298: }

302: PETSC_STATIC_INLINE void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w) {PetscInt d; for (d = 0; d < dim; ++d) w[d] = a*x[d] + y[d];}

306: PETSC_STATIC_INLINE PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d])*y[d]; return sum;}

310: PETSC_STATIC_INLINE PetscReal DMPlex_DotRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*y[d]; return sum;}

314: PETSC_STATIC_INLINE PetscReal DMPlex_NormD_Internal(PetscInt dim, const PetscReal *x) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*x[d]; return PetscSqrtReal(sum);}


319: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
320: {
322: #if defined(PETSC_USE_DEBUG)
323:   {
324:     PetscInt       dof;
326:     *start = *end = 0; /* Silence overzealous compiler warning */
327:     if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
328:     PetscSectionGetOffset(dm->defaultSection, point, start);
329:     PetscSectionGetDof(dm->defaultSection, point, &dof);
330:     *end = *start + dof;
331:   }
332: #else
333:   {
334:     const PetscSection s = dm->defaultSection;
335:     *start = s->atlasOff[point - s->pStart];
336:     *end   = *start + s->atlasDof[point - s->pStart];
337:   }
338: #endif
339:   return(0);
340: }

344: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
345: {
347: #if defined(PETSC_USE_DEBUG)
348:   {
349:     PetscInt       dof;
351:     *start = *end = 0; /* Silence overzealous compiler warning */
352:     if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
353:     PetscSectionGetFieldOffset(dm->defaultSection, point, field, start);
354:     PetscSectionGetFieldDof(dm->defaultSection, point, field, &dof);
355:     *end = *start + dof;
356:   }
357: #else
358:   {
359:     const PetscSection s = dm->defaultSection->field[field];
360:     *start = s->atlasOff[point - s->pStart];
361:     *end   = *start + s->atlasDof[point - s->pStart];
362:   }
363: #endif
364:   return(0);
365: }

369: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
370: {
372: #if defined(PETSC_USE_DEBUG)
373:   {
375:     PetscInt       dof,cdof;
376:     *start = *end = 0; /* Silence overzealous compiler warning */
377:     if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
378:     if (!dm->defaultGlobalSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default global Section. It will be created automatically by DMGetDefaultGlobalSection()");
379:     PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
380:     PetscSectionGetDof(dm->defaultGlobalSection, point, &dof);
381:     PetscSectionGetConstraintDof(dm->defaultGlobalSection, point, &cdof);
382:     *end = *start + dof - cdof + (dof < 0 ? 1 : 0);
383:   }
384: #else
385:   {
386:     const PetscSection s    = dm->defaultGlobalSection;
387:     const PetscInt     dof  = s->atlasDof[point - s->pStart];
388:     const PetscInt     cdof = s->bc ? s->bc->atlasDof[point - s->bc->pStart] : 0;
389:     *start = s->atlasOff[point - s->pStart];
390:     *end   = *start + dof - cdof + (dof < 0 ? 1 : 0);
391:   }
392: #endif
393:   return(0);
394: }

398: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
399: {
401: #if defined(PETSC_USE_DEBUG)
402:   {
403:     PetscInt       loff, lfoff, fdof, fcdof, ffcdof, f;
405:     *start = *end = 0; /* Silence overzealous compiler warning */
406:     if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
407:     if (!dm->defaultGlobalSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default global Section. It will be crated automatically by DMGetDefaultGlobalSection()");
408:     PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
409:     PetscSectionGetOffset(dm->defaultSection, point, &loff);
410:     PetscSectionGetFieldOffset(dm->defaultSection, point, field, &lfoff);
411:     PetscSectionGetFieldDof(dm->defaultSection, point, field, &fdof);
412:     PetscSectionGetFieldConstraintDof(dm->defaultSection, point, field, &fcdof);
413:     *start = *start < 0 ? *start - (lfoff-loff) : *start + lfoff-loff;
414:     for (f = 0; f < field; ++f) {
415:       PetscSectionGetFieldConstraintDof(dm->defaultSection, point, f, &ffcdof);
416:       *start = *start < 0 ? *start + ffcdof : *start - ffcdof;
417:     }
418:     *end   = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
419:   }
420: #else
421:   {
422:     const PetscSection s     = dm->defaultSection;
423:     const PetscSection fs    = dm->defaultSection->field[field];
424:     const PetscSection gs    = dm->defaultGlobalSection;
425:     const PetscInt     loff  = s->atlasOff[point - s->pStart];
426:     const PetscInt     goff  = gs->atlasOff[point - s->pStart];
427:     const PetscInt     lfoff = fs->atlasOff[point - s->pStart];
428:     const PetscInt     fdof  = fs->atlasDof[point - s->pStart];
429:     const PetscInt     fcdof = fs->bc ? fs->bc->atlasDof[point - fs->bc->pStart] : 0;
430:     PetscInt           ffcdof = 0, f;

432:     for (f = 0; f < field; ++f) {
433:       const PetscSection ffs = dm->defaultSection->field[f];
434:       ffcdof += ffs->bc ? ffs->bc->atlasDof[point - ffs->bc->pStart] : 0;
435:     }
436:     *start = goff + (goff < 0 ? loff-lfoff + ffcdof : lfoff-loff - ffcdof);
437:     *end   = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
438:   }
439: #endif
440:   return(0);
441: }

443: #endif /* _PLEXIMPL_H */