Actual source code: dmpleximpl.h
petsc-3.11.4 2019-09-28
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_Distribute;
13: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeCones;
14: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeLabels;
15: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeSF;
16: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeOverlap;
17: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeField;
18: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeData;
19: PETSC_EXTERN PetscLogEvent DMPLEX_Migrate;
20: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolateSF;
21: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalBegin;
22: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalEnd;
23: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalBegin;
24: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalEnd;
25: PETSC_EXTERN PetscLogEvent DMPLEX_Stratify;
26: PETSC_EXTERN PetscLogEvent DMPLEX_Preallocate;
27: PETSC_EXTERN PetscLogEvent DMPLEX_ResidualFEM;
28: PETSC_EXTERN PetscLogEvent DMPLEX_JacobianFEM;
29: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolatorFEM;
30: PETSC_EXTERN PetscLogEvent DMPLEX_InjectorFEM;
31: PETSC_EXTERN PetscLogEvent DMPLEX_IntegralFEM;
32: PETSC_EXTERN PetscLogEvent DMPLEX_CreateGmsh;
34: PETSC_EXTERN PetscBool PetscPartitionerRegisterAllCalled;
35: PETSC_EXTERN PetscErrorCode PetscPartitionerRegisterAll(void);
37: typedef enum {REFINER_NOOP = 0,
38: REFINER_SIMPLEX_1D,
39: REFINER_SIMPLEX_2D,
40: REFINER_HYBRID_SIMPLEX_2D,
41: REFINER_SIMPLEX_TO_HEX_2D,
42: REFINER_HYBRID_SIMPLEX_TO_HEX_2D,
43: REFINER_HEX_2D,
44: REFINER_HYBRID_HEX_2D,
45: REFINER_SIMPLEX_3D,
46: REFINER_HYBRID_SIMPLEX_3D,
47: REFINER_SIMPLEX_TO_HEX_3D,
48: REFINER_HYBRID_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: PetscInt edgeCut; /* The number of edge cut by the partition */
66: PetscReal balance; /* The maximum partition size divided by the minimum size */
67: PetscViewer viewerGraph;
68: PetscViewerFormat formatGraph;
69: PetscBool viewGraph;
70: };
72: typedef struct {
73: PetscInt dummy;
74: } PetscPartitioner_Chaco;
76: typedef struct {
77: PetscInt ptype;
78: PetscReal imbalanceRatio;
79: PetscInt debugFlag;
80: } PetscPartitioner_ParMetis;
82: typedef struct {
83: PetscSection section; /* Sizes for each partition */
84: IS partition; /* Points in each partition */
85: PetscBool random; /* Flag for a random partition */
86: } PetscPartitioner_Shell;
88: typedef struct {
89: PetscInt dummy;
90: } PetscPartitioner_Simple;
92: typedef struct {
93: PetscInt dummy;
94: } PetscPartitioner_Gather;
96: /* Utility struct to store the contents of a Fluent file in memory */
97: typedef struct {
98: int index; /* Type of section */
99: unsigned int zoneID;
100: unsigned int first;
101: unsigned int last;
102: int type;
103: int nd; /* Either ND or element-type */
104: void *data;
105: } FluentSection;
107: struct _PetscGridHash {
108: PetscInt dim;
109: PetscReal lower[3]; /* The lower-left corner */
110: PetscReal upper[3]; /* The upper-right corner */
111: PetscReal extent[3]; /* The box size */
112: PetscReal h[3]; /* The subbox size */
113: PetscInt n[3]; /* The number of subboxes */
114: PetscSection cellSection; /* Offsets for cells in each subbox*/
115: IS cells; /* List of cells in each subbox */
116: DMLabel cellsSparse; /* Sparse storage for cell map */
117: };
119: typedef struct {
120: PetscInt refct;
122: /* Sieve */
123: PetscSection coneSection; /* Layout of cones (inedges for DAG) */
124: PetscInt maxConeSize; /* Cached for fast lookup */
125: PetscInt *cones; /* Cone for each point */
126: 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 */
127: PetscSection supportSection; /* Layout of cones (inedges for DAG) */
128: PetscInt maxSupportSize; /* Cached for fast lookup */
129: PetscInt *supports; /* Cone for each point */
130: PetscBool refinementUniform; /* Flag for uniform cell refinement */
131: PetscReal refinementLimit; /* Maximum volume for refined cell */
132: PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *); /* Function giving the maximum volume for refined cell */
133: PetscInt hybridPointMax[8]; /* Allow segregation of some points, each dimension has a divider (used in VTK output and refinement) */
135: PetscInt *facesTmp; /* Work space for faces operation */
137: /* Hierarchy */
138: PetscBool regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */
140: /* Generation */
141: char *tetgenOpts;
142: char *triangleOpts;
143: PetscPartitioner partitioner;
144: PetscBool partitionBalance; /* Evenly divide partition overlap when distributing */
145: PetscBool remeshBd;
147: /* Submesh */
148: DMLabel subpointMap; /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
150: /* Labels and numbering */
151: PetscObjectState depthState; /* State of depth label, so that we can determine if a user changes it */
152: IS globalVertexNumbers;
153: IS globalCellNumbers;
155: /* Constraints */
156: PetscSection anchorSection; /* maps constrained points to anchor points */
157: IS anchorIS; /* anchors indexed by the above section */
158: PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
159: PetscErrorCode (*computeanchormatrix)(DM,PetscSection,PetscSection,Mat);
161: /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
162: PetscSection parentSection; /* dof == 1 if point has parent */
163: PetscInt *parents; /* point to parent */
164: PetscInt *childIDs; /* point to child ID */
165: PetscSection childSection; /* inverse of parent section */
166: PetscInt *children; /* point to children */
167: DM referenceTree; /* reference tree to which child ID's refer */
168: PetscErrorCode (*getchildsymmetry)(DM,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*);
170: /* MATIS support */
171: PetscSection subdomainSection;
173: /* Adjacency */
174: PetscBool useAnchors; /* Replace constrained points with their anchors in adjacency lists */
175: PetscErrorCode (*useradjacency)(DM,PetscInt,PetscInt*,PetscInt[],void*); /* User callback for adjacency */
176: void *useradjacencyctx; /* User context for callback */
178: /* Projection */
179: PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
181: /* Output */
182: PetscInt vtkCellHeight; /* The height of cells for output, default is 0 */
183: PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */
185: /* Geometry */
186: PetscReal minradius; /* Minimum distance from cell centroid to face */
187: PetscBool useHashLocation; /* Use grid hashing for point location */
188: PetscGridHash lbox; /* Local box for searching */
190: /* Debugging */
191: PetscBool printSetValues;
192: PetscInt printFEM;
193: PetscInt printL2;
194: PetscReal printTol;
195: } DM_Plex;
197: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
198: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
199: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec,PetscViewer);
200: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
201: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
202: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec,PetscViewer);
203: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
204: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
205: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM,PetscViewer);
206: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject,PetscViewer);
207: #if defined(PETSC_HAVE_HDF5)
208: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
209: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
210: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
211: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
212: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
213: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
214: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
215: #endif
217: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *, DM);
218: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
219: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM []);
220: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
221: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM []);
222: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, DMLabel, DM *);
223: PETSC_INTERN PetscErrorCode DMAdaptMetric_Plex(DM, Vec, DMLabel, DM *);
224: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
225: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
226: 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);
227: 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);
228: 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);
229: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
230: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *);
231: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
232: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);
234: PETSC_INTERN PetscErrorCode DMPlexBuildFromCellList_Internal(DM, PetscInt, PetscInt, PetscInt, PetscInt, const int[], PetscBool);
235: PETSC_INTERN PetscErrorCode DMPlexBuildFromCellList_Parallel_Internal(DM, PetscInt, PetscInt, PetscInt, PetscInt, const int[], PetscBool, PetscSF *);
236: PETSC_INTERN PetscErrorCode DMPlexBuildCoordinates_Internal(DM, PetscInt, PetscInt, PetscInt, const double[]);
237: PETSC_INTERN PetscErrorCode DMPlexBuildCoordinates_Parallel_Internal(DM, PetscInt, PetscInt, PetscInt, PetscSF, const PetscReal[]);
238: PETSC_INTERN PetscErrorCode DMPlexLoadLabels_HDF5_Internal(DM, PetscViewer);
239: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
240: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
241: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM, PetscViewer);
242: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
243: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
244: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
245: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
246: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
247: /* TODO Make these INTERN */
248: PETSC_EXTERN PetscErrorCode DMPlexView_ExodusII_Internal(DM, int, PetscInt);
249: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec, int, int);
250: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec, int, int);
251: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec, int, int);
252: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec, int, int);
253: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM,PetscInt,PetscInt,PetscInt*);
254: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
255: PETSC_INTERN PetscErrorCode DMPlexGetFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
256: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,PetscInt,PetscInt,const PetscInt[], PetscInt*,PetscInt*,const PetscInt*[]);
257: PETSC_INTERN PetscErrorCode DMPlexRestoreFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
258: PETSC_INTERN PetscErrorCode DMPlexRefineUniform_Internal(DM,CellRefiner,DM*);
259: PETSC_INTERN PetscErrorCode DMPlexGetCellRefiner_Internal(DM,CellRefiner*);
260: PETSC_INTERN PetscErrorCode CellRefinerGetAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
261: PETSC_INTERN PetscErrorCode CellRefinerRestoreAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
262: PETSC_INTERN PetscErrorCode CellRefinerInCellTest_Internal(CellRefiner, const PetscReal[], PetscBool *);
263: PETSC_INTERN PetscErrorCode DMPlexInvertCell_Internal(PetscInt, PetscInt, PetscInt[]);
264: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], const PetscScalar[], InsertMode);
265: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
266: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
267: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM,DM,const char *,DM*);
268: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM,DM,PetscSF,PetscInt *,Mat);
269: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM,DM,PetscSF,PetscInt *,Mat);
270: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM,PetscSection,PetscInt,PetscInt,const PetscInt[],const PetscInt ***,const PetscScalar[],PetscInt*,PetscInt*,PetscInt*[],PetscScalar*[],PetscInt[],PetscBool);
271: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt []);
272: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt [],PetscBool,PetscInt,PetscInt []);
273: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM,PetscInt,const PetscScalar [],PetscInt,PetscInt *);
274: PETSC_EXTERN PetscErrorCode DMPlexOrientCell_Internal(DM,PetscInt,PetscInt,PetscBool);
275: PETSC_EXTERN PetscErrorCode DMPlexOrientInterface(DM);
277: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
278: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
279: PETSC_INTERN PetscErrorCode DMPlexCreateNumbering_Internal(DM, PetscInt, PetscInt, PetscInt, PetscInt *, PetscSF, IS *);
280: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, DMLabel, DM *);
281: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, DMLabel, DM *);
283: /* invert dihedral symmetry: return a^-1,
284: * using the representation described in
285: * DMPlexGetConeOrientation() */
286: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
287: {
288: return (a <= 0) ? a : (N - a);
289: }
291: /* invert dihedral symmetry: return b * a,
292: * using the representation described in
293: * DMPlexGetConeOrientation() */
294: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
295: {
296: if (!N) return 0;
297: return (a >= 0) ?
298: ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
299: ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
300: }
302: /* swap dihedral symmetries: return b * a^-1,
303: * using the representation described in
304: * DMPlexGetConeOrientation() */
305: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
306: {
307: return DihedralCompose(N,DihedralInvert(N,a),b);
308: }
310: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, IS , PetscReal, Vec, Vec, PetscReal, Vec, void *);
311: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
312: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);
314: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
315: {
316: const PetscReal invDet = 1.0/detJ;
318: invJ[0] = invDet*J[3];
319: invJ[1] = -invDet*J[1];
320: invJ[2] = -invDet*J[2];
321: invJ[3] = invDet*J[0];
322: (void)PetscLogFlops(5.0);
323: }
325: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
326: {
327: const PetscReal invDet = 1.0/detJ;
329: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
330: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
331: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
332: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
333: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
334: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
335: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
336: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
337: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
338: (void)PetscLogFlops(37.0);
339: }
341: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
342: {
343: *detJ = J[0]*J[3] - J[1]*J[2];
344: (void)PetscLogFlops(3.0);
345: }
347: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
348: {
349: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
350: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
351: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
352: (void)PetscLogFlops(12.0);
353: }
355: PETSC_STATIC_INLINE void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
356: {
357: *detJ = PetscRealPart(J[0])*PetscRealPart(J[3]) - PetscRealPart(J[1])*PetscRealPart(J[2]);
358: (void)PetscLogFlops(3.0);
359: }
361: PETSC_STATIC_INLINE void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
362: {
363: *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])) +
364: 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])) +
365: 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])));
366: (void)PetscLogFlops(12.0);
367: }
369: 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];}
371: 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;}
373: 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;}
375: 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);}
377: /* Compare cones of the master and slave face (with the same cone points modulo order), and return relative orientation of the slave. */
378: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Orient_Private(PetscInt coneSize, PetscInt masterConeSize, const PetscInt masterCone[], const PetscInt slaveCone[], PetscInt *start, PetscBool *reverse)
379: {
380: PetscInt i;
383: *start = 0;
384: for (i=0; i<coneSize; i++) {
385: if (slaveCone[i] == masterCone[0]) {
386: *start = i;
387: break;
388: }
389: }
390: if (PetscUnlikely(i==coneSize)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "starting point of master cone not found in slave cone");
391: *reverse = PETSC_FALSE;
392: for (i=0; i<masterConeSize; i++) {if (slaveCone[((*start)+i)%coneSize] != masterCone[i]) break;}
393: if (i == masterConeSize) return(0);
394: *reverse = PETSC_TRUE;
395: for (i=0; i<masterConeSize; i++) {if (slaveCone[(coneSize+(*start)-i)%coneSize] != masterCone[i]) break;}
396: if (i < masterConeSize) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "master and slave cone have non-conforming order of points");
397: return(0);
398: }
400: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Translate_Private(PetscInt ornt, PetscInt *start, PetscBool *reverse)
401: {
403: *reverse = (ornt < 0) ? PETSC_TRUE : PETSC_FALSE;
404: *start = *reverse ? -(ornt+1) : ornt;
405: return(0);
406: }
408: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Combine_Private(PetscInt coneSize, PetscInt origStart, PetscBool origReverse, PetscInt rotateStart, PetscBool rotateReverse, PetscInt *newStart, PetscBool *newReverse)
409: {
411: *newReverse = (origReverse == rotateReverse) ? PETSC_FALSE : PETSC_TRUE;
412: *newStart = rotateReverse ? (coneSize + rotateStart - origStart) : (coneSize + origStart - rotateStart);
413: *newStart %= coneSize;
414: return(0);
415: }
417: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_TranslateBack_Private(PetscInt coneSize, PetscInt start, PetscBool reverse, PetscInt *ornt)
418: {
420: if (coneSize < 3) {
421: /* edges just get flipped if start == 1 regardless direction */
422: *ornt = start ? -2 : 0;
423: } else {
424: *ornt = reverse ? -(start+1) : start;
425: }
426: return(0);
427: }
429: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Permute_Private(PetscInt n, const PetscInt arr[], PetscInt start, PetscBool reverse, PetscInt newarr[])
430: {
431: PetscInt i;
434: if (reverse) {for (i=0; i<n; i++) newarr[i] = arr[(n+start-i)%n];}
435: else {for (i=0; i<n; i++) newarr[i] = arr[(start+i)%n];}
436: return(0);
437: }
439: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM,PetscInt,PetscInt,PetscDualSpace *);
440: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,const PetscInt[],PetscInt[]);
441: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection,PetscInt,PetscInt,PetscInt[],PetscBool,const PetscInt***,PetscInt,PetscInt[]);
443: PETSC_EXTERN PetscErrorCode DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
444: PETSC_EXTERN PetscErrorCode DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
445: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
446: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
447: PETSC_INTERN PetscErrorCode DMCreateSubDomainDM_Plex(DM,DMLabel,PetscInt,IS*,DM*);
449: #endif /* _PLEXIMPL_H */