Actual source code: dmpleximpl.h

petsc-3.13.6 2020-09-29
Report Typos and Errors
  1: #if !defined(_PLEXIMPL_H)
  2: #define _PLEXIMPL_H

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

 10: PETSC_EXTERN PetscLogEvent DMPLEX_Interpolate;
 11: PETSC_EXTERN PetscLogEvent DMPLEX_Partition;
 12: PETSC_EXTERN PetscLogEvent DMPLEX_PartSelf;
 13: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelInvert;
 14: PETSC_EXTERN PetscLogEvent DMPLEX_PartLabelCreateSF;
 15: PETSC_EXTERN PetscLogEvent DMPLEX_PartStratSF;
 16: PETSC_EXTERN PetscLogEvent DMPLEX_CreatePointSF;
 17: PETSC_EXTERN PetscLogEvent DMPLEX_Distribute;
 18: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeCones;
 19: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeLabels;
 20: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeSF;
 21: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeOverlap;
 22: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeField;
 23: PETSC_EXTERN PetscLogEvent DMPLEX_DistributeData;
 24: PETSC_EXTERN PetscLogEvent DMPLEX_Migrate;
 25: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolateSF;
 26: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalBegin;
 27: PETSC_EXTERN PetscLogEvent DMPLEX_GlobalToNaturalEnd;
 28: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalBegin;
 29: PETSC_EXTERN PetscLogEvent DMPLEX_NaturalToGlobalEnd;
 30: PETSC_EXTERN PetscLogEvent DMPLEX_Stratify;
 31: PETSC_EXTERN PetscLogEvent DMPLEX_Symmetrize;
 32: PETSC_EXTERN PetscLogEvent DMPLEX_Preallocate;
 33: PETSC_EXTERN PetscLogEvent DMPLEX_ResidualFEM;
 34: PETSC_EXTERN PetscLogEvent DMPLEX_JacobianFEM;
 35: PETSC_EXTERN PetscLogEvent DMPLEX_InterpolatorFEM;
 36: PETSC_EXTERN PetscLogEvent DMPLEX_InjectorFEM;
 37: PETSC_EXTERN PetscLogEvent DMPLEX_IntegralFEM;
 38: PETSC_EXTERN PetscLogEvent DMPLEX_CreateGmsh;
 39: PETSC_EXTERN PetscLogEvent DMPLEX_RebalanceSharedPoints;
 40: PETSC_EXTERN PetscLogEvent DMPLEX_CreateFromFile;
 41: PETSC_EXTERN PetscLogEvent DMPLEX_CreateFromCellList;
 42: PETSC_EXTERN PetscLogEvent DMPLEX_CreateFromCellList_Coordinates;

 44: PETSC_EXTERN PetscBool      PetscPartitionerRegisterAllCalled;
 45: PETSC_EXTERN PetscErrorCode PetscPartitionerRegisterAll(void);

 47: typedef struct _PetscPartitionerOps *PetscPartitionerOps;
 48: struct _PetscPartitionerOps {
 49:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscPartitioner);
 50:   PetscErrorCode (*setup)(PetscPartitioner);
 51:   PetscErrorCode (*view)(PetscPartitioner,PetscViewer);
 52:   PetscErrorCode (*destroy)(PetscPartitioner);
 53:   PetscErrorCode (*partition)(PetscPartitioner, PetscInt, PetscInt, PetscInt[], PetscInt[], PetscSection, PetscSection, PetscSection, IS *);
 54: };

 56: struct _p_PetscPartitioner {
 57:   PETSCHEADER(struct _PetscPartitionerOps);
 58:   void        *data;            /* Implementation object */
 59:   PetscInt    height;           /* Height of points to partition into non-overlapping subsets */
 60:   PetscInt    edgeCut;          /* The number of edge cut by the partition */
 61:   PetscReal   balance;          /* The maximum partition size divided by the minimum size */
 62:   PetscViewer viewer;
 63:   PetscViewer viewerGraph;
 64:   PetscBool   viewGraph;
 65:   PetscBool   noGraph;          /* if true, the partitioner does not need the connectivity graph, only the number of local vertices */
 66:   PetscBool   usevwgt;          /* if true, the partitioner looks at the local section vertSection to weight the vertices of the graph */
 67: };

 69: typedef struct {
 70:   PetscInt dummy;
 71: } PetscPartitioner_Chaco;

 73: typedef struct {
 74:   MPI_Comm  pcomm;
 75:   PetscInt  ptype;
 76:   PetscReal imbalanceRatio;
 77:   PetscInt  debugFlag;
 78:   PetscInt  randomSeed;
 79: } PetscPartitioner_ParMetis;

 81: typedef struct {
 82:   MPI_Comm  pcomm;
 83:   PetscInt  strategy;
 84:   PetscReal imbalance;
 85: } PetscPartitioner_PTScotch;

 87: static const char *const
 88: PTScotchStrategyList[] = {
 89:   "DEFAULT",
 90:   "QUALITY",
 91:   "SPEED",
 92:   "BALANCE",
 93:   "SAFETY",
 94:   "SCALABILITY",
 95:   "RECURSIVE",
 96:   "REMAP"
 97: };

 99: typedef struct {
100:   PetscSection section;   /* Sizes for each partition */
101:   IS           partition; /* Points in each partition */
102:   PetscBool    random;    /* Flag for a random partition */
103: } PetscPartitioner_Shell;

