Actual source code: dmpleximpl.h
petsc-3.9.4 2018-09-11
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;
13: PETSC_EXTERN PetscLogEvent PETSCPARTITIONER_Partition;
14: PETSC_EXTERN PetscLogEvent DMPLEX_Distribute;
15: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeCones;
16: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeLabels;
17: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeSF;
18: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeOverlap;
19: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeField;
20: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeData;
21: PETSC_EXTERN PetscLogEvent DMPLEX_Migrate;
22: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolateSF;
23: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalBegin;
24: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalEnd;
25: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalBegin;
26: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalEnd;
27: PETSC_EXTERN PetscLogEvent DMPLEX_Stratify;
28: PETSC_EXTERN PetscLogEvent DMPLEX_Preallocate;
29: PETSC_EXTERN PetscLogEvent DMPLEX_ResidualFEM;
30: PETSC_EXTERN PetscLogEvent DMPLEX_JacobianFEM;
31: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolatorFEM;
32: PETSC_EXTERN PetscLogEvent DMPLEX_InjectorFEM;
33: PETSC_EXTERN PetscLogEvent DMPLEX_IntegralFEM;
34: PETSC_EXTERN PetscLogEvent DMPLEX_CreateGmsh;
36: PETSC_EXTERN PetscBool PetscPartitionerRegisterAllCalled;
37: PETSC_EXTERN PetscErrorCode PetscPartitionerRegisterAll(void);
39: typedef enum {REFINER_NOOP = 0,
40: REFINER_SIMPLEX_1D,
41: REFINER_SIMPLEX_2D,
42: REFINER_HYBRID_SIMPLEX_2D,
43: REFINER_SIMPLEX_TO_HEX_2D,
44: REFINER_HEX_2D,
45: REFINER_HYBRID_HEX_2D,
46: REFINER_SIMPLEX_3D,
47: REFINER_HYBRID_SIMPLEX_3D,
48: REFINER_SIMPLEX_TO_HEX_3D,
49: REFINER_HEX_3D,
50: REFINER_HYBRID_HEX_3D} CellRefiner;
52: typedef struct _PetscPartitionerOps *PetscPartitionerOps;
53: struct _PetscPartitionerOps {
54: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscPartitioner);
55: PetscErrorCode (*setup)(PetscPartitioner);
56: PetscErrorCode (*view)(PetscPartitioner,PetscViewer);
57: PetscErrorCode (*destroy)(PetscPartitioner);
58: PetscErrorCode (*partition)(PetscPartitioner, DM, PetscInt, PetscInt, PetscInt[], PetscInt[], PetscSection, IS *);
59: };
61: struct _p_PetscPartitioner {
62: PETSCHEADER(struct _PetscPartitionerOps);
63: void *data; /* Implementation object */
64: PetscInt height; /* Height of points to partition into non-overlapping subsets */
65: };
67: typedef struct {
68: PetscInt dummy;
69: } PetscPartitioner_Chaco;
71: typedef struct {
72: PetscInt ptype;
73: } PetscPartitioner_ParMetis;
75: typedef struct {
76: PetscSection section; /* Sizes for each partition */
77: IS partition; /* Points in each partition */
78: PetscBool random; /* Flag for a random partition */
79: } PetscPartitioner_Shell;
81: typedef struct {
82: PetscInt dummy;
83: } PetscPartitioner_Simple;
85: typedef struct {
86: PetscInt dummy;
87: } PetscPartitioner_Gather;
89: /* Utility struct to store the contents of a Gmsh file in memory */
90: typedef struct {
91: PetscInt dim; /* Entity dimension */
92: PetscInt id; /* Element number */
93: PetscInt numNodes; /* Size of node array */
94: int nodes[8]; /* Node array */
95: PetscInt numTags; /* Size of tag array */
96: int tags[4]; /* Tag array */
97: } GmshElement;
99: /* Utility struct to store the contents of a Fluent file in memory */
100: typedef struct {
101: int index; /* Type of section */
102: unsigned int zoneID;
103: unsigned int first;
104: unsigned int last;
105: int type;
106: int nd; /* Either ND or element-type */
107: void *data;
108: } FluentSection;
110: struct _PetscGridHash {
111: PetscInt dim;
112: PetscReal lower[3]; /* The lower-left corner */
113: PetscReal upper[3]; /* The upper-right corner */
114: PetscReal extent[3]; /* The box size */
115: PetscReal h[3]; /* The subbox size */
116: PetscInt n[3]; /* The number of subboxes */
117: PetscSection cellSection; /* Offsets for cells in each subbox*/
118: IS cells; /* List of cells in each subbox */
119: DMLabel cellsSparse; /* Sparse storage for cell map */
120: };
122: typedef struct {
123: PetscInt refct;
125: /* Sieve */
126: PetscSection coneSection; /* Layout of cones (inedges for DAG) */
127: PetscInt maxConeSize; /* Cached for fast lookup */
128: PetscInt *cones; /* Cone for each point */
129: 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 */
130: PetscSection supportSection; /* Layout of cones (inedges for DAG) */
131: PetscInt maxSupportSize; /* Cached for fast lookup */
132: PetscInt *supports; /* Cone for each point */
133: PetscBool refinementUniform; /* Flag for uniform cell refinement */
134: PetscReal refinementLimit; /* Maximum volume for refined cell */
135: PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *); /* Function giving the maximum volume for refined cell */
136: PetscInt hybridPointMax[8]; /* Allow segregation of some points, each dimension has a divider (used in VTK output and refinement) */
138: PetscInt *facesTmp; /* Work space for faces operation */
140: /* Hierarchy */
141: PetscBool regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */
143: /* Generation */
144: char *tetgenOpts;
145: char *triangleOpts;
146: PetscPartitioner partitioner;
147: PetscBool remeshBd;
149: /* Submesh */
150: DMLabel subpointMap; /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
152: /* Labels and numbering */
153: PetscObjectState depthState; /* State of depth label, so that we can determine if a user changes it */
154: IS globalVertexNumbers;
155: IS globalCellNumbers;
157: /* Constraints */
158: PetscSection anchorSection; /* maps constrained points to anchor points */
159: IS anchorIS; /* anchors indexed by the above section */
160: PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
161: PetscErrorCode (*computeanchormatrix)(DM,PetscSection,PetscSection,Mat);
163: /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
164: PetscSection parentSection; /* dof == 1 if point has parent */
165: PetscInt *parents; /* point to parent */
166: PetscInt *childIDs; /* point to child ID */
167: PetscSection childSection; /* inverse of parent section */
168: PetscInt *children; /* point to children */
169: DM referenceTree; /* reference tree to which child ID's refer */
170: PetscErrorCode (*getchildsymmetry)(DM,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*);
172: /* MATIS support */
173: PetscSection subdomainSection;
175: /* Adjacency */
176: PetscBool useAnchors; /* Replace constrained points with their anchors in adjacency lists */
177: PetscErrorCode (*useradjacency)(DM,PetscInt,PetscInt*,PetscInt[],void*); /* User callback for adjacency */
178: void *useradjacencyctx; /* User context for callback */
180: /* Projection */
181: PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
183: /* Output */
184: PetscInt vtkCellHeight; /* The height of cells for output, default is 0 */
185: PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */
187: /* Geometry */
188: PetscReal minradius; /* Minimum distance from cell centroid to face */
189: PetscBool useHashLocation; /* Use grid hashing for point location */
190: PetscGridHash lbox; /* Local box for searching */
192: /* Debugging */
193: PetscBool printSetValues;
194: PetscInt printFEM;
195: PetscReal printTol;
196: } DM_Plex;
198: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
199: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
200: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec,PetscViewer);
201: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
202: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
203: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec,PetscViewer);
204: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
205: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
206: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM,PetscViewer);
207: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject,PetscViewer);
208: #if defined(PETSC_HAVE_HDF5)
209: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
210: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
211: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
212: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
213: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
214: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
215: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
216: #endif
218: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *, DM);
219: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
220: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM []);
221: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
222: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM []);
223: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, DMLabel, DM *);
224: PETSC_INTERN PetscErrorCode DMAdaptMetric_Plex(DM, Vec, DMLabel, DM *);
225: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
226: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
227: 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);
228: 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);
229: 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);
230: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
231: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *);
232: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
233: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);
235: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
236: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
237: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
238: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
239: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
240: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
241: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
242: /* TODO Make these INTERN */
243: PETSC_EXTERN PetscErrorCode DMPlexView_ExodusII_Internal(DM, int, PetscInt);
244: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec, int, int);
245: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec, int, int);
246: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec, int, int);
247: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec, int, int);
248: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM,PetscInt,PetscInt,PetscInt*);
249: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
250: PETSC_INTERN PetscErrorCode DMPlexGetFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
251: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,PetscInt,PetscInt,const PetscInt[], PetscInt*,PetscInt*,const PetscInt*[]);
252: PETSC_INTERN PetscErrorCode DMPlexRestoreFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
253: PETSC_INTERN PetscErrorCode DMPlexRefineUniform_Internal(DM,CellRefiner,DM*);
254: PETSC_INTERN PetscErrorCode DMPlexGetCellRefiner_Internal(DM,CellRefiner*);
255: PETSC_INTERN PetscErrorCode CellRefinerGetAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
256: PETSC_INTERN PetscErrorCode CellRefinerRestoreAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
257: PETSC_INTERN PetscErrorCode CellRefinerInCellTest_Internal(CellRefiner, const PetscReal[], PetscBool *);
258: PETSC_INTERN PetscErrorCode DMPlexInvertCell_Internal(PetscInt, PetscInt, PetscInt[]);
259: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], const PetscScalar[], InsertMode);
260: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
261: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM,DM,const char *,DM*);
262: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM,DM,PetscSF,PetscInt *,Mat);
263: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM,DM,PetscSF,PetscInt *,Mat);
264: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM,PetscSection,PetscInt,PetscInt,const PetscInt[],const PetscInt ***,const PetscScalar[],PetscInt*,PetscInt*,PetscInt*[],PetscScalar*[],PetscInt[],PetscBool);
265: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt []);
266: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt [],PetscBool,PetscInt,PetscInt []);
267: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM,PetscInt,const PetscScalar [],PetscInt,PetscInt *);
269: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
270: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
271: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, DMLabel, DM *);
272: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, DMLabel, DM *);
274: /* invert dihedral symmetry: return a^-1,
275: * using the representation described in
276: * DMPlexGetConeOrientation() */
277: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
278: {
279: return (a <= 0) ? a : (N - a);
280: }
282: /* invert dihedral symmetry: return b * a,
283: * using the representation described in
284: * DMPlexGetConeOrientation() */
285: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
286: {
287: if (!N) return 0;
288: return (a >= 0) ?
289: ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
290: ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
291: }
293: /* swap dihedral symmetries: return b * a^-1,
294: * using the representation described in
295: * DMPlexGetConeOrientation() */
296: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
297: {
298: return DihedralCompose(N,DihedralInvert(N,a),b);
299: }
301: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, PetscInt, PetscInt, PetscReal, Vec, Vec, PetscReal, Vec, void *);
302: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, PetscInt, PetscInt, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
303: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);
305: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
306: {
307: const PetscReal invDet = 1.0/detJ;
309: invJ[0] = invDet*J[3];
310: invJ[1] = -invDet*J[1];
311: invJ[2] = -invDet*J[2];
312: invJ[3] = invDet*J[0];
313: (void)PetscLogFlops(5.0);
314: }
316: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
317: {
318: const PetscReal invDet = 1.0/detJ;
320: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
321: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
322: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
323: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
324: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
325: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
326: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
327: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
328: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
329: (void)PetscLogFlops(37.0);
330: }
332: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
333: {
334: *detJ = J[0]*J[3] - J[1]*J[2];
335: (void)PetscLogFlops(3.0);
336: }
338: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
339: {
340: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
341: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
342: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
343: (void)PetscLogFlops(12.0);
344: }
346: PETSC_STATIC_INLINE void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
347: {
348: *detJ = PetscRealPart(J[0])*PetscRealPart(J[3]) - PetscRealPart(J[1])*PetscRealPart(J[2]);
349: (void)PetscLogFlops(3.0);
350: }
352: PETSC_STATIC_INLINE void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
353: {
354: *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])) +
355: 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])) +
356: 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])));
357: (void)PetscLogFlops(12.0);
358: }
360: 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];}
362: 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;}
364: 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;}
366: 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);}
368: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
369: {
371: #if defined(PETSC_USE_DEBUG)
372: {
373: PetscInt dof;
375: *start = *end = 0; /* Silence overzealous compiler warning */
376: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
377: PetscSectionGetOffset(dm->defaultSection, point, start);
378: PetscSectionGetDof(dm->defaultSection, point, &dof);
379: *end = *start + dof;
380: }
381: #else
382: {
383: const PetscSection s = dm->defaultSection;
384: *start = s->atlasOff[point - s->pStart];
385: *end = *start + s->atlasDof[point - s->pStart];
386: }
387: #endif
388: return(0);
389: }
391: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
392: {
394: #if defined(PETSC_USE_DEBUG)
395: {
396: PetscInt dof;
398: *start = *end = 0; /* Silence overzealous compiler warning */
399: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
400: PetscSectionGetFieldOffset(dm->defaultSection, point, field, start);
401: PetscSectionGetFieldDof(dm->defaultSection, point, field, &dof);
402: *end = *start + dof;
403: }
404: #else
405: {
406: const PetscSection s = dm->defaultSection->field[field];
407: *start = s->atlasOff[point - s->pStart];
408: *end = *start + s->atlasDof[point - s->pStart];
409: }
410: #endif
411: return(0);
412: }
414: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
415: {
417: #if defined(PETSC_USE_DEBUG)
418: {
420: PetscInt dof,cdof;
421: *start = *end = 0; /* Silence overzealous compiler warning */
422: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
423: 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()");
424: PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
425: PetscSectionGetDof(dm->defaultGlobalSection, point, &dof);
426: PetscSectionGetConstraintDof(dm->defaultGlobalSection, point, &cdof);
427: *end = *start + dof - cdof + (dof < 0 ? 1 : 0);
428: }
429: #else
430: {
431: const PetscSection s = dm->defaultGlobalSection;
432: const PetscInt dof = s->atlasDof[point - s->pStart];
433: const PetscInt cdof = s->bc ? s->bc->atlasDof[point - s->bc->pStart] : 0;
434: *start = s->atlasOff[point - s->pStart];
435: *end = *start + dof - cdof + (dof < 0 ? 1 : 0);
436: }
437: #endif
438: return(0);
439: }
441: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
442: {
444: #if defined(PETSC_USE_DEBUG)
445: {
446: PetscInt loff, lfoff, fdof, fcdof, ffcdof, f;
448: *start = *end = 0; /* Silence overzealous compiler warning */
449: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
450: 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()");
451: PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
452: PetscSectionGetOffset(dm->defaultSection, point, &loff);
453: PetscSectionGetFieldOffset(dm->defaultSection, point, field, &lfoff);
454: PetscSectionGetFieldDof(dm->defaultSection, point, field, &fdof);
455: PetscSectionGetFieldConstraintDof(dm->defaultSection, point, field, &fcdof);
456: *start = *start < 0 ? *start - (lfoff-loff) : *start + lfoff-loff;
457: for (f = 0; f < field; ++f) {
458: PetscSectionGetFieldConstraintDof(dm->defaultSection, point, f, &ffcdof);
459: *start = *start < 0 ? *start + ffcdof : *start - ffcdof;
460: }
461: *end = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
462: }
463: #else
464: {
465: const PetscSection s = dm->defaultSection;
466: const PetscSection fs = dm->defaultSection->field[field];
467: const PetscSection gs = dm->defaultGlobalSection;
468: const PetscInt loff = s->atlasOff[point - s->pStart];
469: const PetscInt goff = gs->atlasOff[point - s->pStart];
470: const PetscInt lfoff = fs->atlasOff[point - s->pStart];
471: const PetscInt fdof = fs->atlasDof[point - s->pStart];
472: const PetscInt fcdof = fs->bc ? fs->bc->atlasDof[point - fs->bc->pStart] : 0;
473: PetscInt ffcdof = 0, f;
475: for (f = 0; f < field; ++f) {
476: const PetscSection ffs = dm->defaultSection->field[f];
477: ffcdof += ffs->bc ? ffs->bc->atlasDof[point - ffs->bc->pStart] : 0;
478: }
479: *start = goff + (goff < 0 ? loff-lfoff + ffcdof : lfoff-loff - ffcdof);
480: *end = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
481: }
482: #endif
483: return(0);
484: }
486: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM,PetscInt,PetscInt,PetscDualSpace *);
487: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,const PetscInt[],PetscInt[]);
488: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection,PetscInt,PetscInt,PetscInt[],PetscBool,const PetscInt***,PetscInt,PetscInt[]);
490: #endif /* _PLEXIMPL_H */