TriangleMesh.h

Go to the documentation of this file.

#ifndef included_AMP_TriangleMesh
#define included_AMP_TriangleMesh
 
#include "AMP/geometry/Geometry.h"
#include "AMP/mesh/Mesh.h"
#include "AMP/mesh/MeshID.h"
#include "AMP/mesh/MeshIterator.h"
 
#include <array>
#include <map>
#include <memory>
#include <vector>
 
 
namespace AMP::Mesh {
 
 
template<uint8_t NG>
class TriangleMeshIterator;
template<uint8_t NG>
class TriangleMeshElement;
 
 
// Class to store vertex data
template<class TYPE, size_t N>
class StoreTriData final
{
public:
    StoreTriData() = default;
    StoreTriData( std::vector<std::array<TYPE, N>> x,
                  std::vector<int> offset,
                  int rank,
                  GeomType type );
    inline int size() const { return d_x.size(); }
    inline int start() const { return d_start; }
    inline int end() const { return d_end; }
    inline int start( int r ) const { return d_offset[r]; }
    inline int end( int r ) const { return d_offset[r + 1]; }
    int rank( int i ) const;
    inline auto data() { return d_x.data(); }
    inline const auto data() const { return d_x.data(); }
    inline auto &offset() const { return d_offset; }
    inline int index( const ElementID &id ) const
    {
        AMP_DEBUG_ASSERT( id.type() == d_type );
        int rank = id.owner_rank();
        return d_offset[rank] + id.local_id();
    }
    ElementID getID( int i ) const;
    int find( const std::array<TYPE, N> &x ) const;
    inline auto &operator[]( int i ) { return d_x[i]; }
    inline auto &operator[]( int i ) const { return d_x[i]; }
    inline auto &operator[]( const ElementID &id ) { return d_x[index( id )]; }
    inline auto &operator[]( const ElementID &id ) const { return d_x[index( id )]; }
    inline auto &operator()() { return d_x; }
    inline const auto &operator()() const { return d_x; }
 
private:
    GeomType d_type = GeomType::Vertex;
    int d_start     = 0;
    int d_end       = 0;
    int d_rank      = 0;
    std::vector<int> d_offset;
    std::vector<std::array<TYPE, N>> d_x;
};
 
 
// Class to store parent data
template<class TYPE>
class StoreCompressedList
{
public:
    inline StoreCompressedList() {}
    inline StoreCompressedList( size_t N )
        : StoreCompressedList( std::vector<std::vector<TYPE>>( N ) )
    {
    }
    inline explicit StoreCompressedList( const std::vector<std::vector<TYPE>> &data )
    {
        size_t Nt = 0;
        for ( size_t i = 0; i < data.size(); i++ )
            Nt += data[i].size();
        d_size.resize( data.size() );
        d_offset.resize( data.size() );
        d_data.resize( Nt );
        for ( size_t i = 0, k = 0; i < data.size(); i++ ) {
            d_size[i]   = data[i].size();
            d_offset[i] = k;
            for ( size_t j = 0; j < d_size[i]; j++, k++ )
                d_data[k] = data[i][j];
        }
    }
    inline const ElementID *begin( size_t i ) const
    {
        if ( i >= d_size.size() )
            return nullptr;
        return &d_data[d_offset[i]];
    }
    inline const ElementID *end( size_t i ) const
    {
        if ( i >= d_size.size() )
            return nullptr;
        return &d_data[d_offset[i]] + d_size[i];
    }
 
private:
    std::vector<size_t> d_size;
    std::vector<size_t> d_offset;
    std::vector<TYPE> d_data;
};
 
 
template<uint8_t NG>
class TriangleMesh : public AMP::Mesh::Mesh
{
public: // Convenience typedefs
    static_assert( NG <= 3, "Not programmed for higher dimensions yet" );
    typedef std::array<double, 3> Point;
    typedef std::array<int, 2> Edge;
    typedef std::array<int, 3> Face;
    typedef std::array<int, NG + 1> TRI;
    using IteratorSet = std::array<TriangleMeshIterator<NG>, NG + 1>;
    typedef StoreCompressedList<ElementID> ElementList;
 
 
public:
    explicit TriangleMesh( int NP,
                           std::vector<Point> vertices,
                           std::vector<TRI> triangles,
                           std::vector<TRI> tri_nab,
                           const AMP_MPI &comm,
                           std::shared_ptr<Geometry::Geometry> geom = {},
                           std::vector<int> block                   = {},
                           int max_gcw                              = 2 );
 
 
    std::string meshClass() const override;
 
 
    std::unique_ptr<Mesh> clone() const override final;
 
 
    bool operator==( const Mesh &mesh ) const override;
 
 
    // Copy/move constructors
    TriangleMesh( const TriangleMesh & );
    TriangleMesh( TriangleMesh && ) = default;
 
