variabledofspace.hpp

Go to the documentation of this file.

#ifndef VARIABLEDOFSPACE_HPP
#define VARIABLEDOFSPACE_HPP
 
#include<mesh.hpp>
 
namespace HArDCore3D
{
 
  class VariableDOFSpace {
  public:
    VariableDOFSpace(
             const Mesh & mesh,
             const Eigen::VectorXi n_local_vertex_dofs,
             const Eigen::VectorXi n_local_edge_dofs,
             const Eigen::VectorXi n_local_face_dofs,
             const Eigen::VectorXi n_local_cell_dofs
             );
 
    VariableDOFSpace(
             const Mesh & mesh,
             size_t n_local_vertex_dofs,
             const Eigen::VectorXi n_local_edge_dofs,
             const Eigen::VectorXi n_local_face_dofs,
             const Eigen::VectorXi n_local_cell_dofs
             );
 
    VariableDOFSpace(
             const Mesh & mesh,
             size_t n_local_vertex_dofs,
             size_t n_local_edge_dofs,
             const Eigen::VectorXi n_local_face_dofs,
             const Eigen::VectorXi n_local_cell_dofs
             );
 
    VariableDOFSpace(
             const Mesh & mesh,
             size_t n_local_vertex_dofs,
             size_t n_local_edge_dofs,
             size_t n_local_face_dofs,
             const Eigen::VectorXi n_local_cell_dofs
             );
 
    VariableDOFSpace(
             const Mesh & mesh,
             size_t n_local_vertex_dofs,
             size_t n_local_edge_dofs,
             size_t n_local_face_dofs,
             size_t n_local_cell_dofs
             );
 
 
    //------------------------------------------------------------------------------
    // Accessors
    //------------------------------------------------------------------------------
    
    const Mesh & mesh() const
    {
      return m_mesh;
    }
 
    inline size_t numLocalDofsVertex(const size_t iV) const
    {
      assert(iV<m_mesh.n_vertices());
      return m_n_local_vertex_dofs[iV];
    }
    
    inline size_t numLocalDofsVertex(const Vertex & V) const
    {
      return numLocalDofsVertex(V.global_index());
    }
 
    inline size_t numLocalDofsEdge(const size_t iE) const
    {
      assert(iE<m_mesh.n_edges());
      return m_n_local_edge_dofs[iE];
    }
    
    inline size_t numLocalDofsEdge(const Edge & E) const
    {
      return numLocalDofsEdge(E.global_index());
    }
 
    inline size_t numLocalDofsFace(const size_t iF) const
    {
      assert(iF<m_mesh.n_faces());
      return m_n_local_face_dofs[iF];
    }
    
    inline size_t numLocalDofsFace(const Face & F) const
    {
      return numLocalDofsFace(F.global_index());
    }
 
    inline size_t numLocalDofsCell(const size_t iT) const
    {
      assert(iT<m_mesh.n_cells());
      return m_n_local_cell_dofs[iT];
    }
 
    inline size_t numLocalDofsCell(const Cell & T) const
    {
      return numLocalDofsCell(T.global_index());
    }
 
    //------------------------------------------------------------------------------
    // Dimensions
    //------------------------------------------------------------------------------
    
   
    inline size_t nDOFs_vertices() const
    {
      return m_n_local_vertex_dofs.sum();
    }
   
    inline size_t nDOFs_edges() const
    {
      return m_n_local_edge_dofs.sum();
    }
 
    inline size_t nDOFs_faces() const
    {
      return m_n_local_face_dofs.sum();
    }
 
    inline size_t nDOFs_cells() const
    {
      return m_n_local_cell_dofs.sum();
    }
 
    inline size_t dimension() const
    {
      return nDOFs_vertices() + nDOFs_edges() + nDOFs_faces() + nDOFs_cells();
    }
 
    //--- Vertices -----//
    inline size_t dimensionVertex(const Vertex & V) const
    {
      return m_n_local_vertex_dofs[V.global_index()];
    }
 
    inline size_t dimensionVertex(size_t iV) const
    {
      assert(iV<m_mesh.n_vertices());
      return dimensionVertex(*m_mesh.vertex(iV));
    }
 
    //--- Edges ----//
    inline size_t dimensionEdge(const Edge & E) const
    {
      return m_n_local_vertex_dofs[E.vertex(0)->global_index()] + m_n_local_vertex_dofs[E.vertex(1)->global_index()]
        + m_n_local_edge_dofs[E.global_index()];
    }
 
    inline size_t dimensionEdge(size_t iE) const
    {
      assert(iE<m_mesh.n_edges());
      return dimensionEdge(*m_mesh.edge(iE));
    }
 
