triangleBoundary3D.h

Go to the documentation of this file.

/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011 FlowKit Sarl
 * Avenue de Chailly 23
 * 1012 Lausanne, Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library 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 Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef INNER_FLOW_BOUNDARY_3D_H
#define INNER_FLOW_BOUNDARY_3D_H

#include "core/globalDefs.h"
#include "offLattice/boundaryShapes3D.h"
#include "offLattice/triangleSet.h"
#include "offLattice/triangularSurfaceMesh.h"
#include "offLattice/offLatticeBoundaryProfiles3D.h"
#include "particles/multiParticleField3D.h"
#include "multiBlock/redistribution3D.h"
#include "multiBlock/multiDataField3D.h"
#include "atomicBlock/atomicContainerBlock3D.h"
#include "multiBlock/multiContainerBlock3D.h"
#include <stack>

namespace plb {

template<typename T>
class DEFscaledMesh {
public:
    DEFscaledMesh(TriangleSet<T> const& triangleSet_);
    DEFscaledMesh (
        TriangleSet<T> const& triangleSet_,
        plint resolution_, plint referenceDirection_,
        plint margin_, plint extraLayer );
    DEFscaledMesh (
        TriangleSet<T> const& triangleSet_,
        plint resolution_, plint referenceDirection_,
        plint margin_, Dot3D location );
    DEFscaledMesh(DEFscaledMesh<T> const& rhs);
    ~DEFscaledMesh();
    DEFscaledMesh<T>& operator=(DEFscaledMesh<T> const& rhs);
    void swap(DEFscaledMesh<T>& rhs);
public: // Mesh usage interface.
    TriangularSurfaceMesh<T>& getMesh();
    TriangularSurfaceMesh<T> const& getMesh() const;
    plint getMargin() const;
    std::vector<Array<T,3> > const& getVertexList() const {
        return vertexList;
    }
    std::vector<plint> const& getEmanatingEdgeList() const {
        return emanatingEdgeList;
    }
    std::vector<typename TriangularSurfaceMesh<T>::Edge> const&
        getEdgeList() const {
            return edgeList;
        }
    Array<T,3> getPhysicalLocation() const {
        return physicalLocation;
    }
    T getDx() const {
        return dx;
    }
    void setPhysicalLocation(Array<T,3> physicalLocation_) {
        physicalLocation = physicalLocation_;
    }
    void setDx(T dx_) {
        dx = dx_;
    }
private:
    void initialize (
        TriangleSet<T> const& triangleSet_, plint resolution_,
        plint referenceDirection_, Dot3D location );
private:
    std::vector<Array<T,3> > vertexList;
    std::vector<plint> emanatingEdgeList;
    std::vector<typename TriangularSurfaceMesh<T>::Edge> edgeList;
    TriangularSurfaceMesh<T>* mesh;
    plint margin;
    Array<T,3> physicalLocation;
    T dx;
};

template<typename T>
struct VertexProperty3D {
    virtual ~VertexProperty3D() { }
    // Wall portions which should not follow the laws of elasticity are rigid.
    virtual bool isRigid() const =0;
    // Inlet/outlet nodes are for example not part of the wall.
    virtual bool isWall() const =0;
    virtual VertexProperty3D<T>* clone() const =0;
};

template<typename T>
struct RigidWallProperty3D : public VertexProperty3D<T> {
    virtual bool isRigid() const { return true; }
    virtual bool isWall() const { return true; }
    virtual RigidWallProperty3D<T>* clone() const {
        return new RigidWallProperty3D<T>(*this);
    }
};

template<typename T>
struct InletOutletProperty3D : public VertexProperty3D<T> {
    virtual bool isRigid() const { return true; }
    virtual bool isWall() const { return false; }
    virtual InletOutletProperty3D<T>* clone() const {
        return new InletOutletProperty3D<T>(*this);
    }
};

template<typename T>
bool isRigid(VertexProperty3D<T> const* property) {
    if (property) {
        return property->isRigid();
    }
    else {
        return false;
    }
}

template<typename T>
bool isWall(VertexProperty3D<T> const* property) {
    if (property) {
        return property->isWall();
    }
    else {
        return true;
    }
}

template<typename T> class TriangleBoundary3D;

template<typename T, class SurfaceData>
class BoundaryProfiles3D {
public:    
    BoundaryProfiles3D();
    ~BoundaryProfiles3D();
    BoundaryProfiles3D(BoundaryProfiles3D<T,SurfaceData> const& rhs);
    BoundaryProfiles3D<T,SurfaceData>& operator=(BoundaryProfiles3D<T,SurfaceData> const& rhs);
    void swap(BoundaryProfiles3D<T,SurfaceData>& rhs);
    void defineProfile(plint tag, BoundaryProfile3D<T,SurfaceData>* profile);
    void resetProfiles(std::map<plint,BoundaryProfile3D<T,SurfaceData>*> profiles_);
    void defineInletOutletTags(TriangleBoundary3D<T> const& boundary, plint sortDirection);
    void setWallProfile(BoundaryProfile3D<T,SurfaceData>* wallProfile_);
    void setInletOutlet( std::vector<BoundaryProfile3D<T,SurfaceData>*> inletOutlets );
    void setInletOutlet( BoundaryProfile3D<T,SurfaceData>* profile1, BoundaryProfile3D<T,SurfaceData>* profile2 );
    void setInletOutlet( BoundaryProfile3D<T,SurfaceData>* profile1, BoundaryProfile3D<T,SurfaceData>* profile2,
                         BoundaryProfile3D<T,SurfaceData>* profile3 );
    void setInletOutlet( BoundaryProfile3D<T,SurfaceData>* profile1, BoundaryProfile3D<T,SurfaceData>* profile2,
                         BoundaryProfile3D<T,SurfaceData>* profile3, BoundaryProfile3D<T,SurfaceData>* profile4 );
    void adjustInletOutlet(TriangleBoundary3D<T> const& boundary, plint sortDirection);
    BoundaryProfile3D<T,SurfaceData> const& getProfile (
            TriangleBoundary3D<T> const& boundary, plint iTriangle ) const;
private:
    void replaceProfile(plint id, BoundaryProfile3D<T,SurfaceData>* newProfile);
    void clearProfiles();
private:
    BoundaryProfile3D<T,SurfaceData>* wallProfile;
    std::map<plint,BoundaryProfile3D<T,SurfaceData>*> profiles;
    std::vector<plint> inletOutletIds;
    std::vector<Array<T,3> > lidNormal;
    std::vector<Array<T,3> > lidCenter;
    std::vector<T> lidRadius;
};

template<typename T>
class TriangleBoundary3D {
public:
    TriangleBoundary3D(DEFscaledMesh<T> const& defMesh, bool automaticCloseHoles=true);
    ~TriangleBoundary3D();
    TriangleBoundary3D(TriangleBoundary3D<T> const& rhs);
    TriangleBoundary3D<T>& operator=(TriangleBoundary3D<T> const& rhs);
    void swap(TriangleBoundary3D<T>& rhs);
public: // Mesh usage interface.
    TriangleBoundary3D<T> const& select (
            plint whichTopology, plint whichVertices ) const;
    TriangleBoundary3D<T> const& pushSelect (
            plint whichTopology, plint whichVertices ) const;
    TriangleBoundary3D<T> const& popSelect() const;
    void getSelection(plint& whichTopology, plint& whichVertices) const;
    TriangularSurfaceMesh<T>& getMesh();
    TriangularSurfaceMesh<T> const& getMesh() const;
    VertexProperty3D<T> const* getVertexProperty(plint iVertex) const;
    bool intersectSegment (
            plint iTriangle, AtomicBlock3D* boundaryArg,
            Array<T,3> const& fromPoint, Array<T,3> const& direction,
            Array<T,3>& locatedPoint, T& distance, Array<T,3>& wallNormal ) const;
    Array<T,3> computeContinuousNormal (
            Array<T,3> const& p, plint iTriangle, bool isAreaWeighted = false ) const;
    void cloneVertexSet(plint whichVertexSet);
    Array<T,3> getPhysicalLocation() const {
        return physicalLocation;
    }
    T getDx() const {
        return dx;
    }
public: // Mesh preparation interface.
    std::vector<typename TriangularSurfaceMesh<T>::Lid> const& getInletOutlet() const;

