Actual source code: dmpleximpl.h
1: #ifndef _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: #if defined(PETSC_HAVE_EXODUSII)
11: #include <exodusII.h>
12: #endif
14: PETSC_EXTERN PetscLogEvent DMPLEX_Interpolate;
15: PETSC_EXTERN PetscLogEvent DMPLEX_Partition;
16: PETSC_EXTERN PetscLogEvent DMPLEX_PartSelf;
17: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelInvert;
18: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelCreateSF;
19: PETSC_EXTERN PetscLogEvent DMPLEX_PartStratSF;
20: PETSC_EXTERN PetscLogEvent DMPLEX_CreatePointSF;
21: PETSC_EXTERN PetscLogEvent DMPLEX_Distribute;
22: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeCones;
23: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeLabels;
24: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeSF;
25: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeOverlap;
26: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeField;
27: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeData;
28: PETSC_EXTERN PetscLogEvent DMPLEX_Migrate;
29: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolateSF;
30: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalBegin;
31: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalEnd;
32: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalBegin;
33: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalEnd;
34: PETSC_EXTERN PetscLogEvent DMPLEX_Stratify;
35: PETSC_EXTERN PetscLogEvent DMPLEX_Symmetrize;
36: PETSC_EXTERN PetscLogEvent DMPLEX_Preallocate;
37: PETSC_EXTERN PetscLogEvent DMPLEX_ResidualFEM;
38: PETSC_EXTERN PetscLogEvent DMPLEX_JacobianFEM;
39: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolatorFEM;
40: PETSC_EXTERN PetscLogEvent DMPLEX_InjectorFEM;
41: PETSC_EXTERN PetscLogEvent DMPLEX_IntegralFEM;
42: PETSC_EXTERN PetscLogEvent DMPLEX_CreateGmsh;
43: PETSC_EXTERN PetscLogEvent DMPLEX_RebalanceSharedPoints;
44: PETSC_EXTERN PetscLogEvent DMPLEX_CreateFromFile;
45: PETSC_EXTERN PetscLogEvent DMPLEX_BuildFromCellList;
46: PETSC_EXTERN PetscLogEvent DMPLEX_BuildCoordinatesFromCellList;
47: PETSC_EXTERN PetscLogEvent DMPLEX_LocatePoints;
48: PETSC_EXTERN PetscLogEvent DMPLEX_TopologyView;
49: PETSC_EXTERN PetscLogEvent DMPLEX_LabelsView;
50: PETSC_EXTERN PetscLogEvent DMPLEX_CoordinatesView;
51: PETSC_EXTERN PetscLogEvent DMPLEX_SectionView;
52: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalVectorView;
53: PETSC_EXTERN PetscLogEvent DMPLEX_LocalVectorView;
54: PETSC_EXTERN PetscLogEvent DMPLEX_TopologyLoad;
55: PETSC_EXTERN PetscLogEvent DMPLEX_LabelsLoad;
56: PETSC_EXTERN PetscLogEvent DMPLEX_CoordinatesLoad;
57: PETSC_EXTERN PetscLogEvent DMPLEX_SectionLoad;
58: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalVectorLoad;
59: PETSC_EXTERN PetscLogEvent DMPLEX_LocalVectorLoad;
60: PETSC_EXTERN PetscLogEvent DMPLEX_MetricEnforceSPD;
61: PETSC_EXTERN PetscLogEvent DMPLEX_MetricNormalize;
62: PETSC_EXTERN PetscLogEvent DMPLEX_MetricAverage;
63: PETSC_EXTERN PetscLogEvent DMPLEX_MetricIntersection;
65: PETSC_EXTERN PetscLogEvent DMPLEX_RebalBuildGraph;
66: PETSC_EXTERN PetscLogEvent DMPLEX_RebalRewriteSF;
67: PETSC_EXTERN PetscLogEvent DMPLEX_RebalGatherGraph;
68: PETSC_EXTERN PetscLogEvent DMPLEX_RebalPartition;
69: PETSC_EXTERN PetscLogEvent DMPLEX_RebalScatterPart;
71: /* Utility struct to store the contents of a Fluent file in memory */
72: typedef struct {
73: int index; /* Type of section */
74: unsigned int zoneID;
75: unsigned int first;
76: unsigned int last;
77: int type;
78: int nd; /* Either ND or element-type */
79: void *data;
80: } FluentSection;
82: struct _PetscGridHash {
83: PetscInt dim;
84: PetscReal lower[3]; /* The lower-left corner */
85: PetscReal upper[3]; /* The upper-right corner */
86: PetscReal extent[3]; /* The box size */
87: PetscReal h[3]; /* The subbox size */
88: PetscInt n[3]; /* The number of subboxes */
89: PetscSection cellSection; /* Offsets for cells in each subbox*/
90: IS cells; /* List of cells in each subbox */
91: DMLabel cellsSparse; /* Sparse storage for cell map */
92: };
94: typedef struct {
95: PetscBool isotropic; /* Is the metric isotropic? */
96: PetscBool uniform; /* Is the metric uniform? */
97: PetscBool restrictAnisotropyFirst; /* Should anisotropy or normalization come first? */
98: PetscBool noInsert; /* Should node insertion/deletion be turned off? */
99: PetscBool noSwap; /* Should facet swapping be turned off? */
100: PetscBool noMove; /* Should node movement be turned off? */
101: PetscBool noSurf; /* Should surface modification be turned off? */
102: PetscReal h_min, h_max; /* Minimum/maximum tolerated metric magnitudes */
103: PetscReal a_max; /* Maximum tolerated anisotropy */
104: PetscReal targetComplexity; /* Target metric complexity */
105: PetscReal p; /* Degree for L-p normalization methods */
106: PetscReal gradationFactor; /* Maximum tolerated length ratio for adjacent edges */
107: PetscReal hausdorffNumber; /* Max. distance between piecewise linear representation of boundary and reconstructed ideal boundary */
108: PetscInt numIter; /* Number of ParMmg mesh adaptation iterations */
109: PetscInt verbosity; /* Level of verbosity for remesher (-1 = no output, 10 = maximum) */
110: } DMPlexMetricCtx;
112: /* Point Numbering in Plex:
114: Points are numbered contiguously by stratum. Strate are organized as follows:
116: First Stratum: Cells [height 0]
117: Second Stratum: Vertices [depth 0]
118: Third Stratum: Faces [height 1]
119: Fourth Stratum: Edges [depth 1]
121: We do this so that the numbering of a cell-vertex mesh does not change after interpolation. Within a given stratum,
122: we allow additional segregation of by cell type.
123: */
124: typedef struct {
125: PetscInt refct;
127: PetscSection coneSection; /* Layout of cones (inedges for DAG) */
128: PetscInt *cones; /* Cone for each point */
129: PetscInt *coneOrientations; /* Orientation of each cone point, means cone traversal should start on point 'o', and if negative start on -(o+1) and go in reverse */
130: PetscSection supportSection; /* Layout of cones (inedges for DAG) */
131: PetscInt *supports; /* Cone for each point */
132: PetscInt *facesTmp; /* Work space for faces operation */
134: /* Transformation */
135: PetscBool refinementUniform; /* Flag for uniform cell refinement */
136: char *transformType; /* Type of transform for uniform cell refinement */
137: PetscReal refinementLimit; /* Maximum volume for refined cell */
138: PetscErrorCode (*refinementFunc)(const PetscReal[], PetscReal *); /* Function giving the maximum volume for refined cell */
140: /* Interpolation */
141: DMPlexInterpolatedFlag interpolated;
142: DMPlexInterpolatedFlag interpolatedCollective;
144: DMPlexReorderDefaultFlag reorderDefault; /* Reorder the DM by default */
146: /* Distribution */
147: PetscBool distDefault; /* Distribute the DM by default */
148: PetscInt overlap; /* Overlap of the partitions as passed to DMPlexDistribute() or DMPlexDistributeOverlap() */
149: PetscInt numOvLabels; /* The number of labels used for candidate overlap points */
150: DMLabel ovLabels[16]; /* Labels used for candidate overlap points */
151: PetscInt ovValues[16]; /* Label values used for candidate overlap points */
152: PetscInt numOvExLabels; /* The number of labels used for exclusion */
153: DMLabel ovExLabels[16]; /* Labels used to exclude points from the overlap */
154: PetscInt ovExValues[16]; /* Label values used to exclude points from the overlap */
155: char *distributionName; /* Name of the specific parallel distribution of the DM */
157: /* Hierarchy */
158: PetscBool regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */
160: /* Generation */
161: char *tetgenOpts;
162: char *triangleOpts;
163: PetscPartitioner partitioner;
164: PetscBool partitionBalance; /* Evenly divide partition overlap when distributing */
165: PetscBool remeshBd;
167: /* Submesh */
168: DMLabel subpointMap; /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
169: IS subpointIS; /* IS holding point number in the enclosing mesh of every point in the submesh chart */
170: PetscObjectState subpointState; /* The state of subpointMap when the subpointIS was last created */
172: /* Labels and numbering */
173: PetscObjectState depthState; /* State of depth label, so that we can determine if a user changes it */
174: PetscObjectState celltypeState; /* State of celltype label, so that we can determine if a user changes it */
175: IS globalVertexNumbers;
176: IS globalCellNumbers;
178: /* Constraints */
179: PetscSection anchorSection; /* maps constrained points to anchor points */
180: IS anchorIS; /* anchors indexed by the above section */
181: PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
182: PetscErrorCode (*computeanchormatrix)(DM, PetscSection, PetscSection, Mat);
184: /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
185: PetscSection parentSection; /* dof == 1 if point has parent */
186: PetscInt *parents; /* point to parent */
187: PetscInt *childIDs; /* point to child ID */
188: PetscSection childSection; /* inverse of parent section */
189: PetscInt *children; /* point to children */
190: DM referenceTree; /* reference tree to which child ID's refer */
191: PetscErrorCode (*getchildsymmetry)(DM, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt *, PetscInt *);
193: /* MATIS support */
194: PetscSection subdomainSection;
196: /* Adjacency */
197: PetscBool useAnchors; /* Replace constrained points with their anchors in adjacency lists */
198: PetscErrorCode (*useradjacency)(DM, PetscInt, PetscInt *, PetscInt[], void *); /* User callback for adjacency */
199: void *useradjacencyctx; /* User context for callback */
201: /* Projection */
202: PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
203: PetscInt activePoint; /* current active point in iteration */
205: /* Output */
206: PetscInt vtkCellHeight; /* The height of cells for output, default is 0 */
207: PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */
209: /* Geometry */
210: PetscBool ignoreModel; /* Ignore the geometry model during refinement */
211: PetscReal minradius; /* Minimum distance from cell centroid to face */
212: PetscBool useHashLocation; /* Use grid hashing for point location */
213: PetscGridHash lbox; /* Local box for searching */
214: void (*coordFunc)(PetscInt, PetscInt, PetscInt, /* Function used to remap newly introduced vertices */
215: 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[]);
217: /* Neighbors */
218: PetscMPIInt *neighbors;
220: /* Metric */
221: DMPlexMetricCtx *metricCtx;
223: /* Debugging */
224: PetscBool printSetValues;
225: PetscInt printFEM;
226: PetscInt printL2;
227: PetscInt printLocate;
228: PetscReal printTol;
229: } DM_Plex;
231: PETSC_INTERN PetscErrorCode DMPlexCopy_Internal(DM, PetscBool, PetscBool, DM);
232: PETSC_INTERN PetscErrorCode DMPlexReplace_Internal(DM, DM *);
234: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM, PetscViewer);
235: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec, PetscViewer);
236: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec, PetscViewer);
237: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec, PetscViewer);
238: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec, PetscViewer);
239: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec, PetscViewer);
240: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec, PetscViewer);
241: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
242: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM, PetscViewer);
243: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject, PetscViewer);
244: #if defined(PETSC_HAVE_HDF5)
245: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
246: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
247: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
248: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
249: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
250: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
251: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
252: PETSC_INTERN PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM, IS, PetscViewer);
253: PETSC_INTERN PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM, PetscViewer);
254: PETSC_INTERN PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM, IS, PetscViewer);
255: PETSC_INTERN PetscErrorCode DMPlexSectionView_HDF5_Internal(DM, PetscViewer, DM);
256: PETSC_INTERN PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
257: PETSC_INTERN PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
258: PETSC_INTERN PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM, PetscViewer, PetscSF *);
259: PETSC_INTERN PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM, PetscViewer, PetscSF);
260: PETSC_INTERN PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM, PetscViewer, PetscSF);
261: PETSC_INTERN PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, PetscSF *, PetscSF *);
262: PETSC_INTERN PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, Vec);
263: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
264: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
265: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM, PetscViewer);
266: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
267: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
268: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
269: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
270: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
272: struct _n_DMPlexStorageVersion {
273: int major, minor, subminor;
274: };
275: typedef struct _n_DMPlexStorageVersion *DMPlexStorageVersion;
277: PETSC_EXTERN PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionReading(PetscViewer, DMPlexStorageVersion *);
278: PETSC_EXTERN PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionWriting(PetscViewer, DMPlexStorageVersion *);
279: #endif
280: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_CGNS(Vec, PetscViewer);
282: PETSC_INTERN PetscErrorCode DMPlexVecGetClosureAtDepth_Internal(DM, PetscSection, Vec, PetscInt, PetscInt, PetscInt *, PetscScalar *[]);
283: PETSC_INTERN PetscErrorCode DMPlexClosurePoints_Private(DM, PetscInt, const PetscInt[], IS *);
284: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(DM, PetscOptionItems *);
285: PETSC_INTERN PetscErrorCode DMSetFromOptions_Overlap_Plex(DM, PetscOptionItems *, PetscInt *);
286: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
287: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM[]);
288: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
289: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM[]);
290: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, Vec, DMLabel, DMLabel, DM *);
291: PETSC_INTERN PetscErrorCode DMExtrude_Plex(DM, PetscInt, DM *);
292: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
293: PETSC_INTERN PetscErrorCode DMPlexInsertTimeDerivativeBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
294: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
295: 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);
296: 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);
297: 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);
298: 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);
299: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
300: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, const PetscReal[], PetscReal *);
301: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
302: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);
304: #if defined(PETSC_HAVE_EXODUSII)
305: PETSC_EXTERN PetscErrorCode DMView_PlexExodusII(DM, PetscViewer);
306: PETSC_INTERN PetscErrorCode VecView_PlexExodusII_Internal(Vec, PetscViewer);
307: PETSC_INTERN PetscErrorCode VecLoad_PlexExodusII_Internal(Vec, PetscViewer);
308: PETSC_INTERN PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec, int, int, int);
309: PETSC_INTERN PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec, int, int, int);
310: PETSC_INTERN PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec, int, int, int);
311: PETSC_INTERN PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec, int, int, int);
312: PETSC_INTERN PetscErrorCode EXOGetVarIndex_Internal(int, ex_entity_type, const char[], int *);
313: #endif
314: PETSC_INTERN PetscErrorCode DMView_PlexCGNS(DM, PetscViewer);
315: PETSC_INTERN PetscErrorCode DMPlexCreateCGNSFromFile_Internal(MPI_Comm, const char *, PetscBool, DM *);
316: PETSC_INTERN PetscErrorCode DMPlexCreateCGNS_Internal(MPI_Comm, PetscInt, PetscBool, DM *);
317: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM, PetscInt, PetscInt, PetscInt *);
318: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM, PetscInt, PetscBool, PetscBool, PetscBool, PetscInt *, PetscInt *[]);
319: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
320: PETSC_INTERN PetscErrorCode DMPlexRestoreRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
321: PETSC_INTERN PetscErrorCode DMPlexComputeCellType_Internal(DM, PetscInt, PetscInt, DMPolytopeType *);
322: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], DMLabel, PetscInt, const PetscScalar[], InsertMode);
323: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
324: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
325: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM, DM, const char *, DM *);
326: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM, DM, PetscSF, PetscInt *, Mat);
327: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM, DM, PetscSF, PetscInt *, Mat);
328: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, const PetscScalar[], PetscInt *, PetscInt *, PetscInt *[], PetscScalar *[], PetscInt[], PetscBool);
329: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection, PetscInt, PetscInt, PetscInt *, PetscBool, PetscInt, PetscInt[]);
330: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection, PetscInt, PetscInt, PetscInt[], PetscBool, PetscInt, PetscInt[]);
331: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM, PetscInt, const PetscScalar[], PetscInt, PetscInt *);
332: /* these two are PETSC_EXTERN just because of src/dm/impls/plex/tests/ex18.c */
333: PETSC_EXTERN PetscErrorCode DMPlexOrientInterface_Internal(DM);
335: /* Applications may use this function */
336: PETSC_EXTERN PetscErrorCode DMPlexCreateNumbering_Plex(DM, PetscInt, PetscInt, PetscInt, PetscInt *, PetscSF, IS *);
338: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
339: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
340: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, Vec, DMLabel, DMLabel, DM *);
341: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, Vec, DMLabel, DMLabel, DM *);
342: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM, Mat *);
344: PETSC_INTERN PetscErrorCode DMPlexGetOverlap_Plex(DM, PetscInt *);
345: PETSC_INTERN PetscErrorCode DMPlexSetOverlap_Plex(DM, DM, PetscInt);
346: PETSC_INTERN PetscErrorCode DMPlexDistributeGetDefault_Plex(DM, PetscBool *);
347: PETSC_INTERN PetscErrorCode DMPlexDistributeSetDefault_Plex(DM, PetscBool);
348: PETSC_INTERN PetscErrorCode DMPlexReorderGetDefault_Plex(DM, DMPlexReorderDefaultFlag *);
349: PETSC_INTERN PetscErrorCode DMPlexReorderSetDefault_Plex(DM, DMPlexReorderDefaultFlag);
351: #if 1
352: static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
353: {
354: return (a <= 0) ? a : (N - a);
355: }
357: static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
358: {
359: if (!N) return 0;
360: return (a >= 0) ? ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) : ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
361: }
363: static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
364: {
365: return DihedralCompose(N, DihedralInvert(N, a), b);
366: }
367: #else
368: /* TODO
369: This is a reimplementation of the tensor dihedral symmetries using the new orientations.
370: These should be turned on when we convert to new-style orientations in p4est.
371: */
372: /* invert dihedral symmetry: return a^-1,
373: * using the representation described in
374: * DMPlexGetConeOrientation() */
375: static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
376: {
377: switch (N) {
378: case 0:
379: return 0;
380: case 2:
381: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, 0, a);
382: case 4:
383: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, 0, a);
384: case 8:
385: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, 0, a);
386: default:
387: SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralInvert()");
388: }
389: return 0;
390: }
392: /* compose dihedral symmetry: return b * a,
393: * using the representation described in
394: * DMPlexGetConeOrientation() */
395: static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
396: {
397: switch (N) {
398: case 0:
399: return 0;
400: case 2:
401: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
402: case 4:
403: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
404: case 8:
405: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
406: default:
407: SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
408: }
409: return 0;
410: }
412: /* swap dihedral symmetries: return b * a^-1,
413: * using the representation described in
414: * DMPlexGetConeOrientation() */
415: static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
416: {
417: switch (N) {
418: case 0:
419: return 0;
420: case 2:
421: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
422: case 4:
423: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
424: case 8:
425: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
426: default:
427: SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
428: }
429: return 0;
430: }
431: #endif
432: PETSC_EXTERN PetscInt DMPolytopeConvertNewOrientation_Internal(DMPolytopeType, PetscInt);
433: PETSC_EXTERN PetscInt DMPolytopeConvertOldOrientation_Internal(DMPolytopeType, PetscInt);
434: PETSC_EXTERN PetscErrorCode DMPlexConvertOldOrientations_Internal(DM);
436: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, PetscFormKey, IS, PetscReal, Vec, Vec, PetscReal, Vec, void *);
437: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Hybrid_Internal(DM, PetscFormKey[], IS, PetscReal, Vec, Vec, PetscReal, Vec, void *);
438: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, PetscFormKey, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
439: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Hybrid_Internal(DM, PetscFormKey[], IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
440: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Action_Internal(DM, PetscFormKey, IS, PetscReal, PetscReal, Vec, Vec, Vec, Vec, void *);
441: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);
443: /* Matvec with A in row-major storage, x and y can be aliased */
444: static inline void DMPlex_Mult2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
445: {
446: const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
447: y[0 * ldx] = A[0] * z[0] + A[1] * z[1];
448: y[1 * ldx] = A[2] * z[0] + A[3] * z[1];
449: (void)PetscLogFlops(6.0);
450: }
451: static inline void DMPlex_Mult3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
452: {
453: const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
454: y[0 * ldx] = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
455: y[1 * ldx] = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
456: y[2 * ldx] = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
457: (void)PetscLogFlops(15.0);
458: }
459: static inline void DMPlex_MultTranspose2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
460: {
461: const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
462: y[0 * ldx] = A[0] * z[0] + A[2] * z[1];
463: y[1 * ldx] = A[1] * z[0] + A[3] * z[1];
464: (void)PetscLogFlops(6.0);
465: }
466: static inline void DMPlex_MultTranspose3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
467: {
468: const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
469: y[0 * ldx] = A[0] * z[0] + A[3] * z[1] + A[6] * z[2];
470: y[1 * ldx] = A[1] * z[0] + A[4] * z[1] + A[7] * z[2];
471: y[2 * ldx] = A[2] * z[0] + A[5] * z[1] + A[8] * z[2];
472: (void)PetscLogFlops(15.0);
473: }
474: static inline void DMPlex_Mult2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
475: {
476: const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
477: y[0 * ldx] = A[0] * z[0] + A[1] * z[1];
478: y[1 * ldx] = A[2] * z[0] + A[3] * z[1];
479: (void)PetscLogFlops(6.0);
480: }
481: static inline void DMPlex_Mult3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
482: {
483: const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
484: y[0 * ldx] = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
485: y[1 * ldx] = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
486: y[2 * ldx] = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
487: (void)PetscLogFlops(15.0);
488: }
489: static inline void DMPlex_MultAdd2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
490: {
491: const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
492: y[0 * ldx] += A[0] * z[0] + A[1] * z[1];
493: y[1 * ldx] += A[2] * z[0] + A[3] * z[1];
494: (void)PetscLogFlops(6.0);
495: }
496: static inline void DMPlex_MultAdd3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
497: {
498: const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
499: y[0 * ldx] += A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
500: y[1 * ldx] += A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
501: y[2 * ldx] += A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
502: (void)PetscLogFlops(15.0);
503: }
504: /*
505: A: packed, row-major m x n array
506: x: length m array
507: y: length n arra
508: ldx: the stride in x and y
510: A[i,j] = A[i * n + j]
511: A^T[j,i] = A[i,j]
512: */
513: static inline void DMPlex_MultTransposeReal_Internal(const PetscReal A[], PetscInt m, PetscInt n, PetscInt ldx, const PetscScalar x[], PetscScalar y[])
514: {
515: PetscScalar z[3];
516: PetscInt i, j;
517: for (i = 0; i < m; ++i) z[i] = x[i * ldx];
518: for (j = 0; j < n; ++j) {
519: const PetscInt l = j * ldx;
520: y[l] = 0;
521: for (i = 0; i < m; ++i) y[l] += A[i * n + j] * z[i];
522: }
523: (void)PetscLogFlops(2 * m * n);
524: }
525: static inline void DMPlex_MultTranspose2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
526: {
527: const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
528: y[0 * ldx] = A[0] * z[0] + A[2] * z[1];
529: y[1 * ldx] = A[1] * z[0] + A[3] * z[1];
530: (void)PetscLogFlops(6.0);
531: }
532: static inline void DMPlex_MultTranspose3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
533: {
534: const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
535: y[0 * ldx] = A[0] * z[0] + A[3] * z[1] + A[6] * z[2];
536: y[1 * ldx] = A[1] * z[0] + A[4] * z[1] + A[7] * z[2];
537: y[2 * ldx] = A[2] * z[0] + A[5] * z[1] + A[8] * z[2];
538: (void)PetscLogFlops(15.0);
539: }
541: static inline void DMPlex_MatMult2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
542: {
543: #define PLEX_DIM__ 2
544: PetscScalar z[PLEX_DIM__];
545: for (PetscInt j = 0; j < n; ++j) {
546: for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
547: DMPlex_Mult2D_Internal(A, 1, z, z);
548: for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
549: }
550: (void)PetscLogFlops(8.0 * n);
551: #undef PLEX_DIM__
552: }
553: static inline void DMPlex_MatMult3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
554: {
555: #define PLEX_DIM__ 3
556: PetscScalar z[PLEX_DIM__];
557: for (PetscInt j = 0; j < n; ++j) {
558: for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
559: DMPlex_Mult3D_Internal(A, 1, z, z);
560: for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
561: }
562: (void)PetscLogFlops(8.0 * n);
563: #undef PLEX_DIM__
564: }
565: static inline void DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
566: {
567: #define PLEX_DIM__ 2
568: PetscScalar z[PLEX_DIM__];
569: for (PetscInt j = 0; j < n; ++j) {
570: for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
571: DMPlex_MultTranspose2D_Internal(A, 1, z, z);
572: for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
573: }
574: (void)PetscLogFlops(8.0 * n);
575: #undef PLEX_DIM__
576: }
577: static inline void DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
578: {
579: #define PLEX_DIM__ 3
580: PetscScalar z[PLEX_DIM__];
581: for (PetscInt j = 0; j < n; ++j) {
582: for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
583: DMPlex_MultTranspose3D_Internal(A, 1, z, z);
584: for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
585: }
586: (void)PetscLogFlops(8.0 * n);
587: #undef PLEX_DIM__
588: }
590: static inline void DMPlex_MatMultLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
591: {
592: for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
593: (void)PetscLogFlops(8.0 * m);
594: }
595: static inline void DMPlex_MatMultLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
596: {
597: for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
598: (void)PetscLogFlops(8.0 * m);
599: }
600: static inline void DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
601: {
602: for (PetscInt j = 0; j < m; ++j) DMPlex_Mult2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
603: (void)PetscLogFlops(8.0 * m);
604: }
605: static inline void DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
606: {
607: for (PetscInt j = 0; j < m; ++j) DMPlex_Mult3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
608: (void)PetscLogFlops(8.0 * m);
609: }
611: static inline void DMPlex_PTAP2DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
612: {
613: PetscScalar out[4];
614: PetscInt i, j, k, l;
615: for (i = 0; i < 2; ++i) {
616: for (j = 0; j < 2; ++j) {
617: out[i * 2 + j] = 0.;
618: for (k = 0; k < 2; ++k) {
619: for (l = 0; l < 2; ++l) out[i * 2 + j] += P[k * 2 + i] * A[k * 2 + l] * P[l * 2 + j];
620: }
621: }
622: }
623: for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
624: (void)PetscLogFlops(48.0);
625: }
626: static inline void DMPlex_PTAP3DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
627: {
628: PetscScalar out[9];
629: PetscInt i, j, k, l;
630: for (i = 0; i < 3; ++i) {
631: for (j = 0; j < 3; ++j) {
632: out[i * 3 + j] = 0.;
633: for (k = 0; k < 3; ++k) {
634: for (l = 0; l < 3; ++l) out[i * 3 + j] += P[k * 3 + i] * A[k * 3 + l] * P[l * 3 + j];
635: }
636: }
637: }
638: for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
639: (void)PetscLogFlops(243.0);
640: }
641: /* TODO Fix for aliasing of A and C */
642: static inline void DMPlex_PTAPReal_Internal(const PetscReal P[], PetscInt m, PetscInt n, const PetscScalar A[], PetscScalar C[])
643: {
644: PetscInt i, j, k, l;
645: for (i = 0; i < n; ++i) {
646: for (j = 0; j < n; ++j) {
647: C[i * n + j] = 0.;
648: for (k = 0; k < m; ++k) {
649: for (l = 0; l < m; ++l) C[i * n + j] += P[k * n + i] * A[k * m + l] * P[l * n + j];
650: }
651: }
652: }
653: (void)PetscLogFlops(243.0);
654: }
656: static inline void DMPlex_Transpose2D_Internal(PetscScalar A[])
657: {
658: PetscScalar tmp;
659: tmp = A[1];
660: A[1] = A[2];
661: A[2] = tmp;
662: }
663: static inline void DMPlex_Transpose3D_Internal(PetscScalar A[])
664: {
665: PetscScalar tmp;
666: tmp = A[1];
667: A[1] = A[3];
668: A[3] = tmp;
669: tmp = A[2];
670: A[2] = A[6];
671: A[6] = tmp;
672: tmp = A[5];
673: A[5] = A[7];
674: A[7] = tmp;
675: }
677: static inline void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
678: {
679: const PetscReal invDet = 1.0 / detJ;
681: invJ[0] = invDet * J[3];
682: invJ[1] = -invDet * J[1];
683: invJ[2] = -invDet * J[2];
684: invJ[3] = invDet * J[0];
685: (void)PetscLogFlops(5.0);
686: }
688: static inline void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
689: {
690: const PetscReal invDet = 1.0 / detJ;
692: invJ[0 * 3 + 0] = invDet * (J[1 * 3 + 1] * J[2 * 3 + 2] - J[1 * 3 + 2] * J[2 * 3 + 1]);
693: invJ[0 * 3 + 1] = invDet * (J[0 * 3 + 2] * J[2 * 3 + 1] - J[0 * 3 + 1] * J[2 * 3 + 2]);
694: invJ[0 * 3 + 2] = invDet * (J[0 * 3 + 1] * J[1 * 3 + 2] - J[0 * 3 + 2] * J[1 * 3 + 1]);
695: invJ[1 * 3 + 0] = invDet * (J[1 * 3 + 2] * J[2 * 3 + 0] - J[1 * 3 + 0] * J[2 * 3 + 2]);
696: invJ[1 * 3 + 1] = invDet * (J[0 * 3 + 0] * J[2 * 3 + 2] - J[0 * 3 + 2] * J[2 * 3 + 0]);
697: invJ[1 * 3 + 2] = invDet * (J[0 * 3 + 2] * J[1 * 3 + 0] - J[0 * 3 + 0] * J[1 * 3 + 2]);
698: invJ[2 * 3 + 0] = invDet * (J[1 * 3 + 0] * J[2 * 3 + 1] - J[1 * 3 + 1] * J[2 * 3 + 0]);
699: invJ[2 * 3 + 1] = invDet * (J[0 * 3 + 1] * J[2 * 3 + 0] - J[0 * 3 + 0] * J[2 * 3 + 1]);
700: invJ[2 * 3 + 2] = invDet * (J[0 * 3 + 0] * J[1 * 3 + 1] - J[0 * 3 + 1] * J[1 * 3 + 0]);
701: (void)PetscLogFlops(37.0);
702: }
704: static inline void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
705: {
706: *detJ = J[0] * J[3] - J[1] * J[2];
707: (void)PetscLogFlops(3.0);
708: }
710: static inline void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
711: {
712: *detJ = (J[0 * 3 + 0] * (J[1 * 3 + 1] * J[2 * 3 + 2] - J[1 * 3 + 2] * J[2 * 3 + 1]) + J[0 * 3 + 1] * (J[1 * 3 + 2] * J[2 * 3 + 0] - J[1 * 3 + 0] * J[2 * 3 + 2]) + J[0 * 3 + 2] * (J[1 * 3 + 0] * J[2 * 3 + 1] - J[1 * 3 + 1] * J[2 * 3 + 0]));
713: (void)PetscLogFlops(12.0);
714: }
716: static inline void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
717: {
718: *detJ = PetscRealPart(J[0]) * PetscRealPart(J[3]) - PetscRealPart(J[1]) * PetscRealPart(J[2]);
719: (void)PetscLogFlops(3.0);
720: }
722: static inline void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
723: {
724: *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])) + 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])) + 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])));
725: (void)PetscLogFlops(12.0);
726: }
728: static inline void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
729: {
730: PetscInt d;
731: for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
732: }
734: static inline PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y)
735: {
736: PetscReal sum = 0.0;
737: PetscInt d;
738: for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d]) * y[d];
739: return sum;
740: }
742: static inline PetscReal DMPlex_DotRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y)
743: {
744: PetscReal sum = 0.0;
745: PetscInt d;
746: for (d = 0; d < dim; ++d) sum += x[d] * y[d];
747: return sum;
748: }
750: static inline PetscReal DMPlex_NormD_Internal(PetscInt dim, const PetscReal *x)
751: {
752: PetscReal sum = 0.0;
753: PetscInt d;
754: for (d = 0; d < dim; ++d) sum += x[d] * x[d];
755: return PetscSqrtReal(sum);
756: }
758: static inline PetscReal DMPlex_DistD_Internal(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
759: {
760: PetscReal sum = 0.0;
761: PetscInt d;
762: for (d = 0; d < dim; ++d) sum += PetscRealPart(PetscConj(x[d] - y[d]) * (x[d] - y[d]));
763: return PetscSqrtReal(sum);
764: }
766: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM, PetscInt, PetscInt, PetscDualSpace *);
767: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection, PetscBool, PetscInt, PetscInt, PetscInt *, PetscBool, const PetscInt[], const PetscInt[], PetscInt[]);
768: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection, PetscBool, PetscInt, PetscInt, PetscInt[], PetscBool, const PetscInt ***, PetscInt, const PetscInt[], PetscInt[]);
769: PETSC_INTERN PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]);
771: PETSC_EXTERN PetscErrorCode DMPlexGetAllCells_Internal(DM, IS *);
772: PETSC_EXTERN PetscErrorCode DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
773: PETSC_EXTERN PetscErrorCode DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
774: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
775: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
776: PETSC_INTERN PetscErrorCode DMCreateSubDomainDM_Plex(DM, DMLabel, PetscInt, IS *, DM *);
777: PETSC_INTERN PetscErrorCode DMPlexBasisTransformPoint_Internal(DM, DM, Vec, PetscInt, PetscBool[], PetscBool, PetscScalar *);
778: PETSC_EXTERN PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM, DM, Vec, PetscInt, PetscBool, PetscInt, PetscScalar *);
779: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscReal *, PetscReal *, void *);
780: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApply_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscScalar *, PetscScalar *, void *);
781: PETSC_INTERN PetscErrorCode DMCreateNeumannOverlap_Plex(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **);
783: /* Functions in the vtable */
784: PETSC_INTERN PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
785: PETSC_INTERN PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat);
786: PETSC_INTERN PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
787: PETSC_INTERN PetscErrorCode DMCreateMassMatrixLumped_Plex(DM dmCoarse, Vec *lm);
788: PETSC_INTERN PetscErrorCode DMCreateLocalSection_Plex(DM dm);
789: PETSC_INTERN PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm);
790: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J);
791: PETSC_INTERN PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
792: PETSC_INTERN PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field);
793: PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
794: PETSC_INTERN PetscErrorCode DMSetUp_Plex(DM dm);
795: PETSC_INTERN PetscErrorCode DMDestroy_Plex(DM dm);
796: PETSC_INTERN PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
797: PETSC_INTERN PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
798: PETSC_INTERN PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm);
799: PETSC_INTERN PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm);
801: #endif /* _PLEXIMPL_H */