Actual source code: dmpleximpl.h

  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: #if defined(PETSC_HAVE_EXODUSII)
 11: #include <exodusII.h>
 12: #endif

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

 49: typedef struct _n_PlexGeneratorFunctionList *PlexGeneratorFunctionList;
 50: struct _n_PlexGeneratorFunctionList {
 51:   PetscErrorCode    (*generate)(DM, PetscBool, DM*);
 52:   PetscErrorCode    (*refine)(DM, PetscReal*, DM*);
 53:   PetscErrorCode    (*adaptlabel)(DM, DMLabel, DM*);
 54:   char              *name;
 55:   PetscInt          dim;
 56:   PlexGeneratorFunctionList next;
 57: };

 59: /* Utility struct to store the contents of a Fluent file in memory */
 60: typedef struct {
 61:   int          index;    /* Type of section */
 62:   unsigned int zoneID;
 63:   unsigned int first;
 64:   unsigned int last;
 65:   int          type;
 66:   int          nd;       /* Either ND or element-type */
 67:   void        *data;
 68: } FluentSection;

 70: struct _PetscGridHash {
 71:   PetscInt     dim;
 72:   PetscReal    lower[3];    /* The lower-left corner */
 73:   PetscReal    upper[3];    /* The upper-right corner */
 74:   PetscReal    extent[3];   /* The box size */
 75:   PetscReal    h[3];        /* The subbox size */
 76:   PetscInt     n[3];        /* The number of subboxes */
 77:   PetscSection cellSection; /* Offsets for cells in each subbox*/
 78:   IS           cells;       /* List of cells in each subbox */
 79:   DMLabel      cellsSparse; /* Sparse storage for cell map */
 80: };

 82: /* Point Numbering in Plex:

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

 86:    First Stratum:  Cells [height 0]
 87:    Second Stratum: Vertices [depth 0]
 88:    Third Stratum:  Faces [height 1]
 89:    Fourth Stratum: Edges [depth 1]

 91:    We do this so that the numbering of a cell-vertex mesh does not change after interpolation. Within a given stratum,
 92:    we allow additional segregation of by cell type.
 93: */
 94: typedef struct {
 95:   PetscInt             refct;

 97:   PetscSection         coneSection;       /* Layout of cones (inedges for DAG) */
 98:   PetscInt             maxConeSize;       /* Cached for fast lookup */
 99:   PetscInt            *cones;             /* Cone for each point */
100:   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 */
101:   PetscSection         supportSection;    /* Layout of cones (inedges for DAG) */
102:   PetscInt             maxSupportSize;    /* Cached for fast lookup */
103:   PetscInt            *supports;          /* Cone for each point */
104:   PetscBool            refinementUniform; /* Flag for uniform cell refinement */
105:   char                *transformType;     /* Type of transform for uniform cell refinement */
106:   PetscReal            refinementLimit;   /* Maximum volume for refined cell */
107:   PetscErrorCode     (*refinementFunc)(const PetscReal [], PetscReal *); /* Function giving the maximum volume for refined cell */
108:   PetscInt             overlap;           /* Overlap of the partitions as passed to DMPlexDistribute() or DMPlexDistributeOverlap() */
109:   DMPlexInterpolatedFlag interpolated;
110:   DMPlexInterpolatedFlag interpolatedCollective;

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

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

117:   /* Generation */
118:   char                *tetgenOpts;
119:   char                *triangleOpts;
120:   PetscPartitioner     partitioner;
121:   PetscBool            partitionBalance;  /* Evenly divide partition overlap when distributing */
122:   PetscBool            remeshBd;

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

129:   /* Labels and numbering */
130:   PetscObjectState     depthState;        /* State of depth label, so that we can determine if a user changes it */
131:   PetscObjectState     celltypeState;     /* State of celltype label, so that we can determine if a user changes it */
132:   IS                   globalVertexNumbers;
133:   IS                   globalCellNumbers;

135:   /* Constraints */
136:   PetscSection         anchorSection;      /* maps constrained points to anchor points */
137:   IS                   anchorIS;           /* anchors indexed by the above section */
138:   PetscErrorCode     (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
139:   PetscErrorCode     (*computeanchormatrix)(DM,PetscSection,PetscSection,Mat);

141:   /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
142:   PetscSection         parentSection;     /* dof == 1 if point has parent */
143:   PetscInt            *parents;           /* point to parent */
144:   PetscInt            *childIDs;          /* point to child ID */
145:   PetscSection         childSection;      /* inverse of parent section */
146:   PetscInt            *children;          /* point to children */
147:   DM                   referenceTree;     /* reference tree to which child ID's refer */
148:   PetscErrorCode      (*getchildsymmetry)(DM,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*);

150:   /* MATIS support */
151:   PetscSection         subdomainSection;

153:   /* Adjacency */
154:   PetscBool            useAnchors;        /* Replace constrained points with their anchors in adjacency lists */
155:   PetscErrorCode      (*useradjacency)(DM,PetscInt,PetscInt*,PetscInt[],void*); /* User callback for adjacency */
156:   void                *useradjacencyctx;  /* User context for callback */

158:   /* Projection */
159:   PetscInt             maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
160:   PetscInt             activePoint;         /* current active point in iteration */

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

166:   /* Geometry */
167:   PetscBool            ignoreModel;       /* Ignore the geometry model during refinement */
168:   PetscReal            minradius;         /* Minimum distance from cell centroid to face */
169:   PetscBool            useHashLocation;   /* Use grid hashing for point location */
170:   PetscGridHash        lbox;              /* Local box for searching */
171:   void               (*coordFunc)(PetscInt, PetscInt, PetscInt, /* Function used to remap newly introduced vertices */
172:                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
173:                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
174:                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);

176:   /* Neighbors */
177:   PetscMPIInt*         neighbors;

179:   /* Debugging */
180:   PetscBool            printSetValues;
181:   PetscInt             printFEM;
182:   PetscInt             printL2;
183:   PetscReal            printTol;
184: } DM_Plex;

186: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
187: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
188: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec,PetscViewer);
189: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
190: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
191: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec,PetscViewer);
192: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
193: PETSC_INTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
194: PETSC_INTERN PetscErrorCode DMPlexView_GLVis(DM,PetscViewer);
195: PETSC_INTERN PetscErrorCode DMSetUpGLVisViewer_Plex(PetscObject,PetscViewer);
196: #if defined(PETSC_HAVE_HDF5)
197: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
198: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
199: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
200: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
201: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
202: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
203: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
204: PETSC_INTERN PetscErrorCode DMPlexTopologyView_HDF5_Internal(DM, IS, PetscViewer);
205: PETSC_INTERN PetscErrorCode DMPlexCoordinatesView_HDF5_Internal(DM, PetscViewer);
206: PETSC_INTERN PetscErrorCode DMPlexLabelsView_HDF5_Internal(DM, IS, PetscViewer);
207: PETSC_INTERN PetscErrorCode DMPlexSectionView_HDF5_Internal(DM, PetscViewer, DM);
208: PETSC_INTERN PetscErrorCode DMPlexGlobalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
209: PETSC_INTERN PetscErrorCode DMPlexLocalVectorView_HDF5_Internal(DM, PetscViewer, DM, Vec);
210: PETSC_INTERN PetscErrorCode DMPlexTopologyLoad_HDF5_Internal(DM, PetscViewer, PetscSF*);
211: PETSC_INTERN PetscErrorCode DMPlexCoordinatesLoad_HDF5_Internal(DM, PetscViewer);
212: PETSC_INTERN PetscErrorCode DMPlexLabelsLoad_HDF5_Internal(DM, PetscViewer);
213: PETSC_INTERN PetscErrorCode DMPlexSectionLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, PetscSF*, PetscSF*);
214: PETSC_INTERN PetscErrorCode DMPlexVecLoad_HDF5_Internal(DM, PetscViewer, DM, PetscSF, Vec);
215: PETSC_INTERN PetscErrorCode DMPlexView_HDF5_Internal(DM, PetscViewer);
216: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Internal(DM, PetscViewer);
217: PETSC_INTERN PetscErrorCode DMPlexLoad_HDF5_Xdmf_Internal(DM, PetscViewer);
218: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Internal(Vec, PetscViewer);
219: PETSC_INTERN PetscErrorCode VecView_Plex_HDF5_Native_Internal(Vec, PetscViewer);
220: PETSC_INTERN PetscErrorCode VecView_Plex_Local_HDF5_Internal(Vec, PetscViewer);
221: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Internal(Vec, PetscViewer);
222: PETSC_INTERN PetscErrorCode VecLoad_Plex_HDF5_Native_Internal(Vec, PetscViewer);
223: #endif

