Actual source code: dmpleximpl.h

  1: #ifndef _PLEXIMPL_H
  2: #define _PLEXIMPL_H

  4: #include <petscmat.h>
  5: #include <petscdmplex.h>
  6: #include <petscdmplextransform.h>
  7: #include <petscbt.h>
  8: #include <petscsf.h>
  9: #include <petsc/private/dmimpl.h>

 11: #if defined(PETSC_HAVE_EXODUSII)
 12:   #include <exodusII.h>
 13: #endif

 15: PETSC_EXTERN PetscLogEvent DMPLEX_Interpolate;
 16: PETSC_EXTERN PetscLogEvent DMPLEX_Partition;
 17: PETSC_EXTERN PetscLogEvent DMPLEX_PartSelf;
 18: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelInvert;
 19: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelCreateSF;
 20: PETSC_EXTERN PetscLogEvent DMPLEX_PartStratSF;
 21: PETSC_EXTERN PetscLogEvent DMPLEX_CreatePointSF;
 22: PETSC_EXTERN PetscLogEvent DMPLEX_Distribute;
 23: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeCones;
 24: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeLabels;
 25: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeSF;
 26: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeOverlap;
 27: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeField;
 28: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeData;
 29: PETSC_EXTERN PetscLogEvent DMPLEX_Migrate;
 30: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolateSF;
 31: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalBegin;
 32: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalEnd;
 33: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalBegin;
 34: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalEnd;
 35: PETSC_EXTERN PetscLogEvent DMPLEX_Stratify;
 36: PETSC_EXTERN PetscLogEvent DMPLEX_Symmetrize;
 37: PETSC_EXTERN PetscLogEvent DMPLEX_Preallocate;
 38: PETSC_EXTERN PetscLogEvent DMPLEX_ResidualFEM;
 39: PETSC_EXTERN PetscLogEvent DMPLEX_JacobianFEM;
 40: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolatorFEM;
 41: PETSC_EXTERN PetscLogEvent DMPLEX_InjectorFEM;
 42: PETSC_EXTERN PetscLogEvent DMPLEX_IntegralFEM;
 43: PETSC_EXTERN PetscLogEvent DMPLEX_CreateGmsh;
 44: PETSC_EXTERN PetscLogEvent DMPLEX_RebalanceSharedPoints;
 45: PETSC_EXTERN PetscLogEvent DMPLEX_CreateFromFile;
 46: PETSC_EXTERN PetscLogEvent DMPLEX_BuildFromCellList;
 47: PETSC_EXTERN PetscLogEvent DMPLEX_BuildCoordinatesFromCellList;
 48: PETSC_EXTERN PetscLogEvent DMPLEX_LocatePoints;
 49: PETSC_EXTERN PetscLogEvent DMPLEX_TopologyView;
 50: PETSC_EXTERN PetscLogEvent DMPLEX_DistributionView;
 51: PETSC_EXTERN PetscLogEvent DMPLEX_LabelsView;
 52: PETSC_EXTERN PetscLogEvent DMPLEX_CoordinatesView;
 53: PETSC_EXTERN PetscLogEvent DMPLEX_SectionView;
 54: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalVectorView;
 55: PETSC_EXTERN PetscLogEvent DMPLEX_LocalVectorView;
 56: PETSC_EXTERN PetscLogEvent DMPLEX_TopologyLoad;
 57: PETSC_EXTERN PetscLogEvent DMPLEX_DistributionLoad;
 58: PETSC_EXTERN PetscLogEvent DMPLEX_LabelsLoad;
 59: PETSC_EXTERN PetscLogEvent DMPLEX_CoordinatesLoad;
 60: PETSC_EXTERN PetscLogEvent DMPLEX_SectionLoad;
 61: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalVectorLoad;
 62: PETSC_EXTERN PetscLogEvent DMPLEX_LocalVectorLoad;
 63: PETSC_EXTERN PetscLogEvent DMPLEX_MetricEnforceSPD;
 64: PETSC_EXTERN PetscLogEvent DMPLEX_MetricNormalize;
 65: PETSC_EXTERN PetscLogEvent DMPLEX_MetricAverage;
 66: PETSC_EXTERN PetscLogEvent DMPLEX_MetricIntersection;
 67: PETSC_EXTERN PetscLogEvent DMPLEX_Generate;

 69: PETSC_EXTERN PetscLogEvent DMPLEX_RebalBuildGraph;
 70: PETSC_EXTERN PetscLogEvent DMPLEX_RebalRewriteSF;
 71: PETSC_EXTERN PetscLogEvent DMPLEX_RebalGatherGraph;
 72: PETSC_EXTERN PetscLogEvent DMPLEX_RebalPartition;
 73: PETSC_EXTERN PetscLogEvent DMPLEX_RebalScatterPart;

 75: /* Utility struct to store the contents of a Fluent file in memory */
 76: typedef struct {
 77:   int          index; /* Type of section */
 78:   unsigned int zoneID;
 79:   unsigned int first;
 80:   unsigned int last;
 81:   int          type;
 82:   int          nd; /* Either ND or element-type */
 83:   void        *data;
 84: } FluentSection;

 86: struct _PetscGridHash {
 87:   PetscInt     dim;
 88:   PetscReal    lower[3];    /* The lower-left corner */
 89:   PetscReal    upper[3];    /* The upper-right corner */
 90:   PetscReal    extent[3];   /* The box size */
 91:   PetscReal    h[3];        /* The subbox size */
 92:   PetscInt     n[3];        /* The number of subboxes */
 93:   PetscSection cellSection; /* Offsets for cells in each subbox*/
 94:   IS           cells;       /* List of cells in each subbox */
 95:   DMLabel      cellsSparse; /* Sparse storage for cell map */
 96: };

 98: typedef struct {
 99:   PetscBool isotropic;               /* Is the metric isotropic? */
100:   PetscBool uniform;                 /* Is the metric uniform? */
101:   PetscBool restrictAnisotropyFirst; /* Should anisotropy or normalization come first? */
102:   PetscBool noInsert;                /* Should node insertion/deletion be turned off? */
103:   PetscBool noSwap;                  /* Should facet swapping be turned off? */
104:   PetscBool noMove;                  /* Should node movement be turned off? */
105:   PetscBool noSurf;                  /* Should surface modification be turned off? */
106:   PetscReal h_min, h_max;            /* Minimum/maximum tolerated metric magnitudes */
107:   PetscReal a_max;                   /* Maximum tolerated anisotropy */
108:   PetscReal targetComplexity;        /* Target metric complexity */
109:   PetscReal p;                       /* Degree for L-p normalization methods */
110:   PetscReal gradationFactor;         /* Maximum tolerated length ratio for adjacent edges */
111:   PetscReal hausdorffNumber;         /* Max. distance between piecewise linear representation of boundary and reconstructed ideal boundary */
112:   PetscInt  numIter;                 /* Number of ParMmg mesh adaptation iterations */
113:   PetscInt  verbosity;               /* Level of verbosity for remesher (-1 = no output, 10 = maximum) */
114: } DMPlexMetricCtx;

