isoSurface.H

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2016-2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.
 
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
Class
    Foam::isoSurface
 
Description
    A surface formed by the iso value.
    After "Regularised Marching Tetrahedra: improved iso-surface extraction",
    G.M. Treece, R.W. Prager and A.H. Gee.
 
    Note:
    - does tets without using cell centres/cell values. Not tested.
    - regularisation can create duplicate triangles/non-manifold surfaces.
    Current handling of those is bit ad-hoc for now and not perfect.
    - regularisation does not do boundary points so as to preserve the
      boundary perfectly.
    - uses geometric merge with fraction of bounding box as distance.
    - triangles can be between two cell centres so constant sampling
      does not make sense.
    - on empty patches behaves like zero gradient.
    - does not do 2D correctly, creates non-flat iso surface.
    - on processor boundaries might two overlapping (identical) triangles
      (one from either side)
 
    The handling on coupled patches is a bit complex. All fields
    (values and coordinates) get rewritten so
    - empty patches get zerogradient (value) and facecentre (coordinates)
    - separated processor patch faces get interpolate (value) and facecentre
      (coordinates). (this is already the default for cyclics)
    - non-separated processor patch faces keep opposite value and cell centre
 
    Now the triangle generation on non-separated processor patch faces
    can use the neighbouring value. Any separated processor face or cyclic
    face gets handled just like any boundary face.
 
SourceFiles
    isoSurface.C
    isoSurfaceTemplates.C
 
\*---------------------------------------------------------------------------*/
 
#ifndef isoSurface_H
#define isoSurface_H
 