105: typedef struct {
106:   PetscInt dummy;
107: } PetscPartitioner_Simple;

109: typedef struct {
110:   PetscInt dummy;
111: } PetscPartitioner_Gather;

113: typedef struct _DMPlexCellRefinerOps *DMPlexCellRefinerOps;
114: struct _DMPlexCellRefinerOps {
115:   PetscErrorCode (*refine)(DMPlexCellRefiner, DMPolytopeType, PetscInt *, DMPolytopeType *[], PetscInt *[], PetscInt *[], PetscInt *[]);
116:   PetscErrorCode (*mapsubcells)(DMPlexCellRefiner, DMPolytopeType, PetscInt, DMPolytopeType, PetscInt, PetscInt, PetscInt *, PetscInt *);
117:   PetscErrorCode (*getaffinetransforms)(DMPlexCellRefiner, DMPolytopeType, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
118:   PetscErrorCode (*getaffinefacetransforms)(DMPlexCellRefiner, DMPolytopeType, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[], PetscReal *[]);
119:   PetscErrorCode (*getcellvertices)(DMPlexCellRefiner, DMPolytopeType, PetscInt *, PetscReal *[]);
120:   PetscErrorCode (*getsubcellvertices)(DMPlexCellRefiner, DMPolytopeType, DMPolytopeType, PetscInt, PetscInt *, PetscInt *[]);
121: };

123: struct _p_DMPlexCellRefiner {
124:   PETSCHEADER(struct _DMPlexCellRefinerOps);
125:   DM        dm;         /* The original DM */
126:   PetscBool setupcalled;
127:   DMPlexCellRefinerType type;
128:   PetscInt *ctOrder;    /* [i] = ct: An array with cell types in depth order */
129:   PetscInt *ctOrderInv; /* [ct] = i: An array with the ordinal numbers for each cell type */
130:   PetscInt *ctStart;    /* The number for the first cell of each polytope type in the original mesh, indexed by cell type */
131:   PetscInt *ctStartNew; /* The number for the first cell of each polytope type in the new mesh, indexed by cell type */
132:   PetscInt *offset;     /* [ct][ctNew]: The offset in the new point numbering of a point of type ctNew produced from an old point of type ct */
133:   PetscFE      *coordFE; /* Finite element for each cell type, used for localized coordinate interpolation */
134:   PetscFEGeom **refGeom; /* Geometry of the reference cell for each cell type */
135: };

137: /* Utility struct to store the contents of a Fluent file in memory */
138: typedef struct {
139:   int          index;    /* Type of section */
140:   unsigned int zoneID;
141:   unsigned int first;
142:   unsigned int last;
143:   int          type;
144:   int          nd;       /* Either ND or element-type */
145:   void        *data;
146: } FluentSection;

148: struct _PetscGridHash {
149:   PetscInt     dim;
150:   PetscReal    lower[3];    /* The lower-left corner */
151:   PetscReal    upper[3];    /* The upper-right corner */
152:   PetscReal    extent[3];   /* The box size */
153:   PetscReal    h[3];        /* The subbox size */
154:   PetscInt     n[3];        /* The number of subboxes */
155:   PetscSection cellSection; /* Offsets for cells in each subbox*/
156:   IS           cells;       /* List of cells in each subbox */
157:   DMLabel      cellsSparse; /* Sparse storage for cell map */
158: };

160: /* Point Numbering in Plex:

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

164:    First Stratum:  Cells [height 0]
165:    Second Stratum: Vertices [depth 0]
166:    Third Stratum:  Faces [height 1]
167:    Fourth Stratum: Edges [depth 1]

169:    We do this so that the numbering of a cell-vertex mesh does not change after interpolation. Within a given stratum,
170:    we allow additional segregation of by cell type.
171: */
172: typedef struct {
173:   PetscInt             refct;

175:   PetscSection         coneSection;       /* Layout of cones (inedges for DAG) */
176:   PetscInt             maxConeSize;       /* Cached for fast lookup */
177:   PetscInt            *cones;             /* Cone for each point */
178:   PetscInt            *coneOrientations;  /* Orientation of each cone point, means cone traveral should start on point 'o', and if negative start on -(o+1) and go in reverse */
179:   PetscSection         supportSection;    /* Layout of cones (inedges for DAG) */
180:   PetscInt             maxSupportSize;    /* Cached for fast lookup */
181:   PetscInt            *supports;          /* Cone for each point */
182:   PetscBool            refinementUniform; /* Flag for uniform cell refinement */
183:   PetscReal            refinementLimit;   /* Maximum volume for refined cell */
184:   PetscErrorCode     (*refinementFunc)(const PetscReal [], PetscReal *); /* Function giving the maximum volume for refined cell */
185:   PetscInt             overlap;           /* Overlap of the partitions as passed to DMPlexDistribute() or DMPlexDistributeOverlap() */
186:   DMPlexInterpolatedFlag interpolated;
187:   DMPlexInterpolatedFlag interpolatedCollective;

189:   PetscInt            *facesTmp;          /* Work space for faces operation */

191:   /* Hierarchy */
192:   DMPlexCellRefinerType cellRefiner;       /* Strategy for refining cells */
193:   PetscBool             regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */

195:   /* Generation */
196:   char                *tetgenOpts;
197:   char                *triangleOpts;
198:   PetscPartitioner     partitioner;
199:   PetscBool            partitionBalance;  /* Evenly divide partition overlap when distributing */
200:   PetscBool            remeshBd;

202:   /* Submesh */
203:   DMLabel              subpointMap;       /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */

205:   /* Labels and numbering */
206:   PetscObjectState     depthState;        /* State of depth label, so that we can determine if a user changes it */
207:   PetscObjectState     celltypeState;     /* State of celltype label, so that we can determine if a user changes it */
208:   IS                   globalVertexNumbers;
209:   IS                   globalCellNumbers;

211:   /* Constraints */
212:   PetscSection         anchorSection;      /* maps constrained points to anchor points */
213:   IS                   anchorIS;           /* anchors indexed by the above section */
214:   PetscErrorCode     (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
215:   PetscErrorCode     (*computeanchormatrix)(DM,PetscSection,PetscSection,Mat);

217:   /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
218:   PetscSection         parentSection;     /* dof == 1 if point has parent */
219:   PetscInt            *parents;           /* point to parent */
220:   PetscInt            *childIDs;          /* point to child ID */
221:   PetscSection         childSection;      /* inverse of parent section */
222:   PetscInt            *children;          /* point to children */
223:   DM                   referenceTree;     /* reference tree to which child ID's refer */
224:   PetscErrorCode      (*getchildsymmetry)(DM,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*);

226:   /* MATIS support */
227:   PetscSection         subdomainSection;

229:   /* Adjacency */
230:   PetscBool            useAnchors;        /* Replace constrained points with their anchors in adjacency lists */
231:   PetscErrorCode      (*useradjacency)(DM,PetscInt,PetscInt*,PetscInt[],void*); /* User callback for adjacency */
232:   void                *useradjacencyctx;  /* User context for callback */

234:   /* Projection */
235:   PetscInt             maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */

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

241:   /* Geometry */
242:   PetscReal            minradius;         /* Minimum distance from cell centroid to face */
243:   PetscBool            useHashLocation;   /* Use grid hashing for point location */
244:   PetscGridHash        lbox;              /* Local box for searching */

246:   /* Neighbors */
247:   PetscMPIInt*         neighbors;

249:   /* Debugging */
250:   PetscBool            printSetValues;
251:   PetscInt             printFEM;
252:   PetscInt             printL2;
253:   PetscReal            printTol;
254: } DM_Plex;

256: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
257: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
258: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec,PetscViewer);
259: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
260: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
261: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec,PetscViewer);
262: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
263: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
264: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM,PetscViewer);
265: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject,PetscViewer);
266: #if defined(PETSC_HAVE_HDF5)
267: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
268: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
269: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
270: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
271: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
272: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
273: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
274: #endif