225: PETSC_INTERN PetscErrorCode DMPlexVecGetClosureAtDepth_Internal(DM, PetscSection, Vec, PetscInt, PetscInt, PetscInt *, PetscScalar *[]);
226: PETSC_INTERN PetscErrorCode DMPlexClosurePoints_Private(DM,PetscInt,const PetscInt[],IS*);
227: PETSC_INTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *, DM);
228: PETSC_INTERN PetscErrorCode DMCoarsen_Plex(DM, MPI_Comm, DM *);
229: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Plex(DM, PetscInt, DM []);
230: PETSC_INTERN PetscErrorCode DMRefine_Plex(DM, MPI_Comm, DM *);
231: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Plex(DM, PetscInt, DM []);
232: PETSC_INTERN PetscErrorCode DMAdaptLabel_Plex(DM, DMLabel, DM *);
233: PETSC_INTERN PetscErrorCode DMAdaptMetric_Plex(DM, Vec, DMLabel, DM *);
234: PETSC_INTERN PetscErrorCode DMPlexInsertBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
235: PETSC_INTERN PetscErrorCode DMPlexInsertTimeDerivativeBoundaryValues_Plex(DM, PetscBool, Vec, PetscReal, Vec, Vec, Vec);
236: PETSC_INTERN PetscErrorCode DMProjectFunctionLocal_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,InsertMode,Vec);
237: 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);
238: 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);
239: 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);
240: 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);
241: PETSC_INTERN PetscErrorCode DMComputeL2Diff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
242: PETSC_INTERN PetscErrorCode DMComputeL2GradientDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[], const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,const PetscReal [],PetscReal *);
243: PETSC_INTERN PetscErrorCode DMComputeL2FieldDiff_Plex(DM,PetscReal,PetscErrorCode(**)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *),void **,Vec,PetscReal *);
244: PETSC_INTERN PetscErrorCode DMLocatePoints_Plex(DM, Vec, DMPointLocationType, PetscSF);