    //--- Faces -----//
    inline size_t dimensionFace(const Face & F) const
    {
      size_t nb_dofs = 0;
      for (Vertex * V : F.get_vertices()){
        nb_dofs += m_n_local_vertex_dofs[V->global_index()];
      }
      for (Edge * E : F.get_edges()){
        nb_dofs += m_n_local_edge_dofs[E->global_index()];
      }
      return nb_dofs + m_n_local_face_dofs[F.global_index()];
    }
    
    inline size_t dimensionFace(size_t iF) const
    {
      return dimensionFace(*m_mesh.face(iF));
    }
 
    //--- Cell -----//
    inline size_t dimensionCell(const Cell & T) const
    {
      size_t nb_dofs = 0;
      for (Vertex * V : T.get_vertices()){
        nb_dofs += m_n_local_vertex_dofs[V->global_index()];
      }
      for (Edge * E : T.get_edges()){
        nb_dofs += m_n_local_edge_dofs[E->global_index()];
      }
      for (Face * F : T.get_faces()){
        nb_dofs += m_n_local_face_dofs[F->global_index()];
      }
      return nb_dofs + m_n_local_cell_dofs[T.global_index()];
    }
    
    inline size_t dimensionCell(size_t iT) const
    {
      return dimensionCell(*m_mesh.cell(iT));
    }
 
    //------------------------------------------------------------------------------
    // Local offsets
    //------------------------------------------------------------------------------
 
    inline size_t localOffset(const Edge & E, const Vertex & V) const
    {
      size_t nb_dofs = 0;
      for (int iV=0; iV<E.index_vertex(&V); iV++){
        nb_dofs += m_n_local_vertex_dofs[E.vertex(iV)->global_index()];
      }
      return nb_dofs;
    }
 
    inline size_t localOffset(const Edge & E) const
    {
      return m_n_local_vertex_dofs[E.vertex(0)->global_index()] + m_n_local_vertex_dofs[E.vertex(1)->global_index()];
    }
    
    inline size_t localOffset(const Face & F, const Vertex & V) const
    {
      size_t nb_dofs = 0;
      for (int iV=0; iV<F.index_vertex(&V); iV++){
        nb_dofs += m_n_local_vertex_dofs[F.vertex(iV)->global_index()];
      }
      return nb_dofs;
    }
 
    inline size_t localOffset(const Face & F, const Edge & E) const
    {
      size_t nb_dofs = 0;
      for (size_t iV=0; iV<F.n_vertices(); iV++){
        nb_dofs += m_n_local_vertex_dofs[F.vertex(iV)->global_index()];
      }
      for (int iE=0; iE<F.index_edge(&E); iE++){
        nb_dofs += m_n_local_edge_dofs[F.edge(iE)->global_index()];
      }
      return nb_dofs;
    }
 
    inline size_t localOffset(const Face & F) const
    {
      size_t nb_dofs = 0;
      for (size_t iV=0; iV<F.n_vertices(); iV++){
        nb_dofs += m_n_local_vertex_dofs[F.vertex(iV)->global_index()];
      }
      for (size_t iE=0; iE<F.n_edges(); iE++){
        nb_dofs += m_n_local_edge_dofs[F.edge(iE)->global_index()];
      }
      return nb_dofs;
    }
 
    inline size_t localOffset(const Cell & T, const Vertex & V) const
    {
      size_t nb_dofs = 0;
      for (int iV=0; iV<T.index_vertex(&V); iV++){
        nb_dofs += m_n_local_vertex_dofs[T.vertex(iV)->global_index()];
      }
      return nb_dofs;
    }
 
    inline size_t localOffset(const Cell & T, const Edge & E) const
    {
      size_t nb_dofs = 0;
      for (size_t iV=0; iV<T.n_vertices(); iV++){
        nb_dofs += m_n_local_vertex_dofs[T.vertex(iV)->global_index()];
      }
      for (int iE=0; iE<T.index_edge(&E); iE++){
        nb_dofs += m_n_local_edge_dofs[T.edge(iE)->global_index()];
      }
      return nb_dofs;
    }
 
    inline size_t localOffset(const Cell & T, const Face & F) const
    {
      size_t nb_dofs = 0;
      for (size_t iV=0; iV<T.n_vertices(); iV++){
        nb_dofs += m_n_local_vertex_dofs[T.vertex(iV)->global_index()];
      }
      for (size_t iE=0; iE<T.n_edges(); iE++){
        nb_dofs += m_n_local_edge_dofs[T.edge(iE)->global_index()];
      }
      for (int iF=0; iF<T.index_face(&F); iF++){
        nb_dofs += m_n_local_face_dofs[T.face(iF)->global_index()];
      }
      return nb_dofs;
    }
 
