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 | Copyright (C) 2011-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Class
25  Foam::isoSurface
26 
27 Description
28  A surface formed by the iso value.
29  After "Regularised Marching Tetrahedra: improved iso-surface extraction",
30  G.M. Treece, R.W. Prager and A.H. Gee.
31 
32  Note:
33  - does tets without using cell centres/cell values. Not tested.
34  - regularisation can create duplicate triangles/non-manifold surfaces.
35  Current handling of those is bit ad-hoc for now and not perfect.
36  - regularisation does not do boundary points so as to preserve the
37  boundary perfectly.
38  - uses geometric merge with fraction of bounding box as distance.
39  - triangles can be between two cell centres so constant sampling
40  does not make sense.
41  - on empty patches behaves like zero gradient.
42  - does not do 2D correctly, creates non-flat iso surface.
43  - on processor boundaries might two overlapping (identical) triangles
44  (one from either side)
45 
46  The handling on coupled patches is a bit complex. All fields
47  (values and coordinates) get rewritten so
48  - empty patches get zerogradient (value) and facecentre (coordinates)
49  - separated processor patch faces get interpolate (value) and facecentre
50  (coordinates). (this is already the default for cyclics)
51  - non-separated processor patch faces keep opposite value and cell centre
52 
53  Now the triangle generation on non-separated processor patch faces
54  can use the neighbouring value. Any separated processor face or cyclic
55  face gets handled just like any boundary face.
56 
57 
58 SourceFiles
59  isoSurface.C
60 
61 \*---------------------------------------------------------------------------*/
62 
63 #ifndef isoSurface_H
64 #define isoSurface_H
65 
66 #include "triSurface.H"
67 #include "labelPair.H"
68 #include "pointIndexHit.H"
69 #include "PackedBoolList.H"
70 #include "volFields.H"
71 #include "slicedVolFields.H"
72 
73 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
74 
75 namespace Foam
76 {
77 
78 class fvMesh;
79 
80 /*---------------------------------------------------------------------------*\
81  Class isoSurface Declaration
82 \*---------------------------------------------------------------------------*/
83 
84 class isoSurface
85 :
86  public triSurface
87 {
88  // Private data
89 
90  enum segmentCutType
91  {
92  NEARFIRST, // intersection close to e.first()
93  NEARSECOND, // ,, e.second()
94  NOTNEAR // no intersection
95  };
96 
97  enum cellCutType
98  {
99  NOTCUT, // not cut
100  SPHERE, // all edges to cell centre cut
101  CUT // normal cut
102  };
103 
104 
105  //- Reference to mesh
106  const fvMesh& mesh_;
107 
108  const scalarField& pVals_;
109 
110  //- Input volScalarField with separated coupled patches rewritten
112 
113  //- Isosurface value
114  const scalar iso_;
115 
116  //- Regularise?
117  const Switch regularise_;
118 
119  //- When to merge points
120  const scalar mergeDistance_;
121 
122 
123  //- Whether face might be cut
124  List<cellCutType> faceCutType_;
125 
126  //- Whether cell might be cut
127  List<cellCutType> cellCutType_;
128 
129  //- Estimated number of cut cells
130  label nCutCells_;
131 
132  //- For every triangle the original cell in mesh
133  labelList meshCells_;
134 
135  //- For every unmerged triangle point the point in the triSurface
136  labelList triPointMergeMap_;
137 
138 
139  // Private Member Functions
140 
141  // Point synchronisation
142 
143  //- Does tensor differ (to within mergeTolerance) from identity
144  bool noTransform(const tensor& tt) const;
145 
146  //- Is patch a collocated (i.e. non-separated) coupled patch?
147  static bool collocatedPatch(const polyPatch&);
148 
149  //- Per face whether is collocated
150  PackedBoolList collocatedFaces(const coupledPolyPatch&) const;
151 
152  //- Synchonise points on all non-separated coupled patches
153  void syncUnseparatedPoints
154  (
155  pointField& collapsedPoint,
156  const point& nullValue
157  ) const;
158 
159 
160  //- Get location of iso value as fraction inbetween s0,s1
161  scalar isoFraction
162  (
163  const scalar s0,
164  const scalar s1
165  ) const;
166 
167  //- Check if any edge of a face is cut
168  bool isEdgeOfFaceCut
169  (
170  const scalarField& pVals,
171  const face& f,
172  const bool ownLower,
173  const bool neiLower
174  ) const;
175 
176  //- Get neighbour value and position.
177  void getNeighbour
178  (
179  const labelList& boundaryRegion,
180  const volVectorField& meshC,
181  const volScalarField& cVals,
182  const label celli,
183  const label facei,
184  scalar& nbrValue,
185  point& nbrPoint
186  ) const;
187 
188  //- Determine for every face/cell whether it (possibly) generates
189  // triangles.
190  void calcCutTypes
191  (
192  const labelList& boundaryRegion,
193  const volVectorField& meshC,
194  const volScalarField& cVals,
195  const scalarField& pVals
196  );
197 
198  static point calcCentre(const triSurface&);
199 
200  //- Determine per cc whether all near cuts can be snapped to single
201  // point.
202  void calcSnappedCc
203  (
204  const labelList& boundaryRegion,
205  const volVectorField& meshC,
206  const volScalarField& cVals,
207  const scalarField& pVals,
208  DynamicList<point>& snappedPoints,
209  labelList& snappedCc
210  ) const;
211 
212  //- Determine per point whether all near cuts can be snapped to single
213  // point.
214  void calcSnappedPoint
215  (
216  const PackedBoolList& isBoundaryPoint,
217  const labelList& boundaryRegion,
218  const volVectorField& meshC,
219  const volScalarField& cVals,
220  const scalarField& pVals,
221  DynamicList<point>& snappedPoints,
222  labelList& snappedPoint
223  ) const;
224 
225 
226  //- Return input field with coupled (and empty) patch values rewritten
227  template<class Type>
230  adaptPatchFields
231  (
233  ) const;
234 
235  //- Generate single point by interpolation or snapping
236  template<class Type>
237  Type generatePoint
238  (
239  const scalar s0,
240  const Type& p0,
241  const bool hasSnap0,
242  const Type& snapP0,
243 
244  const scalar s1,
245  const Type& p1,
246  const bool hasSnap1,
247  const Type& snapP1
248  ) const;
249 
250 
251  //- Note: cannot use simpler isoSurfaceCell::generateTriPoints since
252  // the need here to sometimes pass in remote 'snappedPoints'
253  template<class Type>
254  void generateTriPoints
255  (
256  const scalar s0,
257  const Type& p0,
258  const bool hasSnap0,
259  const Type& snapP0,
260 
261  const scalar s1,
262  const Type& p1,
263  const bool hasSnap1,
264  const Type& snapP1,
265 
266  const scalar s2,
267  const Type& p2,
268  const bool hasSnap2,
269  const Type& snapP2,
270 
271  const scalar s3,
272  const Type& p3,
273  const bool hasSnap3,
274  const Type& snapP3,
275 
277  ) const;
278 
279  template<class Type>
280  label generateFaceTriPoints
281  (
282  const volScalarField& cVals,
283  const scalarField& pVals,
284 
286  const Field<Type>& pCoords,
287 
288  const DynamicList<Type>& snappedPoints,
289  const labelList& snappedCc,
290  const labelList& snappedPoint,
291  const label facei,
292 
293  const scalar neiVal,
294  const Type& neiPt,
295  const bool hasNeiSnap,
296  const Type& neiSnapPt,
297 
298  DynamicList<Type>& triPoints,
299  DynamicList<label>& triMeshCells
300  ) const;
301 
302  template<class Type>
303  void generateTriPoints
304  (
305  const volScalarField& cVals,
306  const scalarField& pVals,
307 
309  const Field<Type>& pCoords,
310 
311  const DynamicList<Type>& snappedPoints,
312  const labelList& snappedCc,
313  const labelList& snappedPoint,
314 
315  DynamicList<Type>& triPoints,
316  DynamicList<label>& triMeshCells
317  ) const;
318 
319  template<class Type>
320  static tmp<Field<Type>> interpolate
321  (
322  const label nPoints,
323  const labelList& triPointMergeMap,
324  const DynamicList<Type>& unmergedValues
325  );
326 
327  triSurface stitchTriPoints
328  (
329  const bool checkDuplicates,
330  const List<point>& triPoints,
331  labelList& triPointReverseMap, // unmerged to merged point
332  labelList& triMap // merged to unmerged triangle
333  ) const;
334 
335  //- Check single triangle for (topological) validity
336  static bool validTri(const triSurface&, const label);
337 
338  static triSurface subsetMesh
339  (
340  const triSurface& s,
341  const labelList& newToOldFaces,
342  labelList& oldToNewPoints,
343  labelList& newToOldPoints
344  );
345 
346 public:
347 
348  //- Declare friendship with isoSurfaceCell to share some functionality
349  friend class isoSurfaceCell;
350 
351 
352  //- Runtime type information
353  TypeName("isoSurface");
354 
355 
356  // Constructors
357 
358  //- Construct from cell values and point values. Uses boundaryField
359  // for boundary values. Holds reference to cellIsoVals and
360  // pointIsoVals.
361  isoSurface
362  (
363  const volScalarField& cellIsoVals,
364  const scalarField& pointIsoVals,
365  const scalar iso,
366  const bool regularise,
367  const scalar mergeTol = 1e-6 // fraction of bounding box
368  );
369 
370 
371  // Member Functions
372 
373  //- For every triangle the original cell in mesh
374  const labelList& meshCells() const
375  {
376  return meshCells_;
377  }
378 
379  //- Interpolates cCoords,pCoords. Uses the references to the original
380  // fields used to create the iso surface.
381  template<class Type>
383  (
385  const Field<Type>& pCoords
386  ) const;
387 
388 };
389 
390 
391 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
392 
393 } // End namespace Foam
394 
395 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
396 
397 #ifdef NoRepository
398  #include "isoSurfaceTemplates.C"
399 #endif
400 
401 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
402 
403 #endif
404 
405 // ************************************************************************* //
Specialization of fvPatchField which creates the underlying fvPatchField as a slice of the given comp...
label nPoints() const
Return number of points supporting patch faces.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const double e
Elementary charge.
Definition: doubleFloat.H:78
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
Definition: Switch.H:60
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:66
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
isoSurface(const volScalarField &cellIsoVals, const scalarField &pointIsoVals, const scalar iso, const bool regularise, const scalar mergeTol=1e-6)
Construct from cell values and point values. Uses boundaryField.
Definition: isoSurface.C:1175
const labelList & meshCells() const
For every triangle the original cell in mesh.
Definition: isoSurface.H:373
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Mesh data needed to do the Finite Volume discretisation.
Definition: volMesh.H:53
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< ' ';}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< ' ';}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< ' ';}gmvFile<< nl;forAll(lagrangianScalarNames, i){ const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
const Field< point > & points() const
Return reference to global points.
A surface formed by the iso value. After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke (http://paulbourke.net/geometry/polygonise) and "Regularised Marching Tetrahedra: improved iso-surface extraction", G.M. Treece, R.W. Prager and A.H. Gee.
labelList f(nPoints)
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
TypeName("isoSurface")
Runtime type information.
A bit-packed bool list.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
A surface formed by the iso value. After "Regularised Marching Tetrahedra: improved iso-surface extra...
Definition: isoSurface.H:83
A class for managing temporary objects.
Definition: PtrList.H:53
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Triangulated surface description with patch information.
Definition: triSurface.H:65
The boundaryRegion persistent data saved as a Map<dictionary>.
Specialization of GeometricField which holds slices of given complete fields in a form that they act ...
Namespace for OpenFOAM.