Actual source code: dmpleximpl.h
petsc-3.7.7 2017-09-25
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 */