    std::vector<typename TriangularSurfaceMesh<T>::Lid> getInletOutlet(plint sortDirection) const;
    template<typename DomainFunctional> plint tagDomain(DomainFunctional functional);
    template<typename DomainFunctional> plint tagDomain(DomainFunctional functional, Array<T,3> normal, T angleTolerance, plint previousTag=-1);
    std::vector<plint> getInletOutletIds(plint sortDirection) const;
    void getLidProperties (
        plint sortDirection, std::vector<Array<T,3> >& normal,
        std::vector<Array<T,3> >& center, std::vector<T>& radius ) const;

    template<typename DomainFunctional>
    plint setVertexProperty (
            VertexProperty3D<T> const& property, DomainFunctional functional );
    plint getMargin() const;
    plint getTag(plint iTriangle) const;
    std::vector<plint> const& getTriangleTags() const { return triangleTagList; }
    std::vector<plint> const& getVertexTags() const { return vertexTagList; }
private:
    plint currentMesh() const;
    void defineMeshes();
    void closeHoles();
    void assignLidVertexProperty();
private:
    void tagInletOutlet (
        std::vector<typename TriangularSurfaceMesh<T>::Lid> const& newLids );
    std::vector< std::vector<Array<T,3> > > vertexLists;
    std::vector<plint> emanatingEdgeLists[2];
    std::vector<typename TriangularSurfaceMesh<T>::Edge> edgeLists[2];
    std::vector<TriangularSurfaceMesh<T> > meshes;
    std::vector<plint> triangleTagList;
    plint currentTagNum;
    std::vector<plint> vertexTagList;
    std::vector<VertexProperty3D<T>*> vertexProperties;
    std::vector<typename TriangularSurfaceMesh<T>::Lid> lids;
    plint margin;
    Array<T,3> physicalLocation;
    T dx;
    mutable std::stack<plint> topology;
    mutable std::stack<plint> vertexSet;
};

template< typename T, class SurfaceData >
class TriangleFlowShape3D : public BoundaryShape3D<T,SurfaceData> {
public:
    TriangleFlowShape3D (
            TriangleBoundary3D<T> const& boundary_,
            BoundaryProfiles3D<T,SurfaceData> const& profiles_ );
    virtual bool isInside(Dot3D const& location) const;
    virtual bool pointOnSurface (
            Array<T,3> const& fromPoint, Array<T,3> const& direction,
            Array<T,3>& locatedPoint, T& distance,
            Array<T,3>& wallNormal, SurfaceData& surfaceData,
            OffBoundary::Type& bdType, plint& id ) const;
    virtual Array<T,3> computeContinuousNormal (
            Array<T,3> const& p, plint id, bool isAreaWeighted = false ) const;
    virtual bool intersectsSurface (
            Array<T,3> const& p1, Array<T,3> const& p2, plint& id ) const;
    virtual plint getTag(plint id) const;
    virtual bool distanceToSurface( Array<T,3> const& point,
                                    T& distance, bool& isBehind ) const;
    virtual TriangleFlowShape3D<T,SurfaceData>* clone() const;