276: PETSC_INTERN PetscErrorCode DMPlexClosurePoints_Private(DM,PetscInt,const PetscInt[],IS*);
277: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *, DM);
278: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
279: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM []);
280: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
281: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM []);
282: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, DMLabel, DM *);
283: PETSC_INTERN PetscErrorCode DMAdaptMetric_Plex(DM, Vec, DMLabel, DM *);
284: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
285: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
286: 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);
287: 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);
288: 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);
289: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
290: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *);
291: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
292: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);

294: PETSC_INTERN PetscErrorCode DMPlexBuildFromCellList_Internal(DM, PetscInt, PetscInt, PetscInt, PetscInt, const int[], PetscBool);
295: PETSC_INTERN PetscErrorCode DMPlexBuildFromCellList_Parallel_Internal(DM, PetscInt, PetscInt, PetscInt, PetscInt, const int[], PetscBool, PetscSF *);
296: PETSC_INTERN PetscErrorCode DMPlexBuildCoordinates_Internal(DM, PetscInt, PetscInt, PetscInt, const double[]);
297: PETSC_INTERN PetscErrorCode DMPlexBuildCoordinates_Parallel_Internal(DM, PetscInt, PetscInt, PetscInt, PetscSF, const PetscReal[]);
298: PETSC_INTERN PetscErrorCode DMPlexLoadLabels_HDF5_Internal(DM, PetscViewer);
299: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
300: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
301: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM, PetscViewer);
302: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
303: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
304: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
305: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
306: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
307: /* TODO Make these INTERN */
308: PETSC_EXTERN PetscErrorCode DMPlexView_ExodusII_Internal(DM, int, PetscInt);
309: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec, int, int);
310: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec, int, int);
311: PETSC_EXTERN PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec, int, int);
312: PETSC_EXTERN PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec, int, int);
313: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM,PetscInt,PetscInt,PetscInt*);
314: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
315: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,DMPolytopeType,const PetscInt[],PetscInt*,const DMPolytopeType*[],const PetscInt*[],const PetscInt*[]);
316: PETSC_INTERN PetscErrorCode DMPlexRestoreRawFaces_Internal(DM,DMPolytopeType,const PetscInt[],PetscInt*,const DMPolytopeType*[],const PetscInt*[],const PetscInt*[]);
317: PETSC_INTERN PetscErrorCode CellRefinerInCellTest_Internal(DMPolytopeType, const PetscReal[], PetscBool *);
318: PETSC_INTERN PetscErrorCode DMPlexComputeCellType_Internal(DM, PetscInt, PetscInt, DMPolytopeType *);
319: PETSC_INTERN PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType, PetscInt *[], PetscInt *[]);
320: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], const PetscScalar[], InsertMode);
321: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
322: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
323: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM,DM,const char *,DM*);
324: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM,DM,PetscSF,PetscInt *,Mat);
325: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM,DM,PetscSF,PetscInt *,Mat);
326: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM,PetscSection,PetscInt,PetscInt,const PetscInt[],const PetscInt ***,const PetscScalar[],PetscInt*,PetscInt*,PetscInt*[],PetscScalar*[],PetscInt[],PetscBool);
327: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt []);
328: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt [],PetscBool,PetscInt,PetscInt []);
329: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM,PetscInt,const PetscScalar [],PetscInt,PetscInt *);
330: /* these two are PETSC_EXTERN just because of src/dm/impls/plex/tests/ex18.c */
331: PETSC_EXTERN PetscErrorCode DMPlexOrientCell_Internal(DM,PetscInt,PetscInt,PetscBool);
332: PETSC_EXTERN PetscErrorCode DMPlexOrientInterface_Internal(DM);

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

