Actual source code: dmpleximpl.h
petsc-3.12.0 2019-09-29
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>
10: PETSC_EXTERN PetscLogEvent DMPLEX_Interpolate;
11: PETSC_EXTERN PetscLogEvent DMPLEX_Partition;
12: PETSC_EXTERN PetscLogEvent DMPLEX_PartSelf;
13: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelInvert;
14: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelCreateSF;
15: PETSC_EXTERN PetscLogEvent DMPLEX_PartStratSF;
16: PETSC_EXTERN PetscLogEvent DMPLEX_CreatePointSF;
17: PETSC_EXTERN PetscLogEvent DMPLEX_Distribute;
18: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeCones;
19: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeLabels;
20: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeSF;
21: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeOverlap;
22: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeField;
23: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeData;
24: PETSC_EXTERN PetscLogEvent DMPLEX_Migrate;
25: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolateSF;
26: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalBegin;
27: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalEnd;
28: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalBegin;
29: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalEnd;
30: PETSC_EXTERN PetscLogEvent DMPLEX_Stratify;
31: PETSC_EXTERN PetscLogEvent DMPLEX_Symmetrize;
32: PETSC_EXTERN PetscLogEvent DMPLEX_Preallocate;
33: PETSC_EXTERN PetscLogEvent DMPLEX_ResidualFEM;
34: PETSC_EXTERN PetscLogEvent DMPLEX_JacobianFEM;
35: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolatorFEM;
36: PETSC_EXTERN PetscLogEvent DMPLEX_InjectorFEM;
37: PETSC_EXTERN PetscLogEvent DMPLEX_IntegralFEM;
38: PETSC_EXTERN PetscLogEvent DMPLEX_CreateGmsh;
39: PETSC_EXTERN PetscLogEvent DMPLEX_RebalanceSharedPoints;
41: PETSC_EXTERN PetscBool PetscPartitionerRegisterAllCalled;
42: PETSC_EXTERN PetscErrorCode PetscPartitionerRegisterAll(void);
44: PETSC_EXTERN const char * const CellRefiners[];
45: typedef enum {REFINER_NOOP = 0,
46: REFINER_SIMPLEX_1D,
47: REFINER_SIMPLEX_2D,
48: REFINER_HYBRID_SIMPLEX_2D,
49: REFINER_SIMPLEX_TO_HEX_2D,
50: REFINER_HYBRID_SIMPLEX_TO_HEX_2D,
51: REFINER_HEX_2D,
52: REFINER_HYBRID_HEX_2D,
53: REFINER_SIMPLEX_3D,
54: REFINER_HYBRID_SIMPLEX_3D,
55: REFINER_SIMPLEX_TO_HEX_3D,
56: REFINER_HYBRID_SIMPLEX_TO_HEX_3D,
57: REFINER_HEX_3D,
58: REFINER_HYBRID_HEX_3D} CellRefiner;
60: typedef struct _PetscPartitionerOps *PetscPartitionerOps;
61: struct _PetscPartitionerOps {
62: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscPartitioner);
63: PetscErrorCode (*setup)(PetscPartitioner);
64: PetscErrorCode (*view)(PetscPartitioner,PetscViewer);
65: PetscErrorCode (*destroy)(PetscPartitioner);
66: PetscErrorCode (*partition)(PetscPartitioner, DM, PetscInt, PetscInt, PetscInt[], PetscInt[], PetscSection, IS *);
67: };
69: struct _p_PetscPartitioner {
70: PETSCHEADER(struct _PetscPartitionerOps);
71: void *data; /* Implementation object */
72: PetscInt height; /* Height of points to partition into non-overlapping subsets */
73: PetscInt edgeCut; /* The number of edge cut by the partition */
74: PetscReal balance; /* The maximum partition size divided by the minimum size */
75: PetscViewer viewerGraph;
76: PetscViewerFormat formatGraph;
77: PetscBool viewGraph;
78: PetscBool noGraph; /* if true, the partitioner does not need the connectivity graph, only the number of local vertices */
79: };
81: typedef struct {
82: PetscInt dummy;
83: } PetscPartitioner_Chaco;
85: typedef struct {
86: PetscInt ptype;
87: PetscReal imbalanceRatio;
88: PetscInt debugFlag;
89: PetscInt randomSeed;
90: } PetscPartitioner_ParMetis;
92: typedef struct {
93: PetscInt strategy;
94: PetscReal imbalance;
95: } PetscPartitioner_PTScotch;
97: static const char *const
98: PTScotchStrategyList[] = {
99: "DEFAULT",
100: "QUALITY",
101: "SPEED",
102: "BALANCE",
103: "SAFETY",
104: "SCALABILITY",
105: "RECURSIVE",
106: "REMAP"
107: };
109: typedef struct {
110: PetscSection section; /* Sizes for each partition */
111: IS partition; /* Points in each partition */
112: PetscBool random; /* Flag for a random partition */
113: } PetscPartitioner_Shell;
115: typedef struct {
116: PetscInt dummy;
117: } PetscPartitioner_Simple;
119: typedef struct {
120: PetscInt dummy;
121: } PetscPartitioner_Gather;
123: /* Utility struct to store the contents of a Fluent file in memory */
124: typedef struct {
125: int index; /* Type of section */
126: unsigned int zoneID;
127: unsigned int first;
128: unsigned int last;
129: int type;
130: int nd; /* Either ND or element-type */
131: void *data;
132: } FluentSection;
134: struct _PetscGridHash {
135: PetscInt dim;
136: PetscReal lower[3]; /* The lower-left corner */
137: PetscReal upper[3]; /* The upper-right corner */
138: PetscReal extent[3]; /* The box size */
139: PetscReal h[3]; /* The subbox size */
140: PetscInt n[3]; /* The number of subboxes */
141: PetscSection cellSection; /* Offsets for cells in each subbox*/
142: IS cells; /* List of cells in each subbox */
143: DMLabel cellsSparse; /* Sparse storage for cell map */
144: };
146: typedef struct {
147: PetscInt refct;
149: /* Sieve */
150: PetscSection coneSection; /* Layout of cones (inedges for DAG) */
151: PetscInt maxConeSize; /* Cached for fast lookup */
152: PetscInt *cones; /* Cone for each point */
153: 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 */
154: PetscSection supportSection; /* Layout of cones (inedges for DAG) */
155: PetscInt maxSupportSize; /* Cached for fast lookup */
156: PetscInt *supports; /* Cone for each point */
157: PetscBool refinementUniform; /* Flag for uniform cell refinement */
158: PetscReal refinementLimit; /* Maximum volume for refined cell */
159: PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *); /* Function giving the maximum volume for refined cell */
160: PetscInt hybridPointMax[8]; /* Allow segregation of some points, each dimension has a divider (used in VTK output and refinement) */
161: PetscInt overlap; /* Overlap of the partitions as passed to DMPlexDistribute() or DMPlexDistributeOverlap() */
163: PetscInt *facesTmp; /* Work space for faces operation */
165: /* Hierarchy */
166: PetscBool regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */
168: /* Generation */
169: char *tetgenOpts;
170: char *triangleOpts;
171: PetscPartitioner partitioner;
172: PetscBool partitionBalance; /* Evenly divide partition overlap when distributing */
173: PetscBool remeshBd;
175: /* Submesh */
176: DMLabel subpointMap; /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
178: /* Labels and numbering */
179: PetscObjectState depthState; /* State of depth label, so that we can determine if a user changes it */
180: IS globalVertexNumbers;
181: IS globalCellNumbers;
183: /* Constraints */
184: PetscSection anchorSection; /* maps constrained points to anchor points */
185: IS anchorIS; /* anchors indexed by the above section */
186: PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
187: PetscErrorCode (*computeanchormatrix)(DM,PetscSection,PetscSection,Mat);
189: /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
190: PetscSection parentSection; /* dof == 1 if point has parent */
191: PetscInt *parents; /* point to parent */
192: PetscInt *childIDs; /* point to child ID */
193: PetscSection childSection; /* inverse of parent section */
194: PetscInt *children; /* point to children */
195: DM referenceTree; /* reference tree to which child ID's refer */
196: PetscErrorCode (*getchildsymmetry)(DM,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*);
198: /* MATIS support */
199: PetscSection subdomainSection;
201: /* Adjacency */
202: PetscBool useAnchors; /* Replace constrained points with their anchors in adjacency lists */
203: PetscErrorCode (*useradjacency)(DM,PetscInt,PetscInt*,PetscInt[],void*); /* User callback for adjacency */
204: void *useradjacencyctx; /* User context for callback */
206: /* Projection */
207: PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
209: /* Output */
210: PetscInt vtkCellHeight; /* The height of cells for output, default is 0 */
211: PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */
213: /* Geometry */
214: PetscReal minradius; /* Minimum distance from cell centroid to face */
215: PetscBool useHashLocation; /* Use grid hashing for point location */
216: PetscGridHash lbox; /* Local box for searching */
218: /* Debugging */
219: PetscBool printSetValues;
220: PetscInt printFEM;
221: PetscInt printL2;
222: PetscReal printTol;
223: } DM_Plex;
225: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
226: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
227: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec,PetscViewer);
228: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
229: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
230: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec,PetscViewer);
231: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
232: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
233: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM,PetscViewer);
234: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject,PetscViewer);
235: #if defined(PETSC_HAVE_HDF5)
236: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
237: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
238: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
239: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
240: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
241: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
242: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
243: #endif
245: PETSC_INTERN PetscErrorCode DMPlexClosurePoints_Private(DM,PetscInt,const PetscInt[],IS*);
246: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *, DM);
247: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
248: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM []);
249: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
250: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM []);
251: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, DMLabel, DM *);
252: PETSC_INTERN PetscErrorCode DMAdaptMetric_Plex(DM, Vec, DMLabel, DM *);
253: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
254: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
255: 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);
256: 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);
257: 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);
258: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
259: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *);
260: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
261: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);
263: PETSC_INTERN PetscErrorCode DMPlexBuildFromCellList_Internal(DM, PetscInt, PetscInt, PetscInt, PetscInt, const int[], PetscBool);
264: PETSC_INTERN PetscErrorCode DMPlexBuildFromCellList_Parallel_Internal(DM, PetscInt, PetscInt, PetscInt, PetscInt, const int[], PetscBool, PetscSF *);
265: PETSC_INTERN PetscErrorCode DMPlexBuildCoordinates_Internal(DM, PetscInt, PetscInt, PetscInt, const double[]);
266: PETSC_INTERN PetscErrorCode DMPlexBuildCoordinates_Parallel_Internal(DM, PetscInt, PetscInt, PetscInt, PetscSF, const PetscReal[]);
267: PETSC_INTERN PetscErrorCode DMPlexLoadLabels_HDF5_Internal(DM, PetscViewer);
268: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
269: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
270: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM, PetscViewer);
271: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
272: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
273: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
274: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
275: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
276: /* TODO Make these INTERN */
277: PETSC_EXTERN PetscErrorCode DMPlexView_ExodusII_Internal(DM, int, PetscInt);
278: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec, int, int);
279: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec, int, int);
280: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec, int, int);
281: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec, int, int);
282: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM,PetscInt,PetscInt,PetscInt*);
283: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
284: PETSC_INTERN PetscErrorCode DMPlexGetFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
285: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,PetscInt,PetscInt,const PetscInt[], PetscInt*,PetscInt*,const PetscInt*[]);
286: PETSC_INTERN PetscErrorCode DMPlexRestoreFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
287: PETSC_INTERN PetscErrorCode DMPlexRefineUniform_Internal(DM,CellRefiner,DM*);
288: PETSC_INTERN PetscErrorCode DMPlexGetCellRefiner_Internal(DM,CellRefiner*);
289: PETSC_INTERN PetscErrorCode CellRefinerGetAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
290: PETSC_INTERN PetscErrorCode CellRefinerGetAffineFaceTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[], PetscReal *[]);
291: PETSC_INTERN PetscErrorCode CellRefinerRestoreAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
292: PETSC_INTERN PetscErrorCode CellRefinerInCellTest_Internal(CellRefiner, const PetscReal[], PetscBool *);
293: PETSC_INTERN PetscErrorCode DMPlexInvertCell_Internal(PetscInt, PetscInt, PetscInt[]);
294: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], const PetscScalar[], InsertMode);
295: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
296: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
297: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM,DM,const char *,DM*);
298: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM,DM,PetscSF,PetscInt *,Mat);
299: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM,DM,PetscSF,PetscInt *,Mat);
300: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM,PetscSection,PetscInt,PetscInt,const PetscInt[],const PetscInt ***,const PetscScalar[],PetscInt*,PetscInt*,PetscInt*[],PetscScalar*[],PetscInt[],PetscBool);
301: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt []);
302: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt [],PetscBool,PetscInt,PetscInt []);
303: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM,PetscInt,const PetscScalar [],PetscInt,PetscInt *);
304: PETSC_EXTERN PetscErrorCode DMPlexOrientCell_Internal(DM,PetscInt,PetscInt,PetscBool);
305: PETSC_EXTERN PetscErrorCode DMPlexOrientInterface(DM);
307: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
308: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
309: PETSC_INTERN PetscErrorCode DMPlexCreateNumbering_Internal(DM, PetscInt, PetscInt, PetscInt, PetscInt *, PetscSF, IS *);
310: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, DMLabel, DM *);
311: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, DMLabel, DM *);
312: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM, Mat*);
314: PETSC_INTERN PetscErrorCode DMPlexGetOverlap_Plex(DM, PetscInt *);
316: /* invert dihedral symmetry: return a^-1,
317: * using the representation described in
318: * DMPlexGetConeOrientation() */
319: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
320: {
321: return (a <= 0) ? a : (N - a);
322: }
324: /* invert dihedral symmetry: return b * a,
325: * using the representation described in
326: * DMPlexGetConeOrientation() */
327: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
328: {
329: if (!N) return 0;
330: return (a >= 0) ?
331: ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
332: ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
333: }
335: /* swap dihedral symmetries: return b * a^-1,
336: * using the representation described in
337: * DMPlexGetConeOrientation() */
338: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
339: {
340: return DihedralCompose(N,DihedralInvert(N,a),b);
341: }
343: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, IS , PetscReal, Vec, Vec, PetscReal, Vec, void *);
344: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
345: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);
347: /* Matvec with A in row-major storage, x and y can be aliased */
348: PETSC_STATIC_INLINE void DMPlex_Mult2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
349: {
350: PetscScalar z[2];
351: z[0] = x[0]; z[1] = x[ldx];
352: y[0] = A[0]*z[0] + A[1]*z[1];
353: y[ldx] = A[2]*z[0] + A[3]*z[1];
354: (void)PetscLogFlops(6.0);
355: }
356: PETSC_STATIC_INLINE void DMPlex_Mult3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
357: {
358: PetscScalar z[3];
359: z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
360: y[0] = A[0]*z[0] + A[1]*z[1] + A[2]*z[2];
361: y[ldx] = A[3]*z[0] + A[4]*z[1] + A[5]*z[2];
362: y[ldx*2] = A[6]*z[0] + A[7]*z[1] + A[8]*z[2];
363: (void)PetscLogFlops(15.0);
364: }
365: PETSC_STATIC_INLINE void DMPlex_MultTranspose2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
366: {
367: PetscScalar z[2];
368: z[0] = x[0]; z[1] = x[ldx];
369: y[0] = A[0]*z[0] + A[2]*z[1];
370: y[ldx] = A[1]*z[0] + A[3]*z[1];
371: (void)PetscLogFlops(6.0);
372: }
373: PETSC_STATIC_INLINE void DMPlex_MultTranspose3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
374: {
375: PetscScalar z[3];
376: z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
377: y[0] = A[0]*z[0] + A[3]*z[1] + A[6]*z[2];
378: y[ldx] = A[1]*z[0] + A[4]*z[1] + A[7]*z[2];
379: y[ldx*2] = A[2]*z[0] + A[5]*z[1] + A[8]*z[2];
380: (void)PetscLogFlops(15.0);
381: }
382: PETSC_STATIC_INLINE void DMPlex_Mult2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
383: {
384: PetscScalar z[2];
385: z[0] = x[0]; z[1] = x[ldx];
386: y[0] = A[0]*z[0] + A[1]*z[1];
387: y[ldx] = A[2]*z[0] + A[3]*z[1];
388: (void)PetscLogFlops(6.0);
389: }
390: PETSC_STATIC_INLINE void DMPlex_Mult3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
391: {
392: PetscScalar z[3];
393: z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
394: y[0] = A[0]*z[0] + A[1]*z[1] + A[2]*z[2];
395: y[ldx] = A[3]*z[0] + A[4]*z[1] + A[5]*z[2];
396: y[ldx*2] = A[6]*z[0] + A[7]*z[1] + A[8]*z[2];
397: (void)PetscLogFlops(15.0);
398: }
399: PETSC_STATIC_INLINE void DMPlex_MultTranspose2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
400: {
401: PetscScalar z[2];
402: z[0] = x[0]; z[1] = x[ldx];
403: y[0] = A[0]*z[0] + A[2]*z[1];
404: y[ldx] = A[1]*z[0] + A[3]*z[1];
405: (void)PetscLogFlops(6.0);
406: }
407: PETSC_STATIC_INLINE void DMPlex_MultTranspose3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
408: {
409: PetscScalar z[3];
410: z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
411: y[0] = A[0]*z[0] + A[3]*z[1] + A[6]*z[2];
412: y[ldx] = A[1]*z[0] + A[4]*z[1] + A[7]*z[2];
413: y[ldx*2] = A[2]*z[0] + A[5]*z[1] + A[8]*z[2];
414: (void)PetscLogFlops(15.0);
415: }
417: PETSC_STATIC_INLINE void DMPlex_MatMult2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
418: {
419: PetscInt j;
420: for (j = 0; j < n; ++j) {
421: PetscScalar z[2];
422: z[0] = B[0+j]; z[1] = B[1*ldb+j];
423: DMPlex_Mult2D_Internal(A, 1, z, z);
424: C[0+j] = z[0]; C[1*ldb+j] = z[1];
425: }
426: (void)PetscLogFlops(8.0*n);
427: }
428: PETSC_STATIC_INLINE void DMPlex_MatMult3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
429: {
430: PetscInt j;
431: for (j = 0; j < n; ++j) {
432: PetscScalar z[3];
433: z[0] = B[0+j]; z[1] = B[1*ldb+j]; z[2] = B[2*ldb+j];
434: DMPlex_Mult3D_Internal(A, 1, z, z);
435: C[0+j] = z[0]; C[1*ldb+j] = z[1]; C[2*ldb+j] = z[2];
436: }
437: (void)PetscLogFlops(8.0*n);
438: }
439: PETSC_STATIC_INLINE void DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
440: {
441: PetscInt j;
442: for (j = 0; j < n; ++j) {
443: PetscScalar z[2];
444: z[0] = B[0+j]; z[1] = B[1*ldb+j];
445: DMPlex_MultTranspose2D_Internal(A, 1, z, z);
446: C[0+j] = z[0]; C[1*ldb+j] = z[1];
447: }
448: (void)PetscLogFlops(8.0*n);
449: }
450: PETSC_STATIC_INLINE void DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
451: {
452: PetscInt j;
453: for (j = 0; j < n; ++j) {
454: PetscScalar z[3];
455: z[0] = B[0+j]; z[1] = B[1*ldb+j]; z[2] = B[2*ldb+j];
456: DMPlex_MultTranspose3D_Internal(A, 1, z, z);
457: C[0+j] = z[0]; C[1*ldb+j] = z[1]; C[2*ldb+j] = z[2];
458: }
459: (void)PetscLogFlops(8.0*n);
460: }
462: PETSC_STATIC_INLINE void DMPlex_MatMultLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
463: {
464: PetscInt j;
465: for (j = 0; j < m; ++j) {
466: DMPlex_MultTranspose2D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
467: }
468: (void)PetscLogFlops(8.0*m);
469: }
470: PETSC_STATIC_INLINE void DMPlex_MatMultLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
471: {
472: PetscInt j;
473: for (j = 0; j < m; ++j) {
474: DMPlex_MultTranspose3D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
475: }
476: (void)PetscLogFlops(8.0*m);
477: }
478: PETSC_STATIC_INLINE void DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
479: {
480: PetscInt j;
481: for (j = 0; j < m; ++j) {
482: DMPlex_Mult2D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
483: }
484: (void)PetscLogFlops(8.0*m);
485: }
486: PETSC_STATIC_INLINE void DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
487: {
488: PetscInt j;
489: for (j = 0; j < m; ++j) {
490: DMPlex_Mult3D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
491: }
492: (void)PetscLogFlops(8.0*m);
493: }
495: PETSC_STATIC_INLINE void DMPlex_Transpose2D_Internal(PetscScalar A[])
496: {
497: PetscScalar tmp;
498: tmp = A[1]; A[1] = A[2]; A[2] = tmp;
499: }
500: PETSC_STATIC_INLINE void DMPlex_Transpose3D_Internal(PetscScalar A[])
501: {
502: PetscScalar tmp;
503: tmp = A[1]; A[1] = A[3]; A[3] = tmp;
504: tmp = A[2]; A[2] = A[6]; A[6] = tmp;
505: tmp = A[5]; A[5] = A[7]; A[7] = tmp;
506: }
508: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
509: {
510: const PetscReal invDet = 1.0/detJ;
512: invJ[0] = invDet*J[3];
513: invJ[1] = -invDet*J[1];
514: invJ[2] = -invDet*J[2];
515: invJ[3] = invDet*J[0];
516: (void)PetscLogFlops(5.0);
517: }
519: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
520: {
521: const PetscReal invDet = 1.0/detJ;
523: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
524: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
525: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
526: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
527: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
528: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
529: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
530: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
531: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
532: (void)PetscLogFlops(37.0);
533: }
535: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
536: {
537: *detJ = J[0]*J[3] - J[1]*J[2];
538: (void)PetscLogFlops(3.0);
539: }
541: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
542: {
543: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
544: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
545: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
546: (void)PetscLogFlops(12.0);
547: }
549: PETSC_STATIC_INLINE void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
550: {
551: *detJ = PetscRealPart(J[0])*PetscRealPart(J[3]) - PetscRealPart(J[1])*PetscRealPart(J[2]);
552: (void)PetscLogFlops(3.0);
553: }
555: PETSC_STATIC_INLINE void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
556: {
557: *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])) +
558: 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])) +
559: 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])));
560: (void)PetscLogFlops(12.0);
561: }
563: 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];}
565: 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;}
567: 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;}
569: 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);}
571: /* Compare cones of the master and slave face (with the same cone points modulo order), and return relative orientation of the slave. */
572: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Orient_Private(PetscInt coneSize, PetscInt masterConeSize, const PetscInt masterCone[], const PetscInt slaveCone[], PetscInt *start, PetscBool *reverse)
573: {
574: PetscInt i;
577: *start = 0;
578: for (i=0; i<coneSize; i++) {
579: if (slaveCone[i] == masterCone[0]) {
580: *start = i;
581: break;
582: }
583: }
584: if (PetscUnlikely(i==coneSize)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "starting point of master cone not found in slave cone");
585: *reverse = PETSC_FALSE;
586: for (i=0; i<masterConeSize; i++) {if (slaveCone[((*start)+i)%coneSize] != masterCone[i]) break;}
587: if (i == masterConeSize) return(0);
588: *reverse = PETSC_TRUE;
589: for (i=0; i<masterConeSize; i++) {if (slaveCone[(coneSize+(*start)-i)%coneSize] != masterCone[i]) break;}
590: if (i < masterConeSize) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "master and slave cone have non-conforming order of points");
591: return(0);
592: }
594: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Translate_Private(PetscInt ornt, PetscInt *start, PetscBool *reverse)
595: {
597: *reverse = (ornt < 0) ? PETSC_TRUE : PETSC_FALSE;
598: *start = *reverse ? -(ornt+1) : ornt;
599: return(0);
600: }
602: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Combine_Private(PetscInt coneSize, PetscInt origStart, PetscBool origReverse, PetscInt rotateStart, PetscBool rotateReverse, PetscInt *newStart, PetscBool *newReverse)
603: {
605: *newReverse = (origReverse == rotateReverse) ? PETSC_FALSE : PETSC_TRUE;
606: *newStart = rotateReverse ? (coneSize + rotateStart - origStart) : (coneSize + origStart - rotateStart);
607: *newStart %= coneSize;
608: return(0);
609: }
611: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_TranslateBack_Private(PetscInt coneSize, PetscInt start, PetscBool reverse, PetscInt *ornt)
612: {
614: if (coneSize < 3) {
615: /* edges just get flipped if start == 1 regardless direction */
616: *ornt = start ? -2 : 0;
617: } else {
618: *ornt = reverse ? -(start+1) : start;
619: }
620: return(0);
621: }
623: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Permute_Private(PetscInt n, const PetscInt arr[], PetscInt start, PetscBool reverse, PetscInt newarr[])
624: {
625: PetscInt i;
628: if (reverse) {for (i=0; i<n; i++) newarr[i] = arr[(n+start-i)%n];}
629: else {for (i=0; i<n; i++) newarr[i] = arr[(start+i)%n];}
630: return(0);
631: }
633: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM,PetscInt,PetscInt,PetscDualSpace *);
634: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,const PetscInt[],const PetscInt[],PetscInt[]);
635: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection,PetscInt,PetscInt,PetscInt[],PetscBool,const PetscInt***,PetscInt,const PetscInt[],PetscInt[]);
636: PETSC_INTERN PetscErrorCode DMPlexGetCompressedClosure(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscSection *, IS *, const PetscInt **);
637: PETSC_INTERN PetscErrorCode DMPlexRestoreCompressedClosure(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscSection *, IS *, const PetscInt **);
639: PETSC_EXTERN PetscErrorCode DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
640: PETSC_EXTERN PetscErrorCode DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
641: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
642: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
643: PETSC_INTERN PetscErrorCode DMCreateSubDomainDM_Plex(DM,DMLabel,PetscInt,IS*,DM*);
644: PETSC_INTERN PetscErrorCode DMPlexBasisTransformPoint_Internal(DM, DM, Vec, PetscInt, PetscBool[], PetscBool, PetscScalar *);
645: PETSC_EXTERN PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM, DM, Vec, PetscInt, PetscBool, PetscInt, PetscScalar *);
646: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscReal *, PetscReal *, void *);
647: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApply_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscScalar *, PetscScalar *, void *);
648: PETSC_INTERN PetscErrorCode DMCreateNeumannOverlap_Plex(DM, IS*, Mat*, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void **);
650: #endif /* _PLEXIMPL_H */