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:     PetscScalar transform[4][4]; // geometric transform
220:     PetscSF     face_sf;         // root(donor faces) <-- leaf(local faces)
221:     // Created eagerly (depends on points)
222:     PetscSF composed_sf; // root(non-periodic global points) <-- leaf(local points)
223:     IS      periodic_points;
224:   } periodic;

226:   /* Projection */
227:   PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
228:   PetscInt activePoint;         /* current active point in iteration */

230:   /* Output */
231:   PetscInt  vtkCellHeight;          /* The height of cells for output, default is 0 */
232:   PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */

234:   /* Geometry */
235:   PetscBool     ignoreModel;                      /* Ignore the geometry model during refinement */
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);

304: struct _n_DMPlexStorageVersion {
305:   int major, minor, subminor;
306: };
307: typedef struct _n_DMPlexStorageVersion *DMPlexStorageVersion;

309: PETSC_EXTERN PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionReading(PetscViewer, DMPlexStorageVersion *);
310: PETSC_EXTERN PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionWriting(PetscViewer, DMPlexStorageVersion *);
311: #endif
312: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_CGNS(Vec, PetscViewer);

314: PETSC_INTERN PetscErrorCode DMPlexVecGetClosure_Internal(DM, PetscSection, PetscBool, Vec, PetscInt, PetscInt *, PetscScalar *[]);
315: PETSC_INTERN PetscErrorCode DMPlexVecGetClosureAtDepth_Internal(DM, PetscSection, Vec, PetscInt, PetscInt, PetscInt *, PetscScalar *[]);
316: PETSC_INTERN PetscErrorCode DMPlexClosurePoints_Private(DM, PetscInt, const PetscInt[], IS *);
317: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(DM, PetscOptionItems *);
318: PETSC_INTERN PetscErrorCode DMSetFromOptions_Overlap_Plex(DM, PetscOptionItems *, PetscInt *);
319: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
320: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM[]);
321: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
322: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM[]);
323: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, Vec, DMLabel, DMLabel, DM *);
324: PETSC_INTERN PetscErrorCode DMExtrude_Plex(DM, PetscInt, DM *);
325: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
326: PETSC_INTERN PetscErrorCode DMPlexInsertTimeDerivativeBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
327: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
328: PETSC_INTERN PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM, PetscReal, DMLabel, PetscInt, const PetscInt[], PetscInt, const PetscInt[], PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
329: PETSC_INTERN PetscErrorCode DMProjectFieldLocal_Plex(DM, PetscReal, Vec, void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode, Vec);
330: PETSC_INTERN PetscErrorCode DMProjectFieldLabelLocal_Plex(DM, PetscReal, DMLabel, PetscInt, const PetscInt[], PetscInt, const PetscInt[], Vec, void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode, Vec);
331: PETSC_INTERN PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM, PetscReal, DMLabel, PetscInt, const PetscInt[], PetscInt, const PetscInt[], Vec, void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode, Vec);
332: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
333: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, const PetscReal[], PetscReal *);
334: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
335: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);