    inline size_t localOffset(const Cell & T) const
    {
      size_t nb_dofs = 0;
      for (size_t iV=0; iV<T.n_vertices(); iV++){
        nb_dofs += m_n_local_vertex_dofs[T.vertex(iV)->global_index()];
      }
      for (size_t iE=0; iE<T.n_edges(); iE++){
        nb_dofs += m_n_local_edge_dofs[T.edge(iE)->global_index()];
      }
      for (size_t iF=0; iF<T.n_faces(); iF++){
        nb_dofs += m_n_local_face_dofs[T.face(iF)->global_index()];
      }
      return nb_dofs;
    }
 
    //------------------------------------------------------------------------------
    // Global offsets
    //------------------------------------------------------------------------------
    
    inline size_t globalOffset(const Vertex & V) const
    {
      return ( m_n_local_vertex_dofs.head(V.global_index()) ).sum();
    }
 
    inline size_t globalOffset(const Edge & E) const
    {
      return m_n_local_vertex_dofs.sum()
          + ( m_n_local_edge_dofs.head(E.global_index()) ).sum();
    }
 
    inline size_t globalOffset(const Face & F) const
    {
      return m_n_local_vertex_dofs.sum()
          + m_n_local_edge_dofs.sum()
          + ( m_n_local_face_dofs.head(F.global_index()) ).sum();
    }
 
    inline size_t globalOffset(const Cell & T) const
    {
      return m_n_local_vertex_dofs.sum()
          + m_n_local_edge_dofs.sum()
          + m_n_local_face_dofs.sum()
          + ( m_n_local_cell_dofs.head(T.global_index()) ).sum();
    }
    
    //------------------------------------------------------------------------------
    // Restrictions
    //------------------------------------------------------------------------------
 
    Eigen::VectorXd restrictEdge(size_t iE, const Eigen::VectorXd & vh) const;
 
    Eigen::VectorXd restrictFace(size_t iF, const Eigen::VectorXd & vh) const;
 
    Eigen::VectorXd restrictCell(size_t iT, const Eigen::VectorXd & vh) const;
 
    inline Eigen::VectorXd restrict(const Edge & E, const Eigen::VectorXd vh) const
    {
      return restrictEdge(E.global_index(), vh);
    }
 
    inline Eigen::VectorXd restrict(const Face & F, const Eigen::VectorXd vh) const
    {
      return restrictFace(F.global_index(), vh);
    }
 
    inline Eigen::VectorXd restrict(const Cell & T, const Eigen::VectorXd vh) const
    {
      return restrictCell(T.global_index(), vh);
    }
    
    //------------------------------------------------------------------------------
    // Extensions
    //------------------------------------------------------------------------------
 
    Eigen::MatrixXd extendOperator(const Cell & T, const Face & F, const Eigen::MatrixXd & opF) const;
 
    Eigen::MatrixXd extendOperator(const Cell & T, const Edge & E, const Eigen::MatrixXd & opE) const;
 
    Eigen::MatrixXd extendOperator(const Face & F, const Edge & E, const Eigen::MatrixXd & opE) const;
 
    void extendOperator(const Cell & T, const Face & F, Eigen::Ref<Eigen::MatrixXd> opT, const Eigen::MatrixXd & opF) const;
 
    void extendOperator(const Cell & T, const Edge & E, Eigen::MatrixXd & opT, const Eigen::MatrixXd & opE) const;
 
    void extendOperator(const Face & F, const Edge & E, Eigen::MatrixXd & opF, const Eigen::MatrixXd & opE) const;
 
    void addInnerProductContribution(const Cell & T, const Face & F, Eigen::MatrixXd & prodT, const Eigen::MatrixXd & prodF) const;
    
    //------------------------------------------------------------------------------
    // Global DOF indices for an element T
    //------------------------------------------------------------------------------
 
    std::vector<size_t> globalDOFIndices(const Cell & T) const;
 
    std::vector<size_t> globalDOFIndices(const Face & F) const;  
    
    private:
    const Mesh & m_mesh;
    Eigen::VectorXi m_n_local_vertex_dofs;
    Eigen::VectorXi m_n_local_edge_dofs;
    Eigen::VectorXi m_n_local_face_dofs;
    Eigen::VectorXi m_n_local_cell_dofs;
 
      
  };
 
} // namespace HArDCore3D
 
#endif