#include "bitSet.H"
#include "volFields.H"
#include "slicedVolFields.H"
#include "isoSurfaceBase.H"
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
namespace Foam
{
 
// Forward declarations
 
class fvMesh;
class plane;
class treeBoundBox;
class triSurface;
 
/*---------------------------------------------------------------------------*\
                       Class isoSurface Declaration
\*---------------------------------------------------------------------------*/
 
class isoSurface
:
    public isoSurfaceBase
{
    // Private data
 
        enum segmentCutType
        {
            NEARFIRST,      // intersection close to e.first()
            NEARSECOND,     //  ,,                   e.second()
            NOTNEAR         // no intersection
        };
 
        enum cellCutType
        {
            NOTCUT,     // not cut
            SPHERE,     // all edges to cell centre cut
            CUT         // normal cut
        };
 
 
        //- Reference to mesh
        const fvMesh& mesh_;
 
        const scalarField& pVals_;
 
        //- Input volScalarField with separated coupled patches rewritten
        autoPtr<slicedVolScalarField> cValsPtr_;
 
        //- Regularise?
        const bool regularise_;
 
        //- When to merge points
        const scalar mergeDistance_;
 
        //- Whether face might be cut
        List<cellCutType> faceCutType_;
 
        //- Whether cell might be cut
        List<cellCutType> cellCutType_;
 
        //- Estimated number of cut cells
        label nCutCells_;
 
        //- For every unmerged triangle point the point in the triSurface
        labelList triPointMergeMap_;
 
        //- triSurface points that have weighted interpolation
        DynamicList<label> interpolatedPoints_;
 
        //- corresponding original, unmerged points
        DynamicList<FixedList<label, 3>> interpolatedOldPoints_;
 
        //- corresponding weights
        DynamicList<FixedList<scalar, 3>> interpolationWeights_;
 
 
    // Private Member Functions
 
        // Point synchronisation
 
            //- Does tensor differ (to within mergeTolerance) from identity
            bool noTransform(const tensor& tt) const;
 
            //- Is patch a collocated (i.e. non-separated) coupled patch?
            static bool collocatedPatch(const polyPatch& pp);
 
            //- Per face whether is collocated
            bitSet collocatedFaces(const coupledPolyPatch&) const;
 
            //- Synchonise points on all non-separated coupled patches
            void syncUnseparatedPoints
            (
                pointField& collapsedPoint,
                const point& nullValue
            ) const;
 
 
        //- Get location of iso value as fraction inbetween s0,s1
        scalar isoFraction
        (
            const scalar s0,
            const scalar s1
        ) const;
 
        //- Check if any edge of a face is cut
        bool isEdgeOfFaceCut
        (
            const scalarField& pVals,
            const face& f,
            const bool ownLower,
            const bool neiLower
        ) const;
 
        //- Get neighbour value and position.
        void getNeighbour
        (
            const labelList& boundaryRegion,
            const volVectorField& meshC,
            const volScalarField& cVals,
            const label celli,
            const label facei,
            scalar& nbrValue,
            point& nbrPoint
        ) const;
 
        //- Determine for every face/cell whether it (possibly) generates
        //  triangles.
        void calcCutTypes
        (
            const labelList& boundaryRegion,
            const volVectorField& meshC,
            const volScalarField& cVals,
            const scalarField& pVals
        );
 
        static point calcCentre(const triSurface&);
 
        //- Determine per cc whether all near cuts can be snapped to single
        //  point.
        void calcSnappedCc
        (
            const labelList& boundaryRegion,
            const volVectorField& meshC,
            const volScalarField& cVals,
            const scalarField& pVals,
            DynamicList<point>& snappedPoints,
            labelList& snappedCc
        ) const;
 
        //- Determine per point whether all near cuts can be snapped to single
        //  point.
        void calcSnappedPoint
        (
            const bitSet& isBoundaryPoint,
            const labelList& boundaryRegion,
            const volVectorField& meshC,
            const volScalarField& cVals,
            const scalarField& pVals,
            DynamicList<point>& snappedPoints,
            labelList& snappedPoint
        ) const;
 
 
        //- Return input field with coupled (and empty) patch values rewritten
        template<class Type>
        tmp<SlicedGeometricField
        <Type, fvPatchField, slicedFvPatchField, volMesh>>
        adaptPatchFields
        (
            const GeometricField<Type, fvPatchField, volMesh>& fld
        ) const;
 
        //- Generate single point by interpolation or snapping
        template<class Type>
        Type generatePoint
        (
            const scalar s0,
            const Type& p0,
            const bool hasSnap0,
            const Type& snapP0,
 
            const scalar s1,
            const Type& p1,
            const bool hasSnap1,
            const Type& snapP1
        ) const;
 
 
        //- Note: cannot use simpler isoSurfaceCell::generateTriPoints since
        //  the need here to sometimes pass in remote 'snappedPoints'
        template<class Type>
        void generateTriPoints
        (
            const scalar s0,
            const Type& p0,
            const bool hasSnap0,
            const Type& snapP0,
 
            const scalar s1,
            const Type& p1,
            const bool hasSnap1,
            const Type& snapP1,
 
            const scalar s2,
            const Type& p2,
            const bool hasSnap2,
            const Type& snapP2,
 
            const scalar s3,
            const Type& p3,
            const bool hasSnap3,
            const Type& snapP3,
 
            DynamicList<Type>& pts
        ) const;
 
        template<class Type>
        label generateFaceTriPoints
        (
            const volScalarField& cVals,
            const scalarField& pVals,
 
            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
            const Field<Type>& pCoords,
 
            const DynamicList<Type>& snappedPoints,
            const labelList& snappedCc,
            const labelList& snappedPoint,
            const label facei,
 
            const scalar neiVal,
            const Type& neiPt,
            const bool hasNeiSnap,
            const Type& neiSnapPt,
 
            DynamicList<Type>& triPoints,
            DynamicList<label>& triMeshCells
        ) const;
 
        template<class Type>
        void generateTriPoints
        (
            const volScalarField& cVals,
            const scalarField& pVals,
 
            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
            const Field<Type>& pCoords,
 
            const DynamicList<Type>& snappedPoints,
            const labelList& snappedCc,
            const labelList& snappedPoint,
 
            DynamicList<Type>& triPoints,
            DynamicList<label>& triMeshCells
        ) const;
 
        template<class Type>
        static tmp<Field<Type>> interpolate
        (
            const label nPoints,
            const labelList& triPointMergeMap,
            const labelList& interpolatedPoints,
            const List<FixedList<label, 3>>& interpolatedOldPoints,
            const List<FixedList<scalar, 3>>& interpolationWeights,
            const DynamicList<Type>& unmergedValues
        );
 
        triSurface stitchTriPoints
        (
            const bool checkDuplicates,
            const List<point>& triPoints,
            labelList& triPointReverseMap,  // unmerged to merged point
            labelList& triMap               // merged to unmerged triangle
        ) const;
 
        //- Trim triangle to planes
        static void trimToPlanes
        (
            const PtrList<plane>& planes,
            const triPointRef& tri,
            DynamicList<point>& newTriPoints
        );
 
        //- Trim all triangles to box
        static void trimToBox
        (
            const treeBoundBox& bb,
            DynamicList<point>& triPoints,
            DynamicList<label>& triMeshCells
        );
 
        //- Trim all triangles to box. Determine interpolation
        //  for existing and new points
        static void trimToBox
        (
            const treeBoundBox& bb,
            DynamicList<point>& triPoints,
            DynamicList<label>& triMap,
            labelList& triPointMap,
            labelList& interpolatedPoints,
            List<FixedList<label, 3>>& interpolatedOldPoints,
            List<FixedList<scalar, 3>>& interpolationWeights
        );
 
        static triSurface subsetMesh
        (
            const triSurface&,
            const labelList& newToOldFaces,
            labelList& oldToNewPoints,
            labelList& newToOldPoints
        );
 
public:
 
    //- Declare friendship to share some functionality
    friend class isoSurfaceCell;
    friend class isoSurfaceTopo;
 
    //- Filtering type
    using isoSurfaceBase::filterType;
 
 
    //- Runtime type information
    TypeName("isoSurface");
 
 
    // Constructors
 
        //- Construct from cell values and point values.
        //  Uses boundaryField for boundary values.
        //  Holds reference to cellIsoVals and pointIsoVals.
        //
        //  \param bounds optional bounding box for trimming
        //  \param mergeTol fraction of mesh bounding box for merging points
        isoSurface
        (
            const volScalarField& cellValues,
            const scalarField& pointValues,
            const scalar iso,
            const isoSurfaceBase::filterType filter,
            const boundBox& bounds = boundBox::invertedBox,
            const scalar mergeTol = 1e-6
        );
 
 
        //- Construct from cell values and point values.
        //  Uses boundaryField for boundary values.
        //  Holds reference to cellIsoVals and pointIsoVals.
        //
        //  \param bounds optional bounding box for trimming
        //  \param mergeTol fraction of mesh bounding box for merging points
        isoSurface
        (
            const volScalarField& cellValues,
            const scalarField& pointValues,
            const scalar iso,
            const bool regularise,
            const boundBox& bounds = boundBox::invertedBox,
            const scalar mergeTol = 1e-6
        );
 
 
    // Member Functions
 
        //- Interpolates cCoords, pCoords.
        //  Uses the references to the original fields used to create the
        //  iso surface.
        template<class Type>
        tmp<Field<Type>> interpolate
        (
            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
            const Field<Type>& pCoords
        ) const;
 
};
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
} // End namespace Foam
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
#ifdef NoRepository
    #include "isoSurfaceTemplates.C"
#endif
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
#endif
 
// ************************************************************************* //