337: #if defined(PETSC_HAVE_EXODUSII)
338: PETSC_INTERN PetscErrorCode DMView_PlexExodusII(DM, PetscViewer);
339: PETSC_INTERN PetscErrorCode VecView_PlexExodusII_Internal(Vec, PetscViewer);
340: PETSC_INTERN PetscErrorCode VecLoad_PlexExodusII_Internal(Vec, PetscViewer);
341: #endif
342: PETSC_INTERN PetscErrorCode DMView_PlexCGNS(DM, PetscViewer);
343: PETSC_INTERN PetscErrorCode DMPlexCreateCGNSFromFile_Internal(MPI_Comm, const char *, PetscBool, DM *);
344: PETSC_INTERN PetscErrorCode DMPlexCreateCGNS_Internal(MPI_Comm, PetscInt, PetscBool, DM *);
345: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM, PetscInt, PetscInt, PetscInt *);
346: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM, PetscInt, PetscBool, PetscBool, PetscBool, PetscInt *, PetscInt *[]);
347: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
348: PETSC_INTERN PetscErrorCode DMPlexRestoreRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
349: PETSC_INTERN PetscErrorCode DMPlexComputeCellType_Internal(DM, PetscInt, PetscInt, DMPolytopeType *);
350: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], DMLabel, PetscInt, const PetscScalar[], InsertMode);
351: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
352: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
353: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM, DM, const char *, DM *);
354: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM, DM, PetscSF, PetscInt *, Mat);
355: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM, DM, PetscSF, PetscInt *, Mat);
356: PETSC_INTERN PetscErrorCode DMPlexAnchorsModifyMat(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, const PetscScalar[], PetscInt *, PetscInt *, PetscInt *[], PetscScalar *[], PetscInt[], PetscBool);
357: PETSC_INTERN PetscErrorCode DMPlexAnchorsModifyMat_Internal(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, PetscInt, PetscInt, const PetscScalar[], PetscInt *, PetscInt *, PetscInt *[], PetscScalar *[], PetscInt[], PetscBool, PetscBool);
358: PETSC_INTERN PetscErrorCode DMPlexAnchorsGetSubMatModification(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, PetscInt *, PetscInt *, PetscInt *[], PetscInt[], PetscScalar *[]);
359: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM, PetscInt, const PetscScalar[], PetscInt, PetscInt *);
360: /* these two are PETSC_EXTERN just because of src/dm/impls/plex/tests/ex18.c */
361: PETSC_EXTERN PetscErrorCode DMPlexOrientInterface_Internal(DM);

363: /* Applications may use this function */
364: PETSC_EXTERN PetscErrorCode DMPlexCreateNumbering_Plex(DM, PetscInt, PetscInt, PetscInt, PetscInt *, PetscSF, IS *);

366: PETSC_INTERN PetscErrorCode DMPlexInterpolateInPlace_Internal(DM);
367: PETSC_INTERN PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Internal(DM, PetscInt, const PetscInt[], const PetscReal[], const PetscReal[], const DMBoundaryType[], PetscBool);
368: PETSC_INTERN PetscErrorCode DMPlexMigrateIsoperiodicFaceSF_Internal(DM, DM, PetscSF);
369: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
370: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
371: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, Vec, DMLabel, DMLabel, DM *);
372: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, Vec, DMLabel, DMLabel, DM *);
373: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM, Mat *);

375: PETSC_INTERN PetscErrorCode DMPlexGetOverlap_Plex(DM, PetscInt *);
376: PETSC_INTERN PetscErrorCode DMPlexSetOverlap_Plex(DM, DM, PetscInt);
377: PETSC_INTERN PetscErrorCode DMPlexDistributeGetDefault_Plex(DM, PetscBool *);
378: PETSC_INTERN PetscErrorCode DMPlexDistributeSetDefault_Plex(DM, PetscBool);
379: PETSC_INTERN PetscErrorCode DMPlexReorderGetDefault_Plex(DM, DMReorderDefaultFlag *);
380: PETSC_INTERN PetscErrorCode DMPlexReorderSetDefault_Plex(DM, DMReorderDefaultFlag);
381: PETSC_INTERN PetscErrorCode DMPlexGetUseCeed_Plex(DM, PetscBool *);
382: PETSC_INTERN PetscErrorCode DMPlexSetUseCeed_Plex(DM, PetscBool);
383: PETSC_INTERN PetscErrorCode DMReorderSectionGetDefault_Plex(DM, DMReorderDefaultFlag *);
384: PETSC_INTERN PetscErrorCode DMReorderSectionSetDefault_Plex(DM, DMReorderDefaultFlag);
385: PETSC_INTERN PetscErrorCode DMReorderSectionGetType_Plex(DM, MatOrderingType *);
386: PETSC_INTERN PetscErrorCode DMReorderSectionSetType_Plex(DM, MatOrderingType);