337: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
338: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
339: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, DMLabel, DM *);
340: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, DMLabel, DM *);
341: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM, Mat*);

343: PETSC_INTERN PetscErrorCode DMPlexGetOverlap_Plex(DM, PetscInt *);

345: /* invert dihedral symmetry: return a^-1,
346:  * using the representation described in
347:  * DMPlexGetConeOrientation() */
348: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
349: {
350:   return (a <= 0) ? a : (N - a);
351: }

353: /* invert dihedral symmetry: return b * a,
354:  * using the representation described in
355:  * DMPlexGetConeOrientation() */
356: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
357: {
358:   if (!N) return 0;
359:   return  (a >= 0) ?
360:          ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
361:          ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
362: }

364: /* swap dihedral symmetries: return b * a^-1,
365:  * using the representation described in
366:  * DMPlexGetConeOrientation() */
367: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
368: {
369:   return DihedralCompose(N,DihedralInvert(N,a),b);
370: }

372: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, IS , PetscReal, Vec, Vec, PetscReal, Vec, void *);
373: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
374: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);

376: /* Matvec with A in row-major storage, x and y can be aliased */
377: PETSC_STATIC_INLINE void DMPlex_Mult2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
378: {
379:   PetscScalar z[2];
380:   z[0] = x[0]; z[1] = x[ldx];
381:   y[0]   = A[0]*z[0] + A[1]*z[1];
382:   y[ldx] = A[2]*z[0] + A[3]*z[1];
383:   (void)PetscLogFlops(6.0);
384: }
385: PETSC_STATIC_INLINE void DMPlex_Mult3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
386: {
387:   PetscScalar z[3];
388:   z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
389:   y[0]     = A[0]*z[0] + A[1]*z[1] + A[2]*z[2];
390:   y[ldx]   = A[3]*z[0] + A[4]*z[1] + A[5]*z[2];
391:   y[ldx*2] = A[6]*z[0] + A[7]*z[1] + A[8]*z[2];
392:   (void)PetscLogFlops(15.0);
393: }
394: PETSC_STATIC_INLINE void DMPlex_MultTranspose2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
395: {
396:   PetscScalar z[2];
397:   z[0] = x[0]; z[1] = x[ldx];
398:   y[0]   = A[0]*z[0] + A[2]*z[1];
399:   y[ldx] = A[1]*z[0] + A[3]*z[1];
400:   (void)PetscLogFlops(6.0);
401: }
402: PETSC_STATIC_INLINE void DMPlex_MultTranspose3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
403: {
404:   PetscScalar z[3];
405:   z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
406:   y[0]     = A[0]*z[0] + A[3]*z[1] + A[6]*z[2];
407:   y[ldx]   = A[1]*z[0] + A[4]*z[1] + A[7]*z[2];
408:   y[ldx*2] = A[2]*z[0] + A[5]*z[1] + A[8]*z[2];
409:   (void)PetscLogFlops(15.0);
410: }
411: PETSC_STATIC_INLINE void DMPlex_Mult2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
412: {
413:   PetscScalar z[2];
414:   z[0] = x[0]; z[1] = x[ldx];
415:   y[0]   = A[0]*z[0] + A[1]*z[1];
416:   y[ldx] = A[2]*z[0] + A[3]*z[1];
417:   (void)PetscLogFlops(6.0);
418: }
419: PETSC_STATIC_INLINE void DMPlex_Mult3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
420: {
421:   PetscScalar z[3];
422:   z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
423:   y[0]     = A[0]*z[0] + A[1]*z[1] + A[2]*z[2];
424:   y[ldx]   = A[3]*z[0] + A[4]*z[1] + A[5]*z[2];
425:   y[ldx*2] = A[6]*z[0] + A[7]*z[1] + A[8]*z[2];
426:   (void)PetscLogFlops(15.0);
427: }
428: PETSC_STATIC_INLINE void DMPlex_MultTranspose2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
429: {
430:   PetscScalar z[2];
431:   z[0] = x[0]; z[1] = x[ldx];
432:   y[0]   = A[0]*z[0] + A[2]*z[1];
433:   y[ldx] = A[1]*z[0] + A[3]*z[1];
434:   (void)PetscLogFlops(6.0);
435: }
436: PETSC_STATIC_INLINE void DMPlex_MultTranspose3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
437: {
438:   PetscScalar z[3];
439:   z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
440:   y[0]     = A[0]*z[0] + A[3]*z[1] + A[6]*z[2];
441:   y[ldx]   = A[1]*z[0] + A[4]*z[1] + A[7]*z[2];
442:   y[ldx*2] = A[2]*z[0] + A[5]*z[1] + A[8]*z[2];
443:   (void)PetscLogFlops(15.0);
444: }