116: /* Point Numbering in Plex:

118:    Points are numbered contiguously by stratum. Strate are organized as follows:

120:    First Stratum:  Cells [height 0]
121:    Second Stratum: Vertices [depth 0]
122:    Third Stratum:  Faces [height 1]
123:    Fourth Stratum: Edges [depth 1]

125:    We do this so that the numbering of a cell-vertex mesh does not change after interpolation. Within a given stratum,
126:    we allow additional segregation of by cell type.
127: */
128: typedef struct {
129:   PetscInt refct;

131:   PetscSection coneSection;      /* Layout of cones (inedges for DAG) */
132:   PetscInt    *cones;            /* Cone for each point */
133:   PetscInt    *coneOrientations; /* Orientation of each cone point, means cone traversal should start on point 'o', and if negative start on -(o+1) and go in reverse */
134:   PetscSection supportSection;   /* Layout of cones (inedges for DAG) */
135:   PetscInt    *supports;         /* Cone for each point */
136:   PetscInt    *facesTmp;         /* Work space for faces operation */

138:   /* Transformation */
139:   DMPlexTransform tr;                                               /* Type of transform used to define an ephemeral mesh */
140:   char           *transformType;                                    /* Type of transform for uniform cell refinement */
141:   PetscBool       refinementUniform;                                /* Flag for uniform cell refinement */
142:   PetscReal       refinementLimit;                                  /* Maximum volume for refined cell */
143:   PetscErrorCode (*refinementFunc)(const PetscReal[], PetscReal *); /* Function giving the maximum volume for refined cell */

145:   /* Interpolation */
146:   DMPlexInterpolatedFlag interpolated;
147:   DMPlexInterpolatedFlag interpolatedCollective;

149:   DMPlexReorderDefaultFlag reorderDefault; /* Reorder the DM by default */

151:   /* Distribution */
152:   PetscBool distDefault;      /* Distribute the DM by default */
153:   PetscInt  overlap;          /* Overlap of the partitions as passed to DMPlexDistribute() or DMPlexDistributeOverlap() */
154:   PetscInt  numOvLabels;      /* The number of labels used for candidate overlap points */
155:   DMLabel   ovLabels[16];     /* Labels used for candidate overlap points */
156:   PetscInt  ovValues[16];     /* Label values used for candidate overlap points */
157:   PetscInt  numOvExLabels;    /* The number of labels used for exclusion */
158:   DMLabel   ovExLabels[16];   /* Labels used to exclude points from the overlap */
159:   PetscInt  ovExValues[16];   /* Label values used to exclude points from the overlap */
160:   char     *distributionName; /* Name of the specific parallel distribution of the DM */

162:   /* Hierarchy */
163:   PetscBool regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */

165:   /* Generation */
166:   char            *tetgenOpts;
167:   char            *triangleOpts;
168:   PetscPartitioner partitioner;
169:   PetscBool        partitionBalance; /* Evenly divide partition overlap when distributing */
170:   PetscBool        remeshBd;

172:   /* Submesh */
173:   DMLabel          subpointMap;   /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
174:   IS               subpointIS;    /* IS holding point number in the enclosing mesh of every point in the submesh chart */
175:   PetscObjectState subpointState; /* The state of subpointMap when the subpointIS was last created */

177:   /* Labels and numbering */
178:   PetscObjectState depthState;    /* State of depth label, so that we can determine if a user changes it */
179:   PetscObjectState celltypeState; /* State of celltype label, so that we can determine if a user changes it */
180:   IS               globalVertexNumbers;
181:   IS               globalCellNumbers;

183:   /* Constraints */
184:   PetscSection anchorSection;          /* maps constrained points to anchor points */
185:   IS           anchorIS;               /* anchors indexed by the above section */
186:   PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
187:   PetscErrorCode (*computeanchormatrix)(DM, PetscSection, PetscSection, Mat);

189:   /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
190:   PetscSection parentSection; /* dof == 1 if point has parent */
191:   PetscInt    *parents;       /* point to parent */
192:   PetscInt    *childIDs;      /* point to child ID */
193:   PetscSection childSection;  /* inverse of parent section */
194:   PetscInt    *children;      /* point to children */
195:   DM           referenceTree; /* reference tree to which child ID's refer */
196:   PetscErrorCode (*getchildsymmetry)(DM, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt, PetscInt *, PetscInt *);

198:   /* MATIS support */
199:   PetscSection subdomainSection;

201:   /* Adjacency */
202:   PetscBool useAnchors;                                                          /* Replace constrained points with their anchors in adjacency lists */
203:   PetscErrorCode (*useradjacency)(DM, PetscInt, PetscInt *, PetscInt[], void *); /* User callback for adjacency */
204:   void *useradjacencyctx;                                                        /* User context for callback */

206:   // Periodicity
207:   struct {
208:     // Specified by the user
209:     PetscScalar transform[4][4]; // geometric transform
210:     PetscSF     face_sf;         // root(donor faces) <-- leaf(local faces)
211:     // Created eagerly (depends on points)
212:     PetscSF composed_sf; // root(non-periodic global points) <-- leaf(local points)
213:     IS      periodic_points;
214:   } periodic;

216:   /* Projection */
217:   PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
218:   PetscInt activePoint;         /* current active point in iteration */

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

224:   /* Geometry */
225:   PetscBool     ignoreModel;                      /* Ignore the geometry model during refinement */
226:   PetscReal     minradius;                        /* Minimum distance from cell centroid to face */
227:   PetscBool     useHashLocation;                  /* Use grid hashing for point location */
228:   PetscGridHash lbox;                             /* Local box for searching */
229:   void (*coordFunc)(PetscInt, PetscInt, PetscInt, /* Function used to remap newly introduced vertices */
230:                     const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);

232:   /* Neighbors */
233:   PetscMPIInt *neighbors;

235:   /* Metric */
236:   DMPlexMetricCtx *metricCtx;

238:   /* Debugging */
239:   PetscBool printSetValues;
240:   PetscInt  printFEM;
241:   PetscInt  printL2;
242:   PetscInt  printLocate;
243:   PetscReal printTol;
244: } DM_Plex;