246: #if defined(PETSC_HAVE_EXODUSII)
247: PETSC_EXTERN PetscErrorCode DMView_PlexExodusII(DM, PetscViewer);
248: PETSC_INTERN PetscErrorCode VecView_PlexExodusII_Internal(Vec, PetscViewer);
249: PETSC_INTERN PetscErrorCode VecLoad_PlexExodusII_Internal(Vec, PetscViewer);
250: PETSC_INTERN PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec, int, int, int);
251: PETSC_INTERN PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec, int, int, int);
252: PETSC_INTERN PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec, int, int, int);
253: PETSC_INTERN PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec, int, int, int);
254: PETSC_INTERN PetscErrorCode EXOGetVarIndex_Internal(int, ex_entity_type, const char[], int*);
255: #endif
256: PETSC_INTERN PetscErrorCode DMPlexVTKGetCellType_Internal(DM,PetscInt,PetscInt,PetscInt*);
257: PETSC_INTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
258: PETSC_INTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,DMPolytopeType,const PetscInt[],PetscInt*,const DMPolytopeType*[],const PetscInt*[],const PetscInt*[]);
259: PETSC_INTERN PetscErrorCode DMPlexRestoreRawFaces_Internal(DM,DMPolytopeType,const PetscInt[],PetscInt*,const DMPolytopeType*[],const PetscInt*[],const PetscInt*[]);
260: PETSC_INTERN PetscErrorCode DMPlexComputeCellType_Internal(DM, PetscInt, PetscInt, DMPolytopeType *);
261: PETSC_INTERN PetscErrorCode DMPlexCreateCellTypeOrder_Internal(DMPolytopeType, PetscInt *[], PetscInt *[]);
262: PETSC_INTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, PetscInt, const PetscInt[], DMLabel, PetscInt, const PetscScalar[], InsertMode);
263: PETSC_INTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
264: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_SetTree(DM, PetscSection, PetscInt[], PetscInt[]);
265: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM,DM,const char *,DM*);
266: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM,DM,PetscSF,PetscInt *,Mat);
267: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM,DM,PetscSF,PetscInt *,Mat);
268: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM,PetscSection,PetscInt,PetscInt,const PetscInt[],const PetscInt ***,const PetscScalar[],PetscInt*,PetscInt*,PetscInt*[],PetscScalar*[],PetscInt[],PetscBool);
269: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt []);
270: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt [],PetscBool,PetscInt,PetscInt []);
271: PETSC_INTERN PetscErrorCode DMPlexLocatePoint_Internal(DM,PetscInt,const PetscScalar [],PetscInt,PetscInt *);
272: /* these two are PETSC_EXTERN just because of src/dm/impls/plex/tests/ex18.c */
273: PETSC_EXTERN PetscErrorCode DMPlexOrientInterface_Internal(DM);

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

