src/sampling/sampledSurface/isoSurface/isoSurface.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::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
 111         autoPtr<slicedVolScalarField> cValsPtr_;
 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         void getNeighbour
 177         (
 178             const labelList& boundaryRegion,
 179             const volVectorField& meshC,
 180             const volScalarField& cVals,
 181             const label cellI,
 182             const label faceI,
 183             scalar& nbrValue,
 184             point& nbrPoint
 185         ) const;
 186
 187         //- Set faceCutType,cellCutType.
 188         void calcCutTypes
 189         (
 190             const labelList& boundaryRegion,
 191             const volVectorField& meshC,
 192             const volScalarField& cVals,
 193             const scalarField& pVals
 194         );
 195
 196         static labelPair findCommonPoints
 197         (
 198             const labelledTri&,
 199             const labelledTri&
 200         );
 201
 202         static point calcCentre(const triSurface&);
 203
 204         static pointIndexHit collapseSurface
 205         (
 206             pointField& localPoints,
 207             DynamicList<labelledTri, 64>& localTris
 208         );
 209
 210         //- Determine per cc whether all near cuts can be snapped to single
 211         //  point.
 212         void calcSnappedCc
 213         (
 214             const labelList& boundaryRegion,
 215             const volVectorField& meshC,
 216             const volScalarField& cVals,
 217             const scalarField& pVals,
 218             DynamicList<point>& snappedPoints,
 219             labelList& snappedCc
 220         ) const;
 221
 222         //- Determine per point whether all near cuts can be snapped to single
 223         //  point.
 224         void calcSnappedPoint
 225         (
 226             const PackedBoolList& isBoundaryPoint,
 227             const labelList& boundaryRegion,
 228             const volVectorField& meshC,
 229             const volScalarField& cVals,
 230             const scalarField& pVals,
 231             DynamicList<point>& snappedPoints,
 232             labelList& snappedPoint
 233         ) const;
 234
 235
 236         //- Return input field with coupled (and empty) patch values rewritten
 237         template<class Type>
 238         tmp<SlicedGeometricField
 239         <Type, fvPatchField, slicedFvPatchField, volMesh> >
 240         adaptPatchFields
 241         (
 242             const GeometricField<Type, fvPatchField, volMesh>& fld
 243         ) const;
 244
 245         //- Generate single point by interpolation or snapping
 246         template<class Type>
 247         Type generatePoint
 248         (
 249             const scalar s0,
 250             const Type& p0,
 251             const bool hasSnap0,
 252             const Type& snapP0,
 253
 254             const scalar s1,
 255             const Type& p1,
 256             const bool hasSnap1,
 257             const Type& snapP1
 258         ) const;
 259
 260         template<class Type>
 261         void generateTriPoints
 262         (
 263             const scalar s0,
 264             const Type& p0,
 265             const bool hasSnap0,
 266             const Type& snapP0,
 267
 268             const scalar s1,
 269             const Type& p1,
 270             const bool hasSnap1,
 271             const Type& snapP1,
 272
 273             const scalar s2,
 274             const Type& p2,
 275             const bool hasSnap2,
 276             const Type& snapP2,
 277
 278             const scalar s3,
 279             const Type& p3,
 280             const bool hasSnap3,
 281             const Type& snapP3,
 282
 283             DynamicList<Type>& points
 284         ) const;
 285
 286         template<class Type>
 287         label generateFaceTriPoints
 288         (
 289             const volScalarField& cVals,
 290             const scalarField& pVals,
 291
 292             const GeometricField<Type, fvPatchField, volMesh>& cCoords,
 293             const Field<Type>& pCoords,
 294
 295             const DynamicList<Type>& snappedPoints,
 296             const labelList& snappedCc,
 297             const labelList& snappedPoint,
 298             const label faceI,
 299
 300             const scalar neiVal,
 301             const Type& neiPt,
 302             const bool hasNeiSnap,
 303             const Type& neiSnapPt,
 304
 305             DynamicList<Type>& triPoints,
 306             DynamicList<label>& triMeshCells
 307         ) const;
 308
 309         template<class Type>
 310         void generateTriPoints
 311         (
 312             const volScalarField& cVals,
 313             const scalarField& pVals,
 314
 315             const GeometricField<Type, fvPatchField, volMesh>& cCoords,
 316             const Field<Type>& pCoords,
 317
 318             const DynamicList<Type>& snappedPoints,
 319             const labelList& snappedCc,
 320             const labelList& snappedPoint,
 321
 322             DynamicList<Type>& triPoints,
 323             DynamicList<label>& triMeshCells
 324         ) const;
 325
 326         triSurface stitchTriPoints
 327         (
 328             const bool checkDuplicates,
 329             const List<point>& triPoints,
 330             labelList& triPointReverseMap,  // unmerged to merged point
 331             labelList& triMap               // merged to unmerged triangle
 332         ) const;
 333
 334         //- Check single triangle for (topological) validity
 335         static bool validTri(const triSurface&, const label);
 336
 337         //- Determine edge-face addressing
 338         void calcAddressing
 339         (
 340             const triSurface& surf,
 341             List<FixedList<label, 3> >& faceEdges,
 342             labelList& edgeFace0,
 343             labelList& edgeFace1,
 344             Map<labelList>& edgeFacesRest
 345         ) const;
 346
 347         //- Determine orientation
 348         static void walkOrientation
 349         (
 350             const triSurface& surf,
 351             const List<FixedList<label, 3> >& faceEdges,
 352             const labelList& edgeFace0,
 353             const labelList& edgeFace1,
 354             const label seedTriI,
 355             labelList& flipState
 356         );
 357
 358         //- Orient surface
 359         static void orientSurface
 360         (
 361             triSurface&,
 362             const List<FixedList<label, 3> >& faceEdges,
 363             const labelList& edgeFace0,
 364             const labelList& edgeFace1,
 365             const Map<labelList>& edgeFacesRest
 366         );
 367
 368         //- Is triangle (given by 3 edges) not fully connected?
 369         static bool danglingTriangle
 370         (
 371             const FixedList<label, 3>& fEdges,
 372             const labelList& edgeFace1
 373         );
 374
 375         //- Mark all non-fully connected triangles
 376         static label markDanglingTriangles
 377         (
 378             const List<FixedList<label, 3> >& faceEdges,
 379             const labelList& edgeFace0,
 380             const labelList& edgeFace1,
 381             const Map<labelList>& edgeFacesRest,
 382             boolList& keepTriangles
 383         );
 384
 385         static triSurface subsetMesh
 386         (
 387             const triSurface& s,
 388             const labelList& newToOldFaces,
 389             labelList& oldToNewPoints,
 390             labelList& newToOldPoints
 391         );
 392
 393 public:
 394
 395     //- Runtime type information
 396     TypeName("isoSurface");
 397
 398
 399     // Constructors
 400
 401         //- Construct from cell values and point values. Uses boundaryField
 402         //  for boundary values. Holds reference to cellIsoVals and
 403         //  pointIsoVals.
 404         isoSurface
 405         (
 406             const volScalarField& cellIsoVals,
 407             const scalarField& pointIsoVals,
 408             const scalar iso,
 409             const bool regularise,
 410             const scalar mergeTol = 1E-6    // fraction of bounding box
 411         );
 412
 413
 414     // Member Functions
 415
 416         //- For every face original cell in mesh
 417         const labelList& meshCells() const
 418         {
 419             return meshCells_;
 420         }
 421
 422         //- For every unmerged triangle point the point in the triSurface
 423         const labelList& triPointMergeMap() const
 424         {
 425             return triPointMergeMap_;
 426         }
 427
 428         //- Interpolates cCoords,pCoords. Uses the references to the original
 429         //  fields used to create the iso surface.
 430         template <class Type>
 431         tmp<Field<Type> > interpolate
 432         (
 433             const GeometricField<Type, fvPatchField, volMesh>& cCoords,
 434             const Field<Type>& pCoords
 435         ) const;
 436
 437 };
 438
 439
 440 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 441
 442 } // End namespace Foam
 443
 444 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 445
 446 #ifdef NoRepository
 447 #   include "isoSurfaceTemplates.C"
 448 #endif
 449
 450 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 451
 452 #endif
 453
 454 // ************************************************************************* //