246: PETSC_INTERN PetscErrorCode DMPlexCopy_Internal(DM, PetscBool, PetscBool, DM);
247: PETSC_INTERN PetscErrorCode DMPlexReplace_Internal(DM, DM *);

249: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM, PetscViewer);
250: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec, PetscViewer);
251: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec, PetscViewer);
252: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec, PetscViewer);
253: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec, PetscViewer);
254: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec, PetscViewer);
255: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec, PetscViewer);
256: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
257: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM, PetscViewer);
258: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject, PetscViewer);
259: #if defined(PETSC_HAVE_HDF5)
260: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
261: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
262: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
263: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
264: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
265: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
266: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
267: PETSC_INTERN PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM, IS, PetscViewer);
268: PETSC_INTERN PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM, PetscViewer);
269: PETSC_INTERN PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM, IS, PetscViewer);
270: PETSC_INTERN PetscErrorCode DMPlexSectionView_HDF5_Internal(DM, PetscViewer, DM);
271: PETSC_INTERN PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
272: PETSC_INTERN PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
273: PETSC_INTERN PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM, PetscViewer, PetscSF *);
274: PETSC_INTERN PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM, PetscViewer, PetscSF);
275: PETSC_INTERN PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM, PetscViewer, PetscSF);
276: PETSC_INTERN PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, PetscSF *, PetscSF *);
277: PETSC_INTERN PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, Vec);
278: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
279: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
280: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM, PetscViewer);
281: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
282: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
283: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
284: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
285: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);

287: struct _n_DMPlexStorageVersion {
288:   int major, minor, subminor;
289: };
290: typedef struct _n_DMPlexStorageVersion *DMPlexStorageVersion;

292: PETSC_EXTERN PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionReading(PetscViewer, DMPlexStorageVersion *);
293: PETSC_EXTERN PetscErrorCode PetscViewerHDF5GetDMPlexStorageVersionWriting(PetscViewer, DMPlexStorageVersion *);
294: #endif
295: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_CGNS(Vec, PetscViewer);