278: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
279: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
280: PETSC_INTERN PetscErrorCode DMPlexRefine_Internal(DM, DMLabel, DM *);
281: PETSC_INTERN PetscErrorCode DMPlexCoarsen_Internal(DM, DMLabel, DM *);
282: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM, Mat*);

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

286: #if 1
287: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
288: {
289:   return (a <= 0) ? a : (N - a);
290: }

292: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
293: {
294:   if (!N) return 0;
295:   return  (a >= 0) ?
296:          ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
297:          ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
298: }

300: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
301: {
302:   return DihedralCompose(N,DihedralInvert(N,a),b);
303: }
304: #else
305: /* TODO
306:    This is a reimplementation of the tensor dihedral symmetries using the new orientations.
307:    These should be turned on when we convert to new-style orientations in p4est.
308: */
309: /* invert dihedral symmetry: return a^-1,
310:  * using the representation described in
311:  * DMPlexGetConeOrientation() */
312: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
313: {
314:   switch (N) {
315:     case 0: return 0;
316:     case 2: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, 0, a);
317:     case 4: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, 0, a);
318:     case 8: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, 0, a);
319:     default: SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralInvert()");
320:   }
321:   return 0;
322: }

324: /* compose dihedral symmetry: return b * a,
325:  * using the representation described in
326:  * DMPlexGetConeOrientation() */
327: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
328: {
329:   switch (N) {
330:     case 0: return 0;
331:     case 2: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
332:     case 4: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
333:     case 8: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
334:     default: SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
335:   }
336:   return 0;
337: }

339: /* swap dihedral symmetries: return b * a^-1,
340:  * using the representation described in
341:  * DMPlexGetConeOrientation() */
342: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
343: {
344:   switch (N) {
345:     case 0: return 0;
346:     case 2: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_SEGMENT, b, a);
347:     case 4: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_QUADRILATERAL, b, a);
348:     case 8: return DMPolytopeTypeComposeOrientationInv(DM_POLYTOPE_HEXAHEDRON, b, a);
349:     default: SETERRABORT(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid celltype for DihedralCompose()");
350:   }
351:   return 0;
352: }
353: #endif
354: PETSC_EXTERN PetscInt       DMPolytopeConvertNewOrientation_Internal(DMPolytopeType, PetscInt);
355: PETSC_EXTERN PetscInt       DMPolytopeConvertOldOrientation_Internal(DMPolytopeType, PetscInt);
356: PETSC_EXTERN PetscErrorCode DMPlexConvertOldOrientations_Internal(DM);

