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) 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.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             //- 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         (
 219             pointField& localPoints,
 220             DynamicList<labelledTri, 64>& localTris
 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>
 251         tmp<SlicedGeometricField
 252         <Type, fvPatchField, slicedFvPatchField, volMesh> >
 253         adaptPatchFields
 254         (
 255             const GeometricField<Type, fvPatchField, volMesh>& fld
 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
 296             DynamicList<Type>& points
 297         ) const;
 298
 299         template<class Type>
 300         label generateFaceTriPoints
 301         (
 302             const volScalarField& cVals,
 303             const scalarField& pVals,
 304
 305             const GeometricField<Type, fvPatchField, volMesh>& cCoords,
 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
 328             const GeometricField<Type, fvPatchField, volMesh>& cCoords,
 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,
 354             List<FixedList<label, 3> >& faceEdges,
 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,
 364             const List<FixedList<label, 3> >& faceEdges,
 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&,
 375             const List<FixedList<label, 3> >& faceEdges,
 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         (
 391             const List<FixedList<label, 3> >& faceEdges,
 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
 436         const labelList& triPointMergeMap() const
 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>
 444         tmp<Field<Type> > interpolate
 445         (
 446             const GeometricField<Type, fvPatchField, volMesh>& cCoords,
 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 // ************************************************************************* //