297: PETSC_INTERN PetscErrorCode DMPlexVecGetClosureAtDepth_Internal(DM, PetscSection, Vec, PetscInt, PetscInt, PetscInt *, PetscScalar *[]);
298: PETSC_INTERN PetscErrorCode DMPlexClosurePoints_Private(DM, PetscInt, const PetscInt[], IS *);
299: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(DM, PetscOptionItems *);
300: PETSC_INTERN PetscErrorCode DMSetFromOptions_Overlap_Plex(DM, PetscOptionItems *, PetscInt *);
301: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
302: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM[]);
303: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
304: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM[]);
305: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, Vec, DMLabel, DMLabel, DM *);
306: PETSC_INTERN PetscErrorCode DMExtrude_Plex(DM, PetscInt, DM *);
307: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
308: PETSC_INTERN PetscErrorCode DMPlexInsertTimeDerivativeBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
309: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
310: PETSC_INTERN PetscErrorCode DMProjectFunctionLabelLocal_Plex(DM, PetscReal, DMLabel, PetscInt, const PetscInt[], PetscInt, const PetscInt[], PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, InsertMode, Vec);
311: PETSC_INTERN PetscErrorCode DMProjectFieldLocal_Plex(DM, PetscReal, Vec, void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode, Vec);
312: PETSC_INTERN PetscErrorCode DMProjectFieldLabelLocal_Plex(DM, PetscReal, DMLabel, PetscInt, const PetscInt[], PetscInt, const PetscInt[], Vec, void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode, Vec);
313: PETSC_INTERN PetscErrorCode DMProjectBdFieldLabelLocal_Plex(DM, PetscReal, DMLabel, PetscInt, const PetscInt[], PetscInt, const PetscInt[], Vec, void (**)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]), InsertMode, Vec);
314: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
315: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, const PetscReal[], PetscReal *);
316: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM, PetscReal, PetscErrorCode (**)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void **, Vec, PetscReal *);
317: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);

319: #if defined(PETSC_HAVE_EXODUSII)
320: PETSC_EXTERN PetscErrorCode DMView_PlexExodusII(DM, PetscViewer);
321: PETSC_INTERN PetscErrorCode VecView_PlexExodusII_Internal(Vec, PetscViewer);
322: PETSC_INTERN PetscErrorCode VecLoad_PlexExodusII_Internal(Vec, PetscViewer);
323: PETSC_INTERN PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec, int, int, int);
324: PETSC_INTERN PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec, int, int, int);
325: PETSC_INTERN PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec, int, int, int);
326: PETSC_INTERN PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec, int, int, int);
327: PETSC_INTERN PetscErrorCode EXOGetVarIndex_Internal(int, ex_entity_type, const char[], int *);
328: #endif
329: PETSC_INTERN PetscErrorCode DMView_PlexCGNS(DM, PetscViewer);
330: PETSC_INTERN PetscErrorCode DMPlexCreateCGNSFromFile_Internal(MPI_Comm, const char *, PetscBool, DM *);
331: PETSC_INTERN PetscErrorCode DMPlexCreateCGNS_Internal(MPI_Comm, PetscInt, PetscBool, DM *);
332: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM, PetscInt, PetscInt, PetscInt *);
333: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM, PetscInt, PetscBool, PetscBool, PetscBool, PetscInt *, PetscInt *[]);
334: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
335: PETSC_INTERN PetscErrorCode DMPlexRestoreRawFaces_Internal(DM, DMPolytopeType, const PetscInt[], PetscInt *, const DMPolytopeType *[], const PetscInt *[], const PetscInt *[]);
336: PETSC_INTERN PetscErrorCode DMPlexComputeCellType_Internal(DM, PetscInt, PetscInt, DMPolytopeType *);
337: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], DMLabel, PetscInt, const PetscScalar[], InsertMode);
338: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
339: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
340: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM, DM, const char *, DM *);
341: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM, DM, PetscSF, PetscInt *, Mat);
342: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM, DM, PetscSF, PetscInt *, Mat);
343: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM, PetscSection, PetscInt, PetscInt, const PetscInt[], const PetscInt ***, const PetscScalar[], PetscInt *, PetscInt *, PetscInt *[], PetscScalar *[], PetscInt[], PetscBool);
344: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection, PetscInt, PetscInt, PetscInt *, PetscBool, PetscInt, PetscInt[]);
345: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection, PetscInt, PetscInt, PetscInt[], PetscBool, PetscInt, PetscInt[]);
346: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM, PetscInt, const PetscScalar[], PetscInt, PetscInt *);
347: /* these two are PETSC_EXTERN just because of src/dm/impls/plex/tests/ex18.c */
348: PETSC_EXTERN PetscErrorCode DMPlexOrientInterface_Internal(DM);

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

353: PETSC_INTERN PetscErrorCode DMPlexInterpolateInPlace_Internal(DM);
354: PETSC_INTERN PetscErrorCode DMPlexCreateBoxMesh_Tensor_SFC_Internal(DM, PetscInt, const PetscInt[], const PetscReal[], const PetscReal[], const DMBoundaryType[], PetscBool);
355: PETSC_INTERN PetscErrorCode DMPlexMigrateIsoperiodicFaceSF_Internal(DM, DM, PetscSF);
356: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
357: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
358: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, Vec, DMLabel, DMLabel, DM *);
359: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, Vec, DMLabel, DMLabel, DM *);
360: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM, Mat *);

