FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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/triSurface.H>
67 #include <OpenFOAM/labelPair.H>
70 #include <finiteVolume/volFields.H>
72 
73 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
74 
75 namespace Foam
76 {
77 
78 class fvMesh;
79 
80 /*---------------------------------------------------------------------------*\
81  Class isoSurface Declaration
82 \*---------------------------------------------------------------------------*/
83 
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  //- Take value at local point patchPointI and assign it to its
153  // correct place in patchValues (for transferral) and sharedValues
154  // (for reduction)
155  void insertPointData
156  (
157  const processorPolyPatch& pp,
158  const Map<label>& meshToShared,
159  const pointField& pointValues,
160  const label patchPointI,
161  pointField& patchValues,
162  pointField& sharedValues
163  ) const;
164 
165  //- Synchonise points on all non-separated coupled patches
166  void syncUnseparatedPoints
167  (
168  pointField& collapsedPoint,
169  const point& nullValue
170  ) const;
171 
172 
173  //- Get location of iso value as fraction inbetween s0,s1
174  scalar isoFraction
175  (
176  const scalar s0,
177  const scalar s1
178  ) const;
179 
180  //- Check if any edge of a face is cut
181  bool isEdgeOfFaceCut
182  (
183  const scalarField& pVals,
184  const face& f,
185  const bool ownLower,
186  const bool neiLower
187  ) const;
188 
189  void getNeighbour
190  (
191  const labelList& boundaryRegion,
192  const volVectorField& meshC,
193  const volScalarField& cVals,
194  const label cellI,
195  const label faceI,
196  scalar& nbrValue,
197  point& nbrPoint
198  ) const;
199 
200  //- Set faceCutType,cellCutType.
201  void calcCutTypes
202  (
203  const labelList& boundaryRegion,
204  const volVectorField& meshC,
205  const volScalarField& cVals,
206  const scalarField& pVals
207  );
208 
209  static labelPair findCommonPoints
210  (
211  const labelledTri&,
212  const labelledTri&
213  );
214 
215  static point calcCentre(const triSurface&);
216 
217  static pointIndexHit collapseSurface
218  (
221  );
222 
223  //- Determine per cc whether all near cuts can be snapped to single
224  // point.
225  void calcSnappedCc
226  (
227  const labelList& boundaryRegion,
228  const volVectorField& meshC,
229  const volScalarField& cVals,
230  const scalarField& pVals,
231  DynamicList<point>& snappedPoints,
232  labelList& snappedCc
233  ) const;
234 
235  //- Determine per point whether all near cuts can be snapped to single
236  // point.
237  void calcSnappedPoint
238  (
239  const PackedBoolList& isBoundaryPoint,
240  const labelList& boundaryRegion,
241  const volVectorField& meshC,
242  const volScalarField& cVals,
243  const scalarField& pVals,
244  DynamicList<point>& snappedPoints,
245  labelList& snappedPoint
246  ) const;
247 
248 
249  //- Return input field with coupled (and empty) patch values rewritten
250  template<class Type>
253  adaptPatchFields
254  (
256  ) const;
257 
258  //- Generate single point by interpolation or snapping
259  template<class Type>
260  Type generatePoint
261  (
262  const scalar s0,
263  const Type& p0,
264  const bool hasSnap0,
265  const Type& snapP0,
266 
267  const scalar s1,
268  const Type& p1,
269  const bool hasSnap1,
270  const Type& snapP1
271  ) const;
272 
273  template<class Type>
274  void generateTriPoints
275  (
276  const scalar s0,
277  const Type& p0,
278  const bool hasSnap0,
279  const Type& snapP0,
280 
281  const scalar s1,
282  const Type& p1,
283  const bool hasSnap1,
284  const Type& snapP1,
285 
286  const scalar s2,
287  const Type& p2,
288  const bool hasSnap2,
289  const Type& snapP2,
290 
291  const scalar s3,
292  const Type& p3,
293  const bool hasSnap3,
294  const Type& snapP3,
295 
297  ) const;
298 
299  template<class Type>
300  label generateFaceTriPoints
301  (
302  const volScalarField& cVals,
303  const scalarField& pVals,
304 
306  const Field<Type>& pCoords,
307 
308  const DynamicList<Type>& snappedPoints,
309  const labelList& snappedCc,
310  const labelList& snappedPoint,
311  const label faceI,
312 
313  const scalar neiVal,
314  const Type& neiPt,
315  const bool hasNeiSnap,
316  const Type& neiSnapPt,
317 
318  DynamicList<Type>& triPoints,
319  DynamicList<label>& triMeshCells
320  ) const;
321 
322  template<class Type>
323  void generateTriPoints
324  (
325  const volScalarField& cVals,
326  const scalarField& pVals,
327 
329  const Field<Type>& pCoords,
330 
331  const DynamicList<Type>& snappedPoints,
332  const labelList& snappedCc,
333  const labelList& snappedPoint,
334 
335  DynamicList<Type>& triPoints,
336  DynamicList<label>& triMeshCells
337  ) const;
338 
339  triSurface stitchTriPoints
340  (
341  const bool checkDuplicates,
342  const List<point>& triPoints,
343  labelList& triPointReverseMap, // unmerged to merged point
344  labelList& triMap // merged to unmerged triangle
345  ) const;
346 
347  //- Check single triangle for (topological) validity
348  static bool validTri(const triSurface&, const label);
349 
350  //- Determine edge-face addressing
351  void calcAddressing
352  (
353  const triSurface& surf,
355  labelList& edgeFace0,
356  labelList& edgeFace1,
357  Map<labelList>& edgeFacesRest
358  ) const;
359 
360  //- Determine orientation
361  static void walkOrientation
362  (
363  const triSurface& surf,
365  const labelList& edgeFace0,
366  const labelList& edgeFace1,
367  const label seedTriI,
368  labelList& flipState
369  );
370 
371  //- Orient surface
372  static void orientSurface
373  (
374  triSurface&,
376  const labelList& edgeFace0,
377  const labelList& edgeFace1,
378  const Map<labelList>& edgeFacesRest
379  );
380 
381  //- Is triangle (given by 3 edges) not fully connected?
382  static bool danglingTriangle
383  (
384  const FixedList<label, 3>& fEdges,
385  const labelList& edgeFace1
386  );
387 
388  //- Mark all non-fully connected triangles
389  static label markDanglingTriangles
390  (
392  const labelList& edgeFace0,
393  const labelList& edgeFace1,
394  const Map<labelList>& edgeFacesRest,
395  boolList& keepTriangles
396  );
397 
398  static triSurface subsetMesh
399  (
400  const triSurface& s,
401  const labelList& newToOldFaces,
402  labelList& oldToNewPoints,
403  labelList& newToOldPoints
404  );
405 
406 public:
407 
408  //- Runtime type information
409  TypeName("isoSurface");
410 
411 
412  // Constructors
413 
414  //- Construct from cell values and point values. Uses boundaryField
415  // for boundary values. Holds reference to cellIsoVals and
416  // pointIsoVals.
417  isoSurface
418  (
419  const volScalarField& cellIsoVals,
420  const scalarField& pointIsoVals,
421  const scalar iso,
422  const bool regularise,
423  const scalar mergeTol = 1E-6 // fraction of bounding box
424  );
425 
426 
427  // Member Functions
428 
429  //- For every face original cell in mesh
430  const labelList& meshCells() const
431  {
432  return meshCells_;
433  }
434 
435  //- For every unmerged triangle point the point in the triSurface
437  {
438  return triPointMergeMap_;
439  }
440 
441  //- Interpolates cCoords,pCoords. Uses the references to the original
442  // fields used to create the iso surface.
443  template <class Type>
445  (
447  const Field<Type>& pCoords
448  ) const;
449 
450 };
451 
452 
453 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
454 
455 } // End namespace Foam
456 
457 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
458 
459 #ifdef NoRepository
460 # include "isoSurfaceTemplates.C"
461 #endif
462 
463 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
464 
465 #endif
466 
467 // ************************ vim: set sw=4 sts=4 et: ************************ //