388: #if 1
389: static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
390: {
391:   return (a <= 0) ? a : (N - a);
392: }

394: static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
395: {
396:   if (!N) return 0;
397:   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));
398: }

400: static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
401: {
402:   return DihedralCompose(N, DihedralInvert(N, a), b);
403: }
404: #else
405: /* TODO
406:    This is a reimplementation of the tensor dihedral symmetries using the new orientations.
407:    These should be turned on when we convert to new-style orientations in p4est.
408: */
409: /* invert dihedral symmetry: return a^-1,
410:  * using the representation described in
411:  * DMPlexGetConeOrientation() */
412: static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
413: {
414:   switch (N) {
415:   case 0:
416:     return 0;
417:   case 2:
418:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, 0, a);
419:   case 4:
420:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, 0, a);
421:   case 8:
422:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, 0, a);
423:   default:
424:     SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralInvert()");
425:   }
426:   return 0;
427: }

429: /* compose dihedral symmetry: return b * a,
430:  * using the representation described in
431:  * DMPlexGetConeOrientation() */
432: static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
433: {
434:   switch (N) {
435:   case 0:
436:     return 0;
437:   case 2:
438:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
439:   case 4:
440:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
441:   case 8:
442:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
443:   default:
444:     SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
445:   }
446:   return 0;
447: }

449: /* swap dihedral symmetries: return b * a^-1,
450:  * using the representation described in
451:  * DMPlexGetConeOrientation() */
452: static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
453: {
454:   switch (N) {
455:   case 0:
456:     return 0;
457:   case 2:
458:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
459:   case 4:
460:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
461:   case 8:
462:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
463:   default:
464:     SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
465:   }
466:   return 0;
467: }
468: #endif
469: PETSC_EXTERN PetscInt       DMPolytopeConvertNewOrientation_Internal(DMPolytopeType, PetscInt);
470: PETSC_EXTERN PetscInt       DMPolytopeConvertOldOrientation_Internal(DMPolytopeType, PetscInt);
471: PETSC_EXTERN PetscErrorCode DMPlexConvertOldOrientations_Internal(DM);

473: PETSC_EXTERN PetscErrorCode DMPlexComputeIntegral_Internal(DM, Vec, PetscInt, PetscInt, PetscScalar *, void *);
474: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, PetscFormKey, IS, PetscReal, Vec, Vec, PetscReal, Vec, void *);
475: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Hybrid_Internal(DM, PetscFormKey[], IS, PetscReal, Vec, Vec, PetscReal, Vec, void *);
476: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, PetscFormKey, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
477: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Hybrid_Internal(DM, PetscFormKey[], IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
478: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Action_Internal(DM, PetscFormKey, IS, PetscReal, PetscReal, Vec, Vec, Vec, Vec, void *);
479: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);