362: PETSC_INTERN PetscErrorCode DMPlexGetOverlap_Plex(DM, PetscInt *);
363: PETSC_INTERN PetscErrorCode DMPlexSetOverlap_Plex(DM, DM, PetscInt);
364: PETSC_INTERN PetscErrorCode DMPlexDistributeGetDefault_Plex(DM, PetscBool *);
365: PETSC_INTERN PetscErrorCode DMPlexDistributeSetDefault_Plex(DM, PetscBool);
366: PETSC_INTERN PetscErrorCode DMPlexReorderGetDefault_Plex(DM, DMPlexReorderDefaultFlag *);
367: PETSC_INTERN PetscErrorCode DMPlexReorderSetDefault_Plex(DM, DMPlexReorderDefaultFlag);

369: #if 1
370: static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
371: {
372:   return (a <= 0) ? a : (N - a);
373: }

375: static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
376: {
377:   if (!N) return 0;
378:   return (a >= 0) ? ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) : ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
379: }

381: static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
382: {
383:   return DihedralCompose(N, DihedralInvert(N, a), b);
384: }
385: #else
386: /* TODO
387:    This is a reimplementation of the tensor dihedral symmetries using the new orientations.
388:    These should be turned on when we convert to new-style orientations in p4est.
389: */
390: /* invert dihedral symmetry: return a^-1,
391:  * using the representation described in
392:  * DMPlexGetConeOrientation() */
393: static inline PetscInt DihedralInvert(PetscInt N, PetscInt a)
394: {
395:   switch (N) {
396:   case 0:
397:     return 0;
398:   case 2:
399:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, 0, a);
400:   case 4:
401:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, 0, a);
402:   case 8:
403:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, 0, a);
404:   default:
405:     SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralInvert()");
406:   }
407:   return 0;
408: }

410: /* compose dihedral symmetry: return b * a,
411:  * using the representation described in
412:  * DMPlexGetConeOrientation() */
413: static inline PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
414: {
415:   switch (N) {
416:   case 0:
417:     return 0;
418:   case 2:
419:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
420:   case 4:
421:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
422:   case 8:
423:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
424:   default:
425:     SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
426:   }
427:   return 0;
428: }

430: /* swap dihedral symmetries: return b * a^-1,
431:  * using the representation described in
432:  * DMPlexGetConeOrientation() */
433: static inline PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
434: {
435:   switch (N) {
436:   case 0:
437:     return 0;
438:   case 2:
439:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
440:   case 4:
441:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
442:   case 8:
443:     return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
444:   default:
445:     SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
446:   }
447:   return 0;
448: }
449: #endif
450: PETSC_EXTERN PetscInt       DMPolytopeConvertNewOrientation_Internal(DMPolytopeType, PetscInt);
451: PETSC_EXTERN PetscInt       DMPolytopeConvertOldOrientation_Internal(DMPolytopeType, PetscInt);
452: PETSC_EXTERN PetscErrorCode DMPlexConvertOldOrientations_Internal(DM);

454: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, PetscFormKey, IS, PetscReal, Vec, Vec, PetscReal, Vec, void *);
455: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Hybrid_Internal(DM, PetscFormKey[], IS, PetscReal, Vec, Vec, PetscReal, Vec, void *);
456: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, PetscFormKey, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
457: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Hybrid_Internal(DM, PetscFormKey[], IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
458: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Action_Internal(DM, PetscFormKey, IS, PetscReal, PetscReal, Vec, Vec, Vec, Vec, void *);
459: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);

