geom.h

Go to the documentation of this file.

#ifndef GEOM_FILE
#define GEOM_FILE
#include <array>
#include <functional>             // function; C++11
#include <iostream>
#include <memory>                  // shared_ptr
#include <string>
#include <vector>
 
class Mesh
{
public:
    explicit Mesh(int ndim, int nvert_e = 0, int ndof_e = 0, int nedge_e = 0);
 
    __attribute__((noinline))
    Mesh(Mesh const &) = default;
 
    Mesh &operator=(Mesh const &) = delete;
 
 
    virtual ~Mesh();
 
    explicit Mesh(std::string const &fname);
 
    void ReadVertexBasedMesh(std::string const &fname);
 
    int Nelems() const
    {
        return _nelem;
    }
 
    int NverticesElements() const
    {
        return _nvert_e;
    }
 
    int NdofsElement() const
    {
        return _ndof_e;
    }
 
    int Nnodes() const
    {
        return _nnode;
    }
 
    int Ndims() const
    {
        return _ndim;
    }
 
    void Resize_Connectivity(int nelem, int nvert_e)
    {
        SetNelem(nelem);               // number of elements
        SetNverticesElement(nvert_e);  // vertices per element
        _ia.resize(nelem * nvert_e);
    }
 
    const std::vector<int>  &GetConnectivity() const
    {
        return _ia;
    }
 
    std::vector<int>  &GetConnectivity()
    {
        return _ia;
    }
 
    void Resize_Coords(int nnodes, int ndim)
    {
        SetNnode(nnodes);       // number of nodes
        SetNdim(ndim);          // space dimension
        _xc.resize(nnodes * ndim);
    }
 
    const std::vector<double> &GetCoords() const
    {
        return _xc;
    }
 
    std::vector<double> &GetCoords()
    {
        return _xc;
    }
 
    void SetValues(std::vector<double> &v, const std::function<double(double, double)> &func) const;
    void SetBoundaryValues(std::vector<double> &v, const std::function<double(double, double)> &func) const;
    void SetDirchletValues(std::vector<double> &v, const std::function<double(double, double)> &func) const;
 
    void Debug() const;
 
    void DebugEdgeBased() const;
 
    virtual std::vector<int> Index_DirichletNodes() const;
    virtual std::vector<int> Index_BoundaryNodes() const;
 
    void Write_ascii_matlab(std::string const &fname, std::vector<double> const &v) const;
 
    void Export_scicomp(std::string const &basename) const;
 
    virtual void Visualize(std::vector<double> const &v) const;
 
    int Nedges() const
    {
        return _nedge;
    }
 
    int NedgesElements() const
    {
        return _nedge_e;
    }
 
    const std::vector<int>  &GetEdgeConnectivity() const
    {
        return _ea;
    }
 
    std::vector<int>  &GetEdgeConnectivity()
    {
        return _ea;
    }
 
    const std::vector<int>  &GetEdges() const
    {
        return _edges;
    }
 
    std::vector<int>  &GetEdges()
    {
        return _edges;
    }
 
    std::vector<std::vector<int>> Node2NodeGraph() const
    {
        //auto v1=Node2NodeGraph_1();
        //auto v2=Node2NodeGraph_2();
        //if ( equal(v1.cbegin(),v1.cend(),v2.begin()) )
        //{
        //std::cout << "\nidentical Versions\n";
        //}
        //else
        //{
        //std::cout << "\nE R R O R   in Versions\n";
        //}
 
        //return Node2NodeGraph_1();
        return Node2NodeGraph_2();        // 2 times faster than version 1
    }
 
    virtual std::vector<int> const &GetFathersOfVertices() const
    {
        return _dummy;
    }
 
    void Del_EdgeConnectivity();
 
protected:
    //public:
    void SetNelem(int nelem)
    {
        _nelem = nelem;
    }
 
    void SetNverticesElement(int nvert)
    {
        _nvert_e = nvert;
    }
 
    void SetNdofsElement(int ndof)
    {
        _ndof_e = ndof;
    }
 
    void SetNnode(int nnode)
    {
        _nnode = nnode;
    }
 
    void SetNdim(int ndim)
    {
        _ndim = ndim;
    }
 
    void SetNedge(int nedge)
    {
        _nedge = nedge;
    }
 
    void ReadVectexBasedMesh(std::string const &fname);
 
    void DeriveEdgeFromVertexBased()
    {
        //DeriveEdgeFromVertexBased_slow();
        //DeriveEdgeFromVertexBased_fast();
        DeriveEdgeFromVertexBased_fast_2();
    }
    void DeriveEdgeFromVertexBased_slow();
    void DeriveEdgeFromVertexBased_fast();
    void DeriveEdgeFromVertexBased_fast_2();
 
 
 
