Actual source code: dmpleximpl.h
petsc-3.8.4 2018-03-24
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 useCone; /* Use cone() first when defining adjacency */
156: PetscBool useClosure; /* Use the transitive closure when defining adjacency */
157: PetscBool useAnchors; /* Replace constrained points with their anchors in adjacency lists */
158: PetscErrorCode (*useradjacency)(DM,PetscInt,PetscInt*,PetscInt[],void*); /* User callback for adjacency */
159: void *useradjacencyctx; /* User context for callback */
161: /* Projection */
162: PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
164: /* Output */
165: PetscInt vtkCellHeight; /* The height of cells for output, default is 0 */
166: PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */
168: /* Geometry */
169: PetscReal minradius; /* Minimum distance from cell centroid to face */
170: PetscBool useHashLocation; /* Use grid hashing for point location */
171: PetscGridHash lbox; /* Local box for searching */
173: /* Debugging */
174: PetscBool printSetValues;
175: PetscInt printFEM;
176: PetscReal printTol;
177: } DM_Plex;
179: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
180: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
181: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec,PetscViewer);
182: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
183: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
184: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec,PetscViewer);
185: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
186: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
187: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM,PetscViewer);
188: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject,PetscViewer);
189: #if defined(PETSC_HAVE_HDF5)
190: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
191: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
192: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
193: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
194: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
195: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
196: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
197: #endif
199: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *, DM);
200: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
201: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM []);
202: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
203: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM []);
204: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, DMLabel, DM *);
205: PETSC_INTERN PetscErrorCode DMAdaptMetric_Plex(DM, Vec, DMLabel, DM *);
206: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
207: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
208: 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);
209: 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);
210: 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);
211: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
212: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *);
213: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
214: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);
216: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
217: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
218: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
219: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
220: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
221: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
222: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
223: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM,PetscInt,PetscInt,PetscInt*);
224: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
225: PETSC_INTERN PetscErrorCode DMPlexGetFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
226: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,PetscInt,PetscInt,const PetscInt[], PetscInt*,PetscInt*,const PetscInt*[]);
227: PETSC_INTERN PetscErrorCode DMPlexRestoreFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
228: PETSC_INTERN PetscErrorCode DMPlexRefineUniform_Internal(DM,CellRefiner,DM*);
229: PETSC_INTERN PetscErrorCode DMPlexGetCellRefiner_Internal(DM,CellRefiner*);
230: PETSC_INTERN PetscErrorCode CellRefinerGetAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
231: PETSC_INTERN PetscErrorCode CellRefinerRestoreAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
232: PETSC_INTERN PetscErrorCode CellRefinerInCellTest_Internal(CellRefiner, const PetscReal[], PetscBool *);
233: PETSC_INTERN PetscErrorCode DMPlexInvertCell_Internal(PetscInt, PetscInt, PetscInt[]);
234: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, const PetscScalar[], InsertMode);
235: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
236: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM,DM,const char *,DM*);
237: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM,DM,PetscSF,PetscInt *,Mat);
238: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM,DM,PetscSF,PetscInt *,Mat);
239: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM,PetscSection,PetscInt,PetscInt,const PetscInt[],const PetscInt ***,const PetscScalar[],PetscInt*,PetscInt*,PetscInt*[],PetscScalar*[],PetscInt[],PetscBool);
240: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt []);
241: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt [],PetscBool,PetscInt,PetscInt []);
242: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM,PetscInt,const PetscScalar [],PetscInt,PetscInt *);
244: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
245: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
246: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, DMLabel, DM *);
247: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, DMLabel, DM *);
248: #if defined(PETSC_HAVE_TRIANGLE)
249: PETSC_INTERN PetscErrorCode DMPlexGenerate_Triangle(DM, PetscBool, DM *);
250: PETSC_INTERN PetscErrorCode DMPlexRefine_Triangle(DM, double[], DM *);
251: #endif
252: #if defined(PETSC_HAVE_TETGEN)
253: PETSC_INTERN PetscErrorCode DMPlexGenerate_Tetgen(DM, PetscBool, DM *);
254: PETSC_INTERN PetscErrorCode DMPlexRefine_Tetgen(DM, double[], DM *);
255: #endif
256: #if defined(PETSC_HAVE_CTETGEN)
257: PETSC_INTERN PetscErrorCode DMPlexGenerate_CTetgen(DM, PetscBool, DM *);
258: PETSC_INTERN PetscErrorCode DMPlexRefine_CTetgen(DM, PetscReal[], DM *);
259: #endif
261: /* invert dihedral symmetry: return a^-1,
262: * using the representation described in
263: * DMPlexGetConeOrientation() */
264: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
265: {
266: return (a <= 0) ? a : (N - a);
267: }
269: /* invert dihedral symmetry: return b * a,
270: * using the representation described in
271: * DMPlexGetConeOrientation() */
272: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
273: {
274: if (!N) return 0;
275: return (a >= 0) ?
276: ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
277: ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
278: }
280: /* swap dihedral symmetries: return b * a^-1,
281: * using the representation described in
282: * DMPlexGetConeOrientation() */
283: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
284: {
285: return DihedralCompose(N,DihedralInvert(N,a),b);
286: }
288: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, PetscInt, PetscInt, PetscReal, Vec, Vec, PetscReal, Vec, void *);
289: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, PetscInt, PetscInt, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
290: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);
292: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
293: {
294: const PetscReal invDet = 1.0/detJ;
296: invJ[0] = invDet*J[3];
297: invJ[1] = -invDet*J[1];
298: invJ[2] = -invDet*J[2];
299: invJ[3] = invDet*J[0];
300: (void)PetscLogFlops(5.0);
301: }
303: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
304: {
305: const PetscReal invDet = 1.0/detJ;
307: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
308: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
309: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
310: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
311: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
312: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
313: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
314: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
315: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
316: (void)PetscLogFlops(37.0);
317: }
319: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, PetscReal J[])
320: {
321: *detJ = J[0]*J[3] - J[1]*J[2];
322: (void)PetscLogFlops(3.0);
323: }
325: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, PetscReal J[])
326: {
327: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
328: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
329: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
330: (void)PetscLogFlops(12.0);
331: }
333: 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];}
335: 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;}
337: 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;}
339: 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);}
341: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
342: {
344: #if defined(PETSC_USE_DEBUG)
345: {
346: PetscInt dof;
348: *start = *end = 0; /* Silence overzealous compiler warning */
349: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
350: PetscSectionGetOffset(dm->defaultSection, point, start);
351: PetscSectionGetDof(dm->defaultSection, point, &dof);
352: *end = *start + dof;
353: }
354: #else
355: {
356: const PetscSection s = dm->defaultSection;
357: *start = s->atlasOff[point - s->pStart];
358: *end = *start + s->atlasDof[point - s->pStart];
359: }
360: #endif
361: return(0);
362: }
364: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
365: {
367: #if defined(PETSC_USE_DEBUG)
368: {
369: PetscInt dof;
371: *start = *end = 0; /* Silence overzealous compiler warning */
372: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
373: PetscSectionGetFieldOffset(dm->defaultSection, point, field, start);
374: PetscSectionGetFieldDof(dm->defaultSection, point, field, &dof);
375: *end = *start + dof;
376: }
377: #else
378: {
379: const PetscSection s = dm->defaultSection->field[field];
380: *start = s->atlasOff[point - s->pStart];
381: *end = *start + s->atlasDof[point - s->pStart];
382: }
383: #endif
384: return(0);
385: }
387: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
388: {
390: #if defined(PETSC_USE_DEBUG)
391: {
393: PetscInt dof,cdof;
394: *start = *end = 0; /* Silence overzealous compiler warning */
395: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
396: 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()");
397: PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
398: PetscSectionGetDof(dm->defaultGlobalSection, point, &dof);
399: PetscSectionGetConstraintDof(dm->defaultGlobalSection, point, &cdof);
400: *end = *start + dof - cdof + (dof < 0 ? 1 : 0);
401: }
402: #else
403: {
404: const PetscSection s = dm->defaultGlobalSection;
405: const PetscInt dof = s->atlasDof[point - s->pStart];
406: const PetscInt cdof = s->bc ? s->bc->atlasDof[point - s->bc->pStart] : 0;
407: *start = s->atlasOff[point - s->pStart];
408: *end = *start + dof - cdof + (dof < 0 ? 1 : 0);
409: }
410: #endif
411: return(0);
412: }
414: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
415: {
417: #if defined(PETSC_USE_DEBUG)
418: {
419: PetscInt loff, lfoff, fdof, fcdof, ffcdof, f;
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 crated automatically by DMGetDefaultGlobalSection()");
424: PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
425: PetscSectionGetOffset(dm->defaultSection, point, &loff);
426: PetscSectionGetFieldOffset(dm->defaultSection, point, field, &lfoff);
427: PetscSectionGetFieldDof(dm->defaultSection, point, field, &fdof);
428: PetscSectionGetFieldConstraintDof(dm->defaultSection, point, field, &fcdof);
429: *start = *start < 0 ? *start - (lfoff-loff) : *start + lfoff-loff;
430: for (f = 0; f < field; ++f) {
431: PetscSectionGetFieldConstraintDof(dm->defaultSection, point, f, &ffcdof);
432: *start = *start < 0 ? *start + ffcdof : *start - ffcdof;
433: }
434: *end = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
435: }
436: #else
437: {
438: const PetscSection s = dm->defaultSection;
439: const PetscSection fs = dm->defaultSection->field[field];
440: const PetscSection gs = dm->defaultGlobalSection;
441: const PetscInt loff = s->atlasOff[point - s->pStart];
442: const PetscInt goff = gs->atlasOff[point - s->pStart];
443: const PetscInt lfoff = fs->atlasOff[point - s->pStart];
444: const PetscInt fdof = fs->atlasDof[point - s->pStart];
445: const PetscInt fcdof = fs->bc ? fs->bc->atlasDof[point - fs->bc->pStart] : 0;
446: PetscInt ffcdof = 0, f;
448: for (f = 0; f < field; ++f) {
449: const PetscSection ffs = dm->defaultSection->field[f];
450: ffcdof += ffs->bc ? ffs->bc->atlasDof[point - ffs->bc->pStart] : 0;
451: }
452: *start = goff + (goff < 0 ? loff-lfoff + ffcdof : lfoff-loff - ffcdof);
453: *end = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
454: }
455: #endif
456: return(0);
457: }
459: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM,PetscInt,PetscInt,PetscDualSpace *);
460: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,const PetscInt[],PetscInt[]);
461: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection,PetscInt,PetscInt,PetscInt[],PetscBool,const PetscInt***,PetscInt,PetscInt[]);
463: #endif /* _PLEXIMPL_H */