461: /* Matvec with A in row-major storage, x and y can be aliased */
462: static inline void DMPlex_Mult2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
463: {
464:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
465:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1];
466:   y[1 * ldx]             = A[2] * z[0] + A[3] * z[1];
467:   (void)PetscLogFlops(6.0);
468: }
469: static inline void DMPlex_Mult3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
470: {
471:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
472:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
473:   y[1 * ldx]             = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
474:   y[2 * ldx]             = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
475:   (void)PetscLogFlops(15.0);
476: }
477: static inline void DMPlex_MultTranspose2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
478: {
479:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
480:   y[0 * ldx]             = A[0] * z[0] + A[2] * z[1];
481:   y[1 * ldx]             = A[1] * z[0] + A[3] * z[1];
482:   (void)PetscLogFlops(6.0);
483: }
484: static inline void DMPlex_MultTranspose3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
485: {
486:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
487:   y[0 * ldx]             = PetscConj(A[0]) * z[0] + PetscConj(A[3]) * z[1] + PetscConj(A[6]) * z[2];
488:   y[1 * ldx]             = PetscConj(A[1]) * z[0] + PetscConj(A[4]) * z[1] + PetscConj(A[7]) * z[2];
489:   y[2 * ldx]             = PetscConj(A[2]) * z[0] + PetscConj(A[5]) * z[1] + PetscConj(A[8]) * z[2];
490:   (void)PetscLogFlops(15.0);
491: }
492: static inline void DMPlex_Mult2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
493: {
494:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
495:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1];
496:   y[1 * ldx]             = A[2] * z[0] + A[3] * z[1];
497:   (void)PetscLogFlops(6.0);
498: }
499: static inline void DMPlex_Mult3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
500: {
501:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
502:   y[0 * ldx]             = A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
503:   y[1 * ldx]             = A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
504:   y[2 * ldx]             = A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
505:   (void)PetscLogFlops(15.0);
506: }
507: static inline void DMPlex_MultAdd2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
508: {
509:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
510:   y[0 * ldx] += A[0] * z[0] + A[1] * z[1];
511:   y[1 * ldx] += A[2] * z[0] + A[3] * z[1];
512:   (void)PetscLogFlops(6.0);
513: }
514: static inline void DMPlex_MultAdd3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
515: {
516:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
517:   y[0 * ldx] += A[0] * z[0] + A[1] * z[1] + A[2] * z[2];
518:   y[1 * ldx] += A[3] * z[0] + A[4] * z[1] + A[5] * z[2];
519:   y[2 * ldx] += A[6] * z[0] + A[7] * z[1] + A[8] * z[2];
520:   (void)PetscLogFlops(15.0);
521: }
522: /*
523:   A: packed, row-major m x n array
524:   x: length m array
525:   y: length n arra
526:   ldx: the stride in x and y

528:   A[i,j] = A[i * n + j]
529:   A^T[j,i] = A[i,j]
530: */
531: static inline void DMPlex_MultTransposeReal_Internal(const PetscReal A[], PetscInt m, PetscInt n, PetscInt ldx, const PetscScalar x[], PetscScalar y[])
532: {
533:   PetscScalar z[3];
534:   PetscInt    i, j;
535:   for (i = 0; i < m; ++i) z[i] = x[i * ldx];
536:   for (j = 0; j < n; ++j) {
537:     const PetscInt l = j * ldx;
538:     y[l]             = 0;
539:     for (i = 0; i < m; ++i) y[l] += A[i * n + j] * z[i];
540:   }
541:   (void)PetscLogFlops(2 * m * n);
542: }
543: static inline void DMPlex_MultTranspose2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
544: {
545:   const PetscScalar z[2] = {x[0 * ldx], x[1 * ldx]};
546:   y[0 * ldx]             = A[0] * z[0] + A[2] * z[1];
547:   y[1 * ldx]             = A[1] * z[0] + A[3] * z[1];
548:   (void)PetscLogFlops(6.0);
549: }
550: static inline void DMPlex_MultTranspose3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
551: {
552:   const PetscScalar z[3] = {x[0 * ldx], x[1 * ldx], x[2 * ldx]};
553:   y[0 * ldx]             = A[0] * z[0] + A[3] * z[1] + A[6] * z[2];
554:   y[1 * ldx]             = A[1] * z[0] + A[4] * z[1] + A[7] * z[2];
555:   y[2 * ldx]             = A[2] * z[0] + A[5] * z[1] + A[8] * z[2];
556:   (void)PetscLogFlops(15.0);
557: }

559: static inline void DMPlex_MatMult2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
560: {
561: #define PLEX_DIM__ 2
562:   PetscScalar z[PLEX_DIM__];
563:   for (PetscInt j = 0; j < n; ++j) {
564:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
565:     DMPlex_Mult2D_Internal(A, 1, z, z);
566:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
567:   }
568:   (void)PetscLogFlops(8.0 * n);
569: #undef PLEX_DIM__
570: }
571: static inline void DMPlex_MatMult3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
572: {
573: #define PLEX_DIM__ 3
574:   PetscScalar z[PLEX_DIM__];
575:   for (PetscInt j = 0; j < n; ++j) {
576:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
577:     DMPlex_Mult3D_Internal(A, 1, z, z);
578:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
579:   }
580:   (void)PetscLogFlops(8.0 * n);
581: #undef PLEX_DIM__
582: }
583: static inline void DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
584: {
585: #define PLEX_DIM__ 2
586:   PetscScalar z[PLEX_DIM__];
587:   for (PetscInt j = 0; j < n; ++j) {
588:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
589:     DMPlex_MultTranspose2D_Internal(A, 1, z, z);
590:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
591:   }
592:   (void)PetscLogFlops(8.0 * n);
593: #undef PLEX_DIM__
594: }
595: static inline void DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
596: {
597: #define PLEX_DIM__ 3
598:   PetscScalar z[PLEX_DIM__];
599:   for (PetscInt j = 0; j < n; ++j) {
600:     for (int d = 0; d < PLEX_DIM__; ++d) z[d] = B[d * ldb + j];
601:     DMPlex_MultTranspose3D_Internal(A, 1, z, z);
602:     for (int d = 0; d < PLEX_DIM__; ++d) C[d * ldb + j] = z[d];
603:   }
604:   (void)PetscLogFlops(8.0 * n);
605: #undef PLEX_DIM__
606: }

608: static inline void DMPlex_MatMultLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
609: {
610:   for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
611:   (void)PetscLogFlops(8.0 * m);
612: }
613: static inline void DMPlex_MatMultLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
614: {
615:   for (PetscInt j = 0; j < m; ++j) DMPlex_MultTranspose3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
616:   (void)PetscLogFlops(8.0 * m);
617: }
618: static inline void DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
619: {
620:   for (PetscInt j = 0; j < m; ++j) DMPlex_Mult2D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
621:   (void)PetscLogFlops(8.0 * m);
622: }
623: static inline void DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
624: {
625:   for (PetscInt j = 0; j < m; ++j) DMPlex_Mult3D_Internal(A, 1, &B[j * ldb], &C[j * ldb]);
626:   (void)PetscLogFlops(8.0 * m);
627: }

