1 #pragma once
2
3 #include <petscmat.h> /*I "petscmat.h" I*/
4 #include <petscdmplex.h> /*I "petscdmplex.h" I*/
5 #include <petscdmplextransform.h>
6 #include <petscbt.h>
7 #include <petscsf.h>
8 #include <petsc/private/dmimpl.h>
9
10 #if defined(PETSC_HAVE_EXODUSII)
11 #include <exodusII.h>
12 #endif
13
14 PETSC_INTERN PetscBool Plexcite;
15 PETSC_INTERN const char PlexCitation[];
16
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;
74
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;
81
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 };
93
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;
111
112 /* Point Numbering in Plex:
113
114 Points are numbered contiguously by stratum. Strate are organized as follows:
115
116 First Stratum: Cells [height 0]
117 Second Stratum: Vertices [depth 0]
118 Third Stratum: Faces [height 1]
119 Fourth Stratum: Edges [depth 1]
120
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;
126
127 PetscSection coneSection; /* Layout of cones (inedges for DAG) */
128 PetscInt *cones; /* Cone for each point */
129 PetscInt *coneOrientations; /* Orientation of each cone point, means cone traversal should start on point 'o', and if negative start on -(o+1) and go in reverse */
130 PetscSection supportSection; /* Layout of cones (inedges for DAG) */
131 PetscInt *supports; /* Cone for each point */
132
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;
136
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 */
143
144 /* Interpolation */
145 DMPlexInterpolatedFlag interpolated;
146 DMPlexInterpolatedFlag interpolatedCollective;
147 PetscBool interpolatePreferTensor; // When different orderings exist, prefer the tensor order
148
149 /* Ordering */
150 DMReorderDefaultFlag reorderDefault; /* Reorder the DM by default */
151
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 */
163
164 /* Hierarchy */
165 PetscBool regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */
166
167 /* Generation */
168 char *tetgenOpts;
169 char *triangleOpts;
170 PetscPartitioner partitioner;
171 PetscBool partitionBalance; /* Evenly divide partition overlap when distributing */
172 PetscBool remeshBd;
173
174 /* Submesh */
175 DMLabel subpointMap; /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
176 IS subpointIS; /* IS holding point number in the enclosing mesh of every point in the submesh chart */
177 PetscObjectState subpointState; /* The state of subpointMap when the subpointIS was last created */
178
179 /* Labels and numbering */
180 PetscObjectState depthState; /* State of depth label, so that we can determine if a user changes it */
181 PetscObjectState celltypeState; /* State of celltype label, so that we can determine if a user changes it */
182 IS globalVertexNumbers;
183 IS globalCellNumbers;
184
185 /* Constraints */
186 PetscSection anchorSection; /* maps constrained points to anchor points */
187 IS anchorIS; /* anchors indexed by the above section */
188 PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
189 PetscErrorCode (*computeanchormatrix)(DM, PetscSection, PetscSection, Mat);
190
191 /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
192 PetscSection parentSection; /* dof == 1 if point has parent */
193 PetscInt *parents; /* point to parent */
194 PetscInt *childIDs; /* point to child ID */
195 PetscSection childSection; /* inverse of parent section */
196 PetscInt *children; /* point to children */
197 DM referenceTree; /* reference tree to which child ID's refer */
198 PetscErrorCode (*getchildsymmetry)(DM, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt *, PetscInt *);
199
200 /* MATIS support */
201 PetscSection subdomainSection;
202
203 /* Adjacency */
204 PetscBool useAnchors; /* Replace constrained points with their anchors in adjacency lists */
205 PetscErrorCode (*useradjacency)(DM, PetscInt, PetscInt *, PetscInt[], void *); /* User callback for adjacency */
206 void *useradjacencyctx; /* User context for callback */
207
208 // Periodicity
209 struct {
210 // Specified by the user
211 PetscInt num_face_sfs; // number of face_sfs
212 PetscSF *face_sfs; // root(donor faces) <-- leaf(local faces)
213 PetscScalar (*transform)[4][4]; // geometric transform
214 // Created eagerly (depends on points)
215 PetscSF composed_sf; // root(non-periodic global points) <-- leaf(local points)
216 IS *periodic_points;
217 } periodic;
218
219 /* Projection */
220 PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
221 PetscInt activePoint; /* current active point in iteration */
222
223 /* Output */
224 PetscInt vtkCellHeight; /* The height of cells for output, default is 0 */
225 PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */
226
227 /* Geometry */
228 PetscReal minradius; /* Minimum distance from cell centroid to face */
229 PetscBool useHashLocation; /* Use grid hashing for point location */
230 PetscGridHash lbox; /* Local box for searching */
231 PetscPointFn *coordFunc; /* Function used to remap newly introduced vertices */
232
233 /* Neighbors */
234 PetscMPIInt *neighbors;
235
236 /* Metric */
237 DMPlexMetricCtx *metricCtx;
238
239 /* FEM */
240 PetscBool useCeed; /* This should convert to a registration system when there are more FEM backends */
241 PetscBool useMatClPerm; /* Use the closure permutation when assembling matrices */
242
243 /* CAD */
244 PetscBool ignoreModel; /* If TRUE, Plex refinement will skip Snap-To-Geometry feature ignoring attached CAD geometry information */
245
246 // Transforms
247 PetscBool saveTransform; // Flag to save transform, if this mesh was produced by a transform
248 DMPlexTransform transform; // Transform producing this mesh (note that this will hold on to the original mesh too)
249
250 /* Debugging */
251 PetscBool printSetValues;
252 PetscInt printAdj;
253 PetscInt printFEM;
254 PetscInt printFVM;
255 PetscInt printL2;
256 PetscInt printLocate;
257 PetscInt printProject;
258 PetscReal printTol;
259 } DM_Plex;
260
261 PETSC_INTERN PetscErrorCode DMPlexCopy_Internal(DM, PetscBool, PetscBool, DM);
262 PETSC_INTERN PetscErrorCode DMPlexReplace_Internal(DM, DM *);
263 PETSC_INTERN PetscErrorCode DMPlexCopyEGADSInfo_Internal(DM, DM);
264
265 PETSC_INTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM, PetscViewer);
266 PETSC_INTERN PetscErrorCode VecView_Plex_Local(Vec, PetscViewer);
267 PETSC_INTERN PetscErrorCode VecView_Plex_Native(Vec, PetscViewer);
268 PETSC_INTERN PetscErrorCode VecView_Plex(Vec, PetscViewer);
269 PETSC_INTERN PetscErrorCode VecLoad_Plex_Local(Vec, PetscViewer);
270 PETSC_INTERN PetscErrorCode VecLoad_Plex_Native(Vec, PetscViewer);
271 PETSC_INTERN PetscErrorCode VecLoad_Plex(Vec, PetscViewer);
272 PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
273 PETSC_INTERN PetscErrorCode DMPlexGetFieldTypes_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscInt **, PetscViewerVTKFieldType **);
274 PETSC_INTERN PetscErrorCode DMPlexRestoreFieldTypes_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscInt **, PetscViewerVTKFieldType **);
275 PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM, PetscViewer);
276 PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject, PetscViewer);
277 #if defined(PETSC_HAVE_HDF5)
278 PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
279 PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
280 PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
281 PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
282 PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
283 PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
284 PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
285 PETSC_INTERN PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM, IS, PetscViewer);
286 PETSC_INTERN PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM, PetscViewer);
287 PETSC_INTERN PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM, IS, PetscViewer);
288 PETSC_INTERN PetscErrorCode DMPlexSectionView_HDF5_Internal(DM, PetscViewer, DM);
289 PETSC_INTERN PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
290 PETSC_INTERN PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
291 PETSC_INTERN PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM, PetscViewer, PetscSF *);
292 PETSC_INTERN PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM, PetscViewer, PetscSF);
293 PETSC_INTERN PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM, PetscViewer, PetscSF);
294 PETSC_INTERN PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, PetscSF *, PetscSF *);
295 PETSC_INTERN PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, Vec);
296 PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
297 PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
298 PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM, PetscViewer);
299 PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
300 PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
301 PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
302 PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
303 PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
304 #endif
305 PETSC_EXTERN PetscErrorCode VecView_Plex_Local_CGNS(Vec, PetscViewer);
306 #if defined(PETSC_HAVE_CGNS)
307 PETSC_EXTERN PetscErrorCode VecLoad_Plex_CGNS_Internal(Vec, PetscViewer);
308 #endif
309
310 PETSC_INTERN PetscErrorCode DMPlexClosurePoints_Private(DM, PetscInt, const PetscInt[], IS *);
311 PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(DM, PetscOptionItems);
312 PETSC_INTERN PetscErrorCode DMSetFromOptions_Overlap_Plex(DM, PetscOptionItems, PetscInt *);
313 PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
314 PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM[]);
315 PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
316 PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM[]);
317 PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, Vec, DMLabel, DMLabel, DM *);
318 PETSC_INTERN PetscErrorCode DMExtrude_Plex(DM, PetscInt, DM *);
319 PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
320 PETSC_INTERN PetscErrorCode DMPlexInsertTimeDerivativeBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
321 PETSC_INTERN PetscErrorCode DMPlexInsertBounds_Plex(DM, PetscBool, PetscReal, Vec);
322 PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
323 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);
324 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);
325 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);
326 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);
327 PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
328 PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, const PetscReal[], PetscReal *);
329 PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
330 PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);
331
332 #if defined(PETSC_HAVE_EXODUSII)
333 PETSC_INTERN PetscErrorCode DMView_PlexExodusII(DM, PetscViewer);
334 PETSC_INTERN PetscErrorCode VecView_PlexExodusII_Internal(Vec, PetscViewer);
335 PETSC_INTERN PetscErrorCode VecLoad_PlexExodusII_Internal(Vec, PetscViewer);
336 #endif
337 PETSC_INTERN PetscErrorCode DMView_PlexCGNS(DM, PetscViewer);
338 PETSC_INTERN PetscErrorCode DMPlexCreateCGNSFromFile_Internal(MPI_Comm, const char *, PetscBool, DM *);
339 PETSC_INTERN PetscErrorCode DMPlexCreateCGNS_Internal_Serial(MPI_Comm, PetscInt, PetscBool, DM *);
340 PETSC_INTERN PetscErrorCode DMPlexCreateCGNS_Internal_Parallel(MPI_Comm, PetscInt, PetscBool, DM *);
341 PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM, PetscInt, PetscInt, PetscInt *);
342 PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM, PetscInt, PetscBool, PetscBool, PetscBool, PetscInt *, PetscInt *[]);
343 PETSC_INTERN PetscErrorCode DMPlexGetMaxAdjacencySize_Internal(DM, PetscBool, PetscInt *);
344 PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
345 PETSC_INTERN PetscErrorCode DMPlexRestoreRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
346 PETSC_INTERN PetscErrorCode DMPlexComputeCellType_Internal(DM, PetscInt, PetscInt, DMPolytopeType *);
347 PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], DMLabel, PetscInt, const PetscScalar[], InsertMode);
348 PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
349 PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
350 PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM, DM, const char *, DM *);
351 PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM, DM, PetscSF, PetscInt *, Mat);
352 PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM, DM, PetscSF, PetscInt *, Mat);
353 PETSC_INTERN PetscErrorCode DMPlexAnchorsModifyMat(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, const PetscScalar[], PetscInt *, PetscInt *, PetscInt *[], PetscScalar *[], PetscInt[], PetscBool);
354 PETSC_INTERN PetscErrorCode DMPlexAnchorsModifyMat_Internal(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, PetscInt, PetscInt, const PetscScalar[], PetscInt *, PetscInt *, PetscInt *[], PetscScalar *[], PetscInt[], PetscBool, PetscBool);
355 PETSC_INTERN PetscErrorCode DMPlexAnchorsGetSubMatModification(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, PetscInt *, PetscInt *, PetscInt *[], PetscInt[], PetscScalar *[]);
356 PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM, PetscInt, const PetscScalar[], PetscInt, PetscInt *);
357 /* this is PETSC_EXTERN just because of src/dm/impls/plex/tests/ex18.c */
358 PETSC_EXTERN PetscErrorCode DMPlexOrientInterface_Internal(DM);
359 PETSC_INTERN PetscErrorCode DMPlexOrientCells_Internal(DM, IS, IS);
360
361 /* Applications may use this function */
362 PETSC_EXTERN PetscErrorCode DMPlexCreateNumbering_Plex(DM, PetscInt, PetscInt, PetscInt, PetscInt *, PetscSF, IS *);
363
364 PETSC_INTERN PetscErrorCode DMPlexInterpolateInPlace_Internal(DM);
365 PETSC_INTERN PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Internal(DM, PetscInt, const PetscInt[], const PetscReal[], const PetscReal[], const DMBoundaryType[], PetscBool);
366 PETSC_INTERN PetscErrorCode DMPlexMigrateIsoperiodicFaceSF_Internal(DM, DM, PetscSF);
367 PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
368 PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, Vec, DMLabel, DMLabel, DM *);
369 PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, Vec, DMLabel, DMLabel, DM *);
370 PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM, Mat *);
371
372 PETSC_INTERN PetscErrorCode DMPlexGetOverlap_Plex(DM, PetscInt *);
373 PETSC_INTERN PetscErrorCode DMPlexSetOverlap_Plex(DM, DM, PetscInt);
374 PETSC_INTERN PetscErrorCode DMPlexDistributeGetDefault_Plex(DM, PetscBool *);
375 PETSC_INTERN PetscErrorCode DMPlexDistributeSetDefault_Plex(DM, PetscBool);
376 PETSC_INTERN PetscErrorCode DMPlexReorderGetDefault_Plex(DM, DMReorderDefaultFlag *);
377 PETSC_INTERN PetscErrorCode DMPlexReorderSetDefault_Plex(DM, DMReorderDefaultFlag);
378 PETSC_INTERN PetscErrorCode DMPlexGetUseCeed_Plex(DM, PetscBool *);
379 PETSC_INTERN PetscErrorCode DMPlexSetUseCeed_Plex(DM, PetscBool);
380 PETSC_INTERN PetscErrorCode DMReorderSectionGetDefault_Plex(DM, DMReorderDefaultFlag *);
381 PETSC_INTERN PetscErrorCode DMReorderSectionSetDefault_Plex(DM, DMReorderDefaultFlag);
382 PETSC_INTERN PetscErrorCode DMReorderSectionGetType_Plex(DM, MatOrderingType *);
383 PETSC_INTERN PetscErrorCode DMReorderSectionSetType_Plex(DM, MatOrderingType);
384
385 #if 1
DihedralInvert(PetscInt N,PetscInt a)386 static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
387 {
388 return (a <= 0) ? a : (N - a);
389 }
390
DihedralCompose(PetscInt N,PetscInt a,PetscInt b)391 static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
392 {
393 if (!N) return 0;
394 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));
395 }
396
DihedralSwap(PetscInt N,PetscInt a,PetscInt b)397 static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
398 {
399 return DihedralCompose(N, DihedralInvert(N, a), b);
400 }
401 #else
402 /* TODO
403 This is a reimplementation of the tensor dihedral symmetries using the new orientations.
404 These should be turned on when we convert to new-style orientations in p4est.
405 */
406 /* invert dihedral symmetry: return a^-1,
407 * using the representation described in
408 * DMPlexGetConeOrientation() */
DihedralInvert(PetscInt N,PetscInt a)409 static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
410 {
411 switch (N) {
412 case 0:
413 return 0;
414 case 2:
415 return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, 0, a);
416 case 4:
417 return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, 0, a);
418 case 8:
419 return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, 0, a);
420 default:
421 SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralInvert()");
422 }
423 return 0;
424 }
425
426 /* compose dihedral symmetry: return b * a,
427 * using the representation described in
428 * DMPlexGetConeOrientation() */
DihedralCompose(PetscInt N,PetscInt a,PetscInt b)429 static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
430 {
431 switch (N) {
432 case 0:
433 return 0;
434 case 2:
435 return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
436 case 4:
437 return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
438 case 8:
439 return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
440 default:
441 SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
442 }
443 return 0;
444 }
445
446 /* swap dihedral symmetries: return b * a^-1,
447 * using the representation described in
448 * DMPlexGetConeOrientation() */
DihedralSwap(PetscInt N,PetscInt a,PetscInt b)449 static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
450 {
451 switch (N) {
452 case 0:
453 return 0;
454 case 2:
455 return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
456 case 4:
457 return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
458 case 8:
459 return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
460 default:
461 SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
462 }
463 return 0;
464 }
465 #endif
466 PETSC_INTERN PetscInt DMPolytopeConvertNewOrientation_Internal(DMPolytopeType, PetscInt);
467 PETSC_INTERN PetscInt DMPolytopeConvertOldOrientation_Internal(DMPolytopeType, PetscInt);
468 PETSC_INTERN PetscErrorCode DMPlexConvertOldOrientations_Internal(DM);
469
470 PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexComputeIntegral_Internal(DM, Vec, PetscInt, PetscInt, PetscScalar *, void *);
471 PETSC_INTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);
472
473 /* Matvec with A in row-major storage, x and y can be aliased */
DMPlex_Mult2D_Internal(const PetscScalar A[],PetscInt ldx,const PetscScalar x[],PetscScalar y[])474 static inline void DMPlex_Mult2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
475 {
476 const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
477 y[0 * ldx] = A[0] * z[0] + A[1] * z[1];
478 y[1 * ldx] = A[2] * z[0] + A[3] * z[1];
479 (void)PetscLogFlops(6.0);
480 }
DMPlex_Mult3D_Internal(const PetscScalar A[],PetscInt ldx,const PetscScalar x[],PetscScalar y[])481 static inline void DMPlex_Mult3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
482 {
483 const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
484 y[0 * ldx] = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
485 y[1 * ldx] = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
486 y[2 * ldx] = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
487 (void)PetscLogFlops(15.0);
488 }
DMPlex_MultTranspose2D_Internal(const PetscScalar A[],PetscInt ldx,const PetscScalar x[],PetscScalar y[])489 static inline void DMPlex_MultTranspose2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
490 {
491 const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
492 y[0 * ldx] = A[0] * z[0] + A[2] * z[1];
493 y[1 * ldx] = A[1] * z[0] + A[3] * z[1];
494 (void)PetscLogFlops(6.0);
495 }
DMPlex_MultTranspose3D_Internal(const PetscScalar A[],PetscInt ldx,const PetscScalar x[],PetscScalar y[])496 static inline void DMPlex_MultTranspose3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
497 {
498 const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
499 y[0 * ldx] = PetscConj(A[0]) * z[0] + PetscConj(A[3]) * z[1] + PetscConj(A[6]) * z[2];
500 y[1 * ldx] = PetscConj(A[1]) * z[0] + PetscConj(A[4]) * z[1] + PetscConj(A[7]) * z[2];
501 y[2 * ldx] = PetscConj(A[2]) * z[0] + PetscConj(A[5]) * z[1] + PetscConj(A[8]) * z[2];
502 (void)PetscLogFlops(15.0);
503 }
DMPlex_Mult2DReal_Internal(const PetscReal A[],PetscInt ldx,const PetscScalar x[],PetscScalar y[])504 static inline void DMPlex_Mult2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
505 {
506 const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
507 y[0 * ldx] = A[0] * z[0] + A[1] * z[1];
508 y[1 * ldx] = A[2] * z[0] + A[3] * z[1];
509 (void)PetscLogFlops(6.0);
510 }
DMPlex_Mult3DReal_Internal(const PetscReal A[],PetscInt ldx,const PetscScalar x[],PetscScalar y[])511 static inline void DMPlex_Mult3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
512 {
513 const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
514 y[0 * ldx] = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
515 y[1 * ldx] = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
516 y[2 * ldx] = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
517 (void)PetscLogFlops(15.0);
518 }
DMPlex_MultAdd2DReal_Internal(const PetscReal A[],PetscInt ldx,const PetscScalar x[],PetscScalar y[])519 static inline void DMPlex_MultAdd2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
520 {
521 const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
522 y[0 * ldx] += A[0] * z[0] + A[1] * z[1];
523 y[1 * ldx] += A[2] * z[0] + A[3] * z[1];
524 (void)PetscLogFlops(6.0);
525 }
DMPlex_MultAdd3DReal_Internal(const PetscReal A[],PetscInt ldx,const PetscScalar x[],PetscScalar y[])526 static inline void DMPlex_MultAdd3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
527 {
528 const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
529 y[0 * ldx] += A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
530 y[1 * ldx] += A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
531 y[2 * ldx] += A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
532 (void)PetscLogFlops(15.0);
533 }
534 /*
535 A: packed, row-major m x n array
536 x: length m array
537 y: length n arra
538 ldx: the stride in x and y
539
540 A[i,j] = A[i * n + j]
541 A^T[j,i] = A[i,j]
542 */
DMPlex_MultTransposeReal_Internal(const PetscReal A[],PetscInt m,PetscInt n,PetscInt ldx,const PetscScalar x[],PetscScalar y[])543 static inline void DMPlex_MultTransposeReal_Internal(const PetscReal A[], PetscInt m, PetscInt n, PetscInt ldx, const PetscScalar x[], PetscScalar y[])
544 {
545 PetscScalar z[3];
546 PetscInt i, j;
547 for (i = 0; i < m; ++i) z[i] = x[i * ldx];
548 for (j = 0; j < n; ++j) {
549 const PetscInt l = j * ldx;
550 y[l] = 0;
551 for (i = 0; i < m; ++i) y[l] += A[i * n + j] * z[i];
552 }
553 (void)PetscLogFlops(2 * m * n);
554 }
DMPlex_MultTranspose2DReal_Internal(const PetscReal A[],PetscInt ldx,const PetscScalar x[],PetscScalar y[])555 static inline void DMPlex_MultTranspose2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
556 {
557 const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
558 y[0 * ldx] = A[0] * z[0] + A[2] * z[1];
559 y[1 * ldx] = A[1] * z[0] + A[3] * z[1];
560 (void)PetscLogFlops(6.0);
561 }
DMPlex_MultTranspose3DReal_Internal(const PetscReal A[],PetscInt ldx,const PetscScalar x[],PetscScalar y[])562 static inline void DMPlex_MultTranspose3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
563 {
564 const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
565 y[0 * ldx] = A[0] * z[0] + A[3] * z[1] + A[6] * z[2];
566 y[1 * ldx] = A[1] * z[0] + A[4] * z[1] + A[7] * z[2];
567 y[2 * ldx] = A[2] * z[0] + A[5] * z[1] + A[8] * z[2];
568 (void)PetscLogFlops(15.0);
569 }
570
DMPlex_MatMult2D_Internal(const PetscScalar A[],PetscInt n,PetscInt ldb,const PetscScalar B[],PetscScalar C[])571 static inline void DMPlex_MatMult2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
572 {
573 #define PLEX_DIM__ 2
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_Mult2D_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 }
DMPlex_MatMult3D_Internal(const PetscScalar A[],PetscInt n,PetscInt ldb,const PetscScalar B[],PetscScalar C[])583 static inline void DMPlex_MatMult3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
584 {
585 #define PLEX_DIM__ 3
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_Mult3D_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 }
DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[],PetscInt n,PetscInt ldb,const PetscScalar B[],PetscScalar C[])595 static inline void DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
596 {
597 #define PLEX_DIM__ 2
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_MultTranspose2D_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 }
DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[],PetscInt n,PetscInt ldb,const PetscScalar B[],PetscScalar C[])607 static inline void DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
608 {
609 #define PLEX_DIM__ 3
610 PetscScalar z[PLEX_DIM__];
611 for (PetscInt j = 0; j < n; ++j) {
612 for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
613 DMPlex_MultTranspose3D_Internal(A, 1, z, z);
614 for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
615 }
616 (void)PetscLogFlops(8.0 * n);
617 #undef PLEX_DIM__
618 }
619
DMPlex_MatMultLeft2D_Internal(const PetscScalar A[],PetscInt m,PetscInt ldb,const PetscScalar B[],PetscScalar C[])620 static inline void DMPlex_MatMultLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
621 {
622 for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
623 (void)PetscLogFlops(8.0 * m);
624 }
DMPlex_MatMultLeft3D_Internal(const PetscScalar A[],PetscInt m,PetscInt ldb,const PetscScalar B[],PetscScalar C[])625 static inline void DMPlex_MatMultLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
626 {
627 for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
628 (void)PetscLogFlops(8.0 * m);
629 }
DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[],PetscInt m,PetscInt ldb,const PetscScalar B[],PetscScalar C[])630 static inline void DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
631 {
632 for (PetscInt j = 0; j < m; ++j) DMPlex_Mult2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
633 (void)PetscLogFlops(8.0 * m);
634 }
DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[],PetscInt m,PetscInt ldb,const PetscScalar B[],PetscScalar C[])635 static inline void DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
636 {
637 for (PetscInt j = 0; j < m; ++j) DMPlex_Mult3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
638 (void)PetscLogFlops(8.0 * m);
639 }
640
DMPlex_PTAP2DReal_Internal(const PetscReal P[],const PetscScalar A[],PetscScalar C[])641 static inline void DMPlex_PTAP2DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
642 {
643 PetscScalar out[4];
644 PetscInt i, j, k, l;
645 for (i = 0; i < 2; ++i) {
646 for (j = 0; j < 2; ++j) {
647 out[i * 2 + j] = 0.;
648 for (k = 0; k < 2; ++k) {
649 for (l = 0; l < 2; ++l) out[i * 2 + j] += P[k * 2 + i] * A[k * 2 + l] * P[l * 2 + j];
650 }
651 }
652 }
653 for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
654 (void)PetscLogFlops(48.0);
655 }
DMPlex_PTAP3DReal_Internal(const PetscReal P[],const PetscScalar A[],PetscScalar C[])656 static inline void DMPlex_PTAP3DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
657 {
658 PetscScalar out[9];
659 PetscInt i, j, k, l;
660 for (i = 0; i < 3; ++i) {
661 for (j = 0; j < 3; ++j) {
662 out[i * 3 + j] = 0.;
663 for (k = 0; k < 3; ++k) {
664 for (l = 0; l < 3; ++l) out[i * 3 + j] += P[k * 3 + i] * A[k * 3 + l] * P[l * 3 + j];
665 }
666 }
667 }
668 for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
669 (void)PetscLogFlops(243.0);
670 }
671 /* TODO Fix for aliasing of A and C */
DMPlex_PTAPReal_Internal(const PetscReal P[],PetscInt m,PetscInt n,const PetscScalar A[],PetscScalar C[])672 static inline void DMPlex_PTAPReal_Internal(const PetscReal P[], PetscInt m, PetscInt n, const PetscScalar A[], PetscScalar C[])
673 {
674 PetscInt i, j, k, l;
675 for (i = 0; i < n; ++i) {
676 for (j = 0; j < n; ++j) {
677 C[i * n + j] = 0.;
678 for (k = 0; k < m; ++k) {
679 for (l = 0; l < m; ++l) C[i * n + j] += P[k * n + i] * A[k * m + l] * P[l * n + j];
680 }
681 }
682 }
683 (void)PetscLogFlops(243.0);
684 }
685
DMPlex_Transpose2D_Internal(PetscScalar A[])686 static inline void DMPlex_Transpose2D_Internal(PetscScalar A[])
687 {
688 PetscScalar tmp;
689 tmp = A[1];
690 A[1] = A[2];
691 A[2] = tmp;
692 }
DMPlex_Transpose3D_Internal(PetscScalar A[])693 static inline void DMPlex_Transpose3D_Internal(PetscScalar A[])
694 {
695 PetscScalar tmp;
696 tmp = A[1];
697 A[1] = A[3];
698 A[3] = tmp;
699 tmp = A[2];
700 A[2] = A[6];
701 A[6] = tmp;
702 tmp = A[5];
703 A[5] = A[7];
704 A[7] = tmp;
705 }
706
DMPlex_Invert2D_Internal(PetscReal invJ[],PetscReal J[],PetscReal detJ)707 static inline void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
708 {
709 // Allow zero volume cells
710 const PetscReal invDet = detJ == 0 ? 1.0 : (PetscReal)1.0 / detJ;
711
712 invJ[0] = invDet * J[3];
713 invJ[1] = -invDet * J[1];
714 invJ[2] = -invDet * J[2];
715 invJ[3] = invDet * J[0];
716 (void)PetscLogFlops(5.0);
717 }
718
DMPlex_Invert3D_Internal(PetscReal invJ[],PetscReal J[],PetscReal detJ)719 static inline void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
720 {
721 // Allow zero volume cells
722 const PetscReal invDet = detJ == 0 ? 1.0 : (PetscReal)1.0 / detJ;
723
724 invJ[0 * 3 + 0] = invDet * (J[1 * 3 + 1] * J[2 * 3 + 2] - J[1 * 3 + 2] * J[2 * 3 + 1]);
725 invJ[0 * 3 + 1] = invDet * (J[0 * 3 + 2] * J[2 * 3 + 1] - J[0 * 3 + 1] * J[2 * 3 + 2]);
726 invJ[0 * 3 + 2] = invDet * (J[0 * 3 + 1] * J[1 * 3 + 2] - J[0 * 3 + 2] * J[1 * 3 + 1]);
727 invJ[1 * 3 + 0] = invDet * (J[1 * 3 + 2] * J[2 * 3 + 0] - J[1 * 3 + 0] * J[2 * 3 + 2]);
728 invJ[1 * 3 + 1] = invDet * (J[0 * 3 + 0] * J[2 * 3 + 2] - J[0 * 3 + 2] * J[2 * 3 + 0]);
729 invJ[1 * 3 + 2] = invDet * (J[0 * 3 + 2] * J[1 * 3 + 0] - J[0 * 3 + 0] * J[1 * 3 + 2]);
730 invJ[2 * 3 + 0] = invDet * (J[1 * 3 + 0] * J[2 * 3 + 1] - J[1 * 3 + 1] * J[2 * 3 + 0]);
731 invJ[2 * 3 + 1] = invDet * (J[0 * 3 + 1] * J[2 * 3 + 0] - J[0 * 3 + 0] * J[2 * 3 + 1]);
732 invJ[2 * 3 + 2] = invDet * (J[0 * 3 + 0] * J[1 * 3 + 1] - J[0 * 3 + 1] * J[1 * 3 + 0]);
733 (void)PetscLogFlops(37.0);
734 }
735
DMPlex_Det2D_Internal(PetscReal * detJ,const PetscReal J[])736 static inline void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
737 {
738 *detJ = J[0] * J[3] - J[1] * J[2];
739 (void)PetscLogFlops(3.0);
740 }
741
DMPlex_Det3D_Internal(PetscReal * detJ,const PetscReal J[])742 static inline void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
743 {
744 *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]));
745 (void)PetscLogFlops(12.0);
746 }
747
DMPlex_Det2D_Scalar_Internal(PetscReal * detJ,const PetscScalar J[])748 static inline void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
749 {
750 *detJ = PetscRealPart(J[0]) * PetscRealPart(J[3]) - PetscRealPart(J[1]) * PetscRealPart(J[2]);
751 (void)PetscLogFlops(3.0);
752 }
753
DMPlex_Det3D_Scalar_Internal(PetscReal * detJ,const PetscScalar J[])754 static inline void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
755 {
756 *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])));
757 (void)PetscLogFlops(12.0);
758 }
759
DMPlex_WaxpyD_Internal(PetscInt dim,PetscReal a,const PetscReal * x,const PetscReal * y,PetscReal * w)760 static inline void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
761 {
762 PetscInt d;
763 for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
764 }
765
DMPlex_DotD_Internal(PetscInt dim,const PetscScalar * x,const PetscReal * y)766 static inline PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y)
767 {
768 PetscReal sum = 0.0;
769 PetscInt d;
770 for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d]) * y[d];
771 return sum;
772 }
773
DMPlex_DotRealD_Internal(PetscInt dim,const PetscReal * x,const PetscReal * y)774 static inline PetscReal DMPlex_DotRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y)
775 {
776 PetscReal sum = 0.0;
777 PetscInt d;
778 for (d = 0; d < dim; ++d) sum += x[d] * y[d];
779 return sum;
780 }
781
DMPlex_NormD_Internal(PetscInt dim,const PetscReal * x)782 static inline PetscReal DMPlex_NormD_Internal(PetscInt dim, const PetscReal *x)
783 {
784 PetscReal sum = 0.0;
785 PetscInt d;
786 for (d = 0; d < dim; ++d) sum += x[d] * x[d];
787 return PetscSqrtReal(sum);
788 }
789
DMPlex_DistD_Internal(PetscInt dim,const PetscScalar * x,const PetscScalar * y)790 static inline PetscReal DMPlex_DistD_Internal(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
791 {
792 PetscReal sum = 0.0;
793 PetscInt d;
794 for (d = 0; d < dim; ++d) sum += PetscRealPart(PetscConj(x[d] - y[d]) * (x[d] - y[d]));
795 return PetscSqrtReal(sum);
796 }
797
DMPlex_DistRealD_Internal(PetscInt dim,const PetscReal * x,const PetscReal * y)798 static inline PetscReal DMPlex_DistRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y)
799 {
800 PetscReal sum = 0.0;
801 PetscInt d;
802 for (d = 0; d < dim; ++d) sum += (x[d] - y[d]) * (x[d] - y[d]);
803 return PetscSqrtReal(sum);
804 }
805
806 PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM, PetscInt, PetscInt, PetscDualSpace *);
807 PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection, PetscBool, PetscInt, PetscInt, PetscInt *, PetscBool, const PetscInt[], const PetscInt[], PetscInt[]);
808 PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection, PetscBool, PetscInt, PetscInt, PetscInt[], PetscBool, const PetscInt ***, PetscInt, const PetscInt[], PetscInt[]);
809 PETSC_INTERN PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]);
810 PETSC_INTERN PetscErrorCode DMPlexMatSetClosure_Internal(DM, PetscSection, PetscSection, PetscBool, Mat, PetscInt, const PetscScalar[], InsertMode);
811
812 PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexGetAllCells_Internal(DM, IS *);
813 PETSC_INTERN PetscErrorCode DMPlexGetAllFaces_Internal(DM, IS *);
814 PETSC_EXTERN PetscErrorCode DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscFEGeomMode, PetscFEGeom **);
815 PETSC_EXTERN PetscErrorCode DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
816 PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
817 PETSC_SINGLE_LIBRARY_VISIBILITY_INTERNAL PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
818 PETSC_INTERN PetscErrorCode DMCreateSubDomainDM_Plex(DM, DMLabel, PetscInt, IS *, DM *);
819 PETSC_INTERN PetscErrorCode DMPlexBasisTransformPoint_Internal(DM, DM, Vec, PetscInt, PetscBool[], PetscBool, PetscScalar *);
820 PETSC_INTERN PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM, DM, Vec, PetscInt, PetscBool, PetscInt, PetscScalar *);
821 PETSC_INTERN PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscReal *, PetscReal *, void *);
822 PETSC_INTERN PetscErrorCode DMPlexBasisTransformApply_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscScalar *, PetscScalar *, void *);
823 PETSC_INTERN PetscErrorCode DMCreateNeumannOverlap_Plex(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **);
824 PETSC_INTERN PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM, PetscInt, PetscInt, DMLabel, PetscBool);
825 PETSC_INTERN PetscErrorCode DMPlexDistributeOverlap_Internal(DM, PetscInt, MPI_Comm, const char *, PetscSF *, DM *);
826
827 PETSC_INTERN PetscErrorCode DMPlexInterpolateFaces_Internal(DM, PetscInt, DM);
828
829 PETSC_INTERN PetscErrorCode DMPlexMarkSubmesh_Interpolated(DM, DMLabel, PetscInt, PetscBool, PetscBool, DMLabel, DM);
830
831 PETSC_INTERN PetscErrorCode DMPeriodicCoordinateSetUp_Internal(DM);
832
833 /* Functions in the vtable */
834 PETSC_INTERN PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
835 PETSC_INTERN PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat);
836 PETSC_INTERN PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
837 PETSC_INTERN PetscErrorCode DMCreateMassMatrixLumped_Plex(DM, Vec *, Vec *);
838 PETSC_INTERN PetscErrorCode DMCreateGradientMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
839 PETSC_INTERN PetscErrorCode DMCreateLocalSection_Plex(DM dm);
840 PETSC_INTERN PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm);
841 PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J);
842 PETSC_INTERN PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
843 PETSC_INTERN PetscErrorCode DMCreateCellCoordinateDM_Plex(DM dm, DM *cdm);
844 PETSC_INTERN PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field);
845 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
846 PETSC_INTERN PetscErrorCode DMSetUp_Plex(DM dm);
847 PETSC_INTERN PetscErrorCode DMDestroy_Plex(DM dm);
848 PETSC_INTERN PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
849 PETSC_INTERN PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
850 PETSC_INTERN PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm);
851 PETSC_INTERN PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm);
852 PETSC_INTERN PetscErrorCode DMCreateDomainDecompositionScatters_Plex(DM, PetscInt, DM *, VecScatter **, VecScatter **, VecScatter **);
853 PETSC_INTERN PetscErrorCode DMCreateDomainDecomposition_Plex(DM, PetscInt *, char ***, IS **, IS **, DM **);
854 PETSC_INTERN PetscErrorCode DMCreateSectionPermutation_Plex(DM dm, IS *permutation, PetscBT *blockStarts);
855
856 // Coordinate mapping functions
857 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[]);
858 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[]);
859 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[]);
860 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[]);
861 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[]);
862 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[]);
863
864 PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADS(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
865
866 // FIXME: DM with Attached CAD Models (STEP, IGES, BRep, EGADS, EGADSlite)
867 PETSC_EXTERN PetscErrorCode DMPlexSnapToGeomModel(DM, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
868
869 // Coordinate <-> Reference mapping functions
870 PETSC_INTERN PetscErrorCode DMPlexCoordinatesToReference_FE(DM, PetscFE, PetscInt, PetscInt, const PetscReal[], PetscReal[], Vec, PetscInt, PetscInt, PetscInt, PetscReal *);
871 PETSC_INTERN PetscErrorCode DMPlexReferenceToCoordinates_FE(DM, PetscFE, PetscInt, PetscInt, const PetscReal[], PetscReal[], Vec, PetscInt, PetscInt);
872