446: PETSC_STATIC_INLINE void DMPlex_MatMult2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
447: {
448:   PetscInt j;
449:   for (j = 0; j < n; ++j) {
450:     PetscScalar z[2];
451:     z[0] = B[0+j]; z[1] = B[1*ldb+j];
452:     DMPlex_Mult2D_Internal(A, 1, z, z);
453:     C[0+j] = z[0]; C[1*ldb+j] = z[1];
454:   }
455:   (void)PetscLogFlops(8.0*n);
456: }
457: PETSC_STATIC_INLINE void DMPlex_MatMult3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
458: {
459:   PetscInt j;
460:   for (j = 0; j < n; ++j) {
461:     PetscScalar z[3];
462:     z[0] = B[0+j]; z[1] = B[1*ldb+j]; z[2] = B[2*ldb+j];
463:     DMPlex_Mult3D_Internal(A, 1, z, z);
464:     C[0+j] = z[0]; C[1*ldb+j] = z[1]; C[2*ldb+j] = z[2];
465:   }
466:   (void)PetscLogFlops(8.0*n);
467: }
468: PETSC_STATIC_INLINE void DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
469: {
470:   PetscInt j;
471:   for (j = 0; j < n; ++j) {
472:     PetscScalar z[2];
473:     z[0] = B[0+j]; z[1] = B[1*ldb+j];
474:     DMPlex_MultTranspose2D_Internal(A, 1, z, z);
475:     C[0+j] = z[0]; C[1*ldb+j] = z[1];
476:   }
477:   (void)PetscLogFlops(8.0*n);
478: }
479: PETSC_STATIC_INLINE void DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
480: {
481:   PetscInt j;
482:   for (j = 0; j < n; ++j) {
483:     PetscScalar z[3];
484:     z[0] = B[0+j]; z[1] = B[1*ldb+j]; z[2] = B[2*ldb+j];
485:     DMPlex_MultTranspose3D_Internal(A, 1, z, z);
486:     C[0+j] = z[0]; C[1*ldb+j] = z[1]; C[2*ldb+j] = z[2];
487:   }
488:   (void)PetscLogFlops(8.0*n);
489: }

