Actual source code: dmpleximpl.h
petsc-master 2020-10-19
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;
40: PETSC_EXTERN PetscLogEvent DMPLEX_CreateFromFile;
41: PETSC_EXTERN PetscLogEvent DMPLEX_BuildFromCellList;
42: PETSC_EXTERN PetscLogEvent DMPLEX_BuildCoordinatesFromCellList;
43: PETSC_EXTERN PetscLogEvent DMPLEX_LocatePoints;
45: typedef struct _DMPlexCellRefinerOps *DMPlexCellRefinerOps;
46: struct _DMPlexCellRefinerOps {
47: PetscErrorCode (*refine)(DMPlexCellRefiner, DMPolytopeType, PetscInt *, DMPolytopeType *[], PetscInt *[], PetscInt *[], PetscInt *[]);
48: PetscErrorCode (*mapsubcells)(DMPlexCellRefiner, DMPolytopeType, PetscInt, DMPolytopeType, PetscInt, PetscInt, PetscInt *, PetscInt *);
49: PetscErrorCode (*getaffinetransforms)(DMPlexCellRefiner, DMPolytopeType, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
50: PetscErrorCode (*getaffinefacetransforms)(DMPlexCellRefiner, DMPolytopeType, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[], PetscReal *[]);
51: PetscErrorCode (*getcellvertices)(DMPlexCellRefiner, DMPolytopeType, PetscInt *, PetscReal *[]);
52: PetscErrorCode (*getsubcellvertices)(DMPlexCellRefiner, DMPolytopeType, DMPolytopeType, PetscInt, PetscInt *, PetscInt *[]);
53: PetscErrorCode (*mapcoords)(DMPlexCellRefiner, DMPolytopeType, DMPolytopeType, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
54: PetscErrorCode (*setup)(DMPlexCellRefiner);
55: PetscErrorCode (*destroy)(DMPlexCellRefiner);
56: };
58: struct _p_DMPlexCellRefiner {
59: PETSCHEADER(struct _DMPlexCellRefinerOps);
60: DM dm; /* The original DM */
61: PetscBool setupcalled;
62: DMPlexCellRefinerType type;
63: PetscInt *ctOrder; /* [i] = ct: An array with cell types in depth order */
64: PetscInt *ctOrderInv; /* [ct] = i: An array with the ordinal numbers for each cell type */
65: PetscInt *ctStart; /* The number for the first cell of each polytope type in the original mesh, indexed by cell type */
66: PetscInt *ctStartNew; /* The number for the first cell of each polytope type in the new mesh, indexed by cell type */
67: PetscInt *offset; /* [ct][ctNew]: The offset in the new point numbering of a point of type ctNew produced from an old point of type ct */
68: PetscFE *coordFE; /* Finite element for each cell type, used for localized coordinate interpolation */
69: PetscFEGeom **refGeom; /* Geometry of the reference cell for each cell type */
70: void *data; /* refiner private data */
71: };
73: /* Utility struct to store the contents of a Fluent file in memory */
74: typedef struct {
75: int index; /* Type of section */
76: unsigned int zoneID;
77: unsigned int first;
78: unsigned int last;
79: int type;
80: int nd; /* Either ND or element-type */
81: void *data;
82: } FluentSection;
84: struct _PetscGridHash {
85: PetscInt dim;
86: PetscReal lower[3]; /* The lower-left corner */
87: PetscReal upper[3]; /* The upper-right corner */
88: PetscReal extent[3]; /* The box size */
89: PetscReal h[3]; /* The subbox size */
90: PetscInt n[3]; /* The number of subboxes */
91: PetscSection cellSection; /* Offsets for cells in each subbox*/
92: IS cells; /* List of cells in each subbox */
93: DMLabel cellsSparse; /* Sparse storage for cell map */
94: };
96: /* Point Numbering in Plex:
98: Points are numbered contiguously by stratum. Strate are organized as follows:
100: First Stratum: Cells [height 0]
101: Second Stratum: Vertices [depth 0]
102: Third Stratum: Faces [height 1]
103: Fourth Stratum: Edges [depth 1]
105: We do this so that the numbering of a cell-vertex mesh does not change after interpolation. Within a given stratum,
106: we allow additional segregation of by cell type.
107: */
108: typedef struct {
109: PetscInt refct;
111: PetscSection coneSection; /* Layout of cones (inedges for DAG) */
112: PetscInt maxConeSize; /* Cached for fast lookup */
113: PetscInt *cones; /* Cone for each point */
114: 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 */
115: PetscSection supportSection; /* Layout of cones (inedges for DAG) */
116: PetscInt maxSupportSize; /* Cached for fast lookup */
117: PetscInt *supports; /* Cone for each point */
118: PetscBool refinementUniform; /* Flag for uniform cell refinement */
119: PetscReal refinementLimit; /* Maximum volume for refined cell */
120: PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *); /* Function giving the maximum volume for refined cell */
121: PetscInt overlap; /* Overlap of the partitions as passed to DMPlexDistribute() or DMPlexDistributeOverlap() */
122: DMPlexInterpolatedFlag interpolated;
123: DMPlexInterpolatedFlag interpolatedCollective;
125: PetscInt *facesTmp; /* Work space for faces operation */
127: /* Hierarchy */
128: DMPlexCellRefinerType cellRefiner; /* Strategy for refining cells */
129: PetscBool regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */
131: /* Generation */
132: char *tetgenOpts;
133: char *triangleOpts;
134: PetscPartitioner partitioner;
135: PetscBool partitionBalance; /* Evenly divide partition overlap when distributing */
136: PetscBool remeshBd;
138: /* Submesh */
139: DMLabel subpointMap; /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
140: IS subpointIS; /* IS holding point number in the enclosing mesh of every point in the submesh chart */
141: PetscObjectState subpointState; /* The state of subpointMap when the subpointIS was last created */
143: /* Labels and numbering */
144: PetscObjectState depthState; /* State of depth label, so that we can determine if a user changes it */
145: PetscObjectState celltypeState; /* State of celltype label, so that we can determine if a user changes it */
146: IS globalVertexNumbers;
147: IS globalCellNumbers;
149: /* Constraints */
150: PetscSection anchorSection; /* maps constrained points to anchor points */
151: IS anchorIS; /* anchors indexed by the above section */
152: PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
153: PetscErrorCode (*computeanchormatrix)(DM,PetscSection,PetscSection,Mat);
155: /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
156: PetscSection parentSection; /* dof == 1 if point has parent */
157: PetscInt *parents; /* point to parent */
158: PetscInt *childIDs; /* point to child ID */
159: PetscSection childSection; /* inverse of parent section */
160: PetscInt *children; /* point to children */
161: DM referenceTree; /* reference tree to which child ID's refer */
162: PetscErrorCode (*getchildsymmetry)(DM,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*);
164: /* MATIS support */
165: PetscSection subdomainSection;
167: /* Adjacency */
168: PetscBool useAnchors; /* Replace constrained points with their anchors in adjacency lists */
169: PetscErrorCode (*useradjacency)(DM,PetscInt,PetscInt*,PetscInt[],void*); /* User callback for adjacency */
170: void *useradjacencyctx; /* User context for callback */
172: /* Projection */
173: PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
174: PetscInt activePoint; /* current active point in iteration */
176: /* Output */
177: PetscInt vtkCellHeight; /* The height of cells for output, default is 0 */
178: PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */
180: /* Geometry */
181: PetscReal minradius; /* Minimum distance from cell centroid to face */
182: PetscBool useHashLocation; /* Use grid hashing for point location */
183: PetscGridHash lbox; /* Local box for searching */
184: void (*coordFunc)(PetscInt, PetscInt, PetscInt, /* Function used to remap newly introduced vertices */
185: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
186: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
187: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);
189: /* Neighbors */
190: PetscMPIInt* neighbors;
192: /* Debugging */
193: PetscBool printSetValues;
194: PetscInt printFEM;
195: PetscInt printL2;
196: PetscReal printTol;
197: } DM_Plex;
199: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
200: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
201: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec,PetscViewer);
202: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
203: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
204: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec,PetscViewer);
205: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
206: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
207: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM,PetscViewer);
208: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject,PetscViewer);
209: #if defined(PETSC_HAVE_HDF5)
210: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
211: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
212: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
213: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
214: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
215: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
216: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
217: #endif
219: PETSC_INTERN PetscErrorCode DMPlexVecGetClosureAtDepth_Internal(DM, PetscSection, Vec, PetscInt, PetscInt, PetscInt *, PetscScalar *[]);
220: PETSC_INTERN PetscErrorCode DMPlexClosurePoints_Private(DM,PetscInt,const PetscInt[],IS*);
221: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *, DM);
222: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
223: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM []);
224: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
225: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM []);
226: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, DMLabel, DM *);
227: PETSC_INTERN PetscErrorCode DMAdaptMetric_Plex(DM, Vec, DMLabel, DM *);
228: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
229: PETSC_INTERN PetscErrorCode DMPlexInsertTimeDerivativeBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
230: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
231: 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);
232: 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);
233: 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);
234: PETSC_INTERN PetscErrorCode DMProjectBdFieldLabelLocal_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[],const PetscReal[],PetscInt,const PetscScalar[],PetscScalar[]),InsertMode,Vec);
235: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
236: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *);
237: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
238: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);
240: PETSC_INTERN PetscErrorCode DMPlexLoadLabels_HDF5_Internal(DM, PetscViewer);
241: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
242: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
243: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM, PetscViewer);
244: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
245: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
246: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
247: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
248: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
249: /* TODO Make these INTERN */
250: PETSC_EXTERN PetscErrorCode DMPlexView_ExodusII_Internal(DM, int, PetscInt);
251: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec, int, int);
252: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec, int, int);
253: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec, int, int);
254: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec, int, int);
255: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM,PetscInt,PetscInt,PetscInt*);
256: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
257: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,DMPolytopeType,const PetscInt[],PetscInt*,const DMPolytopeType*[],const PetscInt*[],const PetscInt*[]);
258: PETSC_INTERN PetscErrorCode DMPlexRestoreRawFaces_Internal(DM,DMPolytopeType,const PetscInt[],PetscInt*,const DMPolytopeType*[],const PetscInt*[],const PetscInt*[]);
259: PETSC_INTERN PetscErrorCode CellRefinerInCellTest_Internal(DMPolytopeType, const PetscReal[], PetscBool *);
260: PETSC_INTERN PetscErrorCode DMPlexComputeCellType_Internal(DM, PetscInt, PetscInt, DMPolytopeType *);
261: PETSC_INTERN PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType, PetscInt *[], PetscInt *[]);
262: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], DMLabel, PetscInt, const PetscScalar[], InsertMode);
263: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
264: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
265: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM,DM,const char *,DM*);
266: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM,DM,PetscSF,PetscInt *,Mat);
267: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM,DM,PetscSF,PetscInt *,Mat);
268: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM,PetscSection,PetscInt,PetscInt,const PetscInt[],const PetscInt ***,const PetscScalar[],PetscInt*,PetscInt*,PetscInt*[],PetscScalar*[],PetscInt[],PetscBool);
269: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt []);
270: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt [],PetscBool,PetscInt,PetscInt []);
271: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM,PetscInt,const PetscScalar [],PetscInt,PetscInt *);
272: /* these two are PETSC_EXTERN just because of src/dm/impls/plex/tests/ex18.c */
273: PETSC_EXTERN PetscErrorCode DMPlexOrientCell_Internal(DM,PetscInt,PetscInt,PetscBool);
274: PETSC_EXTERN PetscErrorCode DMPlexOrientInterface_Internal(DM);
276: /* Applications may use this function */
277: PETSC_EXTERN PetscErrorCode DMPlexCreateNumbering_Plex(DM, PetscInt, PetscInt, PetscInt, PetscInt *, PetscSF, IS *);
279: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
280: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
281: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, DMLabel, DM *);
282: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, DMLabel, DM *);
283: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM, Mat*);
285: PETSC_INTERN PetscErrorCode DMPlexGetOverlap_Plex(DM, PetscInt *);
287: /* invert dihedral symmetry: return a^-1,
288: * using the representation described in
289: * DMPlexGetConeOrientation() */
290: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
291: {
292: return (a <= 0) ? a : (N - a);
293: }
295: /* invert dihedral symmetry: return b * a,
296: * using the representation described in
297: * DMPlexGetConeOrientation() */
298: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
299: {
300: if (!N) return 0;
301: return (a >= 0) ?
302: ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
303: ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
304: }
306: /* swap dihedral symmetries: return b * a^-1,
307: * using the representation described in
308: * DMPlexGetConeOrientation() */
309: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
310: {
311: return DihedralCompose(N,DihedralInvert(N,a),b);
312: }
314: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, IS , PetscReal, Vec, Vec, PetscReal, Vec, void *);
315: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Hybrid_Internal(DM, IS , PetscReal, Vec, Vec, PetscReal, Vec, void *);
316: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
317: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Hybrid_Internal(DM, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
318: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);
320: /* Matvec with A in row-major storage, x and y can be aliased */
321: PETSC_STATIC_INLINE void DMPlex_Mult2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
322: {
323: PetscScalar z[2];
324: z[0] = x[0]; z[1] = x[ldx];
325: y[0] = A[0]*z[0] + A[1]*z[1];
326: y[ldx] = A[2]*z[0] + A[3]*z[1];
327: (void)PetscLogFlops(6.0);
328: }
329: PETSC_STATIC_INLINE void DMPlex_Mult3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
330: {
331: PetscScalar z[3];
332: z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
333: y[0] = A[0]*z[0] + A[1]*z[1] + A[2]*z[2];
334: y[ldx] = A[3]*z[0] + A[4]*z[1] + A[5]*z[2];
335: y[ldx*2] = A[6]*z[0] + A[7]*z[1] + A[8]*z[2];
336: (void)PetscLogFlops(15.0);
337: }
338: PETSC_STATIC_INLINE void DMPlex_MultTranspose2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
339: {
340: PetscScalar z[2];
341: z[0] = x[0]; z[1] = x[ldx];
342: y[0] = A[0]*z[0] + A[2]*z[1];
343: y[ldx] = A[1]*z[0] + A[3]*z[1];
344: (void)PetscLogFlops(6.0);
345: }
346: PETSC_STATIC_INLINE void DMPlex_MultTranspose3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
347: {
348: PetscScalar z[3];
349: z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
350: y[0] = A[0]*z[0] + A[3]*z[1] + A[6]*z[2];
351: y[ldx] = A[1]*z[0] + A[4]*z[1] + A[7]*z[2];
352: y[ldx*2] = A[2]*z[0] + A[5]*z[1] + A[8]*z[2];
353: (void)PetscLogFlops(15.0);
354: }
355: PETSC_STATIC_INLINE void DMPlex_Mult2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
356: {
357: PetscScalar z[2];
358: z[0] = x[0]; z[1] = x[ldx];
359: y[0] = A[0]*z[0] + A[1]*z[1];
360: y[ldx] = A[2]*z[0] + A[3]*z[1];
361: (void)PetscLogFlops(6.0);
362: }
363: PETSC_STATIC_INLINE void DMPlex_Mult3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
364: {
365: PetscScalar z[3];
366: z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
367: y[0] = A[0]*z[0] + A[1]*z[1] + A[2]*z[2];
368: y[ldx] = A[3]*z[0] + A[4]*z[1] + A[5]*z[2];
369: y[ldx*2] = A[6]*z[0] + A[7]*z[1] + A[8]*z[2];
370: (void)PetscLogFlops(15.0);
371: }
372: PETSC_STATIC_INLINE void DMPlex_MultTransposeReal_Internal(const PetscReal A[], PetscInt m, PetscInt n, PetscInt ldx, const PetscScalar x[], PetscScalar y[])
373: {
374: PetscScalar z[3];
375: PetscInt i, j;
376: for (i = 0; i < m; ++i) z[i] = x[i*ldx];
377: for (j = 0; j < n; ++j) {
378: const PetscInt l = j*ldx;
379: y[l] = 0;
380: for (i = 0; i < m; ++i) {
381: y[l] += A[j*n+i]*z[i];
382: }
383: }
384: (void)PetscLogFlops(2*m*n);
385: }
386: PETSC_STATIC_INLINE void DMPlex_MultTranspose2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
387: {
388: PetscScalar z[2];
389: z[0] = x[0]; z[1] = x[ldx];
390: y[0] = A[0]*z[0] + A[2]*z[1];
391: y[ldx] = A[1]*z[0] + A[3]*z[1];
392: (void)PetscLogFlops(6.0);
393: }
394: PETSC_STATIC_INLINE void DMPlex_MultTranspose3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
395: {
396: PetscScalar z[3];
397: z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
398: y[0] = A[0]*z[0] + A[3]*z[1] + A[6]*z[2];
399: y[ldx] = A[1]*z[0] + A[4]*z[1] + A[7]*z[2];
400: y[ldx*2] = A[2]*z[0] + A[5]*z[1] + A[8]*z[2];
401: (void)PetscLogFlops(15.0);
402: }
404: PETSC_STATIC_INLINE void DMPlex_MatMult2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
405: {
406: PetscInt j;
407: for (j = 0; j < n; ++j) {
408: PetscScalar z[2];
409: z[0] = B[0+j]; z[1] = B[1*ldb+j];
410: DMPlex_Mult2D_Internal(A, 1, z, z);
411: C[0+j] = z[0]; C[1*ldb+j] = z[1];
412: }
413: (void)PetscLogFlops(8.0*n);
414: }
415: PETSC_STATIC_INLINE void DMPlex_MatMult3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
416: {
417: PetscInt j;
418: for (j = 0; j < n; ++j) {
419: PetscScalar z[3];
420: z[0] = B[0+j]; z[1] = B[1*ldb+j]; z[2] = B[2*ldb+j];
421: DMPlex_Mult3D_Internal(A, 1, z, z);
422: C[0+j] = z[0]; C[1*ldb+j] = z[1]; C[2*ldb+j] = z[2];
423: }
424: (void)PetscLogFlops(8.0*n);
425: }
426: PETSC_STATIC_INLINE void DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
427: {
428: PetscInt j;
429: for (j = 0; j < n; ++j) {
430: PetscScalar z[2];
431: z[0] = B[0+j]; z[1] = B[1*ldb+j];
432: DMPlex_MultTranspose2D_Internal(A, 1, z, z);
433: C[0+j] = z[0]; C[1*ldb+j] = z[1];
434: }
435: (void)PetscLogFlops(8.0*n);
436: }
437: PETSC_STATIC_INLINE void DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
438: {
439: PetscInt j;
440: for (j = 0; j < n; ++j) {
441: PetscScalar z[3];
442: z[0] = B[0+j]; z[1] = B[1*ldb+j]; z[2] = B[2*ldb+j];
443: DMPlex_MultTranspose3D_Internal(A, 1, z, z);
444: C[0+j] = z[0]; C[1*ldb+j] = z[1]; C[2*ldb+j] = z[2];
445: }
446: (void)PetscLogFlops(8.0*n);
447: }
449: PETSC_STATIC_INLINE void DMPlex_MatMultLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
450: {
451: PetscInt j;
452: for (j = 0; j < m; ++j) {
453: DMPlex_MultTranspose2D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
454: }
455: (void)PetscLogFlops(8.0*m);
456: }
457: PETSC_STATIC_INLINE void DMPlex_MatMultLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
458: {
459: PetscInt j;
460: for (j = 0; j < m; ++j) {
461: DMPlex_MultTranspose3D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
462: }
463: (void)PetscLogFlops(8.0*m);
464: }
465: PETSC_STATIC_INLINE void DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
466: {
467: PetscInt j;
468: for (j = 0; j < m; ++j) {
469: DMPlex_Mult2D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
470: }
471: (void)PetscLogFlops(8.0*m);
472: }
473: PETSC_STATIC_INLINE void DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
474: {
475: PetscInt j;
476: for (j = 0; j < m; ++j) {
477: DMPlex_Mult3D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
478: }
479: (void)PetscLogFlops(8.0*m);
480: }
482: PETSC_STATIC_INLINE void DMPlex_Transpose2D_Internal(PetscScalar A[])
483: {
484: PetscScalar tmp;
485: tmp = A[1]; A[1] = A[2]; A[2] = tmp;
486: }
487: PETSC_STATIC_INLINE void DMPlex_Transpose3D_Internal(PetscScalar A[])
488: {
489: PetscScalar tmp;
490: tmp = A[1]; A[1] = A[3]; A[3] = tmp;
491: tmp = A[2]; A[2] = A[6]; A[6] = tmp;
492: tmp = A[5]; A[5] = A[7]; A[7] = tmp;
493: }
495: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
496: {
497: const PetscReal invDet = 1.0/detJ;
499: invJ[0] = invDet*J[3];
500: invJ[1] = -invDet*J[1];
501: invJ[2] = -invDet*J[2];
502: invJ[3] = invDet*J[0];
503: (void)PetscLogFlops(5.0);
504: }
506: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
507: {
508: const PetscReal invDet = 1.0/detJ;
510: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
511: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
512: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
513: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
514: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
515: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
516: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
517: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
518: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
519: (void)PetscLogFlops(37.0);
520: }
522: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
523: {
524: *detJ = J[0]*J[3] - J[1]*J[2];
525: (void)PetscLogFlops(3.0);
526: }
528: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
529: {
530: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
531: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
532: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
533: (void)PetscLogFlops(12.0);
534: }
536: PETSC_STATIC_INLINE void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
537: {
538: *detJ = PetscRealPart(J[0])*PetscRealPart(J[3]) - PetscRealPart(J[1])*PetscRealPart(J[2]);
539: (void)PetscLogFlops(3.0);
540: }
542: PETSC_STATIC_INLINE void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
543: {
544: *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])) +
545: 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])) +
546: 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])));
547: (void)PetscLogFlops(12.0);
548: }
550: 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];}
552: 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;}
554: 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;}
556: 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);}
558: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Translate_Private(PetscInt ornt, PetscInt *start, PetscBool *reverse)
559: {
561: *reverse = (ornt < 0) ? PETSC_TRUE : PETSC_FALSE;
562: *start = *reverse ? -(ornt+1) : ornt;
563: return(0);
564: }
566: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Combine_Private(PetscInt coneSize, PetscInt origStart, PetscBool origReverse, PetscInt rotateStart, PetscBool rotateReverse, PetscInt *newStart, PetscBool *newReverse)
567: {
569: *newReverse = (origReverse == rotateReverse) ? PETSC_FALSE : PETSC_TRUE;
570: *newStart = rotateReverse ? (coneSize + rotateStart - origStart) : (coneSize + origStart - rotateStart);
571: *newStart %= coneSize;
572: return(0);
573: }
575: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_TranslateBack_Private(PetscInt coneSize, PetscInt start, PetscBool reverse, PetscInt *ornt)
576: {
578: if (coneSize < 3) {
579: /* edges just get flipped if start == 1 regardless direction */
580: *ornt = start ? -2 : 0;
581: } else {
582: *ornt = reverse ? -(start+1) : start;
583: }
584: return(0);
585: }
587: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Permute_Private(PetscInt n, const PetscInt arr[], PetscInt start, PetscBool reverse, PetscInt newarr[])
588: {
589: PetscInt i;
592: if (reverse) {for (i=0; i<n; i++) newarr[i] = arr[(n+start-i)%n];}
593: else {for (i=0; i<n; i++) newarr[i] = arr[(start+i)%n];}
594: return(0);
595: }
597: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM,PetscInt,PetscInt,PetscDualSpace *);
598: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection,PetscBool,PetscInt,PetscInt,PetscInt *,PetscBool,const PetscInt[],const PetscInt[],PetscInt[]);
599: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection,PetscBool,PetscInt,PetscInt,PetscInt[],PetscBool,const PetscInt***,PetscInt,const PetscInt[],PetscInt[]);
600: PETSC_INTERN PetscErrorCode DMPlexGetCompressedClosure(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscSection *, IS *, const PetscInt **);
601: PETSC_INTERN PetscErrorCode DMPlexRestoreCompressedClosure(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscSection *, IS *, const PetscInt **);
603: PETSC_EXTERN PetscErrorCode DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
604: PETSC_EXTERN PetscErrorCode DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
605: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
606: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
607: PETSC_INTERN PetscErrorCode DMCreateSubDomainDM_Plex(DM,DMLabel,PetscInt,IS*,DM*);
608: PETSC_INTERN PetscErrorCode DMPlexBasisTransformPoint_Internal(DM, DM, Vec, PetscInt, PetscBool[], PetscBool, PetscScalar *);
609: PETSC_EXTERN PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM, DM, Vec, PetscInt, PetscBool, PetscInt, PetscScalar *);
610: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscReal *, PetscReal *, void *);
611: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApply_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscScalar *, PetscScalar *, void *);
612: PETSC_INTERN PetscErrorCode DMCreateNeumannOverlap_Plex(DM, IS*, Mat*, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void **);
614: #endif /* _PLEXIMPL_H */