Actual source code: dmpleximpl.h
1: #pragma once
3: #include <petscmat.h>
4: #include <petscdmplex.h>
5: #include <petscdmplextransform.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_INTERN PetscBool Plexcite;
15: PETSC_INTERN const char PlexCitation[];
17: PETSC_EXTERN PetscLogEvent DMPLEX_Interpolate;
18: PETSC_EXTERN PetscLogEvent DMPLEX_Partition;
19: PETSC_EXTERN PetscLogEvent DMPLEX_PartSelf;
20: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelInvert;
21: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelCreateSF;
22: PETSC_EXTERN PetscLogEvent DMPLEX_PartStratSF;
23: PETSC_EXTERN PetscLogEvent DMPLEX_CreatePointSF;
24: PETSC_EXTERN PetscLogEvent DMPLEX_Distribute;
25: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeMultistage;
26: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeCones;
27: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeLabels;
28: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeSF;
29: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeOverlap;
30: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeField;
31: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeData;
32: PETSC_EXTERN PetscLogEvent DMPLEX_Migrate;
33: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolateSF;
34: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalBegin;
35: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalEnd;
36: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalBegin;
37: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalEnd;
38: PETSC_EXTERN PetscLogEvent DMPLEX_Stratify;
39: PETSC_EXTERN PetscLogEvent DMPLEX_Symmetrize;
40: PETSC_EXTERN PetscLogEvent DMPLEX_Preallocate;
41: PETSC_EXTERN PetscLogEvent DMPLEX_ResidualFEM;
42: PETSC_EXTERN PetscLogEvent DMPLEX_JacobianFEM;
43: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolatorFEM;
44: PETSC_EXTERN PetscLogEvent DMPLEX_InjectorFEM;
45: PETSC_EXTERN PetscLogEvent DMPLEX_IntegralFEM;
46: PETSC_EXTERN PetscLogEvent DMPLEX_CreateGmsh;
47: PETSC_EXTERN PetscLogEvent DMPLEX_CreateBoxSFC;
48: PETSC_EXTERN PetscLogEvent DMPLEX_RebalanceSharedPoints;
49: PETSC_EXTERN PetscLogEvent DMPLEX_CreateFromFile;
50: PETSC_EXTERN PetscLogEvent DMPLEX_CreateFromOptions;
51: PETSC_EXTERN PetscLogEvent DMPLEX_BuildFromCellList;
52: PETSC_EXTERN PetscLogEvent DMPLEX_BuildCoordinatesFromCellList;
53: PETSC_EXTERN PetscLogEvent DMPLEX_LocatePoints;
54: PETSC_EXTERN PetscLogEvent DMPLEX_TopologyView;
55: PETSC_EXTERN PetscLogEvent DMPLEX_DistributionView;
56: PETSC_EXTERN PetscLogEvent DMPLEX_LabelsView;
57: PETSC_EXTERN PetscLogEvent DMPLEX_CoordinatesView;
58: PETSC_EXTERN PetscLogEvent DMPLEX_SectionView;
59: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalVectorView;
60: PETSC_EXTERN PetscLogEvent DMPLEX_LocalVectorView;
61: PETSC_EXTERN PetscLogEvent DMPLEX_TopologyLoad;
62: PETSC_EXTERN PetscLogEvent DMPLEX_DistributionLoad;
63: PETSC_EXTERN PetscLogEvent DMPLEX_LabelsLoad;
64: PETSC_EXTERN PetscLogEvent DMPLEX_CoordinatesLoad;
65: PETSC_EXTERN PetscLogEvent DMPLEX_SectionLoad;
66: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalVectorLoad;
67: PETSC_EXTERN PetscLogEvent DMPLEX_LocalVectorLoad;
68: PETSC_EXTERN PetscLogEvent DMPLEX_MetricEnforceSPD;
69: PETSC_EXTERN PetscLogEvent DMPLEX_MetricNormalize;
70: PETSC_EXTERN PetscLogEvent DMPLEX_MetricAverage;
71: PETSC_EXTERN PetscLogEvent DMPLEX_MetricIntersection;
72: PETSC_EXTERN PetscLogEvent DMPLEX_Generate;
73: PETSC_EXTERN PetscLogEvent DMPLEX_GetLocalOffsets;
75: PETSC_EXTERN PetscLogEvent DMPLEX_RebalBuildGraph;
76: PETSC_EXTERN PetscLogEvent DMPLEX_RebalRewriteSF;
77: PETSC_EXTERN PetscLogEvent DMPLEX_RebalGatherGraph;
78: PETSC_EXTERN PetscLogEvent DMPLEX_RebalPartition;
79: PETSC_EXTERN PetscLogEvent DMPLEX_RebalScatterPart;
80: PETSC_EXTERN PetscLogEvent DMPLEX_Uninterpolate;
82: struct _n_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 */
133: struct { // DMPolytopeType is an enum (usually size 4), but this needs frequent access
134: uint8_t value_as_uint8; // in a struct to guard for stronger typing
135: } *cellTypes;
137: /* Transformation */
138: DMPlexTransform tr; /* Type of transform used to define an ephemeral mesh */
139: char *transformType; /* Type of transform for uniform cell refinement */
140: PetscBool refinementUniform; /* Flag for uniform cell refinement */
141: PetscReal refinementLimit; /* Maximum volume for refined cell */
142: PetscErrorCode (*refinementFunc)(const PetscReal[], PetscReal *); /* Function giving the maximum volume for refined cell */
144: /* Interpolation */
145: DMPlexInterpolatedFlag interpolated;
146: DMPlexInterpolatedFlag interpolatedCollective;
147: PetscBool interpolatePreferTensor; // When different orderings exist, prefer the tensor order
149: /* Ordering */
150: DMReorderDefaultFlag reorderDefault; /* Reorder the DM by default */
152: /* Distribution */
153: PetscBool distDefault; /* Distribute the DM by default */
154: PetscInt overlap; /* Overlap of the partitions as passed to DMPlexDistribute() or DMPlexDistributeOverlap() */
155: PetscInt numOvLabels; /* The number of labels used for candidate overlap points */
156: DMLabel ovLabels[16]; /* Labels used for candidate overlap points */
157: PetscInt ovValues[16]; /* Label values used for candidate overlap points */
158: PetscInt numOvExLabels; /* The number of labels used for exclusion */
159: DMLabel ovExLabels[16]; /* Labels used to exclude points from the overlap */
160: PetscInt ovExValues[16]; /* Label values used to exclude points from the overlap */
161: char *distributionName; /* Name of the specific parallel distribution of the DM */
162: MPI_Comm nonempty_comm; /* Communicator used for visualization when some processes do not have cells */
164: /* Hierarchy */
165: PetscBool regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */
167: /* Generation */
168: char *tetgenOpts; // COmmand line options for Tetgen
169: PetscReal tetgenRadiusEdgeBound; // Maximum tetgen radius-edge ratio
170: PetscReal tetgenDihedralBound; // Minimum tetgen dihedral angle
171: char *triangleOpts; // Comand line options for Triangle
172: PetscReal triangleAngBound; // Minimum triangle angle
173: PetscPartitioner partitioner; // Partitioner object
174: PetscBool partitionBalance; // Evenly divide partition overlap when distributing
175: PetscBool remeshBd;
177: /* Submesh */
178: DMLabel subpointMap; /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
179: IS subpointIS; /* IS holding point number in the enclosing mesh of every point in the submesh chart */
180: PetscObjectState subpointState; /* The state of subpointMap when the subpointIS was last created */
182: /* Labels and numbering */
183: PetscObjectState depthState; /* State of depth label, so that we can determine if a user changes it */
184: PetscObjectState celltypeState; /* State of celltype label, so that we can determine if a user changes it */
185: IS globalVertexNumbers;
186: IS globalCellNumbers;
188: /* Constraints */
189: PetscSection anchorSection; /* maps constrained points to anchor points */
190: IS anchorIS; /* anchors indexed by the above section */
191: PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
192: PetscErrorCode (*computeanchormatrix)(DM, PetscSection, PetscSection, Mat);
194: /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
195: PetscSection parentSection; /* dof == 1 if point has parent */
196: PetscInt *parents; /* point to parent */
197: PetscInt *childIDs; /* point to child ID */
198: PetscSection childSection; /* inverse of parent section */
199: PetscInt *children; /* point to children */
200: DM referenceTree; /* reference tree to which child ID's refer */
201: PetscErrorCode (*getchildsymmetry)(DM, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt *, PetscInt *);
203: /* MATIS support */
204: PetscSection subdomainSection;
206: /* Adjacency */
207: PetscBool useAnchors; /* Replace constrained points with their anchors in adjacency lists */
208: PetscErrorCode (*useradjacency)(DM, PetscInt, PetscInt *, PetscInt[], void *); /* User callback for adjacency */
209: void *useradjacencyctx; /* User context for callback */
211: // Periodicity
212: struct {
213: // Specified by the user
214: PetscInt num_face_sfs; // number of face_sfs
215: PetscSF *face_sfs; // root(donor faces) <-- leaf(local faces)
216: PetscScalar (*transform)[4][4]; // geometric transform
217: // Created eagerly (depends on points)
218: PetscSF composed_sf; // root(non-periodic global points) <-- leaf(local points)
219: IS *periodic_points;
220: } periodic;
222: /* Projection */
223: PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
224: PetscInt activePoint; /* current active point in iteration */
226: /* Output */
227: PetscInt vtkCellHeight; /* The height of cells for output, default is 0 */
228: PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */
230: /* Geometry */
231: PetscReal minradius; /* Minimum distance from cell centroid to face */
232: PetscBool useHashLocation; /* Use grid hashing for point location */
233: PetscGridHash lbox; /* Local box for searching */
234: PetscPointFn *coordFunc; /* Function used to remap newly introduced vertices */
236: /* Neighbors */
237: PetscMPIInt *neighbors;
239: /* Metric */
240: DMPlexMetricCtx *metricCtx;
242: /* FEM */
243: PetscBool useCeed; /* This should convert to a registration system when there are more FEM backends */
244: PetscBool useMatClPerm; /* Use the closure permutation when assembling matrices */
246: /* CAD */
247: PetscBool ignoreModel; /* If TRUE, Plex refinement will skip Snap-To-Geometry feature ignoring attached CAD geometry information */
249: // Transforms
250: PetscBool saveTransform; // Flag to save transform, if this mesh was produced by a transform
251: DMPlexTransform transform; // Transform producing this mesh (note that this will hold on to the original mesh too)
253: /* Debugging */
254: PetscBool printSetValues;
255: PetscInt printAdj;
256: PetscInt printFEM;
257: PetscInt printFVM;
258: PetscInt printL2;
259: PetscInt printOrient;
260: PetscInt printLocate;
261: PetscInt printProject;
262: PetscInt printCohesive;
263: PetscReal printTol;
264: } DM_Plex;
266: PETSC_INTERN PetscErrorCode DMPlexCopy_Internal(DM, PetscBool, PetscBool, DM);
267: PETSC_INTERN PetscErrorCode DMPlexReplace_Internal(DM, DM *);
268: PETSC_INTERN PetscErrorCode DMPlexCopyEGADSInfo_Internal(DM, DM);
270: PETSC_INTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM, PetscViewer);
271: PETSC_INTERN PetscErrorCode VecView_Plex_Local(Vec, PetscViewer);
272: PETSC_INTERN PetscErrorCode VecView_Plex_Native(Vec, PetscViewer);
273: PETSC_INTERN PetscErrorCode VecView_Plex(Vec, PetscViewer);
274: PETSC_INTERN PetscErrorCode VecLoad_Plex_Local(Vec, PetscViewer);
275: PETSC_INTERN PetscErrorCode VecLoad_Plex_Native(Vec, PetscViewer);
276: PETSC_INTERN PetscErrorCode VecLoad_Plex(Vec, PetscViewer);
277: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
278: PETSC_INTERN PetscErrorCode DMPlexGetFieldTypes_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscInt **, PetscViewerVTKFieldType **);
279: PETSC_INTERN PetscErrorCode DMPlexRestoreFieldTypes_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscInt **, PetscViewerVTKFieldType **);
280: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM, PetscViewer);
281: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject, PetscViewer);
282: #if defined(PETSC_HAVE_HDF5)
283: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
284: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
285: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
286: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
287: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
288: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
289: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
290: PETSC_INTERN PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM, IS, PetscViewer);
291: PETSC_INTERN PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM, PetscViewer);
292: PETSC_INTERN PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM, IS, PetscViewer);
293: PETSC_INTERN PetscErrorCode DMPlexSectionView_HDF5_Internal(DM, PetscViewer, DM);
294: PETSC_INTERN PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
295: PETSC_INTERN PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
296: PETSC_INTERN PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM, PetscViewer, PetscSF *);
297: PETSC_INTERN PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM, PetscViewer, PetscSF);
298: PETSC_INTERN PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM, PetscViewer, PetscSF);
299: PETSC_INTERN PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, PetscSF *, PetscSF *);
300: PETSC_INTERN PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, Vec);
301: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
302: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
303: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM, PetscViewer);
304: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
305: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
306: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
307: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
308: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
309: #endif
310: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_CGNS(Vec, PetscViewer);
311: #if defined(PETSC_HAVE_CGNS)
312: PETSC_EXTERN PetscErrorCode VecLoad_Plex_CGNS_Internal(Vec, PetscViewer);
313: #endif
315: PETSC_INTERN PetscErrorCode DMPlexClosurePoints_Private(DM, PetscInt, const PetscInt[], IS *);
316: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(DM, PetscOptionItems);
317: PETSC_INTERN PetscErrorCode DMSetFromOptions_Overlap_Plex(DM, PetscOptionItems, PetscInt *);
318: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
319: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM[]);
320: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
321: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM[]);
322: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, Vec, DMLabel, DMLabel, DM *);
323: PETSC_INTERN PetscErrorCode DMExtrude_Plex(DM, PetscInt, DM *);
324: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
325: PETSC_INTERN PetscErrorCode DMPlexInsertTimeDerivativeBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
326: PETSC_INTERN PetscErrorCode DMPlexInsertBounds_Plex(DM, PetscBool, PetscReal, Vec);
327: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
328: 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);
329: 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);
330: 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);
331: 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);
332: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
333: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, const PetscReal[], PetscReal *);
334: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
335: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);
337: #if defined(PETSC_HAVE_EXODUSII)
338: PETSC_INTERN PetscErrorCode DMView_PlexExodusII(DM, PetscViewer);
339: PETSC_INTERN PetscErrorCode VecView_PlexExodusII_Internal(Vec, PetscViewer);
340: PETSC_INTERN PetscErrorCode VecLoad_PlexExodusII_Internal(Vec, PetscViewer);
341: #endif
342: PETSC_INTERN PetscErrorCode DMView_PlexCGNS(DM, PetscViewer);
343: PETSC_INTERN PetscErrorCode DMPlexCreateCGNSFromFile_Internal(MPI_Comm, const char *, PetscBool, DM *);
344: PETSC_INTERN PetscErrorCode DMPlexCreateCGNS_Internal_Serial(MPI_Comm, PetscInt, PetscBool, DM *);
345: PETSC_INTERN PetscErrorCode DMPlexCreateCGNS_Internal_Parallel(MPI_Comm, PetscInt, PetscBool, DM *);
346: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM, PetscInt, PetscInt, PetscInt *);
347: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM, PetscInt, PetscBool, PetscBool, PetscBool, PetscInt *, PetscInt *[]);
348: PETSC_INTERN PetscErrorCode DMPlexGetMaxAdjacencySize_Internal(DM, PetscBool, PetscInt *);
349: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
350: PETSC_INTERN PetscErrorCode DMPlexRestoreRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
351: PETSC_INTERN PetscErrorCode DMPlexComputeCellType_Internal(DM, PetscInt, PetscInt, DMPolytopeType *);
352: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], DMLabel, PetscInt, const PetscScalar[], InsertMode);
353: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
354: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
355: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM, DM, const char *, DM *);
356: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM, DM, PetscSF, PetscInt *, Mat);
357: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM, DM, PetscSF, PetscInt *, Mat);
358: PETSC_INTERN PetscErrorCode DMPlexAnchorsModifyMat(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, const PetscScalar[], PetscInt *, PetscInt *, PetscInt *[], PetscScalar *[], PetscInt[], PetscBool);
359: PETSC_INTERN PetscErrorCode DMPlexAnchorsModifyMat_Internal(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, PetscInt, PetscInt, const PetscScalar[], PetscInt *, PetscInt *, PetscInt *[], PetscScalar *[], PetscInt[], PetscBool, PetscBool);
360: PETSC_INTERN PetscErrorCode DMPlexAnchorsGetSubMatModification(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, PetscInt *, PetscInt *, PetscInt *[], PetscInt[], PetscScalar *[]);
361: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM, PetscInt, const PetscScalar[], PetscInt, PetscInt *);
362: /* this is PETSC_EXTERN just because of src/dm/impls/plex/tests/ex18.c */
363: PETSC_EXTERN PetscErrorCode DMPlexOrientInterface_Internal(DM);
364: PETSC_INTERN PetscErrorCode DMPlexOrientCells_Internal(DM, IS, IS);
366: /* Applications may use this function */
367: PETSC_EXTERN PetscErrorCode DMPlexCreateNumbering_Plex(DM, PetscInt, PetscInt, PetscInt, PetscInt *, PetscSF, IS *);
369: PETSC_INTERN PetscErrorCode DMPlexInterpolateInPlace_Internal(DM);
370: PETSC_INTERN PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Internal(DM, PetscInt, const PetscInt[], const PetscReal[], const PetscReal[], const DMBoundaryType[], PetscBool);
371: PETSC_INTERN PetscErrorCode DMPlexMigrateIsoperiodicFaceSF_Internal(DM, DM, PetscSF);
372: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
373: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, Vec, DMLabel, DMLabel, DM *);
374: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, Vec, DMLabel, DMLabel, DM *);
375: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM, Mat *);
377: PETSC_INTERN PetscErrorCode DMPlexGetOverlap_Plex(DM, PetscInt *);
378: PETSC_INTERN PetscErrorCode DMPlexSetOverlap_Plex(DM, DM, PetscInt);
379: PETSC_INTERN PetscErrorCode DMPlexDistributeGetDefault_Plex(DM, PetscBool *);
380: PETSC_INTERN PetscErrorCode DMPlexDistributeSetDefault_Plex(DM, PetscBool);
381: PETSC_INTERN PetscErrorCode DMPlexReorderGetDefault_Plex(DM, DMReorderDefaultFlag *);
382: PETSC_INTERN PetscErrorCode DMPlexReorderSetDefault_Plex(DM, DMReorderDefaultFlag);
383: PETSC_INTERN PetscErrorCode DMPlexGetUseCeed_Plex(DM, PetscBool *);
384: PETSC_INTERN PetscErrorCode DMPlexSetUseCeed_Plex(DM, PetscBool);
385: PETSC_INTERN PetscErrorCode DMReorderSectionGetDefault_Plex(DM, DMReorderDefaultFlag *);
386: PETSC_INTERN PetscErrorCode DMReorderSectionSetDefault_Plex(DM, DMReorderDefaultFlag);
387: PETSC_INTERN PetscErrorCode DMReorderSectionGetType_Plex(DM, MatOrderingType *);
388: PETSC_INTERN PetscErrorCode DMReorderSectionSetType_Plex(DM, MatOrderingType);
390: #if 1
391: static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
392: {
393: return (a <= 0) ? a : (N - a);
394: }
396: static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
397: {
398: if (!N) return 0;
399: 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));
400: }
402: static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
403: {
404: return DihedralCompose(N, DihedralInvert(N, a), b);
405: }
406: #else
407: /* TODO
408: This is a reimplementation of the tensor dihedral symmetries using the new orientations.
409: These should be turned on when we convert to new-style orientations in p4est.
410: */
411: /* invert dihedral symmetry: return a^-1,
412: * using the representation described in
413: * DMPlexGetConeOrientation() */
414: static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
415: {
416: switch (N) {
417: case 0:
418: return 0;
419: case 2:
420: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, 0, a);
421: case 4:
422: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, 0, a);
423: case 8:
424: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, 0, a);
425: default:
426: SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralInvert()");
427: }
428: return 0;
429: }
431: /* compose dihedral symmetry: return b * a,
432: * using the representation described in
433: * DMPlexGetConeOrientation() */
434: static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
435: {
436: switch (N) {
437: case 0:
438: return 0;
439: case 2:
440: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
441: case 4:
442: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
443: case 8:
444: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
445: default:
446: SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
447: }
448: return 0;
449: }
451: /* swap dihedral symmetries: return b * a^-1,
452: * using the representation described in
453: * DMPlexGetConeOrientation() */
454: static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
455: {
456: switch (N) {
457: case 0:
458: return 0;
459: case 2:
460: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
461: case 4:
462: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
463: case 8:
464: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
465: default:
466: SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
467: }
468: return 0;
469: }
470: #endif
471: PETSC_INTERN PetscInt DMPolytopeConvertNewOrientation_Internal(DMPolytopeType, PetscInt);
472: PETSC_INTERN PetscInt DMPolytopeConvertOldOrientation_Internal(DMPolytopeType, PetscInt);
473: PETSC_INTERN PetscErrorCode DMPlexConvertOldOrientations_Internal(DM);
475: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexComputeIntegral_Internal(DM, Vec, PetscInt, PetscInt, PetscScalar *, void *);
476: PETSC_INTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);
478: /* Matvec with A in row-major storage, x and y can be aliased */
479: static inline void DMPlex_Mult2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
480: {
481: const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
482: y[0 * ldx] = A[0] * z[0] + A[1] * z[1];
483: y[1 * ldx] = A[2] * z[0] + A[3] * z[1];
484: (void)PetscLogFlops(6.0);
485: }
486: static inline void DMPlex_Mult3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
487: {
488: const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
489: y[0 * ldx] = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
490: y[1 * ldx] = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
491: y[2 * ldx] = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
492: (void)PetscLogFlops(15.0);
493: }
494: static inline void DMPlex_MultTranspose2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
495: {
496: const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
497: y[0 * ldx] = A[0] * z[0] + A[2] * z[1];
498: y[1 * ldx] = A[1] * z[0] + A[3] * z[1];
499: (void)PetscLogFlops(6.0);
500: }
501: static inline void DMPlex_MultTranspose3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
502: {
503: const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
504: y[0 * ldx] = PetscConj(A[0]) * z[0] + PetscConj(A[3]) * z[1] + PetscConj(A[6]) * z[2];
505: y[1 * ldx] = PetscConj(A[1]) * z[0] + PetscConj(A[4]) * z[1] + PetscConj(A[7]) * z[2];
506: y[2 * ldx] = PetscConj(A[2]) * z[0] + PetscConj(A[5]) * z[1] + PetscConj(A[8]) * z[2];
507: (void)PetscLogFlops(15.0);
508: }
509: static inline void DMPlex_Mult2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
510: {
511: const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
512: y[0 * ldx] = A[0] * z[0] + A[1] * z[1];
513: y[1 * ldx] = A[2] * z[0] + A[3] * z[1];
514: (void)PetscLogFlops(6.0);
515: }
516: static inline void DMPlex_Mult3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
517: {
518: const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
519: y[0 * ldx] = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
520: y[1 * ldx] = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
521: y[2 * ldx] = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
522: (void)PetscLogFlops(15.0);
523: }
524: static inline void DMPlex_MultAdd2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
525: {
526: const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
527: y[0 * ldx] += A[0] * z[0] + A[1] * z[1];
528: y[1 * ldx] += A[2] * z[0] + A[3] * z[1];
529: (void)PetscLogFlops(6.0);
530: }
531: static inline void DMPlex_MultAdd3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
532: {
533: const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
534: y[0 * ldx] += A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
535: y[1 * ldx] += A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
536: y[2 * ldx] += A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
537: (void)PetscLogFlops(15.0);
538: }
539: /*
540: A: packed, row-major m x n array
541: x: length m array
542: y: length n arra
543: ldx: the stride in x and y
545: A[i,j] = A[i * n + j]
546: A^T[j,i] = A[i,j]
547: */
548: static inline void DMPlex_MultTransposeReal_Internal(const PetscReal A[], PetscInt m, PetscInt n, PetscInt ldx, const PetscScalar x[], PetscScalar y[])
549: {
550: PetscScalar z[3];
551: PetscInt i, j;
552: for (i = 0; i < m; ++i) z[i] = x[i * ldx];
553: for (j = 0; j < n; ++j) {
554: const PetscInt l = j * ldx;
555: y[l] = 0;
556: for (i = 0; i < m; ++i) y[l] += A[i * n + j] * z[i];
557: }
558: (void)PetscLogFlops(2 * m * n);
559: }
560: static inline void DMPlex_MultTranspose2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
561: {
562: const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
563: y[0 * ldx] = A[0] * z[0] + A[2] * z[1];
564: y[1 * ldx] = A[1] * z[0] + A[3] * z[1];
565: (void)PetscLogFlops(6.0);
566: }
567: static inline void DMPlex_MultTranspose3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
568: {
569: const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
570: y[0 * ldx] = A[0] * z[0] + A[3] * z[1] + A[6] * z[2];
571: y[1 * ldx] = A[1] * z[0] + A[4] * z[1] + A[7] * z[2];
572: y[2 * ldx] = A[2] * z[0] + A[5] * z[1] + A[8] * z[2];
573: (void)PetscLogFlops(15.0);
574: }
576: static inline void DMPlex_MatMult2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
577: {
578: #define PLEX_DIM__ 2
579: PetscScalar z[PLEX_DIM__];
580: for (PetscInt j = 0; j < n; ++j) {
581: for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
582: DMPlex_Mult2D_Internal(A, 1, z, z);
583: for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
584: }
585: (void)PetscLogFlops(8.0 * n);
586: #undef PLEX_DIM__
587: }
588: static inline void DMPlex_MatMult3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
589: {
590: #define PLEX_DIM__ 3
591: PetscScalar z[PLEX_DIM__];
592: for (PetscInt j = 0; j < n; ++j) {
593: for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
594: DMPlex_Mult3D_Internal(A, 1, z, z);
595: for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
596: }
597: (void)PetscLogFlops(8.0 * n);
598: #undef PLEX_DIM__
599: }
600: static inline void DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
601: {
602: #define PLEX_DIM__ 2
603: PetscScalar z[PLEX_DIM__];
604: for (PetscInt j = 0; j < n; ++j) {
605: for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
606: DMPlex_MultTranspose2D_Internal(A, 1, z, z);
607: for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
608: }
609: (void)PetscLogFlops(8.0 * n);
610: #undef PLEX_DIM__
611: }
612: static inline void DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
613: {
614: #define PLEX_DIM__ 3
615: PetscScalar z[PLEX_DIM__];
616: for (PetscInt j = 0; j < n; ++j) {
617: for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
618: DMPlex_MultTranspose3D_Internal(A, 1, z, z);
619: for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
620: }
621: (void)PetscLogFlops(8.0 * n);
622: #undef PLEX_DIM__
623: }
625: static inline void DMPlex_MatMultLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
626: {
627: for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
628: (void)PetscLogFlops(8.0 * m);
629: }
630: static inline void DMPlex_MatMultLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
631: {
632: for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
633: (void)PetscLogFlops(8.0 * m);
634: }
635: static inline void DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
636: {
637: for (PetscInt j = 0; j < m; ++j) DMPlex_Mult2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
638: (void)PetscLogFlops(8.0 * m);
639: }
640: static inline void DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
641: {
642: for (PetscInt j = 0; j < m; ++j) DMPlex_Mult3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
643: (void)PetscLogFlops(8.0 * m);
644: }
646: static inline void DMPlex_PTAP2DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
647: {
648: PetscScalar out[4];
649: PetscInt i, j, k, l;
650: for (i = 0; i < 2; ++i) {
651: for (j = 0; j < 2; ++j) {
652: out[i * 2 + j] = 0.;
653: for (k = 0; k < 2; ++k) {
654: for (l = 0; l < 2; ++l) out[i * 2 + j] += P[k * 2 + i] * A[k * 2 + l] * P[l * 2 + j];
655: }
656: }
657: }
658: for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
659: (void)PetscLogFlops(48.0);
660: }
661: static inline void DMPlex_PTAP3DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
662: {
663: PetscScalar out[9];
664: PetscInt i, j, k, l;
665: for (i = 0; i < 3; ++i) {
666: for (j = 0; j < 3; ++j) {
667: out[i * 3 + j] = 0.;
668: for (k = 0; k < 3; ++k) {
669: for (l = 0; l < 3; ++l) out[i * 3 + j] += P[k * 3 + i] * A[k * 3 + l] * P[l * 3 + j];
670: }
671: }
672: }
673: for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
674: (void)PetscLogFlops(243.0);
675: }
676: /* TODO Fix for aliasing of A and C */
677: static inline void DMPlex_PTAPReal_Internal(const PetscReal P[], PetscInt m, PetscInt n, const PetscScalar A[], PetscScalar C[])
678: {
679: PetscInt i, j, k, l;
680: for (i = 0; i < n; ++i) {
681: for (j = 0; j < n; ++j) {
682: C[i * n + j] = 0.;
683: for (k = 0; k < m; ++k) {
684: for (l = 0; l < m; ++l) C[i * n + j] += P[k * n + i] * A[k * m + l] * P[l * n + j];
685: }
686: }
687: }
688: (void)PetscLogFlops(243.0);
689: }
691: static inline void DMPlex_Transpose2D_Internal(PetscScalar A[])
692: {
693: PetscScalar tmp;
694: tmp = A[1];
695: A[1] = A[2];
696: A[2] = tmp;
697: }
698: static inline void DMPlex_Transpose3D_Internal(PetscScalar A[])
699: {
700: PetscScalar tmp;
701: tmp = A[1];
702: A[1] = A[3];
703: A[3] = tmp;
704: tmp = A[2];
705: A[2] = A[6];
706: A[6] = tmp;
707: tmp = A[5];
708: A[5] = A[7];
709: A[7] = tmp;
710: }
712: static inline void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
713: {
714: // Allow zero volume cells
715: const PetscReal invDet = detJ == 0 ? 1.0 : (PetscReal)1.0 / detJ;
717: invJ[0] = invDet * J[3];
718: invJ[1] = -invDet * J[1];
719: invJ[2] = -invDet * J[2];
720: invJ[3] = invDet * J[0];
721: (void)PetscLogFlops(5.0);
722: }
724: static inline void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
725: {
726: // Allow zero volume cells
727: const PetscReal invDet = detJ == 0 ? 1.0 : (PetscReal)1.0 / detJ;
729: invJ[0 * 3 + 0] = invDet * (J[1 * 3 + 1] * J[2 * 3 + 2] - J[1 * 3 + 2] * J[2 * 3 + 1]);
730: invJ[0 * 3 + 1] = invDet * (J[0 * 3 + 2] * J[2 * 3 + 1] - J[0 * 3 + 1] * J[2 * 3 + 2]);
731: invJ[0 * 3 + 2] = invDet * (J[0 * 3 + 1] * J[1 * 3 + 2] - J[0 * 3 + 2] * J[1 * 3 + 1]);
732: invJ[1 * 3 + 0] = invDet * (J[1 * 3 + 2] * J[2 * 3 + 0] - J[1 * 3 + 0] * J[2 * 3 + 2]);
733: invJ[1 * 3 + 1] = invDet * (J[0 * 3 + 0] * J[2 * 3 + 2] - J[0 * 3 + 2] * J[2 * 3 + 0]);
734: invJ[1 * 3 + 2] = invDet * (J[0 * 3 + 2] * J[1 * 3 + 0] - J[0 * 3 + 0] * J[1 * 3 + 2]);
735: invJ[2 * 3 + 0] = invDet * (J[1 * 3 + 0] * J[2 * 3 + 1] - J[1 * 3 + 1] * J[2 * 3 + 0]);
736: invJ[2 * 3 + 1] = invDet * (J[0 * 3 + 1] * J[2 * 3 + 0] - J[0 * 3 + 0] * J[2 * 3 + 1]);
737: invJ[2 * 3 + 2] = invDet * (J[0 * 3 + 0] * J[1 * 3 + 1] - J[0 * 3 + 1] * J[1 * 3 + 0]);
738: (void)PetscLogFlops(37.0);
739: }
741: static inline void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
742: {
743: *detJ = J[0] * J[3] - J[1] * J[2];
744: (void)PetscLogFlops(3.0);
745: }
747: static inline void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
748: {
749: *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]));
750: (void)PetscLogFlops(12.0);
751: }
753: static inline void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
754: {
755: *detJ = PetscRealPart(J[0]) * PetscRealPart(J[3]) - PetscRealPart(J[1]) * PetscRealPart(J[2]);
756: (void)PetscLogFlops(3.0);
757: }
759: static inline void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
760: {
761: *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])));
762: (void)PetscLogFlops(12.0);
763: }
765: static inline void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
766: {
767: PetscInt d;
768: for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
769: }
771: static inline PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y)
772: {
773: PetscReal sum = 0.0;
774: PetscInt d;
775: for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d]) * y[d];
776: return sum;
777: }
779: static inline PetscReal DMPlex_DotRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y)
780: {
781: PetscReal sum = 0.0;
782: PetscInt d;
783: for (d = 0; d < dim; ++d) sum += x[d] * y[d];
784: return sum;
785: }
787: static inline PetscReal DMPlex_NormD_Internal(PetscInt dim, const PetscReal *x)
788: {
789: PetscReal sum = 0.0;
790: PetscInt d;
791: for (d = 0; d < dim; ++d) sum += x[d] * x[d];
792: return PetscSqrtReal(sum);
793: }
795: static inline PetscReal DMPlex_DistD_Internal(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
796: {
797: PetscReal sum = 0.0;
798: PetscInt d;
799: for (d = 0; d < dim; ++d) sum += PetscRealPart(PetscConj(x[d] - y[d]) * (x[d] - y[d]));
800: return PetscSqrtReal(sum);
801: }
803: static inline PetscReal DMPlex_DistRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y)
804: {
805: PetscReal sum = 0.0;
806: PetscInt d;
807: for (d = 0; d < dim; ++d) sum += (x[d] - y[d]) * (x[d] - y[d]);
808: return PetscSqrtReal(sum);
809: }
811: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM, PetscInt, PetscInt, PetscDualSpace *);
812: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection, PetscBool, PetscInt, PetscInt, PetscInt *, PetscBool, const PetscInt[], const PetscInt[], PetscInt[]);
813: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection, PetscBool, PetscInt, PetscInt, PetscInt[], PetscBool, const PetscInt ***, PetscInt, const PetscInt[], PetscInt[]);
814: PETSC_INTERN PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]);
815: PETSC_INTERN PetscErrorCode DMPlexMatSetClosure_Internal(DM, PetscSection, PetscSection, PetscBool, Mat, PetscInt, const PetscScalar[], InsertMode);
817: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexGetAllCells_Internal(DM, IS *);
818: PETSC_INTERN PetscErrorCode DMPlexGetAllFaces_Internal(DM, IS *);
819: PETSC_EXTERN PetscErrorCode DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscFEGeomMode, PetscFEGeom **);
820: PETSC_EXTERN PetscErrorCode DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
821: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
822: PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
823: PETSC_INTERN PetscErrorCode DMCreateSubDomainDM_Plex(DM, DMLabel, PetscInt, IS *, DM *);
824: PETSC_INTERN PetscErrorCode DMPlexBasisTransformPoint_Internal(DM, DM, Vec, PetscInt, PetscBool[], PetscBool, PetscScalar *);
825: PETSC_INTERN PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM, DM, Vec, PetscInt, PetscBool, PetscInt, PetscScalar *);
826: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscReal *, PetscReal *, void *);
827: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApply_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscScalar *, PetscScalar *, void *);
828: PETSC_INTERN PetscErrorCode DMCreateNeumannOverlap_Plex(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **);
829: PETSC_INTERN PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM, PetscInt, PetscInt, DMLabel, PetscBool);
830: PETSC_INTERN PetscErrorCode DMPlexDistributeOverlap_Internal(DM, PetscInt, MPI_Comm, const char *, PetscSF *, DM *);
832: PETSC_INTERN PetscErrorCode DMPlexInterpolateFaces_Internal(DM, PetscInt, DM);
834: PETSC_INTERN PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM, DMLabel, PetscInt, PetscBool, PetscBool, DMLabel, DM);
836: PETSC_INTERN PetscErrorCode DMPeriodicCoordinateSetUp_Internal(DM);
838: /* Functions in the vtable */
839: PETSC_INTERN PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
840: PETSC_INTERN PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat);
841: PETSC_INTERN PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
842: PETSC_INTERN PetscErrorCode DMCreateMassMatrixLumped_Plex(DM, Vec *, Vec *);
843: PETSC_INTERN PetscErrorCode DMCreateGradientMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
844: PETSC_INTERN PetscErrorCode DMCreateLocalSection_Plex(DM dm);
845: PETSC_INTERN PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm);
846: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J);
847: PETSC_INTERN PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
848: PETSC_INTERN PetscErrorCode DMCreateCellCoordinateDM_Plex(DM dm, DM *cdm);
849: PETSC_INTERN PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field);
850: PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
851: PETSC_INTERN PetscErrorCode DMSetUp_Plex(DM dm);
852: PETSC_INTERN PetscErrorCode DMDestroy_Plex(DM dm);
853: PETSC_INTERN PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
854: PETSC_INTERN PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
855: PETSC_INTERN PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm);
856: PETSC_INTERN PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm);
857: PETSC_INTERN PetscErrorCode DMCreateDomainDecompositionScatters_Plex(DM, PetscInt, DM *, VecScatter **, VecScatter **, VecScatter **);
858: PETSC_INTERN PetscErrorCode DMCreateDomainDecomposition_Plex(DM, PetscInt *, char ***, IS **, IS **, DM **);
859: PETSC_INTERN PetscErrorCode DMCreateSectionPermutation_Plex(DM dm, IS *permutation, PetscBT *blockStarts);
860: PETSC_INTERN PetscErrorCode DMPlexCellUnsplitVertices_Internal(DM, PetscInt, DMPolytopeType, PetscInt *);
862: // Coordinate mapping functions
863: PETSC_INTERN void coordMap_identity(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[]);
864: PETSC_INTERN void coordMap_rotate(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[]);
865: PETSC_INTERN void coordMap_shear(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[]);
866: PETSC_INTERN void coordMap_flare(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[]);
867: PETSC_INTERN void coordMap_annulus(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[]);
868: PETSC_INTERN void coordMap_shell(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[]);
869: PETSC_INTERN void coordMap_sinusoid(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[]);
870: PETSC_INTERN void coordMap_torus(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[]);
872: PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADS(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
874: // FIXME: DM with Attached CAD Models (STEP, IGES, BRep, EGADS, EGADSlite)
875: PETSC_EXTERN PetscErrorCode DMPlexSnapToGeomModel(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
877: // Coordinate <-> Reference mapping functions
878: PETSC_INTERN PetscErrorCode DMPlexCoordinatesToReference_FE(DM, PetscFE, PetscInt, PetscInt, const PetscReal[], PetscReal[], Vec, PetscInt, PetscInt, PetscInt, PetscReal *);
879: PETSC_INTERN PetscErrorCode DMPlexReferenceToCoordinates_FE(DM, PetscFE, PetscInt, PetscInt, const PetscReal[], PetscReal[], Vec, PetscInt, PetscInt);