xdivdiv.hpp

Go to the documentation of this file.

#ifndef XDIVDIV_HPP
#define XDIVDIV_HPP
 
#include <globaldofspace.hpp>
#include <platescore.hpp>
#include <integralweight.hpp>
 
namespace HArDCore2D
{
 
  class XDivDiv : public GlobalDOFSpace
  {
  public:
    typedef std::function<Eigen::Matrix2d(const Eigen::Vector2d &)> FunctionType;
    typedef std::function<double(const Eigen::Vector2d &, const Edge &)> EdgeFunctionType;
    typedef std::function<Eigen::Matrix2d(const Eigen::Matrix2d &)> ConstitutiveLawType;
 
    struct LocalOperators
    {
      LocalOperators(
             const Eigen::MatrixXd & _divdiv,     
                     const Eigen::MatrixXd & _divdiv_rhs, 
             const Eigen::MatrixXd & _potential   
             )
        : divdiv(_divdiv),
              divdiv_rhs(_divdiv_rhs),
          potential(_potential)
      {
        // Do nothing
      }
 
      Eigen::MatrixXd divdiv;
      Eigen::MatrixXd divdiv_rhs;
      Eigen::MatrixXd potential;
    };
    
    struct SymmetricMatrixBasisVertex {
      typedef MatrixRd FunctionValue;
 
      SymmetricMatrixBasisVertex();      
 
      inline FunctionValue function(size_t i) const
      {
        assert(i >= 0 && i <= 2);
        return m_basis[i];
      }
 
      inline size_t size() const {
        return 3;
      }
      
    private:
      std::array<MatrixRd, 3> m_basis;
    };
 
    XDivDiv(const PlatesCore & plates_core, bool use_threads = true, std::ostream & output = std::cout);
 
    const Mesh & mesh() const
    {
      return m_plates_core.mesh();
    }
    
    const size_t & degree() const
    {
      return m_plates_core.degree();
    }
 
    Eigen::VectorXd interpolate(
                                const FunctionType & tau,      
                                const EdgeFunctionType & Dtau, 
                                const int deg_quad = -1        
                                ) const;
 
    Eigen::MatrixXd computeL2Product(
                                     const size_t iT, 
                                     const ConstitutiveLawType & A = [](const Eigen::Matrix2d tau)->Eigen::Matrix2d { return tau; }, 
                                     const double & penalty_factor = 1. 
                                     ) const;
 
 
    inline const LocalOperators & cellOperators(size_t iT) const
    {
      return *m_cell_operators[iT];
    }
 
    inline const Eigen::MatrixXd & edgePotential(size_t iE) const
    {
      return m_edge_potentials[iE];
    }
 
    inline const Eigen::MatrixXd & edgePotential(const Edge & E) const
    {
      return m_edge_potentials[E.global_index()];
    }
 
    inline const PlatesCore::CellBases & cellBases(size_t iT) const
    {
      return m_plates_core.cellBases(iT);
    }
 
    inline const PlatesCore::CellBases & cellBases(const Cell & T) const
    {
      return m_plates_core.cellBases(T.global_index());
    }
        
    inline const PlatesCore::EdgeBases & edgeBases(size_t iE) const
    {
      return m_plates_core.edgeBases(iE);
    }
 
    inline const PlatesCore::EdgeBases & edgeBases(const Edge & E) const
    {
      return m_plates_core.edgeBases(E.global_index());
    }
    
    inline const SymmetricMatrixBasisVertex & SymBasisVertex() const
    {
      return m_symm_vertex;
    }
    
  private:
    LocalOperators _compute_cell_divdiv_potential(size_t iT);
    Eigen::MatrixXd _compute_edge_potential(size_t iE);
 
    const PlatesCore & m_plates_core;
    bool m_use_threads;
    std::ostream & m_output;
 
    // Basis for the space of symmetric tensor-valued fields on a vertex
    SymmetricMatrixBasisVertex m_symm_vertex;
 
    // Containers for local operators
    std::vector<std::unique_ptr<LocalOperators> > m_cell_operators;
    std::vector<Eigen::MatrixXd> m_edge_potentials;
  };
  
} // end of namespace HArDCore2D
 
#endif