Actual source code: dmpleximpl.h

petsc-3.9.1 2018-04-29
Report Typos and Errors
  1: #if !defined(_PLEXIMPL_H)
  2: #define _PLEXIMPL_H

  4:  #include <petscmat.h>
  5:  #include <petscdmplex.h>
  6:  #include <petscbt.h>
  7:  #include <petscsf.h>
  8:  #include <petsc/private/dmimpl.h>
  9:  #include <petsc/private/isimpl.h>
 10:  #include <petsc/private/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_InterpolateSF, 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_INTERN 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_SIMPLEX_TO_HEX_2D,
 23:               REFINER_HEX_2D,
 24:               REFINER_HYBRID_HEX_2D,
 25:               REFINER_SIMPLEX_3D,
 26:               REFINER_HYBRID_SIMPLEX_3D,
 27:               REFINER_SIMPLEX_TO_HEX_3D,
 28:               REFINER_HEX_3D,
 29:               REFINER_HYBRID_HEX_3D} CellRefiner;

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

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

 46: typedef struct {
 47:   PetscInt dummy;
 48: } PetscPartitioner_Chaco;

 50: typedef struct {
 51:   PetscInt dummy;
 52: } PetscPartitioner_ParMetis;

 54: typedef struct {
 55:   PetscSection section;   /* Sizes for each partition */
 56:   IS           partition; /* Points in each partition */
 57:   PetscBool    random;    /* Flag for a random partition */
 58: } PetscPartitioner_Shell;

 60: typedef struct {
 61:   PetscInt dummy;
 62: } PetscPartitioner_Simple;

 64: typedef struct {
 65:   PetscInt dummy;
 66: } PetscPartitioner_Gather;

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

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

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

101: typedef struct {
102:   PetscInt             refct;

104:   /* Sieve */
105:   PetscSection         coneSection;       /* Layout of cones (inedges for DAG) */
106:   PetscInt             maxConeSize;       /* Cached for fast lookup */
107:   PetscInt            *cones;             /* Cone for each point */
108:   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 */
109:   PetscSection         supportSection;    /* Layout of cones (inedges for DAG) */
110:   PetscInt             maxSupportSize;    /* Cached for fast lookup */
111:   PetscInt            *supports;          /* Cone for each point */
112:   PetscBool            refinementUniform; /* Flag for uniform cell refinement */
113:   PetscReal            refinementLimit;   /* Maximum volume for refined cell */
114:   PetscErrorCode     (*refinementFunc)(const PetscReal [], PetscReal *); /* Function giving the maximum volume for refined cell */
115:   PetscInt             hybridPointMax[8]; /* Allow segregation of some points, each dimension has a divider (used in VTK output and refinement) */

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

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

122:   /* Generation */
123:   char                *tetgenOpts;
124:   char                *triangleOpts;
125:   PetscPartitioner     partitioner;
126:   PetscBool            remeshBd;

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

131:   /* Labels and numbering */
132:   PetscObjectState     depthState;        /* State of depth label, so that we can determine if a user changes it */
133:   IS                   globalVertexNumbers;
134:   IS                   globalCellNumbers;

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

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

151:   /* MATIS support */
152:   PetscSection         subdomainSection;

154:   /* Adjacency */
155:   PetscBool            useAnchors;        /* Replace constrained points with their anchors in adjacency lists */
156:   PetscErrorCode      (*useradjacency)(DM,PetscInt,PetscInt*,PetscInt[],void*); /* User callback for adjacency */
157:   void                *useradjacencyctx;  /* User context for callback */

159:   /* Projection */
160:   PetscInt             maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */

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

166:   /* Geometry */
167:   PetscReal            minradius;         /* Minimum distance from cell centroid to face */
168:   PetscBool            useHashLocation;   /* Use grid hashing for point location */
169:   PetscGridHash        lbox;              /* Local box for searching */

171:   /* Debugging */
172:   PetscBool            printSetValues;
173:   PetscInt             printFEM;
174:   PetscReal            printTol;
175: } DM_Plex;

177: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
178: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
179: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec,PetscViewer);
180: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
181: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
182: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec,PetscViewer);
183: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
184: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
185: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM,PetscViewer);
186: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject,PetscViewer);
187: #if defined(PETSC_HAVE_HDF5)
188: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
189: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
190: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
191: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
192: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
193: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
194: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
195: #endif

197: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *, DM);
198: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
199: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM []);
200: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
201: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM []);
202: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, DMLabel, DM *);
203: PETSC_INTERN PetscErrorCode DMAdaptMetric_Plex(DM, Vec, DMLabel, DM *);
204: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
205: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
206: PETSC_INTERN PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],PetscInt,const PetscInt[],PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
207: PETSC_INTERN PetscErrorCode DMProjectFieldLocal_Plex(DM,PetscReal,Vec,void (**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscInt,const PetscScalar[],PetscScalar[]),InsertMode,Vec);
208: PETSC_INTERN PetscErrorCode DMProjectFieldLabelLocal_Plex(DM,PetscReal,DMLabel,PetscInt,const PetscInt[],PetscInt,const PetscInt[],Vec,void (**)(PetscInt,PetscInt,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],const PetscInt[],const PetscInt[],const PetscScalar[],const PetscScalar[],const PetscScalar[],PetscReal,const PetscReal[],PetscInt,const PetscScalar[],PetscScalar[]),InsertMode,Vec);
209: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
210: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *);
211: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
212: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);