    virtual TriangleFlowShape3D<T,SurfaceData>*
                clone(std::vector<AtomicBlock3D*> args) const;
private:
    TriangleBoundary3D<T> const& boundary;
    BoundaryProfiles3D<T,SurfaceData> const& profiles;
    ScalarField3D<int>* voxelFlags;
    AtomicContainerBlock3D* hashContainer;
    AtomicBlock3D* boundaryArg;
};

template<typename T>
class VoxelizedDomain3D {
public:
    VoxelizedDomain3D(TriangleBoundary3D<T> const& boundary_,
                      int flowType_, plint extraLayer_, plint borderWidth_,
                      plint envelopeWidth_, plint blockSize_,
                      plint gridLevel_=0, bool dynamicMesh_ = false);
    VoxelizedDomain3D(TriangleBoundary3D<T> const& boundary_,
                      int flowType_, Box3D const& boundingBox, plint borderWidth_,
                      plint envelopeWidth_, plint blockSize_,
                      plint gridLevel_=0, bool dynamicMesh_ = false);
    VoxelizedDomain3D(VoxelizedDomain3D<T> const& rhs);
    ~VoxelizedDomain3D();
    MultiScalarField3D<int>& getVoxelMatrix();
    MultiScalarField3D<int> const& getVoxelMatrix() const;
    MultiContainerBlock3D& getTriangleHash();
    MultiBlockManagement3D const& getMultiBlockManagement() const;
    template<class ParticleFieldT>
    void adjustVoxelization(MultiParticleField3D<ParticleFieldT>& particles, bool dynamicMesh);
    void reparallelize(MultiBlockRedistribute3D const& redistribute);
    TriangleBoundary3D<T> const& getBoundary() const { return boundary; }
    int getFlowType() const { return flowType; }
private:
    VoxelizedDomain3D<T>& operator=(VoxelizedDomain3D<T> const& rhs) { }
    void createSparseVoxelMatrix (
        MultiScalarField3D<int>& fullVoxelMatrix,
        plint blockSize_, plint envelopeWidth_ );
    void computeSparseVoxelMatrix (
            MultiScalarField3D<int>& fullVoxelMatrix,
            plint blockSize, plint envelopeWidth );
    void extendEnvelopeWidth (
            MultiScalarField3D<int>& fullVoxelMatrix, plint envelopeWidth );
    void createTriangleHash();
    template<class ParticleFieldT>
    void reCreateTriangleHash(MultiParticleField3D<ParticleFieldT>& particles);
    void computeOuterMask();
private:
    int flowType;
    plint borderWidth;
    TriangleBoundary3D<T> const& boundary;
    MultiScalarField3D<int>* voxelMatrix;
    MultiContainerBlock3D* triangleHash;
};

template<typename T>
void addLayer(MultiScalarField3D<T>& matrix, Box3D const& domain,
              T previousLayer);

template<typename T>
class AddLayerFunctional3D : public BoxProcessingFunctional3D_S<T> {
public:
    AddLayerFunctional3D(T previousLayer_);
    virtual void process(Box3D domain, ScalarField3D<T>& voxels);
    virtual AddLayerFunctional3D<T>* clone() const;
    virtual void getTypeOfModification(std::vector<modif::ModifT>& modified) const;
    virtual BlockDomain::DomainT appliesTo() const;
private:
    T previousLayer;
};

}  // namespace plb

#endif  // INNER_FLOW_BOUNDARY_3D_H