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