214: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
215: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
216: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
217: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
218: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
219: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
220: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
221: /* TODO Make these INTERN */
222: PETSC_EXTERN PetscErrorCode DMPlexView_ExodusII_Internal(DM, int, PetscInt);
223: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec, int, int);
224: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec, int, int);
225: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec, int, int);
226: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec, int, int);
227: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM,PetscInt,PetscInt,PetscInt*);
228: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
229: PETSC_INTERN PetscErrorCode DMPlexGetFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
230: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,PetscInt,PetscInt,const PetscInt[], PetscInt*,PetscInt*,const PetscInt*[]);
231: PETSC_INTERN PetscErrorCode DMPlexRestoreFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
232: PETSC_INTERN PetscErrorCode DMPlexRefineUniform_Internal(DM,CellRefiner,DM*);
233: PETSC_INTERN PetscErrorCode DMPlexGetCellRefiner_Internal(DM,CellRefiner*);
234: PETSC_INTERN PetscErrorCode CellRefinerGetAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
235: PETSC_INTERN PetscErrorCode CellRefinerRestoreAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
236: PETSC_INTERN PetscErrorCode CellRefinerInCellTest_Internal(CellRefiner, const PetscReal[], PetscBool *);
237: PETSC_INTERN PetscErrorCode DMPlexInvertCell_Internal(PetscInt, PetscInt, PetscInt[]);
238: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], const PetscScalar[], InsertMode);
239: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
240: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM,DM,const char *,DM*);
241: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM,DM,PetscSF,PetscInt *,Mat);
242: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM,DM,PetscSF,PetscInt *,Mat);
243: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM,PetscSection,PetscInt,PetscInt,const PetscInt[],const PetscInt ***,const PetscScalar[],PetscInt*,PetscInt*,PetscInt*[],PetscScalar*[],PetscInt[],PetscBool);
244: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt []);
245: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt [],PetscBool,PetscInt,PetscInt []);
246: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM,PetscInt,const PetscScalar [],PetscInt,PetscInt *);

248: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
249: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
250: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, DMLabel, DM *);
251: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, DMLabel, DM *);

253: /* invert dihedral symmetry: return a^-1,
254:  * using the representation described in
255:  * DMPlexGetConeOrientation() */
256: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
257: {
258:   return (a <= 0) ? a : (N - a);
259: }

261: /* invert dihedral symmetry: return b * a,
262:  * using the representation described in
263:  * DMPlexGetConeOrientation() */
264: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
265: {
266:   if (!N) return 0;
267:   return  (a >= 0) ?
268:          ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
269:          ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
270: }

272: /* swap dihedral symmetries: return b * a^-1,
273:  * using the representation described in
274:  * DMPlexGetConeOrientation() */
275: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
276: {
277:   return DihedralCompose(N,DihedralInvert(N,a),b);
278: }

280: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, PetscInt, PetscInt, PetscReal, Vec, Vec, PetscReal, Vec, void *);
281: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, PetscInt, PetscInt, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
282: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);

284: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
285: {
286:   const PetscReal invDet = 1.0/detJ;

288:   invJ[0] =  invDet*J[3];
289:   invJ[1] = -invDet*J[1];
290:   invJ[2] = -invDet*J[2];
291:   invJ[3] =  invDet*J[0];
292:   (void)PetscLogFlops(5.0);
293: }

295: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
296: {
297:   const PetscReal invDet = 1.0/detJ;

299:   invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
300:   invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
301:   invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
302:   invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
303:   invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
304:   invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
305:   invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
306:   invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
307:   invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
308:   (void)PetscLogFlops(37.0);
309: }

311: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
312: {
313:   *detJ = J[0]*J[3] - J[1]*J[2];
314:   (void)PetscLogFlops(3.0);
315: }

317: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
318: {
319:   *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
320:            J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
321:            J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
322:   (void)PetscLogFlops(12.0);
323: }

325: PETSC_STATIC_INLINE void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
326: {
327:   *detJ = PetscRealPart(J[0])*PetscRealPart(J[3]) - PetscRealPart(J[1])*PetscRealPart(J[2]);
328:   (void)PetscLogFlops(3.0);
329: }