    TriangleMesh &operator=( const TriangleMesh & ) = delete;
    TriangleMesh &operator=( TriangleMesh && ) = default;
 
 
    virtual ~TriangleMesh();
 
 
    /* Return the number of local element of the given type
     * \param type   Geometric type
     */
    size_t numLocalElements( const GeomType type ) const override final;
 
 
    /* Return the global number of elements of the given type
     * Note: depending on the mesh this routine may require global communication across the mesh.
     * \param type   Geometric type
     */
    size_t numGlobalElements( const GeomType type ) const override final;
 
 
    /* Return the number of ghost elements of the given type on the current processor
     * \param type   Geometric type
     */
    size_t numGhostElements( const GeomType type, const int gcw ) const override final;
 
 
    MeshIterator getIterator( const GeomType type, const int gcw = 0 ) const override final;
 
 
    virtual MeshIterator getSurfaceIterator( const GeomType type,
                                             const int gcw = 0 ) const override final;
 
 
    std::vector<int> getBoundaryIDs() const override final;
 
 
    virtual MeshIterator getBoundaryIDIterator( const GeomType type,
                                                const int id,
                                                const int gcw = 0 ) const override final;
 
    std::vector<int> getBlockIDs() const override final;
 
 
    virtual MeshIterator
    getBlockIDIterator( const GeomType type, const int id, const int gcw = 0 ) const override final;
 
 
    MeshElementPtr getElement( const MeshElementID &id ) const override final;
 
 
    ElementListPtr getElementParents( const MeshElement &elem,
                                      const GeomType type ) const override final;
 
    Movable isMeshMovable() const override { return Movable::Deform; };
 
 
    uint64_t positionHash() const override;
 
 
    void displaceMesh( const std::vector<double> &x ) override;
 
 
    void displaceMesh( std::shared_ptr<const AMP::LinearAlgebra::Vector> x ) override;
 
 
    void writeRestart( int64_t fid ) const override;
 
 
protected:
    TriangleMesh() = default;
    explicit TriangleMesh( std::shared_ptr<const MeshParameters> );
    void initialize();
    void initializeIterators();
    void initializeBoundingBox();
    std::vector<IteratorSet> createBlockIterators( int block );
    void createSurfaceIterators();
 
 
public:
    // Create an iterator from a list
    TriangleMeshIterator<NG> createIterator( GeomType type, int gcw ) const;
 
    // Return the IDs of the elements composing the current element
    void getElementsIDs( const ElementID &id, const GeomType type, ElementID *IDs ) const;
    void getVertexCoord( const ElementID &id, std::array<double, 3> *x ) const;
 
    // Return the IDs of the neighboring elements
    void getNeighborIDs( const ElementID &id, std::vector<ElementID> &IDs ) const;
 
    // Return the IDs of the parent elements
    std::pair<const ElementID *, const ElementID *> getElementParents( const ElementID &id,
                                                                       const GeomType type ) const;
 
    // Return a new element (user must delete)
    MeshElement *getElement2( const MeshElementID &id ) const;
 
 
    // Check if the element is on the given boundry, block, etc
    bool isOnSurface( const ElementID &elemID ) const;
    bool isOnBoundary( const ElementID &elemID, int id ) const;
    bool isInBlock( const ElementID &elemID, int id ) const;
    static bool inIterator( const ElementID &id, const TriangleMeshIterator<NG> *it );
 
    template<uint8_t TYPE>
    std::array<int, TYPE + 1> getElem( const ElementID &id ) const;
    template<uint8_t TYPE>
    ElementID getID( const std::array<int, TYPE + 1> &id ) const;
 
 
private:
    void loadBalance( const std::vector<Point> &vertices,
                      const std::vector<TRI> &tri,
                      const std::vector<TRI> &tri_nab,
                      const std::vector<int> &block );
    void buildChildren();
    ElementList computeNodeParents( int parentType );
    ElementList getParents( int childType, int parentType );
 
 
private:
    std::array<size_t, 4> d_N_global = { 0 };  
    StoreTriData<double, 3> d_vertex;          
    StoreTriData<int, NG + 1> d_globalTri;     
    std::vector<TRI> d_globalNab;              
    std::vector<int> d_blockID;                
    std::vector<std::vector<int>> d_remoteTri; 
    StoreTriData<int, 2> d_childEdge;          
    StoreTriData<int, 3> d_childFace;          
    ElementList d_parents[NG + 1][NG + 1];     
    std::vector<int> d_block_ids;              
    std::vector<int> d_boundary_ids;           
    std::vector<bool> d_isSurface[NG + 1];     
    std::vector<IteratorSet> d_iterators;      
    std::vector<IteratorSet> d_surface_it;     
    std::vector<std::vector<IteratorSet>> d_boundary_it; 
    std::vector<std::vector<IteratorSet>> d_block_it;    
    uint64_t d_pos_hash = 0;                             
};
 
 
} // namespace AMP::Mesh
 
 
#endif