629: static inline void DMPlex_PTAP2DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
630: {
631:   PetscScalar out[4];
632:   PetscInt    i, j, k, l;
633:   for (i = 0; i < 2; ++i) {
634:     for (j = 0; j < 2; ++j) {
635:       out[i * 2 + j] = 0.;
636:       for (k = 0; k < 2; ++k) {
637:         for (l = 0; l < 2; ++l) out[i * 2 + j] += P[k * 2 + i] * A[k * 2 + l] * P[l * 2 + j];
638:       }
639:     }
640:   }
641:   for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
642:   (void)PetscLogFlops(48.0);
643: }
644: static inline void DMPlex_PTAP3DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
645: {
646:   PetscScalar out[9];
647:   PetscInt    i, j, k, l;
648:   for (i = 0; i < 3; ++i) {
649:     for (j = 0; j < 3; ++j) {
650:       out[i * 3 + j] = 0.;
651:       for (k = 0; k < 3; ++k) {
652:         for (l = 0; l < 3; ++l) out[i * 3 + j] += P[k * 3 + i] * A[k * 3 + l] * P[l * 3 + j];
653:       }
654:     }
655:   }
656:   for (i = 0; i < 2 * 2; ++i) C[i] = out[i];
657:   (void)PetscLogFlops(243.0);
658: }
659: /* TODO Fix for aliasing of A and C */
660: static inline void DMPlex_PTAPReal_Internal(const PetscReal P[], PetscInt m, PetscInt n, const PetscScalar A[], PetscScalar C[])
661: {
662:   PetscInt i, j, k, l;
663:   for (i = 0; i < n; ++i) {
664:     for (j = 0; j < n; ++j) {
665:       C[i * n + j] = 0.;
666:       for (k = 0; k < m; ++k) {
667:         for (l = 0; l < m; ++l) C[i * n + j] += P[k * n + i] * A[k * m + l] * P[l * n + j];
668:       }
669:     }
670:   }
671:   (void)PetscLogFlops(243.0);
672: }

674: static inline void DMPlex_Transpose2D_Internal(PetscScalar A[])
675: {
676:   PetscScalar tmp;
677:   tmp  = A[1];
678:   A[1] = A[2];
679:   A[2] = tmp;
680: }
681: static inline void DMPlex_Transpose3D_Internal(PetscScalar A[])
682: {
683:   PetscScalar tmp;
684:   tmp  = A[1];
685:   A[1] = A[3];
686:   A[3] = tmp;
687:   tmp  = A[2];
688:   A[2] = A[6];
689:   A[6] = tmp;
690:   tmp  = A[5];
691:   A[5] = A[7];
692:   A[7] = tmp;
693: }

695: static inline void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
696: {
697:   // Allow zero volume cells
698:   const PetscReal invDet = detJ == 0 ? 1.0 : 1.0 / detJ;

700:   invJ[0] = invDet * J[3];
701:   invJ[1] = -invDet * J[1];
702:   invJ[2] = -invDet * J[2];
703:   invJ[3] = invDet * J[0];
704:   (void)PetscLogFlops(5.0);
705: }

707: static inline void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
708: {
709:   // Allow zero volume cells
710:   const PetscReal invDet = detJ == 0 ? 1.0 : 1.0 / detJ;

712:   invJ[0 * 3 + 0] = invDet * (J[1 * 3 + 1] * J[2 * 3 + 2] - J[1 * 3 + 2] * J[2 * 3 + 1]);
713:   invJ[0 * 3 + 1] = invDet * (J[0 * 3 + 2] * J[2 * 3 + 1] - J[0 * 3 + 1] * J[2 * 3 + 2]);
714:   invJ[0 * 3 + 2] = invDet * (J[0 * 3 + 1] * J[1 * 3 + 2] - J[0 * 3 + 2] * J[1 * 3 + 1]);
715:   invJ[1 * 3 + 0] = invDet * (J[1 * 3 + 2] * J[2 * 3 + 0] - J[1 * 3 + 0] * J[2 * 3 + 2]);
716:   invJ[1 * 3 + 1] = invDet * (J[0 * 3 + 0] * J[2 * 3 + 2] - J[0 * 3 + 2] * J[2 * 3 + 0]);
717:   invJ[1 * 3 + 2] = invDet * (J[0 * 3 + 2] * J[1 * 3 + 0] - J[0 * 3 + 0] * J[1 * 3 + 2]);
718:   invJ[2 * 3 + 0] = invDet * (J[1 * 3 + 0] * J[2 * 3 + 1] - J[1 * 3 + 1] * J[2 * 3 + 0]);
719:   invJ[2 * 3 + 1] = invDet * (J[0 * 3 + 1] * J[2 * 3 + 0] - J[0 * 3 + 0] * J[2 * 3 + 1]);
720:   invJ[2 * 3 + 2] = invDet * (J[0 * 3 + 0] * J[1 * 3 + 1] - J[0 * 3 + 1] * J[1 * 3 + 0]);
721:   (void)PetscLogFlops(37.0);
722: }

724: static inline void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
725: {
726:   *detJ = J[0] * J[3] - J[1] * J[2];
727:   (void)PetscLogFlops(3.0);
728: }

730: static inline void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
731: {
732:   *detJ = (J[0 * 3 + 0] * (J[1 * 3 + 1] * J[2 * 3 + 2] - J[1 * 3 + 2] * J[2 * 3 + 1]) + J[0 * 3 + 1] * (J[1 * 3 + 2] * J[2 * 3 + 0] - J[1 * 3 + 0] * J[2 * 3 + 2]) + J[0 * 3 + 2] * (J[1 * 3 + 0] * J[2 * 3 + 1] - J[1 * 3 + 1] * J[2 * 3 + 0]));
733:   (void)PetscLogFlops(12.0);
734: }

