src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-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::isoSurfaceCell
  26
  27 Description
  28     A surface formed by the iso value.
  29     After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke
  30     and
  31     "Regularised Marching Tetrahedra: improved iso-surface extraction",
  32     G.M. Treece, R.W. Prager and A.H. Gee.
  33
  34     See isoSurface. This is a variant which does tetrahedrisation from
  35     triangulation of face to cell centre instead of edge of face to two
  36     neighbouring cell centres. This gives much lower quality triangles
  37     but they are local to a cell.
  38
  39 SourceFiles
  40     isoSurfaceCell.C
  41
  42 \*---------------------------------------------------------------------------*/
  43
  44 #ifndef isoSurfaceCell_H
  45 #define isoSurfaceCell_H
  46
  47 #include "triSurface.H"
  48 #include "labelPair.H"
  49 #include "pointIndexHit.H"
  50 #include "PackedBoolList.H"
  51
  52 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  53
  54 namespace Foam
  55 {
  56
  57 class polyMesh;
  58
  59 /*---------------------------------------------------------------------------*\
  60                        Class isoSurfaceCell Declaration
  61 \*---------------------------------------------------------------------------*/
  62
  63 class isoSurfaceCell
  64 :
  65     public triSurface
  66 {
  67     // Private data
  68
  69         enum segmentCutType
  70         {
  71             NEARFIRST,      // intersection close to e.first()
  72             NEARSECOND,     //  ,,                   e.second()
  73             NOTNEAR         // no intersection
  74         };
  75
  76         enum cellCutType
  77         {
  78             NOTCUT,     // not cut
  79             SPHERE,     // all edges to cell centre cut
  80             CUT         // normal cut
  81         };
  82
  83
  84         //- Reference to mesh
  85         const polyMesh& mesh_;
  86
  87         const scalarField& cVals_;
  88
  89         const scalarField& pVals_;
  90
  91         //- isoSurfaceCell value
  92         const scalar iso_;
  93
  94         //- When to merge points
  95         const scalar mergeDistance_;
  96
  97         //- Whether cell might be cut
  98         List<cellCutType> cellCutType_;
  99
 100         //- Estimated number of cut cells
 101         label nCutCells_;
 102
 103         //- For every triangle the original cell in mesh
 104         labelList meshCells_;
 105
 106         //- For every unmerged triangle point the point in the triSurface
 107         labelList triPointMergeMap_;
 108
 109
 110     // Private Member Functions
 111
 112         //- Get location of iso value as fraction inbetween s0,s1
 113         scalar isoFraction
 114         (
 115             const scalar s0,
 116             const scalar s1
 117         ) const;
 118
 119         //List<triFace> triangulate(const face& f) const;
 120
 121         //- Determine whether cell is cut
 122         cellCutType calcCutType
 123         (
 124             const PackedBoolList&,
 125             const scalarField& cellValues,
 126             const scalarField& pointValues,
 127             const label
 128         ) const;
 129
 130         void calcCutTypes
 131         (
 132             const PackedBoolList&,
 133             const scalarField& cellValues,
 134             const scalarField& pointValues
 135         );
 136
 137         static labelPair findCommonPoints
 138         (
 139             const labelledTri&,
 140             const labelledTri&
 141         );
 142
 143         static point calcCentre(const triSurface&);
 144
 145         static pointIndexHit collapseSurface
 146         (
 147             pointField& localPoints,
 148             DynamicList<labelledTri, 64>& localTris
 149         );
 150
 151         //- Determine per cc whether all near cuts can be snapped to single
 152         //  point.
 153         void calcSnappedCc
 154         (
 155             const PackedBoolList& isTet,
 156             const scalarField& cVals,
 157             const scalarField& pVals,
 158             DynamicList<point>& triPoints,
 159             labelList& snappedCc
 160         ) const;
 161
 162         //- Generate triangles for face connected to pointI
 163         void genPointTris
 164         (
 165             const scalarField& cellValues,
 166             const scalarField& pointValues,
 167             const label pointI,
 168             const face& f,
 169             const label cellI,
 170             DynamicList<point, 64>& localTriPoints
 171         ) const;
 172
 173         //- Generate triangles for tet connected to pointI
 174         void genPointTris
 175         (
 176             const scalarField& pointValues,
 177             const label pointI,
 178             const label cellI,
 179             DynamicList<point, 64>& localTriPoints
 180         ) const;
 181
 182         //- Determine per point whether all near cuts can be snapped to single
 183         //  point.
 184         void calcSnappedPoint
 185         (
 186             const PackedBoolList& isBoundaryPoint,
 187             const PackedBoolList& isTet,
 188             const scalarField& cVals,
 189             const scalarField& pVals,
 190             DynamicList<point>& triPoints,
 191             labelList& snappedPoint
 192         ) const;
 193
 194         //- Generate single point by interpolation or snapping
 195         template<class Type>
 196         Type generatePoint
 197         (
 198             const DynamicList<Type>& snappedPoints,
 199             const scalar s0,
 200             const Type& p0,
 201             const label p0Index,
 202             const scalar s1,
 203             const Type& p1,
 204             const label p1Index
 205         ) const;
 206
 207         template<class Type>
 208         void generateTriPoints
 209         (
 210             const DynamicList<Type>& snapped,
 211             const scalar s0,
 212             const Type& p0,
 213             const label p0Index,
 214             const scalar s1,
 215             const Type& p1,
 216             const label p1Index,
 217             const scalar s2,
 218             const Type& p2,
 219             const label p2Index,
 220             const scalar s3,
 221             const Type& p3,
 222             const label p3Index,
 223             DynamicList<Type>& points
 224         ) const;
 225
 226         template<class Type>
 227         void generateTriPoints
 228         (
 229             const scalarField& cVals,
 230             const scalarField& pVals,
 231
 232             const Field<Type>& cCoords,
 233             const Field<Type>& pCoords,
 234
 235             const DynamicList<Type>& snappedPoints,
 236             const labelList& snappedCc,
 237             const labelList& snappedPoint,
 238
 239             DynamicList<Type>& triPoints,
 240             DynamicList<label>& triMeshCells
 241         ) const;
 242
 243         triSurface stitchTriPoints
 244         (
 245             const bool checkDuplicates,
 246             const List<point>& triPoints,
 247             labelList& triPointReverseMap,  // unmerged to merged point
 248             labelList& triMap               // merged to unmerged triangle
 249         ) const;
 250
 251         //- Check single triangle for (topological) validity
 252         static bool validTri(const triSurface&, const label);
 253
 254         //- Determine edge-face addressing
 255         void calcAddressing
 256         (
 257             const triSurface& surf,
 258             List<FixedList<label, 3> >& faceEdges,
 259             labelList& edgeFace0,
 260             labelList& edgeFace1,
 261             Map<labelList>& edgeFacesRest
 262         ) const;
 263
 264         //- Determine orientation
 265         static void walkOrientation
 266         (
 267             const triSurface& surf,
 268             const List<FixedList<label, 3> >& faceEdges,
 269             const labelList& edgeFace0,
 270             const labelList& edgeFace1,
 271             const label seedTriI,
 272             labelList& flipState
 273         );
 274
 275         //- Orient surface
 276         static void orientSurface
 277         (
 278             triSurface&,
 279             const List<FixedList<label, 3> >& faceEdges,
 280             const labelList& edgeFace0,
 281             const labelList& edgeFace1,
 282             const Map<labelList>& edgeFacesRest
 283         );
 284
 285         //- Is triangle (given by 3 edges) not fully connected?
 286         static bool danglingTriangle
 287         (
 288             const FixedList<label, 3>& fEdges,
 289             const labelList& edgeFace1
 290         );
 291
 292         //- Mark all non-fully connected triangles
 293         static label markDanglingTriangles
 294         (
 295             const List<FixedList<label, 3> >& faceEdges,
 296             const labelList& edgeFace0,
 297             const labelList& edgeFace1,
 298             const Map<labelList>& edgeFacesRest,
 299             boolList& keepTriangles
 300         );
 301
 302         static triSurface subsetMesh
 303         (
 304             const triSurface& s,
 305             const labelList& newToOldFaces,
 306             labelList& oldToNewPoints,
 307             labelList& newToOldPoints
 308         );
 309
 310         //- Combine all triangles inside a cell into a minimal triangulation
 311         void combineCellTriangles();
 312
 313 public:
 314
 315     //- Runtime type information
 316     TypeName("isoSurfaceCell");
 317
 318
 319     // Constructors
 320
 321         //- Construct from dictionary
 322         isoSurfaceCell
 323         (
 324             const polyMesh& mesh,
 325             const scalarField& cellValues,
 326             const scalarField& pointValues,
 327             const scalar iso,
 328             const bool regularise,
 329             const scalar mergeTol = 1E-6    // fraction of bounding box
 330         );
 331
 332
 333     // Member Functions
 334
 335         //- For every face original cell in mesh
 336         const labelList& meshCells() const
 337         {
 338             return meshCells_;
 339         }
 340
 341         //- For every unmerged triangle point the point in the triSurface
 342         const labelList triPointMergeMap() const
 343         {
 344             return triPointMergeMap_;
 345         }
 346
 347
 348         //- Interpolates cCoords,pCoords. Takes the original fields
 349         //  used to create the iso surface.
 350         template <class Type>
 351         tmp<Field<Type> > interpolate
 352         (
 353             const scalarField& cVals,
 354             const scalarField& pVals,
 355             const Field<Type>& cCoords,
 356             const Field<Type>& pCoords
 357         ) const;
 358
 359         //- Interpolates cCoords,pCoords.
 360         template <class Type>
 361         tmp<Field<Type> > interpolate
 362         (
 363             const Field<Type>& cCoords,
 364             const Field<Type>& pCoords
 365         ) const;
 366 };
 367
 368
 369 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 370
 371 } // End namespace Foam
 372
 373 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 374
 375 #ifdef NoRepository
 376 #   include "isoSurfaceCellTemplates.C"
 377 #endif
 378
 379 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 380
 381 #endif
 382
 383 // ************************************************************************* //