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