    void DeriveVertexFromEdgeBased();
 
    int Nnbedges() const
    {
        return static_cast<int>(_bedges.size());
    }
 
    virtual bool Check_array_dimensions() const;
 
    void PermuteVertices_EdgeBased(std::vector<int> const &old2new);
 
private:
    std::vector<std::vector<int>> Node2NodeGraph_1() const;  // is correct
 
    std::vector<std::vector<int>> Node2NodeGraph_2() const;  // is correct
 
    //private:
protected:
    int _nelem;         
    int _nvert_e;       
    int _ndof_e;        
    int _nnode;         
    int _ndim;          
    std::vector<int> _ia;    
    std::vector<double> _xc; 
 
protected:
    // B.C.
    std::vector<int> _bedges;     
// 2020-01-08
    std::vector<int> _sdedges;    
 
    //private:
protected:
    // edge based connectivity
    int _nedge;              
    int _nedge_e;            
    std::vector<int> _edges; 
    std::vector<int> _ea;    
    // B.C.
    std::vector<int> _ebedges; 
 
private:
    const std::vector<int> _dummy; 
 
};
 
 
// *********************************************************************
 
class RefinedMesh: public Mesh
{
public:
    //explicit RefinedMesh(Mesh const &cmesh, std::vector<bool> const &ibref = std::vector<bool>(0));
    RefinedMesh(Mesh const &cmesh, std::vector<bool> const &ibref);
    //RefinedMesh(Mesh const &cmesh, std::vector<bool> const &ibref);
 
    explicit RefinedMesh(Mesh const &cmesh)
        : RefinedMesh(cmesh, std::vector<bool>(0))
    {}
 
 
    RefinedMesh(RefinedMesh const &) = delete;
    //RefinedMesh(RefinedMesh const&&) = delete;
 
    RefinedMesh &operator=(RefinedMesh const &) = delete;
    //RefinedMesh& operator=(RefinedMesh const&&) = delete;
 
    virtual ~RefinedMesh() override;
 
    Mesh RefineElements(std::vector<bool> const &ibref);
 
    void RefineAllElements(int nref = 1);
 
    std::vector<int> const &GetFathersOfVertices() const override
    {
        return _vfathers;
    }
 
protected:
    bool Check_array_dimensions() const override;
 
    void PermuteVertices_EdgeBased(std::vector<int> const &old2new);
 
 
private:
    //Mesh const              & _cmesh; //!< coarse mesh
    std::vector<bool> const   _ibref; 
    int                       _nref;  
    std::vector<int>          _vfathers; 
 
};
 
// *********************************************************************
 
class gMesh_Hierarchy
{
public:
    gMesh_Hierarchy(Mesh const &cmesh, int nlevel);
 
    size_t size() const
    {
        return _gmesh.size();
    }
 
    Mesh const &operator[](int lev) const
    {
        return *_gmesh.at(lev);
    }
 
    Mesh const &finest() const
    {
        return *_gmesh.back();
    }
 
    Mesh const &coarsest() const
    {
        return *_gmesh.front();
    }
 
private:
    std::vector<std::shared_ptr<Mesh>> _gmesh; 
 
};
 
 
 
// *********************************************************************
class Mesh_2d_3_square: public Mesh
{
public:
    Mesh_2d_3_square(int nx, int ny, int myid = 0, int procx = 1, int procy = 1);
 
    ~Mesh_2d_3_square() override;
 
    void SetU(std::vector<double> &u) const;
 
    void SetF(std::vector<double> &f) const;
 
    std::vector<int> Index_DirichletNodes() const override;
    std::vector<int> Index_BoundaryNodes() const override;
 
    void SaveVectorP(std::string const &name, std::vector<double> const &u) const;
 
    // here will still need to implement in the class
    //  GetBound(), AddBound()
    //  or better a generalized way with indices and their appropriate ranks for MPI communication
 
private:
    void GetCoordsInRectangle(int nx, int ny, double xl, double xr, double yb, double yt,
                              double xc[]);
    void GetConnectivityInRectangle(int nx, int ny, int ia[]);
 
private:
    int _myid;          
    int _procx;         
    int _procy;         
    std::array<int, 4> _neigh; 
    int _color;         
 
    double _xl;         
    double _xr;         
    double _yb;         
    double _yt;         
    int    _nx;         
    int    _ny;         
};
 
// *********************************************************************
 
 
 
 
#endif