358: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, PetscFormKey, IS, PetscReal, Vec, Vec, PetscReal, Vec, void *);
359: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Hybrid_Internal(DM, PetscFormKey[], IS, PetscReal, Vec, Vec, PetscReal, Vec, void *);
360: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, PetscFormKey, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
361: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Hybrid_Internal(DM, PetscFormKey[], IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
362: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Action_Internal(DM, PetscFormKey, IS, PetscReal, PetscReal, Vec, Vec, Vec, Vec, void *);
363: PETSC_EXTERN PetscErrorCode DMPlexReconstructGradients_Internal(DM, PetscFV, PetscInt, PetscInt, Vec, Vec, Vec, Vec);

365: /* Matvec with A in row-major storage, x and y can be aliased */
366: PETSC_STATIC_INLINE void DMPlex_Mult2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
367: {
368:   PetscScalar z[2];
369:   z[0] = x[0]; z[1] = x[ldx];
370:   y[0]   = A[0]*z[0] + A[1]*z[1];
371:   y[ldx] = A[2]*z[0] + A[3]*z[1];
372:   (void)PetscLogFlops(6.0);
373: }
374: PETSC_STATIC_INLINE void DMPlex_Mult3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
375: {
376:   PetscScalar z[3];
377:   z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
378:   y[0]     = A[0]*z[0] + A[1]*z[1] + A[2]*z[2];
379:   y[ldx]   = A[3]*z[0] + A[4]*z[1] + A[5]*z[2];
380:   y[ldx*2] = A[6]*z[0] + A[7]*z[1] + A[8]*z[2];
381:   (void)PetscLogFlops(15.0);
382: }
383: PETSC_STATIC_INLINE void DMPlex_MultTranspose2D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
384: {
385:   PetscScalar z[2];
386:   z[0] = x[0]; z[1] = x[ldx];
387:   y[0]   = A[0]*z[0] + A[2]*z[1];
388:   y[ldx] = A[1]*z[0] + A[3]*z[1];
389:   (void)PetscLogFlops(6.0);
390: }
391: PETSC_STATIC_INLINE void DMPlex_MultTranspose3D_Internal(const PetscScalar A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
392: {
393:   PetscScalar z[3];
394:   z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
395:   y[0]     = A[0]*z[0] + A[3]*z[1] + A[6]*z[2];
396:   y[ldx]   = A[1]*z[0] + A[4]*z[1] + A[7]*z[2];
397:   y[ldx*2] = A[2]*z[0] + A[5]*z[1] + A[8]*z[2];
398:   (void)PetscLogFlops(15.0);
399: }
400: PETSC_STATIC_INLINE void DMPlex_Mult2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
401: {
402:   PetscScalar z[2];
403:   z[0] = x[0]; z[1] = x[ldx];
404:   y[0]   = A[0]*z[0] + A[1]*z[1];
405:   y[ldx] = A[2]*z[0] + A[3]*z[1];
406:   (void)PetscLogFlops(6.0);
407: }
408: PETSC_STATIC_INLINE void DMPlex_Mult3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
409: {
410:   PetscScalar z[3];
411:   z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
412:   y[0]     = A[0]*z[0] + A[1]*z[1] + A[2]*z[2];
413:   y[ldx]   = A[3]*z[0] + A[4]*z[1] + A[5]*z[2];
414:   y[ldx*2] = A[6]*z[0] + A[7]*z[1] + A[8]*z[2];
415:   (void)PetscLogFlops(15.0);
416: }
417: PETSC_STATIC_INLINE void DMPlex_MultAdd2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
418: {
419:   PetscScalar z[2];
420:   z[0] = x[0]; z[1] = x[ldx];
421:   y[0]   += A[0]*z[0] + A[1]*z[1];
422:   y[ldx] += A[2]*z[0] + A[3]*z[1];
423:   (void)PetscLogFlops(6.0);
424: }
425: PETSC_STATIC_INLINE void DMPlex_MultAdd3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
426: {
427:   PetscScalar z[3];
428:   z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
429:   y[0]     += A[0]*z[0] + A[1]*z[1] + A[2]*z[2];
430:   y[ldx]   += A[3]*z[0] + A[4]*z[1] + A[5]*z[2];
431:   y[ldx*2] += A[6]*z[0] + A[7]*z[1] + A[8]*z[2];
432:   (void)PetscLogFlops(15.0);
433: }
434: PETSC_STATIC_INLINE void DMPlex_MultTransposeReal_Internal(const PetscReal A[], PetscInt m, PetscInt n, PetscInt ldx, const PetscScalar x[], PetscScalar y[])
435: {
436:   PetscScalar z[3];
437:   PetscInt    i, j;
438:   for (i = 0; i < m; ++i) z[i] = x[i*ldx];
439:   for (j = 0; j < n; ++j) {
440:     const PetscInt l = j*ldx;
441:     y[l] = 0;
442:     for (i = 0; i < m; ++i) {
443:       y[l] += A[j*n+i]*z[i];
444:     }
445:   }
446:   (void)PetscLogFlops(2*m*n);
447: }
448: PETSC_STATIC_INLINE void DMPlex_MultTranspose2DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
449: {
450:   PetscScalar z[2];
451:   z[0] = x[0]; z[1] = x[ldx];
452:   y[0]   = A[0]*z[0] + A[2]*z[1];
453:   y[ldx] = A[1]*z[0] + A[3]*z[1];
454:   (void)PetscLogFlops(6.0);
455: }
456: PETSC_STATIC_INLINE void DMPlex_MultTranspose3DReal_Internal(const PetscReal A[], PetscInt ldx, const PetscScalar x[], PetscScalar y[])
457: {
458:   PetscScalar z[3];
459:   z[0] = x[0]; z[1] = x[ldx]; z[2] = x[ldx*2];
460:   y[0]     = A[0]*z[0] + A[3]*z[1] + A[6]*z[2];
461:   y[ldx]   = A[1]*z[0] + A[4]*z[1] + A[7]*z[2];
462:   y[ldx*2] = A[2]*z[0] + A[5]*z[1] + A[8]*z[2];
463:   (void)PetscLogFlops(15.0);
464: }

466: PETSC_STATIC_INLINE void DMPlex_MatMult2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
467: {
468:   PetscInt j;
469:   for (j = 0; j < n; ++j) {
470:     PetscScalar z[2];
471:     z[0] = B[0+j]; z[1] = B[1*ldb+j];
472:     DMPlex_Mult2D_Internal(A, 1, z, z);
473:     C[0+j] = z[0]; C[1*ldb+j] = z[1];
474:   }
475:   (void)PetscLogFlops(8.0*n);
476: }
477: PETSC_STATIC_INLINE void DMPlex_MatMult3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
478: {
479:   PetscInt j;
480:   for (j = 0; j < n; ++j) {
481:     PetscScalar z[3];
482:     z[0] = B[0+j]; z[1] = B[1*ldb+j]; z[2] = B[2*ldb+j];
483:     DMPlex_Mult3D_Internal(A, 1, z, z);
484:     C[0+j] = z[0]; C[1*ldb+j] = z[1]; C[2*ldb+j] = z[2];
485:   }
486:   (void)PetscLogFlops(8.0*n);
487: }
488: PETSC_STATIC_INLINE void DMPlex_MatMultTranspose2D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
489: {
490:   PetscInt j;
491:   for (j = 0; j < n; ++j) {
492:     PetscScalar z[2];
493:     z[0] = B[0+j]; z[1] = B[1*ldb+j];
494:     DMPlex_MultTranspose2D_Internal(A, 1, z, z);
495:     C[0+j] = z[0]; C[1*ldb+j] = z[1];
496:   }
497:   (void)PetscLogFlops(8.0*n);
498: }
499: PETSC_STATIC_INLINE void DMPlex_MatMultTranspose3D_Internal(const PetscScalar A[], PetscInt n, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
500: {
501:   PetscInt j;
502:   for (j = 0; j < n; ++j) {
503:     PetscScalar z[3];
504:     z[0] = B[0+j]; z[1] = B[1*ldb+j]; z[2] = B[2*ldb+j];
505:     DMPlex_MultTranspose3D_Internal(A, 1, z, z);
506:     C[0+j] = z[0]; C[1*ldb+j] = z[1]; C[2*ldb+j] = z[2];
507:   }
508:   (void)PetscLogFlops(8.0*n);
509: }

511: PETSC_STATIC_INLINE void DMPlex_MatMultLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
512: {
513:   PetscInt j;
514:   for (j = 0; j < m; ++j) {
515:     DMPlex_MultTranspose2D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
516:   }
517:   (void)PetscLogFlops(8.0*m);
518: }
519: PETSC_STATIC_INLINE void DMPlex_MatMultLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
520: {
521:   PetscInt j;
522:   for (j = 0; j < m; ++j) {
523:     DMPlex_MultTranspose3D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
524:   }
525:   (void)PetscLogFlops(8.0*m);
526: }
527: PETSC_STATIC_INLINE void DMPlex_MatMultTransposeLeft2D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
528: {
529:   PetscInt j;
530:   for (j = 0; j < m; ++j) {
531:     DMPlex_Mult2D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
532:   }
533:   (void)PetscLogFlops(8.0*m);
534: }
535: PETSC_STATIC_INLINE void DMPlex_MatMultTransposeLeft3D_Internal(const PetscScalar A[], PetscInt m, PetscInt ldb, const PetscScalar B[], PetscScalar C[])
536: {
537:   PetscInt j;
538:   for (j = 0; j < m; ++j) {
539:     DMPlex_Mult3D_Internal(A, 1, &B[j*ldb], &C[j*ldb]);
540:   }
541:   (void)PetscLogFlops(8.0*m);
542: }

544: PETSC_STATIC_INLINE void DMPlex_PTAP2DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
545: {
546:   PetscScalar out[4];
547:   PetscInt i, j, k, l;
548:   for (i = 0; i < 2; ++i) {
549:     for (j = 0; j < 2; ++j) {
550:       out[i*2+j] = 0.;
551:       for (k = 0; k < 2; ++k) {
552:         for (l = 0; l < 2; ++l) {
553:           out[i*2+j] += P[k*2+i]*A[k*2+l]*P[l*2+j];
554:         }
555:       }
556:     }
557:   }
558:   for (i = 0; i < 2*2; ++i) C[i] = out[i];
559:   (void)PetscLogFlops(48.0);
560: }
561: PETSC_STATIC_INLINE void DMPlex_PTAP3DReal_Internal(const PetscReal P[], const PetscScalar A[], PetscScalar C[])
562: {
563:   PetscScalar out[9];
564:   PetscInt i, j, k, l;
565:   for (i = 0; i < 3; ++i) {
566:     for (j = 0; j < 3; ++j) {
567:       out[i*3+j] = 0.;
568:       for (k = 0; k < 3; ++k) {
569:         for (l = 0; l < 3; ++l) {
570:           out[i*3+j] += P[k*3+i]*A[k*3+l]*P[l*3+j];
571:         }
572:       }
573:     }
574:   }
575:   for (i = 0; i < 2*2; ++i) C[i] = out[i];
576:   (void)PetscLogFlops(243.0);
577: }
578: /* TODO Fix for aliasing of A and C */
579: PETSC_STATIC_INLINE void DMPlex_PTAPReal_Internal(const PetscReal P[], PetscInt m, PetscInt n, const PetscScalar A[], PetscScalar C[])
580: {
581:   PetscInt i, j, k, l;
582:   for (i = 0; i < n; ++i) {
583:     for (j = 0; j < n; ++j) {
584:       C[i*n+j] = 0.;
585:       for (k = 0; k < m; ++k) {
586:         for (l = 0; l < m; ++l) {
587:           C[i*n+j] += P[k*n+i]*A[k*m+l]*P[l*n+j];
588:         }
589:       }
590:     }
591:   }
592:   (void)PetscLogFlops(243.0);
593: }

595: PETSC_STATIC_INLINE void DMPlex_Transpose2D_Internal(PetscScalar A[])
596: {
597:   PetscScalar tmp;
598:   tmp = A[1]; A[1] = A[2]; A[2] = tmp;
599: }
600: PETSC_STATIC_INLINE void DMPlex_Transpose3D_Internal(PetscScalar A[])
601: {
602:   PetscScalar tmp;
603:   tmp = A[1]; A[1] = A[3]; A[3] = tmp;
604:   tmp = A[2]; A[2] = A[6]; A[6] = tmp;
605:   tmp = A[5]; A[5] = A[7]; A[7] = tmp;
606: }

608: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
609: {
610:   const PetscReal invDet = 1.0/detJ;

612:   invJ[0] =  invDet*J[3];
613:   invJ[1] = -invDet*J[1];
614:   invJ[2] = -invDet*J[2];
615:   invJ[3] =  invDet*J[0];
616:   (void)PetscLogFlops(5.0);
617: }

619: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
620: {
621:   const PetscReal invDet = 1.0/detJ;

623:   invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
624:   invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
625:   invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
626:   invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
627:   invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
628:   invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
629:   invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
630:   invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
631:   invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
632:   (void)PetscLogFlops(37.0);
633: }

635: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, const PetscReal J[])
636: {
637:   *detJ = J[0]*J[3] - J[1]*J[2];
638:   (void)PetscLogFlops(3.0);
639: }

641: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, const PetscReal J[])
642: {
643:   *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
644:            J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
645:            J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
646:   (void)PetscLogFlops(12.0);
647: }

649: PETSC_STATIC_INLINE void DMPlex_Det2D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
650: {
651:   *detJ = PetscRealPart(J[0])*PetscRealPart(J[3]) - PetscRealPart(J[1])*PetscRealPart(J[2]);
652:   (void)PetscLogFlops(3.0);
653: }

655: PETSC_STATIC_INLINE void DMPlex_Det3D_Scalar_Internal(PetscReal *detJ, const PetscScalar J[])
656: {
657:   *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])) +
658:            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])) +
659:            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])));
660:   (void)PetscLogFlops(12.0);
661: }