491: PETSC_STATIC_INLINE void DMPlex_MatMultLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
492: {
493:   PetscInt j;
494:   for (j = 0; j < m; ++j) {
495:     DMPlex_MultTranspose2D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
496:   }
497:   (void)PetscLogFlops(8.0*m);
498: }
499: PETSC_STATIC_INLINE void DMPlex_MatMultLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
500: {
501:   PetscInt j;
502:   for (j = 0; j < m; ++j) {
503:     DMPlex_MultTranspose3D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
504:   }
505:   (void)PetscLogFlops(8.0*m);
506: }
507: PETSC_STATIC_INLINE void DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
508: {
509:   PetscInt j;
510:   for (j = 0; j < m; ++j) {
511:     DMPlex_Mult2D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
512:   }
513:   (void)PetscLogFlops(8.0*m);
514: }
515: PETSC_STATIC_INLINE void DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
516: {
517:   PetscInt j;
518:   for (j = 0; j < m; ++j) {
519:     DMPlex_Mult3D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
520:   }
521:   (void)PetscLogFlops(8.0*m);
522: }

524: PETSC_STATIC_INLINE void DMPlex_Transpose2D_Internal(PetscScalar A[])
525: {
526:   PetscScalar tmp;
527:   tmp = A[1]; A[1] = A[2]; A[2] = tmp;
528: }
529: PETSC_STATIC_INLINE void DMPlex_Transpose3D_Internal(PetscScalar A[])
530: {
531:   PetscScalar tmp;
532:   tmp = A[1]; A[1] = A[3]; A[3] = tmp;
533:   tmp = A[2]; A[2] = A[6]; A[6] = tmp;
534:   tmp = A[5]; A[5] = A[7]; A[7] = tmp;
535: }