481: /* Matvec with A in row-major storage, x and y can be aliased */
482: static inline void DMPlex_Mult2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
483: {
484:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
485:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1];
486:   y[1 * ldx]             = A[2] * z[0] + A[3] * z[1];
487:   (void)PetscLogFlops(6.0);
488: }
489: static inline void DMPlex_Mult3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
490: {
491:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
492:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
493:   y[1 * ldx]             = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
494:   y[2 * ldx]             = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
495:   (void)PetscLogFlops(15.0);
496: }
497: static inline void DMPlex_MultTranspose2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
498: {
499:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
500:   y[0 * ldx]             = A[0] * z[0] + A[2] * z[1];
501:   y[1 * ldx]             = A[1] * z[0] + A[3] * z[1];
502:   (void)PetscLogFlops(6.0);
503: }
504: static inline void DMPlex_MultTranspose3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
505: {
506:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
507:   y[0 * ldx]             = PetscConj(A[0]) * z[0] + PetscConj(A[3]) * z[1] + PetscConj(A[6]) * z[2];
508:   y[1 * ldx]             = PetscConj(A[1]) * z[0] + PetscConj(A[4]) * z[1] + PetscConj(A[7]) * z[2];
509:   y[2 * ldx]             = PetscConj(A[2]) * z[0] + PetscConj(A[5]) * z[1] + PetscConj(A[8]) * z[2];
510:   (void)PetscLogFlops(15.0);
511: }
512: static inline void DMPlex_Mult2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
513: {
514:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
515:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1];
516:   y[1 * ldx]             = A[2] * z[0] + A[3] * z[1];
517:   (void)PetscLogFlops(6.0);
518: }
519: static inline void DMPlex_Mult3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
520: {
521:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
522:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
523:   y[1 * ldx]             = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
524:   y[2 * ldx]             = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
525:   (void)PetscLogFlops(15.0);
526: }
527: static inline void DMPlex_MultAdd2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
528: {
529:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
530:   y[0 * ldx] += A[0] * z[0] + A[1] * z[1];
531:   y[1 * ldx] += A[2] * z[0] + A[3] * z[1];
532:   (void)PetscLogFlops(6.0);
533: }
534: static inline void DMPlex_MultAdd3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
535: {
536:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
537:   y[0 * ldx] += A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
538:   y[1 * ldx] += A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
539:   y[2 * ldx] += A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
540:   (void)PetscLogFlops(15.0);
541: }
542: /*
543:   A: packed, row-major m x n array
544:   x: length m array
545:   y: length n arra
546:   ldx: the stride in x and y

548:   A[i,j] = A[i * n + j]
549:   A^T[j,i] = A[i,j]
550: */
551: static inline void DMPlex_MultTransposeReal_Internal(const PetscReal A[], PetscInt m, PetscInt n, PetscInt ldx, const PetscScalar x[], PetscScalar y[])
552: {
553:   PetscScalar z[3];
554:   PetscInt    i, j;
555:   for (i = 0; i < m; ++i) z[i] = x[i * ldx];
556:   for (j = 0; j < n; ++j) {
557:     const PetscInt l = j * ldx;
558:     y[l]             = 0;
559:     for (i = 0; i < m; ++i) y[l] += A[i * n + j] * z[i];
560:   }
561:   (void)PetscLogFlops(2 * m * n);
562: }
563: static inline void DMPlex_MultTranspose2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
564: {
565:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
566:   y[0 * ldx]             = A[0] * z[0] + A[2] * z[1];
567:   y[1 * ldx]             = A[1] * z[0] + A[3] * z[1];
568:   (void)PetscLogFlops(6.0);
569: }
570: static inline void DMPlex_MultTranspose3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
571: {
572:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
573:   y[0 * ldx]             = A[0] * z[0] + A[3] * z[1] + A[6] * z[2];
574:   y[1 * ldx]             = A[1] * z[0] + A[4] * z[1] + A[7] * z[2];
575:   y[2 * ldx]             = A[2] * z[0] + A[5] * z[1] + A[8] * z[2];
576:   (void)PetscLogFlops(15.0);
577: }