331: PETSC_STATIC_INLINE void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
332: {
333:   *detJ = (PetscRealPart(J[0*3+0])*(PetscRealPart(J[1*3+1])*PetscRealPart(J[2*3+2]) - PetscRealPart(J[1*3+2])*PetscRealPart(J[2*3+1])) +
334:            PetscRealPart(J[0*3+1])*(PetscRealPart(J[1*3+2])*PetscRealPart(J[2*3+0]) - PetscRealPart(J[1*3+0])*PetscRealPart(J[2*3+2])) +
335:            PetscRealPart(J[0*3+2])*(PetscRealPart(J[1*3+0])*PetscRealPart(J[2*3+1]) - PetscRealPart(J[1*3+1])*PetscRealPart(J[2*3+0])));
336:   (void)PetscLogFlops(12.0);
337: }

339: 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];}

341: 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;}

343: 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;}

345: 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);}

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

370: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
371: {
373: #if defined(PETSC_USE_DEBUG)
374:   {
375:     PetscInt       dof;
377:     *start = *end = 0; /* Silence overzealous compiler warning */
378:     if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
379:     PetscSectionGetFieldOffset(dm->defaultSection, point, field, start);
380:     PetscSectionGetFieldDof(dm->defaultSection, point, field, &dof);
381:     *end = *start + dof;
382:   }
383: #else
384:   {
385:     const PetscSection s = dm->defaultSection->field[field];
386:     *start = s->atlasOff[point - s->pStart];
387:     *end   = *start + s->atlasDof[point - s->pStart];
388:   }
389: #endif
390:   return(0);
391: }

393: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
394: {
396: #if defined(PETSC_USE_DEBUG)
397:   {
399:     PetscInt       dof,cdof;
400:     *start = *end = 0; /* Silence overzealous compiler warning */
401:     if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
402:     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()");
403:     PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
404:     PetscSectionGetDof(dm->defaultGlobalSection, point, &dof);
405:     PetscSectionGetConstraintDof(dm->defaultGlobalSection, point, &cdof);
406:     *end = *start + dof - cdof + (dof < 0 ? 1 : 0);
407:   }
408: #else
409:   {
410:     const PetscSection s    = dm->defaultGlobalSection;
411:     const PetscInt     dof  = s->atlasDof[point - s->pStart];
412:     const PetscInt     cdof = s->bc ? s->bc->atlasDof[point - s->bc->pStart] : 0;
413:     *start = s->atlasOff[point - s->pStart];
414:     *end   = *start + dof - cdof + (dof < 0 ? 1 : 0);
415:   }
416: #endif
417:   return(0);
418: }

420: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
421: {
423: #if defined(PETSC_USE_DEBUG)
424:   {
425:     PetscInt       loff, lfoff, fdof, fcdof, ffcdof, f;
427:     *start = *end = 0; /* Silence overzealous compiler warning */
428:     if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
429:     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()");
430:     PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
431:     PetscSectionGetOffset(dm->defaultSection, point, &loff);
432:     PetscSectionGetFieldOffset(dm->defaultSection, point, field, &lfoff);
433:     PetscSectionGetFieldDof(dm->defaultSection, point, field, &fdof);
434:     PetscSectionGetFieldConstraintDof(dm->defaultSection, point, field, &fcdof);
435:     *start = *start < 0 ? *start - (lfoff-loff) : *start + lfoff-loff;
436:     for (f = 0; f < field; ++f) {
437:       PetscSectionGetFieldConstraintDof(dm->defaultSection, point, f, &ffcdof);
438:       *start = *start < 0 ? *start + ffcdof : *start - ffcdof;
439:     }
440:     *end   = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
441:   }
442: #else
443:   {
444:     const PetscSection s     = dm->defaultSection;
445:     const PetscSection fs    = dm->defaultSection->field[field];
446:     const PetscSection gs    = dm->defaultGlobalSection;
447:     const PetscInt     loff  = s->atlasOff[point - s->pStart];
448:     const PetscInt     goff  = gs->atlasOff[point - s->pStart];
449:     const PetscInt     lfoff = fs->atlasOff[point - s->pStart];
450:     const PetscInt     fdof  = fs->atlasDof[point - s->pStart];
451:     const PetscInt     fcdof = fs->bc ? fs->bc->atlasDof[point - fs->bc->pStart] : 0;
452:     PetscInt           ffcdof = 0, f;

454:     for (f = 0; f < field; ++f) {
455:       const PetscSection ffs = dm->defaultSection->field[f];
456:       ffcdof += ffs->bc ? ffs->bc->atlasDof[point - ffs->bc->pStart] : 0;
457:     }
458:     *start = goff + (goff < 0 ? loff-lfoff + ffcdof : lfoff-loff - ffcdof);
459:     *end   = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
460:   }
461: #endif
462:   return(0);
463: }

465: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM,PetscInt,PetscInt,PetscDualSpace *);
466: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,const PetscInt[],PetscInt[]);
467: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection,PetscInt,PetscInt,PetscInt[],PetscBool,const PetscInt***,PetscInt,PetscInt[]);

469: #endif /* _PLEXIMPL_H */