537: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
538: {
539:   const PetscReal invDet = 1.0/detJ;

541:   invJ[0] =  invDet*J[3];
542:   invJ[1] = -invDet*J[1];
543:   invJ[2] = -invDet*J[2];
544:   invJ[3] =  invDet*J[0];
545:   (void)PetscLogFlops(5.0);
546: }

548: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
549: {
550:   const PetscReal invDet = 1.0/detJ;

552:   invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
553:   invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
554:   invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
555:   invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
556:   invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
557:   invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
558:   invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
559:   invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
560:   invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
561:   (void)PetscLogFlops(37.0);
562: }

564: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
565: {
566:   *detJ = J[0]*J[3] - J[1]*J[2];
567:   (void)PetscLogFlops(3.0);
568: }

570: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
571: {
572:   *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
573:            J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
574:            J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
575:   (void)PetscLogFlops(12.0);
576: }

578: PETSC_STATIC_INLINE void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
579: {
580:   *detJ = PetscRealPart(J[0])*PetscRealPart(J[3]) - PetscRealPart(J[1])*PetscRealPart(J[2]);
581:   (void)PetscLogFlops(3.0);
582: }

584: PETSC_STATIC_INLINE void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
585: {
586:   *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])) +
587:            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])) +
588:            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])));
589:   (void)PetscLogFlops(12.0);
590: }

592: PETSC_STATIC_INLINE void DMPlex_WaxpyD_Internal(PetscInt dim, PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w) {PetscInt d; for (d = 0; d < dim; ++d) w[d] = a*x[d] + y[d];}