663: 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];}

665: 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;}

667: 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;}

669: 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);}

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

673: PETSC_INTERN PetscErrorCode DMPlexGetPointDualSpaceFEM(DM,PetscInt,PetscInt,PetscDualSpace *);
674: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPoint_Internal(PetscSection,PetscBool,PetscInt,PetscInt,PetscInt *,PetscBool,const PetscInt[],const PetscInt[],PetscInt[]);
675: PETSC_INTERN PetscErrorCode DMPlexGetIndicesPointFields_Internal(PetscSection,PetscBool,PetscInt,PetscInt,PetscInt[],PetscBool,const PetscInt***,PetscInt,const PetscInt[],PetscInt[]);
676: PETSC_INTERN PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]);

678: PETSC_EXTERN PetscErrorCode DMPlexGetAllCells_Internal(DM, IS *);
679: PETSC_EXTERN PetscErrorCode DMSNESGetFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
680: PETSC_EXTERN PetscErrorCode DMSNESRestoreFEGeom(DMField, IS, PetscQuadrature, PetscBool, PetscFEGeom **);
681: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Patch_Internal(DM, PetscSection, IS, PetscReal, Vec, Vec, Vec, void *);
682: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Patch_Internal(DM, PetscSection, PetscSection, IS, PetscReal, PetscReal, Vec, Vec, Mat, Mat, void *);
683: PETSC_INTERN PetscErrorCode DMCreateSubDomainDM_Plex(DM,DMLabel,PetscInt,IS*,DM*);
684: PETSC_INTERN PetscErrorCode DMPlexBasisTransformPoint_Internal(DM, DM, Vec, PetscInt, PetscBool[], PetscBool, PetscScalar *);
685: PETSC_EXTERN PetscErrorCode DMPlexBasisTransformPointTensor_Internal(DM, DM, Vec, PetscInt, PetscBool, PetscInt, PetscScalar *);
686: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApplyReal_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscReal *, PetscReal *, void *);
687: PETSC_INTERN PetscErrorCode DMPlexBasisTransformApply_Internal(DM, const PetscReal[], PetscBool, PetscInt, const PetscScalar *, PetscScalar *, void *);
688: PETSC_INTERN PetscErrorCode DMCreateNeumannOverlap_Plex(DM, IS*, Mat*, PetscErrorCode (**)(Mat, PetscReal, Vec, Vec, PetscReal, IS, void*), void **);

690: /* Functions in the vtable */
691: PETSC_INTERN PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
692: PETSC_INTERN PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat);
693: PETSC_INTERN PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
694: PETSC_INTERN PetscErrorCode DMCreateLocalSection_Plex(DM dm);
695: PETSC_INTERN PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm);
696: PETSC_INTERN PetscErrorCode DMCreateMatrix_Plex(DM dm,  Mat *J);
697: PETSC_INTERN PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
698: PETSC_INTERN PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field);
699: PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
700: PETSC_INTERN PetscErrorCode DMSetUp_Plex(DM dm);
701: PETSC_INTERN PetscErrorCode DMDestroy_Plex(DM dm);
702: PETSC_INTERN PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
703: PETSC_INTERN PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
704: PETSC_INTERN PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm);
705: PETSC_INTERN PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm);

707: #endif /* _PLEXIMPL_H */