Actual source code: dmpleximpl.h
petsc-3.6.3 2015-12-03
1: #if !defined(_PLEXIMPL_H)
2: #define _PLEXIMPL_H
4: #include <petscmat.h> /*I "petscmat.h" I*/
5: #include <petscdmplex.h> /*I "petscdmplex.h" I*/
6: #include <petscbt.h>
7: #include <petscsf.h>
8: #include <petsc/private/dmimpl.h>
9: #include <petsc/private/isimpl.h> /* for inline access to atlasOff */
10: #include <../src/sys/utils/hash.h>
12: PETSC_EXTERN PetscLogEvent DMPLEX_Interpolate, PETSCPARTITIONER_Partition, DMPLEX_Distribute, DMPLEX_DistributeCones, DMPLEX_DistributeLabels, DMPLEX_DistributeSF, DMPLEX_DistributeOverlap, DMPLEX_DistributeField, DMPLEX_DistributeData, DMPLEX_Migrate, DMPLEX_Stratify, DMPLEX_Preallocate, DMPLEX_ResidualFEM, DMPLEX_JacobianFEM, DMPLEX_InterpolatorFEM, DMPLEX_InjectorFEM, DMPLEX_IntegralFEM, DMPLEX_CreateGmsh;
14: PETSC_EXTERN PetscBool PetscPartitionerRegisterAllCalled;
15: PETSC_EXTERN PetscErrorCode PetscPartitionerRegisterAll(void);
16: PETSC_EXTERN PetscErrorCode PetscPartitionerSetTypeFromOptions_Internal(PetscPartitioner);
18: typedef enum {REFINER_NOOP = 0,
19: REFINER_SIMPLEX_1D,
20: REFINER_SIMPLEX_2D,
21: REFINER_HYBRID_SIMPLEX_2D,
22: REFINER_HEX_2D,
23: REFINER_HYBRID_HEX_2D,
24: REFINER_SIMPLEX_3D,
25: REFINER_HYBRID_SIMPLEX_3D,
26: REFINER_HEX_3D,
27: REFINER_HYBRID_HEX_3D} CellRefiner;
29: typedef struct _PetscPartitionerOps *PetscPartitionerOps;
30: struct _PetscPartitionerOps {
31: PetscErrorCode (*setfromoptions)(PetscPartitioner);
32: PetscErrorCode (*setup)(PetscPartitioner);
33: PetscErrorCode (*view)(PetscPartitioner,PetscViewer);
34: PetscErrorCode (*destroy)(PetscPartitioner);
35: PetscErrorCode (*partition)(PetscPartitioner, DM, PetscInt, PetscInt, PetscInt[], PetscInt[], PetscSection, IS *);
36: };
38: struct _p_PetscPartitioner {
39: PETSCHEADER(struct _PetscPartitionerOps);
40: void *data; /* Implementation object */
41: PetscInt height; /* Height of points to partition into non-overlapping subsets */
42: };
44: typedef struct {
45: PetscInt dummy;
46: } PetscPartitioner_Chaco;
48: typedef struct {
49: PetscInt dummy;
50: } PetscPartitioner_ParMetis;
52: typedef struct {
53: PetscSection section; /* Sizes for each partition */
54: IS partition; /* Points in each partition */
55: } PetscPartitioner_Shell;
57: typedef struct {
58: PetscInt dummy;
59: } PetscPartitioner_Simple;
61: /* This is an integer map, in addition it is also a container class
62: Design points:
63: - Low storage is the most important design point
64: - We want flexible insertion and deletion
65: - We can live with O(log) query, but we need O(1) iteration over strata
66: */
67: struct _n_DMLabel {
68: PetscInt refct;
69: char *name; /* Label name */
70: PetscInt numStrata; /* Number of integer values */
71: PetscInt *stratumValues; /* Value of each stratum */
72: /* Basic sorted array storage */
73: PetscBool *arrayValid; /* The array storage is valid (no additions need to be merged in) */
74: PetscInt *stratumSizes; /* Size of each stratum */
75: PetscInt **points; /* Points for each stratum, always sorted */
76: /* Hashtable for fast insertion */
77: PetscHashI *ht; /* Hash table for fast insertion */
78: /* Index for fast search */
79: PetscInt pStart, pEnd; /* Bounds for index lookup */
80: PetscBT bt; /* A bit-wise index */
81: };
83: struct _n_PlexLabel {
84: DMLabel label;
85: PetscBool output;
86: struct _n_PlexLabel *next;
87: };
88: typedef struct _n_PlexLabel *PlexLabel;
90: /* Utility struct to store the contents of a Gmsh file in memory */
91: typedef struct {
92: PetscInt dim; /* Entity dimension */
93: PetscInt id; /* Element number */
94: PetscInt numNodes; /* Size of node array */
95: int nodes[8]; /* Node array */
96: PetscInt numTags; /* Size of tag array */
97: int tags[4]; /* Tag array */
98: } GmshElement;
101: /* Utility struct to store the contents of a Fluent file in memory */
102: typedef struct {
103: int index; /* Type of section */
104: int zoneID;
105: int first;
106: int last;
107: int type;
108: int nd; /* Either ND or element-type */
109: void *data;
110: } FluentSection;
112: struct _n_Boundary {
113: const char *name;
114: const char *labelname;
115: DMLabel label;
116: PetscBool essential;
117: PetscInt field;
118: PetscInt numcomps;
119: PetscInt *comps;
120: void (*func)();
121: PetscInt numids;
122: PetscInt *ids;
123: void *ctx;
124: DMBoundary next;
125: };
127: typedef struct {
128: PetscInt refct;
130: /* Sieve */
131: PetscSection coneSection; /* Layout of cones (inedges for DAG) */
132: PetscInt maxConeSize; /* Cached for fast lookup */
133: PetscInt *cones; /* Cone for each point */
134: 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 */
135: PetscSection supportSection; /* Layout of cones (inedges for DAG) */
136: PetscInt maxSupportSize; /* Cached for fast lookup */
137: PetscInt *supports; /* Cone for each point */
138: PetscBool refinementUniform; /* Flag for uniform cell refinement */
139: PetscReal refinementLimit; /* Maximum volume for refined cell */
140: PetscInt hybridPointMax[8]; /* Allow segregation of some points, each dimension has a divider (used in VTK output and refinement) */
142: PetscInt *facesTmp; /* Work space for faces operation */
144: /* Hierarchy */
145: DM coarseMesh; /* This mesh was obtained from coarse mesh using DMRefineHierarchy() */
147: /* Generation */
148: char *tetgenOpts;
149: char *triangleOpts;
150: PetscPartitioner partitioner;
152: /* Submesh */
153: DMLabel subpointMap; /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
155: /* Labels and numbering */
156: PlexLabel labels; /* Linked list of labels */
157: DMLabel depthLabel;
158: IS globalVertexNumbers;
159: IS globalCellNumbers;
161: /* Constraints */
162: PetscSection anchorSection; /* maps constrained points to anchor points */
163: IS anchorIS; /* anchors indexed by the above section */
164: PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
165: PetscErrorCode (*computeanchormatrix)(DM,PetscSection,PetscSection,Mat);
167: /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
168: PetscSection parentSection; /* dof == 1 if point has parent */
169: PetscInt *parents; /* point to parent */
170: PetscInt *childIDs; /* point to child ID */
171: PetscSection childSection; /* inverse of parent section */
172: PetscInt *children; /* point to children */
173: DM referenceTree; /* reference tree to which child ID's refer */
174: PetscErrorCode (*getchildsymmetry)(DM,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*);
176: /* Adjacency */
177: PetscBool useCone; /* Use cone() first when defining adjacency */
178: PetscBool useClosure; /* Use the transitive closure when defining adjacency */
179: PetscBool useAnchors; /* Replace constrained points with their anchors in adjacency lists */
181: /* Projection */
182: PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
184: /* Output */
185: PetscInt vtkCellHeight; /* The height of cells for output, default is 0 */
186: PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */
188: /* Problem definition */
189: DMBoundary boundary; /* List of boundary conditions */
191: /* Geometry */
192: PetscReal minradius; /* Minimum distance from cell centroid to face */
195: /* Debugging */
196: PetscBool printSetValues;
197: PetscInt printFEM;
198: PetscReal printTol;
199: } DM_Plex;
201: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
202: PETSC_EXTERN PetscErrorCode DMPlexVTKGetCellType(DM,PetscInt,PetscInt,PetscInt*);
203: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
204: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
205: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
206: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
207: PETSC_EXTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
208: #if defined(PETSC_HAVE_HDF5)
209: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
210: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
211: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
212: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
213: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
214: #endif
216: PETSC_EXTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
217: PETSC_EXTERN PetscErrorCode DMPlexGetFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
218: PETSC_EXTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,PetscInt,PetscInt,const PetscInt[], PetscInt*,PetscInt*,const PetscInt*[]);
219: PETSC_EXTERN PetscErrorCode DMPlexRestoreFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
220: PETSC_EXTERN PetscErrorCode DMPlexRefineUniform_Internal(DM,CellRefiner,DM*);
221: PETSC_EXTERN PetscErrorCode DMPlexGetCellRefiner_Internal(DM,CellRefiner*);
222: PETSC_EXTERN PetscErrorCode CellRefinerGetAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
223: PETSC_EXTERN PetscErrorCode CellRefinerRestoreAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
224: PETSC_EXTERN PetscErrorCode CellRefinerInCellTest_Internal(CellRefiner, const PetscReal[], PetscBool *);
225: PETSC_EXTERN PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer, PetscInt, PetscBool, PetscBool, GmshElement **);
226: PETSC_EXTERN PetscErrorCode DMPlexInvertCell_Internal(PetscInt, PetscInt, PetscInt[]);
227: PETSC_EXTERN PetscErrorCode DMPlexLocalizeCoordinate_Internal(DM, PetscInt, const PetscScalar[], const PetscScalar[], PetscScalar[]);
228: PETSC_EXTERN PetscErrorCode DMPlexLocalizeCoordinateReal_Internal(DM, PetscInt, const PetscReal[], const PetscReal[], PetscReal[]);
229: PETSC_EXTERN PetscErrorCode DMPlexLocalizeAddCoordinate_Internal(DM, PetscInt, const PetscScalar[], const PetscScalar[], PetscScalar[]);
230: PETSC_EXTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, const PetscScalar[], InsertMode);
231: PETSC_EXTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
235: /* invert dihedral symmetry: return a^-1,
236: * using the representation described in
237: * DMPlexGetConeOrientation() */
238: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
239: {
240: return (a <= 0) ? a : (N - a);
241: }
245: /* invert dihedral symmetry: return b * a,
246: * using the representation described in
247: * DMPlexGetConeOrientation() */
248: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
249: {
250: if (!N) return 0;
251: return (a >= 0) ?
252: ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
253: ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
254: }
258: /* swap dihedral symmetries: return b * a^-1,
259: * using the representation described in
260: * DMPlexGetConeOrientation() */
261: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
262: {
263: return DihedralCompose(N,DihedralInvert(N,a),b);
264: }
266: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, PetscInt, PetscInt, PetscReal, Vec, Vec, Vec, void *);
267: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, PetscInt, PetscInt, PetscReal, PetscReal, Vec, Vec, Mat, Mat,void *);
271: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
272: {
273: const PetscReal invDet = 1.0/detJ;
275: invJ[0] = invDet*J[3];
276: invJ[1] = -invDet*J[1];
277: invJ[2] = -invDet*J[2];
278: invJ[3] = invDet*J[0];
279: PetscLogFlops(5.0);
280: }
284: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
285: {
286: const PetscReal invDet = 1.0/detJ;
288: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
289: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
290: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
291: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
292: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
293: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
294: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
295: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
296: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
297: PetscLogFlops(37.0);
298: }
302: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, PetscReal J[])
303: {
304: *detJ = J[0]*J[3] - J[1]*J[2];
305: PetscLogFlops(3.0);
306: }
310: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, PetscReal J[])
311: {
312: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
313: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
314: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
315: PetscLogFlops(12.0);
316: }
320: 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];}
324: 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;}
328: 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;}
332: 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);}
337: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
338: {
340: #if defined(PETSC_USE_DEBUG)
341: {
342: PetscInt dof;
344: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
345: PetscSectionGetOffset(dm->defaultSection, point, start);
346: PetscSectionGetDof(dm->defaultSection, point, &dof);
347: *end = *start + dof;
348: }
349: #else
350: {
351: const PetscSection s = dm->defaultSection;
352: *start = s->atlasOff[point - s->pStart];
353: *end = *start + s->atlasDof[point - s->pStart];
354: }
355: #endif
356: return(0);
357: }
361: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
362: {
364: #if defined(PETSC_USE_DEBUG)
365: {
366: PetscInt dof;
368: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
369: PetscSectionGetFieldOffset(dm->defaultSection, point, field, start);
370: PetscSectionGetFieldDof(dm->defaultSection, point, field, &dof);
371: *end = *start + dof;
372: }
373: #else
374: {
375: const PetscSection s = dm->defaultSection->field[field];
376: *start = s->atlasOff[point - s->pStart];
377: *end = *start + s->atlasDof[point - s->pStart];
378: }
379: #endif
380: return(0);
381: }
385: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
386: {
388: #if defined(PETSC_USE_DEBUG)
389: {
391: PetscInt dof,cdof;
392: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
393: if (!dm->defaultGlobalSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default global Section. It will be crated automatically by DMGetDefaultGlobalSection()");
394: PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
395: PetscSectionGetDof(dm->defaultGlobalSection, point, &dof);
396: PetscSectionGetConstraintDof(dm->defaultGlobalSection, point, &cdof);
397: *end = *start + dof-cdof;
398: }
399: #else
400: {
401: const PetscSection s = dm->defaultGlobalSection;
402: const PetscInt dof = s->atlasDof[point - s->pStart];
403: const PetscInt cdof = s->bc ? s->bc->atlasDof[point - s->bc->pStart] : 0;
404: *start = s->atlasOff[point - s->pStart];
405: *end = *start + dof-cdof;
406: }
407: #endif
408: return(0);
409: }
413: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
414: {
416: #if defined(PETSC_USE_DEBUG)
417: {
418: PetscInt loff, lfoff, fdof, fcdof, ffcdof, f;
421: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
422: if (!dm->defaultGlobalSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default global Section. It will be crated automatically by DMGetDefaultGlobalSection()");
423: PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
424: PetscSectionGetOffset(dm->defaultSection, point, &loff);
425: PetscSectionGetFieldOffset(dm->defaultSection, point, field, &lfoff);
426: PetscSectionGetFieldDof(dm->defaultSection, point, field, &fdof);
427: PetscSectionGetFieldConstraintDof(dm->defaultSection, point, field, &fcdof);
428: *start = *start < 0 ? *start - (lfoff-loff) : *start + lfoff-loff;
429: for (f = 0; f < field; ++f) {
430: PetscSectionGetFieldConstraintDof(dm->defaultSection, point, f, &ffcdof);
431: *start = *start < 0 ? *start + ffcdof : *start - ffcdof;
432: }
433: *end = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
434: }
435: #else
436: {
437: const PetscSection s = dm->defaultSection;
438: const PetscSection fs = dm->defaultSection->field[field];
439: const PetscSection gs = dm->defaultGlobalSection;
440: const PetscInt loff = s->atlasOff[point - s->pStart];
441: const PetscInt goff = gs->atlasOff[point - s->pStart];
442: const PetscInt lfoff = fs->atlasOff[point - s->pStart];
443: const PetscInt fdof = fs->atlasDof[point - s->pStart];
444: const PetscInt fcdof = fs->bc ? fs->bc->atlasDof[point - fs->bc->pStart] : 0;
445: PetscInt ffcdof = 0, f;
447: for (f = 0; f < field; ++f) {
448: const PetscSection ffs = dm->defaultSection->field[f];
449: ffcdof += ffs->bc ? ffs->bc->atlasDof[point - ffs->bc->pStart] : 0;
450: }
451: *start = goff + (goff < 0 ? loff-lfoff + ffcdof : lfoff-loff - ffcdof);
452: *end = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
453: }
454: #endif
455: return(0);
456: }
458: #endif /* _PLEXIMPL_H */