579: static inline void DMPlex_MatMult2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
580: {
581: #define PLEX_DIM__ 2
582:   PetscScalar z[PLEX_DIM__];
583:   for (PetscInt j = 0; j < n; ++j) {
584:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
585:     DMPlex_Mult2D_Internal(A, 1, z, z);
586:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
587:   }
588:   (void)PetscLogFlops(8.0 * n);
589: #undef PLEX_DIM__
590: }
591: static inline void DMPlex_MatMult3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
592: {
593: #define PLEX_DIM__ 3
594:   PetscScalar z[PLEX_DIM__];
595:   for (PetscInt j = 0; j < n; ++j) {
596:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
597:     DMPlex_Mult3D_Internal(A, 1, z, z);
598:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
599:   }
600:   (void)PetscLogFlops(8.0 * n);
601: #undef PLEX_DIM__
602: }
603: static inline void DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
604: {
605: #define PLEX_DIM__ 2
606:   PetscScalar z[PLEX_DIM__];
607:   for (PetscInt j = 0; j < n; ++j) {
608:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
609:     DMPlex_MultTranspose2D_Internal(A, 1, z, z);
610:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
611:   }
612:   (void)PetscLogFlops(8.0 * n);
613: #undef PLEX_DIM__
614: }
615: static inline void DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
616: {
617: #define PLEX_DIM__ 3
618:   PetscScalar z[PLEX_DIM__];
619:   for (PetscInt j = 0; j < n; ++j) {
620:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
621:     DMPlex_MultTranspose3D_Internal(A, 1, z, z);
622:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
623:   }
624:   (void)PetscLogFlops(8.0 * n);
625: #undef PLEX_DIM__
626: }

628: static inline void DMPlex_MatMultLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
629: {
630:   for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
631:   (void)PetscLogFlops(8.0 * m);
632: }
633: static inline void DMPlex_MatMultLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
634: {
635:   for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
636:   (void)PetscLogFlops(8.0 * m);
637: }
638: static inline void DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
639: {
640:   for (PetscInt j = 0; j < m; ++j) DMPlex_Mult2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
641:   (void)PetscLogFlops(8.0 * m);
642: }
643: static inline void DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
644: {
645:   for (PetscInt j = 0; j < m; ++j) DMPlex_Mult3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
646:   (void)PetscLogFlops(8.0 * m);
647: }

649: static inline void DMPlex_PTAP2DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
650: {
651:   PetscScalar out[4];
652:   PetscInt    i, j, k, l;
653:   for (i = 0; i < 2; ++i) {
654:     for (j = 0; j < 2; ++j) {
655:       out[i * 2 + j] = 0.;
656:       for (k = 0; k < 2; ++k) {
657:         for (l = 0; l < 2; ++l) out[i * 2 + j] += P[k * 2 + i] * A[k * 2 + l] * P[l * 2 + j];
658:       }
659:     }
660:   }
661:   for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
662:   (void)PetscLogFlops(48.0);
663: }
664: static inline void DMPlex_PTAP3DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
665: {
666:   PetscScalar out[9];
667:   PetscInt    i, j, k, l;
668:   for (i = 0; i < 3; ++i) {
669:     for (j = 0; j < 3; ++j) {
670:       out[i * 3 + j] = 0.;
671:       for (k = 0; k < 3; ++k) {
672:         for (l = 0; l < 3; ++l) out[i * 3 + j] += P[k * 3 + i] * A[k * 3 + l] * P[l * 3 + j];
673:       }
674:     }
675:   }
676:   for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
677:   (void)PetscLogFlops(243.0);
678: }
679: /* TODO Fix for aliasing of A and C */
680: static inline void DMPlex_PTAPReal_Internal(const PetscReal P[], PetscInt m, PetscInt n, const PetscScalar A[], PetscScalar C[])
681: {
682:   PetscInt i, j, k, l;
683:   for (i = 0; i < n; ++i) {
684:     for (j = 0; j < n; ++j) {
685:       C[i * n + j] = 0.;
686:       for (k = 0; k < m; ++k) {
687:         for (l = 0; l < m; ++l) C[i * n + j] += P[k * n + i] * A[k * m + l] * P[l * n + j];
688:       }
689:     }
690:   }
691:   (void)PetscLogFlops(243.0);
692: }

694: static inline void DMPlex_Transpose2D_Internal(PetscScalar A[])
695: {
696:   PetscScalar tmp;
697:   tmp  = A[1];
698:   A[1] = A[2];
699:   A[2] = tmp;
700: }
701: static inline void DMPlex_Transpose3D_Internal(PetscScalar A[])
702: {
703:   PetscScalar tmp;
704:   tmp  = A[1];
705:   A[1] = A[3];
706:   A[3] = tmp;
707:   tmp  = A[2];
708:   A[2] = A[6];
709:   A[6] = tmp;
710:   tmp  = A[5];
711:   A[5] = A[7];
712:   A[7] = tmp;
713: }

