formula/matrixElementFormula.hh

Go to the documentation of this file.

#ifndef matrixelementformula_hh
#define matrixelementformula_hh

#include <vector>
#include "toolbox/sharedPointer.hh"
#include "elementFormula.hh"

namespace concepts {


  // ************************************************** MatrixElementFormula **

  template<typename F, int dim, typename G = typename Realtype<F>::type>
  class MatrixElementFormula : public ElementFormula< Mapping<F, dim>, G >
  {
  public:
    typedef Mapping<F, dim> value_type;

    MatrixElementFormula<F,dim,G>() : formulas_() 
    {
      for(int i=0; i < dim*dim; ++i)
        formulas_.push_back(ElementFormulaContainer<F>());
    }


    MatrixElementFormula<F,dim,G>(const std::vector< ElementFormulaContainer<F> >& formulas) 
     : formulas_(formulas)
    { 
      conceptsAssert(formulas_.size() == dim * dim, concepts::Assertion());
    }

    Mapping<F, dim> operator() (const ElementWithCell<G>& elm, const Real p,
                                const Real t = 0.0) const
    {
      Mapping<F, dim> ret;
      for (int i=0; i < dim; ++i) {
        for (int j=0; j < dim; ++j) 
          ret(i,j) = 
            !formulas_[i*dim + j].empty() ? formulas_[i*dim + j](elm, p, t) : F(0);
      }
      return ret;
    }

    Mapping<F, dim> operator() (const ElementWithCell<G>& elm, const Real2d& p,
                                const Real t = 0.0) const
    {
      Mapping<F, dim> ret;
      for (int i=0; i < dim; ++i) {
        for (int j=0; j < dim; ++j)
          ret(i,j) = 
            !formulas_[i*dim + j].empty() ? formulas_[i*dim + j](elm, p, t) : F(0);
      }
      return ret;
    }

    Mapping<F, dim> operator() (const ElementWithCell<G>& elm, const Real3d& p,
                                const Real t = 0.0) const
    {
      Mapping<F, dim> ret;
      for (int i=0; i < dim; ++i) {
        for (int j=0; j < dim; ++j)
          ret(i,j) = 
            !formulas_[i*dim + j].empty() ? formulas_[i*dim + j](elm, p, t) : F(0);
      }
      return ret;
    }

    virtual MatrixElementFormula<F,dim,G>* clone() const
    {
      return new MatrixElementFormula<F, dim, G>(*this);
    }

#if 0
    MatrixElementFormula<F,dim,G>( const MatrixElementFormula<F,dim,G>& other) 
      : formulas_(other.formulas_)
    { }
#endif

    virtual ~MatrixElementFormula() { }

    const std::vector< ElementFormulaContainer<F> >& getFormulas() 
    {
      return formulas_;
    }

    ElementFormulaContainer<F>& operator()(int i, int j)
    {
      conceptsAssert(i < dim*dim, concepts::Assertion());

      return formulas_[i*dim +j];
    }
  protected:
    std::vector< ElementFormulaContainer<F> > formulas_;
  };


} // namespace concepts


#endif // matrixelementformula_hh