594: PETSC_STATIC_INLINE PetscReal DMPlex_DotD_Internal(PetscInt dim, const PetscScalar *x, const PetscReal *y) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += PetscRealPart(x[d])*y[d]; return sum;}

596: PETSC_STATIC_INLINE PetscReal DMPlex_DotRealD_Internal(PetscInt dim, const PetscReal *x, const PetscReal *y) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*y[d]; return sum;}

598: PETSC_STATIC_INLINE PetscReal DMPlex_NormD_Internal(PetscInt dim, const PetscReal *x) {PetscReal sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*x[d]; return PetscSqrtReal(sum);}

600: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Translate_Private(PetscInt ornt, PetscInt *start, PetscBool *reverse)
601: {
603:   *reverse = (ornt < 0) ? PETSC_TRUE : PETSC_FALSE;
604:   *start = *reverse ? -(ornt+1) : ornt;
605:   return(0);
606: }

608: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Combine_Private(PetscInt coneSize, PetscInt origStart, PetscBool origReverse, PetscInt rotateStart, PetscBool rotateReverse, PetscInt *newStart, PetscBool *newReverse)
609: {
611:   *newReverse = (origReverse == rotateReverse) ? PETSC_FALSE : PETSC_TRUE;
612:   *newStart = rotateReverse ? (coneSize + rotateStart - origStart) : (coneSize + origStart - rotateStart);
613:   *newStart %= coneSize;
614:   return(0);
615: }

617: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_TranslateBack_Private(PetscInt coneSize, PetscInt start, PetscBool reverse, PetscInt *ornt)
618: {
620:   if (coneSize < 3) {
621:     /* edges just get flipped if start == 1 regardless direction */
622:     *ornt = start ? -2 : 0;
623:   } else {
624:     *ornt = reverse ? -(start+1) : start;
625:   }
626:   return(0);
627: }

629: PETSC_STATIC_INLINE PetscErrorCode DMPlexFixFaceOrientations_Permute_Private(PetscInt n, const PetscInt arr[], PetscInt start, PetscBool reverse, PetscInt newarr[])
630: {
631:   PetscInt i;

634:   if (reverse) {for (i=0; i<n; i++) newarr[i] = arr[(n+start-i)%n];}
635:   else         {for (i=0; i<n; i++) newarr[i] = arr[(start+i)%n];}
636:   return(0);
637: }

639: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM,PetscInt,PetscInt,PetscDualSpace *);
640: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection,PetscBool,PetscInt,PetscInt,PetscInt *,PetscBool,const PetscInt[],const PetscInt[],PetscInt[]);
641: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection,PetscBool,PetscInt,PetscInt,PetscInt[],PetscBool,const PetscInt***,PetscInt,const PetscInt[],PetscInt[]);
642: PETSC_INTERN PetscErrorCode DMPlexGetCompressedClosure(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscSection *, IS *, const PetscInt **);
643: PETSC_INTERN PetscErrorCode DMPlexRestoreCompressedClosure(DM, PetscSection, PetscInt, PetscInt *, PetscInt **, PetscSection *, IS *, const PetscInt **);

645: PETSC_EXTERN PetscErrorCode DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
646: PETSC_EXTERN PetscErrorCode DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
647: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
648: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
649: PETSC_INTERN PetscErrorCode DMCreateSubDomainDM_Plex(DM,DMLabel,PetscInt,IS*,DM*);
650: PETSC_INTERN PetscErrorCode DMPlexBasisTransformPoint_Internal(DM, DM, Vec, PetscInt, PetscBool[], PetscBool, PetscScalar *);
651: PETSC_EXTERN PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM, DM, Vec, PetscInt, PetscBool, PetscInt, PetscScalar *);
652: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscReal *, PetscReal *, void *);
653: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApply_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscScalar *, PetscScalar *, void *);
654: PETSC_INTERN PetscErrorCode DMCreateNeumannOverlap_Plex(DM, IS*, Mat*, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void **);

656: #endif /* _PLEXIMPL_H */