isoSurface.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2016-2019 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 Class
28  Foam::isoSurface
29 
30 Description
31  A surface formed by the iso value.
32  After "Regularised Marching Tetrahedra: improved iso-surface extraction",
33  G.M. Treece, R.W. Prager and A.H. Gee.
34 
35  Note:
36  - does tets without using cell centres/cell values. Not tested.
37  - regularisation can create duplicate triangles/non-manifold surfaces.
38  Current handling of those is bit ad-hoc for now and not perfect.
39  - regularisation does not do boundary points so as to preserve the
40  boundary perfectly.
41  - uses geometric merge with fraction of bounding box as distance.
42  - triangles can be between two cell centres so constant sampling
43  does not make sense.
44  - on empty patches behaves like zero gradient.
45  - does not do 2D correctly, creates non-flat iso surface.
46  - on processor boundaries might two overlapping (identical) triangles
47  (one from either side)
48 
49  The handling on coupled patches is a bit complex. All fields
50  (values and coordinates) get rewritten so
51  - empty patches get zerogradient (value) and facecentre (coordinates)
52  - separated processor patch faces get interpolate (value) and facecentre
53  (coordinates). (this is already the default for cyclics)
54  - non-separated processor patch faces keep opposite value and cell centre
55 
56  Now the triangle generation on non-separated processor patch faces
57  can use the neighbouring value. Any separated processor face or cyclic
58  face gets handled just like any boundary face.
59 
60 SourceFiles
61  isoSurface.C
62  isoSurfaceTemplates.C
63 
64 \*---------------------------------------------------------------------------*/
65 
66 #ifndef isoSurface_H
67 #define isoSurface_H
68 
69 #include "bitSet.H"
70 #include "volFields.H"
71 #include "slicedVolFields.H"
72 #include "isoSurfaceBase.H"
73 
74 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
75 
76 namespace Foam
77 {
78 
79 // Forward declarations
80 
81 class fvMesh;
82 class plane;
83 class treeBoundBox;
84 class triSurface;
85 
86 /*---------------------------------------------------------------------------*\
87  Class isoSurface Declaration
88 \*---------------------------------------------------------------------------*/
89 
90 class isoSurface
91 :
92  public isoSurfaceBase
93 {
94  // Private data
95 
96  enum segmentCutType
97  {
98  NEARFIRST, // intersection close to e.first()
99  NEARSECOND, // ,, e.second()
100  NOTNEAR // no intersection
101  };
102 
103  enum cellCutType
104  {
105  NOTCUT, // not cut
106  SPHERE, // all edges to cell centre cut
107  CUT // normal cut
108  };
109 
110 
111  //- Reference to mesh
112  const fvMesh& mesh_;
113 
114  const scalarField& pVals_;
115 
116  //- Input volScalarField with separated coupled patches rewritten
118 
119  //- Regularise?
120  const bool regularise_;
121 
122  //- When to merge points
123  const scalar mergeDistance_;
124 
125  //- Whether face might be cut
126  List<cellCutType> faceCutType_;
127 
128  //- Whether cell might be cut
129  List<cellCutType> cellCutType_;
130 
131  //- Estimated number of cut cells
132  label nCutCells_;
133 
134  //- For every unmerged triangle point the point in the triSurface
135  labelList triPointMergeMap_;
136 
137  //- triSurface points that have weighted interpolation
138  DynamicList<label> interpolatedPoints_;
139 
140  //- corresponding original, unmerged points
141  DynamicList<FixedList<label, 3>> interpolatedOldPoints_;
142 
143  //- corresponding weights
144  DynamicList<FixedList<scalar, 3>> interpolationWeights_;
145 
146 
147  // Private Member Functions
148 
149  // Point synchronisation
150 
151  //- Does tensor differ (to within mergeTolerance) from identity
152  bool noTransform(const tensor& tt) const;
153 
154  //- Is patch a collocated (i.e. non-separated) coupled patch?
155  static bool collocatedPatch(const polyPatch& pp);
156 
157  //- Per face whether is collocated
158  bitSet collocatedFaces(const coupledPolyPatch&) const;
159 
160  //- Synchronise points on all non-separated coupled patches
161  void syncUnseparatedPoints
162  (
163  pointField& collapsedPoint,
164  const point& nullValue
165  ) const;
166 
167 
168  //- Get location of iso value as fraction inbetween s0,s1
169  scalar isoFraction
170  (
171  const scalar s0,
172  const scalar s1
173  ) const;
174 
175  //- Check if any edge of a face is cut
176  bool isEdgeOfFaceCut
177  (
178  const scalarField& pVals,
179  const face& f,
180  const bool ownLower,
181  const bool neiLower
182  ) const;
183 
184  //- Get neighbour value and position.
185  void getNeighbour
186  (
187  const labelList& boundaryRegion,
188  const volVectorField& meshC,
189  const volScalarField& cVals,
190  const label celli,
191  const label facei,
192  scalar& nbrValue,
193  point& nbrPoint
194  ) const;
195 
196  //- Determine for every face/cell whether it (possibly) generates
197  // triangles.
198  void calcCutTypes
199  (
200  const labelList& boundaryRegion,
201  const volVectorField& meshC,
202  const volScalarField& cVals,
203  const scalarField& pVals
204  );
205 
206  static point calcCentre(const triSurface&);
207 
208  //- Determine per cc whether all near cuts can be snapped to single
209  // point.
210  void calcSnappedCc
211  (
212  const labelList& boundaryRegion,
213  const volVectorField& meshC,
214  const volScalarField& cVals,
215  const scalarField& pVals,
216  DynamicList<point>& snappedPoints,
217  labelList& snappedCc
218  ) const;
219 
220  //- Determine per point whether all near cuts can be snapped to single
221  // point.
222  void calcSnappedPoint
223  (
224  const bitSet& isBoundaryPoint,
225  const labelList& boundaryRegion,
226  const volVectorField& meshC,
227  const volScalarField& cVals,
228  const scalarField& pVals,
229  DynamicList<point>& snappedPoints,
230  labelList& snappedPoint
231  ) const;
232 
233 
234  //- Return input field with coupled (and empty) patch values rewritten
235  template<class Type>
238  adaptPatchFields
239  (
241  ) const;
242 
243  //- Generate single point by interpolation or snapping
244  template<class Type>
245  Type generatePoint
246  (
247  const scalar s0,
248  const Type& p0,
249  const bool hasSnap0,
250  const Type& snapP0,
251 
252  const scalar s1,
253  const Type& p1,
254  const bool hasSnap1,
255  const Type& snapP1
256  ) const;
257 
258 
259  //- Note: cannot use simpler isoSurfaceCell::generateTriPoints since
260  // the need here to sometimes pass in remote 'snappedPoints'
261  template<class Type>
262  void generateTriPoints
263  (
264  const scalar s0,
265  const Type& p0,
266  const bool hasSnap0,
267  const Type& snapP0,
268 
269  const scalar s1,
270  const Type& p1,
271  const bool hasSnap1,
272  const Type& snapP1,
273 
274  const scalar s2,
275  const Type& p2,
276  const bool hasSnap2,
277  const Type& snapP2,
278 
279  const scalar s3,
280  const Type& p3,
281  const bool hasSnap3,
282  const Type& snapP3,
283 
284  DynamicList<Type>& pts
285  ) const;
286 
287  template<class Type>
288  label generateFaceTriPoints
289  (
290  const volScalarField& cVals,
291  const scalarField& pVals,
292 
294  const Field<Type>& pCoords,
295 
296  const DynamicList<Type>& snappedPoints,
297  const labelList& snappedCc,
298  const labelList& snappedPoint,
299  const label facei,
300 
301  const scalar neiVal,
302  const Type& neiPt,
303  const bool hasNeiSnap,
304  const Type& neiSnapPt,
305 
307  DynamicList<label>& triMeshCells
308  ) const;
309 
310  template<class Type>
311  void generateTriPoints
312  (
313  const volScalarField& cVals,
314  const scalarField& pVals,
315 
317  const Field<Type>& pCoords,
318 
319  const DynamicList<Type>& snappedPoints,
320  const labelList& snappedCc,
321  const labelList& snappedPoint,
322 
324  DynamicList<label>& triMeshCells
325  ) const;
326 
327  template<class Type>
328  static tmp<Field<Type>> interpolate
329  (
330  const label nPoints,
331  const labelList& triPointMergeMap,
332  const labelList& interpolatedPoints,
333  const List<FixedList<label, 3>>& interpolatedOldPoints,
335  const DynamicList<Type>& unmergedValues
336  );
337 
338  triSurface stitchTriPoints
339  (
340  const bool checkDuplicates,
341  const List<point>& triPoints,
342  labelList& triPointReverseMap, // unmerged to merged point
343  labelList& triMap // merged to unmerged triangle
344  ) const;
345 
346  //- Trim triangle to planes
347  static void trimToPlanes
348  (
349  const PtrList<plane>& planes,
350  const triPointRef& tri,
351  DynamicList<point>& newTriPoints
352  );
353 
354  //- Trim all triangles to box
355  static void trimToBox
356  (
357  const treeBoundBox& bb,
359  DynamicList<label>& triMeshCells
360  );
361 
362  //- Trim all triangles to box. Determine interpolation
363  // for existing and new points
364  static void trimToBox
365  (
366  const treeBoundBox& bb,
368  DynamicList<label>& triMap,
369  labelList& triPointMap,
370  labelList& interpolatedPoints,
371  List<FixedList<label, 3>>& interpolatedOldPoints,
373  );
374 
375  static triSurface subsetMesh
376  (
377  const triSurface&,
378  const labelList& newToOldFaces,
379  labelList& oldToNewPoints,
380  labelList& newToOldPoints
381  );
382 
383 public:
384 
385  //- Declare friendship to share some functionality
386  friend class isoSurfaceCell;
387  friend class isoSurfaceTopo;
388 
389  //- Filtering type
391 
392 
393  //- Runtime type information
394  TypeName("isoSurface");
395 
396 
397  // Constructors
398 
399  //- Construct from cell values and point values.
400  // Uses boundaryField for boundary values.
401  // Holds reference to cellIsoVals and pointIsoVals.
402  //
403  // \param bounds optional bounding box for trimming
404  // \param mergeTol fraction of mesh bounding box for merging points
405  isoSurface
406  (
407  const volScalarField& cellValues,
408  const scalarField& pointValues,
409  const scalar iso,
410  const isoSurfaceBase::filterType filter,
411  const boundBox& bounds = boundBox::invertedBox,
412  const scalar mergeTol = 1e-6
413  );
414 
415 
416  //- Construct from cell values and point values.
417  // Uses boundaryField for boundary values.
418  // Holds reference to cellIsoVals and pointIsoVals.
419  //
420  // \param bounds optional bounding box for trimming
421  // \param mergeTol fraction of mesh bounding box for merging points
422  isoSurface
423  (
424  const volScalarField& cellValues,
425  const scalarField& pointValues,
426  const scalar iso,
427  const bool regularise,
428  const boundBox& bounds = boundBox::invertedBox,
429  const scalar mergeTol = 1e-6
430  );
431 
432 
433  // Member Functions
434 
435  //- Interpolates cCoords, pCoords.
436  // Uses the references to the original fields used to create the
437  // iso surface.
438  template<class Type>
440  (
442  const Field<Type>& pCoords
443  ) const;
444 
445 };
446 
447 
448 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
449 
450 } // End namespace Foam
451 
452 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
453 
454 #ifdef NoRepository
455  #include "isoSurfaceTemplates.C"
456 #endif
457 
458 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
459 
460 #endif
461 
462 // ************************************************************************* //
Foam::fvPatchField
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: volSurfaceMapping.H:50
volFields.H
Foam::Tensor< scalar >
Foam::isoSurface::TypeName
TypeName("isoSurface")
Runtime type information.
Foam::isoSurfaceBase
Low-level components common to various iso-surface algorithms.
Definition: isoSurfaceBase.H:54
Foam::bitSet
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:64
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:59
isoSurfaceTemplates.C
Foam::DynamicList< label >
slicedVolFields.H
Foam::volMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: volMesh.H:50
Foam::boundBox::invertedBox
static const boundBox invertedBox
A large inverted boundBox: min/max == +/- ROOTVGREAT.
Definition: boundBox.H:86
Foam::treeBoundBox
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:87
Foam::boundaryRegion
The boundaryRegion persistent data saved as a Map<dictionary>.
Definition: boundaryRegion.H:72
Foam::coupledPolyPatch
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
Definition: coupledPolyPatch.H:53
bitSet.H
isoSurfaceBase.H
Foam::isoSurfaceCell
A surface formed by the iso value. After "Polygonising A Scalar Field Using Tetrahedrons",...
Definition: isoSurfaceCell.H:69
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
Foam::triangle
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:62
Foam::triPoints
Triangle storage. Null constructable (unfortunately triangle<point, point> is not)
Definition: triPoints.H:52
Foam::isoSurface
A surface formed by the iso value. After "Regularised Marching Tetrahedra: improved iso-surface extra...
Definition: isoSurface.H:89
Foam::Field< scalar >
Foam::triSurface
Triangulated surface description with patch information.
Definition: triSurface.H:76
Foam::polyPatch
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:67
Foam::PtrList
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: List.H:62
fld
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< ' ';}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< ' ';}gmvFile<< nl;for(const word &name :lagrangianScalarNames){ IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputLagrangian.H:23
Foam::interpolationWeights
Abstract base class for interpolating in 1D.
Definition: interpolationWeights.H:58
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:84
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::slicedFvPatchField
Specialization of fvPatchField which creates the underlying fvPatchField as a slice of the given comp...
Definition: slicedFvPatchField.H:64
Foam::isoSurface::isoSurface
isoSurface(const volScalarField &cellValues, const scalarField &pointValues, const scalar iso, const isoSurfaceBase::filterType filter, const boundBox &bounds=boundBox::invertedBox, const scalar mergeTol=1e-6)
Construct from cell values and point values.
Definition: isoSurface.C:1339
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: HashPtrTable.H:53
f
labelList f(nPoints)
Foam::Vector< scalar >
Foam::List< cellCutType >
Foam::FixedList< label, 3 >
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::boundBox
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:63
Foam::SlicedGeometricField
Specialization of GeometricField which holds slices of given complete fields in a form that they act ...
Definition: slicedSurfaceFieldsFwd.H:58
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:72
p0
const volScalarField & p0
Definition: EEqn.H:36
Foam::isoSurfaceTopo
Marching tet iso surface algorithm with optional filtering to keep only points originating from mesh ...
Definition: isoSurfaceTopo.H:59
Foam::GeometricField< vector, fvPatchField, volMesh >
Foam::fac::interpolate
static tmp< GeometricField< Type, faePatchField, edgeMesh > > interpolate(const GeometricField< Type, faPatchField, areaMesh > &tvf, const edgeScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Foam::isoSurfaceBase::filterType
filterType
The filtering (regularization) to apply.
Definition: isoSurfaceBase.H:87