Actual source code: dmpleximpl.h
petsc-3.7.6 2017-04-24
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_GlobalToNaturalBegin, DMPLEX_GlobalToNaturalEnd, DMPLEX_NaturalToGlobalBegin, DMPLEX_NaturalToGlobalEnd, 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: typedef struct {
62: PetscInt dummy;
63: } PetscPartitioner_Gather;
65: /* Utility struct to store the contents of a Gmsh file in memory */
66: typedef struct {
67: PetscInt dim; /* Entity dimension */
68: PetscInt id; /* Element number */
69: PetscInt numNodes; /* Size of node array */
70: int nodes[8]; /* Node array */
71: PetscInt numTags; /* Size of tag array */
72: int tags[4]; /* Tag array */
73: } GmshElement;
75: /* Utility struct to store the contents of a Fluent file in memory */
76: typedef struct {
77: int index; /* Type of section */
78: int zoneID;
79: int first;
80: int last;
81: int type;
82: int nd; /* Either ND or element-type */
83: void *data;
84: } FluentSection;
86: struct _PetscGridHash {
87: PetscInt dim;
88: PetscReal lower[3]; /* The lower-left corner */
89: PetscReal upper[3]; /* The upper-right corner */
90: PetscReal extent[3]; /* The box size */
91: PetscReal h[3]; /* The subbox size */
92: PetscInt n[3]; /* The number of subboxes */
93: PetscSection cellSection; /* Offsets for cells in each subbox*/
94: IS cells; /* List of cells in each subbox */
95: DMLabel cellsSparse; /* Sparse storage for cell map */
96: };
98: typedef struct {
99: PetscInt refct;
101: /* Sieve */
102: PetscSection coneSection; /* Layout of cones (inedges for DAG) */
103: PetscInt maxConeSize; /* Cached for fast lookup */
104: PetscInt *cones; /* Cone for each point */
105: 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 */
106: PetscSection supportSection; /* Layout of cones (inedges for DAG) */
107: PetscInt maxSupportSize; /* Cached for fast lookup */
108: PetscInt *supports; /* Cone for each point */
109: PetscBool refinementUniform; /* Flag for uniform cell refinement */
110: PetscReal refinementLimit; /* Maximum volume for refined cell */
111: PetscErrorCode (*refinementFunc)(const PetscReal [], PetscReal *); /* Function giving the maximum volume for refined cell */
112: PetscInt hybridPointMax[8]; /* Allow segregation of some points, each dimension has a divider (used in VTK output and refinement) */
114: PetscInt *facesTmp; /* Work space for faces operation */
116: /* Hierarchy */
117: PetscBool regularRefinement; /* This flag signals that we are a regular refinement of coarseMesh */
119: /* Generation */
120: char *tetgenOpts;
121: char *triangleOpts;
122: PetscPartitioner partitioner;
124: /* Submesh */
125: DMLabel subpointMap; /* Label each original mesh point in the submesh with its depth, subpoint are the implicit numbering */
127: /* Labels and numbering */
128: PetscObjectState depthState; /* State of depth label, so that we can determine if a user changes it */
129: IS globalVertexNumbers;
130: IS globalCellNumbers;
132: /* Constraints */
133: PetscSection anchorSection; /* maps constrained points to anchor points */
134: IS anchorIS; /* anchors indexed by the above section */
135: PetscErrorCode (*createanchors)(DM); /* automatically compute anchors (probably from tree constraints) */
136: PetscErrorCode (*computeanchormatrix)(DM,PetscSection,PetscSection,Mat);
138: /* Tree: automatically construct constraints for hierarchically non-conforming meshes */
139: PetscSection parentSection; /* dof == 1 if point has parent */
140: PetscInt *parents; /* point to parent */
141: PetscInt *childIDs; /* point to child ID */
142: PetscSection childSection; /* inverse of parent section */
143: PetscInt *children; /* point to children */
144: DM referenceTree; /* reference tree to which child ID's refer */
145: PetscErrorCode (*getchildsymmetry)(DM,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt*,PetscInt*);
147: /* Adjacency */
148: PetscBool useCone; /* Use cone() first when defining adjacency */
149: PetscBool useClosure; /* Use the transitive closure when defining adjacency */
150: PetscBool useAnchors; /* Replace constrained points with their anchors in adjacency lists */
152: /* Projection */
153: PetscInt maxProjectionHeight; /* maximum height of cells used in DMPlexProject functions */
155: /* Output */
156: PetscInt vtkCellHeight; /* The height of cells for output, default is 0 */
157: PetscReal scale[NUM_PETSC_UNITS]; /* The scale for each SI unit */
159: /* Geometry */
160: PetscReal minradius; /* Minimum distance from cell centroid to face */
161: PetscBool useHashLocation; /* Use grid hashing for point location */
162: PetscGridHash lbox; /* Local box for searching */
164: /* Debugging */
165: PetscBool printSetValues;
166: PetscInt printFEM;
167: PetscReal printTol;
168: } DM_Plex;
170: PETSC_EXTERN PetscErrorCode DMPlexVTKWriteAll_VTU(DM,PetscViewer);
171: PETSC_EXTERN PetscErrorCode DMPlexVTKGetCellType(DM,PetscInt,PetscInt,PetscInt*);
172: PETSC_EXTERN PetscErrorCode VecView_Plex_Local(Vec,PetscViewer);
173: PETSC_EXTERN PetscErrorCode VecView_Plex_Native(Vec,PetscViewer);
174: PETSC_EXTERN PetscErrorCode VecView_Plex(Vec,PetscViewer);
175: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Local(Vec,PetscViewer);
176: PETSC_EXTERN PetscErrorCode VecLoad_Plex_Native(Vec,PetscViewer);
177: PETSC_EXTERN PetscErrorCode VecLoad_Plex(Vec,PetscViewer);
178: PETSC_EXTERN PetscErrorCode DMPlexGetFieldType_Internal(DM, PetscSection, PetscInt, PetscInt *, PetscInt *, PetscViewerVTKFieldType *);
179: #if defined(PETSC_HAVE_HDF5)
180: PETSC_EXTERN PetscErrorCode VecView_Plex_Local_HDF5(Vec, PetscViewer);
181: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5(Vec, PetscViewer);
182: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5(Vec, PetscViewer);
183: PETSC_EXTERN PetscErrorCode VecView_Plex_HDF5_Native(Vec, PetscViewer);
184: PETSC_EXTERN PetscErrorCode VecLoad_Plex_HDF5_Native(Vec, PetscViewer);
185: PETSC_EXTERN PetscErrorCode DMPlexView_HDF5(DM, PetscViewer);
186: PETSC_EXTERN PetscErrorCode DMPlexLoad_HDF5(DM, PetscViewer);
187: #endif
189: PETSC_EXTERN PetscErrorCode DMPlexGetAdjacency_Internal(DM,PetscInt,PetscBool,PetscBool,PetscBool,PetscInt*,PetscInt*[]);
190: PETSC_EXTERN PetscErrorCode DMPlexGetFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
191: PETSC_EXTERN PetscErrorCode DMPlexGetRawFaces_Internal(DM,PetscInt,PetscInt,const PetscInt[], PetscInt*,PetscInt*,const PetscInt*[]);
192: PETSC_EXTERN PetscErrorCode DMPlexRestoreFaces_Internal(DM,PetscInt,PetscInt,PetscInt*,PetscInt*,const PetscInt*[]);
193: PETSC_EXTERN PetscErrorCode DMPlexRefineUniform_Internal(DM,CellRefiner,DM*);
194: PETSC_EXTERN PetscErrorCode DMPlexGetCellRefiner_Internal(DM,CellRefiner*);
195: PETSC_EXTERN PetscErrorCode CellRefinerGetAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
196: PETSC_EXTERN PetscErrorCode CellRefinerRestoreAffineTransforms_Internal(CellRefiner, PetscInt *, PetscReal *[], PetscReal *[], PetscReal *[]);
197: PETSC_EXTERN PetscErrorCode CellRefinerInCellTest_Internal(CellRefiner, const PetscReal[], PetscBool *);
198: PETSC_EXTERN PetscErrorCode DMPlexCreateGmsh_ReadElement(PetscViewer, PetscInt, PetscBool, PetscBool, GmshElement **);
199: PETSC_EXTERN PetscErrorCode DMPlexInvertCell_Internal(PetscInt, PetscInt, PetscInt[]);
200: PETSC_EXTERN PetscErrorCode DMPlexVecSetFieldClosure_Internal(DM, PetscSection, Vec, PetscBool[], PetscInt, const PetscScalar[], InsertMode);
201: PETSC_EXTERN PetscErrorCode DMPlexProjectConstraints_Internal(DM, Vec, Vec);
202: PETSC_EXTERN PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems*,DM);
203: PETSC_EXTERN PetscErrorCode DMPlexLabelComplete_Internal(DM,DMLabel,PetscBool);
204: PETSC_EXTERN PetscErrorCode DMPlexCreateReferenceTree_Union(DM,DM,const char *,DM*);
205: PETSC_EXTERN PetscErrorCode DMPlexComputeInterpolatorTree(DM,DM,PetscSF,PetscInt *,Mat);
206: PETSC_EXTERN PetscErrorCode DMPlexComputeInjectorTree(DM,DM,PetscSF,PetscInt *,Mat);
207: PETSC_EXTERN PetscErrorCode DMPlexAnchorsModifyMat(DM,PetscSection,PetscInt,PetscInt,const PetscInt*,const PetscScalar*,PetscInt*,PetscInt*,PetscInt**,PetscScalar**,PetscInt*,PetscBool);
208: PETSC_EXTERN PetscErrorCode indicesPoint_private(PetscSection,PetscInt,PetscInt,PetscInt *,PetscBool,PetscInt,PetscInt []);
209: PETSC_EXTERN PetscErrorCode indicesPointFields_private(PetscSection,PetscInt,PetscInt,PetscInt [],PetscBool,PetscInt,PetscInt []);
210: PETSC_EXTERN PetscErrorCode DMPlexLocatePoint_Internal(DM,PetscInt,const PetscScalar [],PetscInt,PetscInt *);
212: PETSC_INTERN PetscErrorCode DMPlexCreateCellNumbering_Internal(DM, PetscBool, IS *);
213: PETSC_INTERN PetscErrorCode DMPlexCreateVertexNumbering_Internal(DM, PetscBool, IS *);
217: /* invert dihedral symmetry: return a^-1,
218: * using the representation described in
219: * DMPlexGetConeOrientation() */
220: PETSC_STATIC_INLINE PetscInt DihedralInvert(PetscInt N, PetscInt a)
221: {
222: return (a <= 0) ? a : (N - a);
223: }
227: /* invert dihedral symmetry: return b * a,
228: * using the representation described in
229: * DMPlexGetConeOrientation() */
230: PETSC_STATIC_INLINE PetscInt DihedralCompose(PetscInt N, PetscInt a, PetscInt b)
231: {
232: if (!N) return 0;
233: return (a >= 0) ?
234: ((b >= 0) ? ((a + b) % N) : -(((a - b - 1) % N) + 1)) :
235: ((b >= 0) ? -(((N - b - a - 1) % N) + 1) : ((N + b - a) % N));
236: }
240: /* swap dihedral symmetries: return b * a^-1,
241: * using the representation described in
242: * DMPlexGetConeOrientation() */
243: PETSC_STATIC_INLINE PetscInt DihedralSwap(PetscInt N, PetscInt a, PetscInt b)
244: {
245: return DihedralCompose(N,DihedralInvert(N,a),b);
246: }
248: PETSC_EXTERN PetscErrorCode DMPlexComputeResidual_Internal(DM, PetscInt, PetscInt, PetscReal, Vec, Vec, Vec, void *);
249: PETSC_EXTERN PetscErrorCode DMPlexComputeJacobian_Internal(DM, PetscInt, PetscInt, PetscReal, PetscReal, Vec, Vec, Mat, Mat,void *);
253: PETSC_STATIC_INLINE void DMPlex_Invert2D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
254: {
255: const PetscReal invDet = 1.0/detJ;
257: invJ[0] = invDet*J[3];
258: invJ[1] = -invDet*J[1];
259: invJ[2] = -invDet*J[2];
260: invJ[3] = invDet*J[0];
261: (void)PetscLogFlops(5.0);
262: }
266: PETSC_STATIC_INLINE void DMPlex_Invert3D_Internal(PetscReal invJ[], PetscReal J[], PetscReal detJ)
267: {
268: const PetscReal invDet = 1.0/detJ;
270: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
271: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
272: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
273: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
274: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
275: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
276: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
277: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
278: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
279: (void)PetscLogFlops(37.0);
280: }
284: PETSC_STATIC_INLINE void DMPlex_Det2D_Internal(PetscReal *detJ, PetscReal J[])
285: {
286: *detJ = J[0]*J[3] - J[1]*J[2];
287: (void)PetscLogFlops(3.0);
288: }
292: PETSC_STATIC_INLINE void DMPlex_Det3D_Internal(PetscReal *detJ, PetscReal J[])
293: {
294: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
295: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
296: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
297: (void)PetscLogFlops(12.0);
298: }
302: 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];}
306: 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;}
310: 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;}
314: 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);}
319: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
320: {
322: #if defined(PETSC_USE_DEBUG)
323: {
324: PetscInt dof;
326: *start = *end = 0; /* Silence overzealous compiler warning */
327: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
328: PetscSectionGetOffset(dm->defaultSection, point, start);
329: PetscSectionGetDof(dm->defaultSection, point, &dof);
330: *end = *start + dof;
331: }
332: #else
333: {
334: const PetscSection s = dm->defaultSection;
335: *start = s->atlasOff[point - s->pStart];
336: *end = *start + s->atlasDof[point - s->pStart];
337: }
338: #endif
339: return(0);
340: }
344: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetLocalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
345: {
347: #if defined(PETSC_USE_DEBUG)
348: {
349: PetscInt dof;
351: *start = *end = 0; /* Silence overzealous compiler warning */
352: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
353: PetscSectionGetFieldOffset(dm->defaultSection, point, field, start);
354: PetscSectionGetFieldDof(dm->defaultSection, point, field, &dof);
355: *end = *start + dof;
356: }
357: #else
358: {
359: const PetscSection s = dm->defaultSection->field[field];
360: *start = s->atlasOff[point - s->pStart];
361: *end = *start + s->atlasDof[point - s->pStart];
362: }
363: #endif
364: return(0);
365: }
369: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalOffset_Private(DM dm, PetscInt point, PetscInt *start, PetscInt *end)
370: {
372: #if defined(PETSC_USE_DEBUG)
373: {
375: PetscInt dof,cdof;
376: *start = *end = 0; /* Silence overzealous compiler warning */
377: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
378: if (!dm->defaultGlobalSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default global Section. It will be created automatically by DMGetDefaultGlobalSection()");
379: PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
380: PetscSectionGetDof(dm->defaultGlobalSection, point, &dof);
381: PetscSectionGetConstraintDof(dm->defaultGlobalSection, point, &cdof);
382: *end = *start + dof - cdof + (dof < 0 ? 1 : 0);
383: }
384: #else
385: {
386: const PetscSection s = dm->defaultGlobalSection;
387: const PetscInt dof = s->atlasDof[point - s->pStart];
388: const PetscInt cdof = s->bc ? s->bc->atlasDof[point - s->bc->pStart] : 0;
389: *start = s->atlasOff[point - s->pStart];
390: *end = *start + dof - cdof + (dof < 0 ? 1 : 0);
391: }
392: #endif
393: return(0);
394: }
398: PETSC_STATIC_INLINE PetscErrorCode DMPlexGetGlobalFieldOffset_Private(DM dm, PetscInt point, PetscInt field, PetscInt *start, PetscInt *end)
399: {
401: #if defined(PETSC_USE_DEBUG)
402: {
403: PetscInt loff, lfoff, fdof, fcdof, ffcdof, f;
405: *start = *end = 0; /* Silence overzealous compiler warning */
406: if (!dm->defaultSection) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "DM must have a default Section, see DMSetDefaultSection()");
407: 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()");
408: PetscSectionGetOffset(dm->defaultGlobalSection, point, start);
409: PetscSectionGetOffset(dm->defaultSection, point, &loff);
410: PetscSectionGetFieldOffset(dm->defaultSection, point, field, &lfoff);
411: PetscSectionGetFieldDof(dm->defaultSection, point, field, &fdof);
412: PetscSectionGetFieldConstraintDof(dm->defaultSection, point, field, &fcdof);
413: *start = *start < 0 ? *start - (lfoff-loff) : *start + lfoff-loff;
414: for (f = 0; f < field; ++f) {
415: PetscSectionGetFieldConstraintDof(dm->defaultSection, point, f, &ffcdof);
416: *start = *start < 0 ? *start + ffcdof : *start - ffcdof;
417: }
418: *end = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
419: }
420: #else
421: {
422: const PetscSection s = dm->defaultSection;
423: const PetscSection fs = dm->defaultSection->field[field];
424: const PetscSection gs = dm->defaultGlobalSection;
425: const PetscInt loff = s->atlasOff[point - s->pStart];
426: const PetscInt goff = gs->atlasOff[point - s->pStart];
427: const PetscInt lfoff = fs->atlasOff[point - s->pStart];
428: const PetscInt fdof = fs->atlasDof[point - s->pStart];
429: const PetscInt fcdof = fs->bc ? fs->bc->atlasDof[point - fs->bc->pStart] : 0;
430: PetscInt ffcdof = 0, f;
432: for (f = 0; f < field; ++f) {
433: const PetscSection ffs = dm->defaultSection->field[f];
434: ffcdof += ffs->bc ? ffs->bc->atlasDof[point - ffs->bc->pStart] : 0;
435: }
436: *start = goff + (goff < 0 ? loff-lfoff + ffcdof : lfoff-loff - ffcdof);
437: *end = *start < 0 ? *start - (fdof-fcdof) : *start + fdof-fcdof;
438: }
439: #endif
440: return(0);
441: }
443: #endif /* _PLEXIMPL_H */