715: static inline void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
716: {
717:   // Allow zero volume cells
718:   const PetscReal invDet = detJ == 0 ? 1.0 : 1.0 / detJ;

720:   invJ[0] = invDet * J[3];
721:   invJ[1] = -invDet * J[1];
722:   invJ[2] = -invDet * J[2];
723:   invJ[3] = invDet * J[0];
724:   (void)PetscLogFlops(5.0);
725: }

727: static inline void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
728: {
729:   // Allow zero volume cells
730:   const PetscReal invDet = detJ == 0 ? 1.0 : 1.0 / detJ;

732:   invJ[0 * 3 + 0] = invDet * (J[1 * 3 + 1] * J[2 * 3 + 2] - J[1 * 3 + 2] * J[2 * 3 + 1]);
733:   invJ[0 * 3 + 1] = invDet * (J[0 * 3 + 2] * J[2 * 3 + 1] - J[0 * 3 + 1] * J[2 * 3 + 2]);
734:   invJ[0 * 3 + 2] = invDet * (J[0 * 3 + 1] * J[1 * 3 + 2] - J[0 * 3 + 2] * J[1 * 3 + 1]);
735:   invJ[1 * 3 + 0] = invDet * (J[1 * 3 + 2] * J[2 * 3 + 0] - J[1 * 3 + 0] * J[2 * 3 + 2]);
736:   invJ[1 * 3 + 1] = invDet * (J[0 * 3 + 0] * J[2 * 3 + 2] - J[0 * 3 + 2] * J[2 * 3 + 0]);
737:   invJ[1 * 3 + 2] = invDet * (J[0 * 3 + 2] * J[1 * 3 + 0] - J[0 * 3 + 0] * J[1 * 3 + 2]);
738:   invJ[2 * 3 + 0] = invDet * (J[1 * 3 + 0] * J[2 * 3 + 1] - J[1 * 3 + 1] * J[2 * 3 + 0]);
739:   invJ[2 * 3 + 1] = invDet * (J[0 * 3 + 1] * J[2 * 3 + 0] - J[0 * 3 + 0] * J[2 * 3 + 1]);
740:   invJ[2 * 3 + 2] = invDet * (J[0 * 3 + 0] * J[1 * 3 + 1] - J[0 * 3 + 1] * J[1 * 3 + 0]);
741:   (void)PetscLogFlops(37.0);
742: }

744: static inline void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
745: {
746:   *detJ = J[0] * J[3] - J[1] * J[2];
747:   (void)PetscLogFlops(3.0);
748: }

750: static inline void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
751: {
752:   *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]));
753:   (void)PetscLogFlops(12.0);
754: }

756: static inline void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
757: {
758:   *detJ = PetscRealPart(J[0]) * PetscRealPart(J[3]) - PetscRealPart(J[1]) * PetscRealPart(J[2]);
759:   (void)PetscLogFlops(3.0);
760: }

762: static inline void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
763: {
764:   *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])));
765:   (void)PetscLogFlops(12.0);
766: }

768: static inline void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
769: {
770:   PetscInt d;
771:   for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
772: }

774: static inline PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y)
775: {
776:   PetscReal sum = 0.0;
777:   PetscInt  d;
778:   for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d]) * y[d];
779:   return sum;
780: }

782: static inline PetscReal DMPlex_DotRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y)
783: {
784:   PetscReal sum = 0.0;
785:   PetscInt  d;
786:   for (d = 0; d < dim; ++d) sum += x[d] * y[d];
787:   return sum;
788: }

790: static inline PetscReal DMPlex_NormD_Internal(PetscInt dim, const PetscReal *x)
791: {
792:   PetscReal sum = 0.0;
793:   PetscInt  d;
794:   for (d = 0; d < dim; ++d) sum += x[d] * x[d];
795:   return PetscSqrtReal(sum);
796: }