736: static inline void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
737: {
738:   *detJ = PetscRealPart(J[0]) * PetscRealPart(J[3]) - PetscRealPart(J[1]) * PetscRealPart(J[2]);
739:   (void)PetscLogFlops(3.0);
740: }

742: static inline void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
743: {
744:   *detJ = (PetscRealPart(J[0 * 3 + 0]) * (PetscRealPart(J[1 * 3 + 1]) * PetscRealPart(J[2 * 3 + 2]) - PetscRealPart(J[1 * 3 + 2]) * PetscRealPart(J[2 * 3 + 1])) + PetscRealPart(J[0 * 3 + 1]) * (PetscRealPart(J[1 * 3 + 2]) * PetscRealPart(J[2 * 3 + 0]) - PetscRealPart(J[1 * 3 + 0]) * PetscRealPart(J[2 * 3 + 2])) + PetscRealPart(J[0 * 3 + 2]) * (PetscRealPart(J[1 * 3 + 0]) * PetscRealPart(J[2 * 3 + 1]) - PetscRealPart(J[1 * 3 + 1]) * PetscRealPart(J[2 * 3 + 0])));
745:   (void)PetscLogFlops(12.0);
746: }

748: static inline void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
749: {
750:   PetscInt d;
751:   for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
752: }

754: static inline PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y)
755: {
756:   PetscReal sum = 0.0;
757:   PetscInt  d;
758:   for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d]) * y[d];
759:   return sum;
760: }

762: static inline PetscReal DMPlex_DotRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y)
763: {
764:   PetscReal sum = 0.0;
765:   PetscInt  d;
766:   for (d = 0; d < dim; ++d) sum += x[d] * y[d];
767:   return sum;
768: }

770: static inline PetscReal DMPlex_NormD_Internal(PetscInt dim, const PetscReal *x)
771: {
772:   PetscReal sum = 0.0;
773:   PetscInt  d;
774:   for (d = 0; d < dim; ++d) sum += x[d] * x[d];
775:   return PetscSqrtReal(sum);
776: }

778: static inline PetscReal DMPlex_DistD_Internal(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
779: {
780:   PetscReal sum = 0.0;
781:   PetscInt  d;
782:   for (d = 0; d < dim; ++d) sum += PetscRealPart(PetscConj(x[d] - y[d]) * (x[d] - y[d]));
783:   return PetscSqrtReal(sum);
784: }

786: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM, PetscInt, PetscInt, PetscDualSpace *);
787: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection, PetscBool, PetscInt, PetscInt, PetscInt *, PetscBool, const PetscInt[], const PetscInt[], PetscInt[]);
788: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection, PetscBool, PetscInt, PetscInt, PetscInt[], PetscBool, const PetscInt ***, PetscInt, const PetscInt[], PetscInt[]);
789: PETSC_INTERN PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]);
790: PETSC_INTERN PetscErrorCode DMPlexVecGetOrientedClosure_Internal(DM, PetscSection, Vec, PetscInt, PetscInt, PetscInt *, PetscScalar *[]);

792: PETSC_EXTERN PetscErrorCode DMPlexGetAllCells_Internal(DM, IS *);
793: PETSC_EXTERN PetscErrorCode DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
794: PETSC_EXTERN PetscErrorCode DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
795: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
796: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
797: PETSC_INTERN PetscErrorCode DMCreateSubDomainDM_Plex(DM, DMLabel, PetscInt, IS *, DM *);
798: PETSC_INTERN PetscErrorCode DMPlexBasisTransformPoint_Internal(DM, DM, Vec, PetscInt, PetscBool[], PetscBool, PetscScalar *);
799: PETSC_EXTERN PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM, DM, Vec, PetscInt, PetscBool, PetscInt, PetscScalar *);
800: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscReal *, PetscReal *, void *);
801: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApply_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscScalar *, PetscScalar *, void *);
802: PETSC_INTERN PetscErrorCode DMCreateNeumannOverlap_Plex(DM, IS *, Mat *, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void *), void **);

804: PETSC_INTERN PetscErrorCode DMPeriodicCoordinateSetUp_Internal(DM);

806: /* Functions in the vtable */
807: PETSC_INTERN PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
808: PETSC_INTERN PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat);
809: PETSC_INTERN PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
810: PETSC_INTERN PetscErrorCode DMCreateMassMatrixLumped_Plex(DM dmCoarse, Vec *lm);
811: PETSC_INTERN PetscErrorCode DMCreateLocalSection_Plex(DM dm);
812: PETSC_INTERN PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm);
813: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J);
814: PETSC_INTERN PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
815: PETSC_INTERN PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field);
816: PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
817: PETSC_INTERN PetscErrorCode DMSetUp_Plex(DM dm);
818: PETSC_INTERN PetscErrorCode DMDestroy_Plex(DM dm);
819: PETSC_INTERN PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
820: PETSC_INTERN PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
821: PETSC_INTERN PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm);
822: PETSC_INTERN PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm);

824: #endif /* _PLEXIMPL_H */