798: static inline PetscReal DMPlex_DistD_Internal(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
799: {
800:   PetscReal sum = 0.0;
801:   PetscInt  d;
802:   for (d = 0; d < dim; ++d) sum += PetscRealPart(PetscConj(x[d] - y[d]) * (x[d] - y[d]));
803:   return PetscSqrtReal(sum);
804: }

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 DMPlexVecGetOrientedClosure_Internal(DM, PetscSection, PetscBool, Vec, PetscInt, PetscInt, PetscInt *, PetscScalar *[]);
811: PETSC_INTERN PetscErrorCode DMPlexMatSetClosure_Internal(DM, PetscSection, PetscSection, PetscBool, Mat, PetscInt, const PetscScalar[], InsertMode);

813: PETSC_EXTERN PetscErrorCode DMPlexGetAllCells_Internal(DM, IS *);
814: PETSC_EXTERN PetscErrorCode DMPlexGetAllFaces_Internal(DM, IS *);
815: PETSC_EXTERN PetscErrorCode DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
816: PETSC_EXTERN PetscErrorCode DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
817: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
818: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
819: PETSC_INTERN PetscErrorCode DMCreateSubDomainDM_Plex(DM, DMLabel, PetscInt, IS *, DM *);
820: PETSC_INTERN PetscErrorCode DMPlexBasisTransformPoint_Internal(DM, DM, Vec, PetscInt, PetscBool[], PetscBool, PetscScalar *);
821: PETSC_EXTERN PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM, DM, Vec, PetscInt, PetscBool, PetscInt, PetscScalar *);
822: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscReal *, PetscReal *, void *);
823: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApply_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscScalar *, PetscScalar *, void *);
824: PETSC_INTERN PetscErrorCode DMCreateNeumannOverlap_Plex(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **);
825: PETSC_INTERN PetscErrorCode DMPlexMarkBoundaryFaces_Internal(DM, PetscInt, PetscInt, DMLabel, PetscBool);
826: PETSC_INTERN PetscErrorCode DMPlexDistributeOverlap_Internal(DM, PetscInt, MPI_Comm, const char *, PetscSF *, DM *);

828: PETSC_INTERN PetscErrorCode DMPeriodicCoordinateSetUp_Internal(DM);

830: /* Functions in the vtable */
831: PETSC_INTERN PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
832: PETSC_INTERN PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat);
833: PETSC_INTERN PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
834: PETSC_INTERN PetscErrorCode DMCreateMassMatrixLumped_Plex(DM dmCoarse, Vec *lm);
835: PETSC_INTERN PetscErrorCode DMCreateLocalSection_Plex(DM dm);
836: PETSC_INTERN PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm);
837: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J);
838: PETSC_INTERN PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
839: PETSC_INTERN PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field);
840: PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
841: PETSC_INTERN PetscErrorCode DMSetUp_Plex(DM dm);
842: PETSC_INTERN PetscErrorCode DMDestroy_Plex(DM dm);
843: PETSC_INTERN PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
844: PETSC_INTERN PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
845: PETSC_INTERN PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm);
846: PETSC_INTERN PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm);
847: PETSC_INTERN PetscErrorCode DMCreateDomainDecompositionScatters_Plex(DM, PetscInt, DM *, VecScatter **, VecScatter **, VecScatter **);
848: PETSC_INTERN PetscErrorCode DMCreateDomainDecomposition_Plex(DM, PetscInt *, char ***, IS **, IS **, DM **);
849: PETSC_INTERN PetscErrorCode DMCreateSectionPermutation_Plex(DM dm, IS *permutation, PetscBT *blockStarts);

851: // Coordinate mapping functions
852: 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[]);
853: 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[]);
854: 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[]);
855: 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[]);
856: 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[]);