formula.hh

Go to the documentation of this file.

 
#ifndef geomFormula_hh
#define geomFormula_hh
 
#include <functional>
#include "basics/debug.hh"
#include "toolbox/dynArray.hh"
#include "toolbox/hashMap.hh"
#include "toolbox/set.hh"
#include "formula/piecewiseFormula.hh"
#include "formula/elementFormulaContainer.hh"
#include "basics/cloneable.hh"
#include "basics/operations.hh"
#include "cell.hh"
#include "mesh.hh"
 
namespace concepts {
 
  // forward declaration
  template<class F>
  class Formula;          // declared in formula/formula.hh
 
  // ***************************************************** PiecewiseFormula0 **
 
  template<typename F>
  class PiecewiseFormula0 : public PiecewiseFormulaBase<F> {
  public:
    PiecewiseFormula0(const Formula<F>& formula) : formula_(formula.clone()) {}
    virtual ~PiecewiseFormula0() {}
    virtual PiecewiseFormula0<F>* clone() const {
      return new PiecewiseFormula0(*formula_);
    }
    virtual F operator()(const Connector& cell, const Real p,
       const Real t = 0.0) const {
      return (*formula_)(p,t); }
    virtual F operator()(const Connector& cell, const Real2d& p,
       const Real t = 0.0) const {
      return (*formula_)(p,t); }
    virtual F operator()(const Connector& cell, const Real3d& p,
       const Real t = 0.0) const {
      return (*formula_)(p,t); }
  protected:
    virtual std::ostream& info(std::ostream& os) const  {
      return os << concepts::typeOf(*this)<<"(" << *formula_ << ")"; }
  private:
    std::unique_ptr<Formula<F> > formula_;
  };
 
  // ************************************************* PiecewiseConstFormula **
 
  template<typename F = Real>
  class PiecewiseConstFormula : public PiecewiseFormulaBase<F> {
  public:
    PiecewiseConstFormula(const F dflt_value = 0);
    PiecewiseConstFormula(const HashMap<F> formula, const F dflt_value = 0);
    PiecewiseConstFormula(const std::map<int,F> formula, const F dflt_value = 0);
    virtual ~PiecewiseConstFormula();
    virtual PiecewiseConstFormula<F>* clone() const;
 
    inline virtual F operator()(const Connector& cell, const Real p,
        const Real t = 0.0) const
      { return (*this)[cell.attrib()]; }
 
    inline virtual F operator()(const Connector& cell, const Real2d& p,
        const Real t = 0.0) const
      { return (*this)[cell.attrib()]; }
 
/*    inline virtual F operator()(const Element<Real>& elm_, const Real2d& p,
        const Real t = 0.0) const
      { 
        //FIXME: this conversion should not be neccesary
        const ElementWithCell<Real>& elm = static_cast<const ElementWithCell<Real>&> (elm_);
        const Connector& cell = elm.cell().connector();
        return operator[](cell.attrib()); 
      }*/
 
    inline virtual F operator()(const Connector& cell, const Real3d& p,
        const Real t = 0.0) const
      { return (*this)[cell.attrib()]; }
    virtual const F& operator[](const uint attrib) const;
    virtual F& operator[](const uint attrib);
    // Multiplies all values with this factor
    template<class H>
    PiecewiseConstFormula<F>& operator*=(const H& factor);
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    HashMap<F> formula_;
  };
 
  template<typename F>
  template<class H>
  PiecewiseConstFormula<F>&
  PiecewiseConstFormula<F>::operator*=(const H& factor){
    typename HashMap<F>::iterator i = formula_.begin();
    while(i != formula_.end())
      i++->second *= factor;
    return *this;
  }
 
  // ******************************************* PiecewiseConstImportFormula **
 
  template<typename F = Real>
  class PiecewiseConstImportFormula : public PiecewiseConstFormula<F> {
  public:
    PiecewiseConstImportFormula(const std::string material,
        const uint idxStart = 1,
        const F dflt_value = 0);
  private:
    void read_(std::istringstream& iss);
  };
 
  // ***************************************** PiecewiseConstDynArrayFormula **
 
  template<typename F = Real>
  class PiecewiseConstDynArrayFormula : public PiecewiseFormulaBase<F> {
  public:
    PiecewiseConstDynArrayFormula(const DynArray<F>& formula);
    ~PiecewiseConstDynArrayFormula();
    virtual PiecewiseConstDynArrayFormula<F>* clone() const;
    inline virtual F operator()(const Connector& cell, const Real p,
        const Real t = 0.0) const
      { return (*formula_)[cell.attrib()]; }
    inline virtual F operator()(const Connector& cell, const Real2d& p,
        const Real t = 0.0) const
      { return (*formula_)[cell.attrib()]; }
    inline virtual F operator()(const Connector& cell, const Real3d& p,
        const Real t = 0.0) const
      { return (*formula_)[cell.attrib()]; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const DynArray<F>* formula_;
  };
 
  // ****************************************************** PiecewiseFormula **
 
  template<typename F>
  class PiecewiseFormula : public PiecewiseFormulaBase<F> {
  public:
    PiecewiseFormula(const F dflt_value = 0);
    virtual ~PiecewiseFormula();
    virtual PiecewiseFormula<F>* clone() const;
    virtual F operator()(const Connector& cell, const Real p,
                         const Real t = 0.0) const;
    virtual F operator()(const Connector& cell, const Real2d& p,
                         const Real t = 0.0) const;
    virtual F operator()(const Connector& cell, const Real3d& p,
                         const Real t = 0.0) const;
    void set(const uint attrib, const Formula<F>& formula);
    void set(const uint attrib, const F a);
    template<class H>
    void multiply(const H& factor);
    template<class H>
    void multiply(const uint attrib, const H& factor);
  protected:
    class FormulaFactor : public OutputOperator {
    public:
      FormulaFactor() : formula_(nullptr) {}
      FormulaFactor(const FormulaFactor& frm) 
        : formula_(nullptr), factor_(frm.factor_) { 
        formula_.reset(frm.formula_.get());
      }
      std::unique_ptr<Formula<F> > formula_;
      typename Datatype<F>::type factor_;
    protected:
      virtual std::ostream& info(std::ostream& os) const { return os;}
    };
    virtual std::ostream& info(std::ostream& os) const;
    HashMap<FormulaFactor> formula_;
  };
 
  template<typename F>
  template<class H>
  void PiecewiseFormula<F>::multiply(const H& factor){
    typename HashMap<FormulaFactor>::iterator i = formula_.begin();
    while(i != formula_.end())
      i++->second.factor_ *= factor;
  }
 
  template<typename F>
  template<class H>
  void PiecewiseFormula<F>::multiply(const uint attrib, const H& factor){
    typename HashMap<FormulaFactor>::iterator i = formula_.find(attrib);
    if (i != formula_.end())
      i->second.factor_ *= factor;
  }
 
 
  // *********************************************** PiecewiseElementFormula **
 
  template<class F, class G = typename Realtype<F>::type>
  class PiecewiseElementFormula : public ElementFormula<F, G> {
  public:
    PiecewiseElementFormula(const F defaultV = 0)
    : defaultV_(defaultV)
    {
    }
 
    virtual ~PiecewiseElementFormula() { }
 
    virtual PiecewiseElementFormula<F, G>* clone() const { 
      return new PiecewiseElementFormula(*this); 
    }
 
    virtual F operator() (const ElementWithCell< G > &elm, const Real3d &p, 
                          const Real t=0.0) const
    {
      typename ElemFormulaMap::const_iterator 
        i = formulas_.find(elm.cell().connector().attrib());
      if (i != formulas_.end()) 
        return i->second.operator()(elm, p, t);
 
      return defaultV_;
    }
 
    virtual F operator() (const ElementWithCell< G > &elm, const Real2d &p, 
                          const Real t=0.0) const
    {
      typename ElemFormulaMap::const_iterator 
        i = formulas_.find(elm.cell().connector().attrib());
      if (i != formulas_.end()) 
        return i->second.operator()(elm, p, t);
 
      return defaultV_;
    }
 
    virtual F operator() (const ElementWithCell< G > &elm, const Real p, 
                          const Real t=0.0) const
    {
      typename ElemFormulaMap::const_iterator 
        i = formulas_.find(elm.cell().connector().attrib());
      if (i != formulas_.end()) 
        return i->second.operator()(elm, p, t);
 
      return defaultV_;
    }
 
    void set(const uint attrib, 
             const ElementFormulaContainer<F,G> formula) {
      formulas_[attrib] = formula;
      //std::cout <<  "atrib: " << attrib << ", " << *formulas[attrib] << std::endl;
    }
 
    void set(concepts::Set<uint> attribs,
             const ElementFormulaContainer<F,G> formula) 
    {
      for (concepts::Set<uint>::const_iterator it = attribs.begin();
          it != attribs.end(); ++it)
      {
        set(*it, formula);
      }
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << concepts::typeOf(*this)<<"( " << defaultV_ << ", "
                << formulas_ << ")" << std::endl;
    }
  public:
    typedef 
    HashMap< ElementFormulaContainer<F, G> > ElemFormulaMap;
    ElemFormulaMap formulas_;
    F defaultV_;
  };
 
 
  // ************************************************************* Operation **
 
  template<class F>
  class Operation : public Cloneable, public OutputOperator {
  public:
    virtual F operator()(const F val) const = 0;
    virtual Operation<F>* clone() const = 0;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };
 
  // **************************************************************** OpMult **
 
  template<class F>
  class OpMult : public Operation<F> {
  public:
    OpMult(const F factor) : factor_(factor) {}
    virtual OpMult<F>* clone() const { return new OpMult<F>(factor_); };
    virtual F operator()(const F val) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const F factor_;
  };
 
  // ***************************************************************** OpAdd **
 
  template<class F>
  class OpAdd : public Operation<F> {
  public:
    OpAdd(const F summand) : summand_(summand) {}
    virtual OpAdd<F>* clone() const { return new OpAdd<F>(summand_); };
    virtual F operator()(const F val) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const F summand_;
  };
 
  // ************************************************************** OpRecipr **
 
  template<class F>
  class OpRecipr : public Operation<F> {
  public:
    virtual F operator()(const F val) const;
    virtual OpRecipr<F>* clone() const { return new OpRecipr<F>; };
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };
 
  // *************************************************** PiecewiseFormulaFun **
 
  template<typename F, typename G = F>
  class PiecewiseFormulaFun : public PiecewiseFormulaBase<F> {
  public:
    PiecewiseFormulaFun(PiecewiseFormulaBase<G>& formula);
    virtual ~PiecewiseFormulaFun();
    virtual PiecewiseFormulaFun<F,G>* clone() const;
    virtual F operator()(const Connector& cell, const Real p,
       const Real t = 0.0) const;
    virtual F operator()(const Connector& cell, const Real2d& p,
       const Real t = 0.0) const;
    virtual F operator()(const Connector& cell, const Real3d& p,
       const Real t = 0.0) const;
 
    PiecewiseFormulaFun<F,G>& operator*=(const F n);
    PiecewiseFormulaFun<F,G>& operator+=(const F n);
    PiecewiseFormulaFun<F,G>& reciprocal();
    void clear();
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    PiecewiseFormulaBase<G>& formula_;
    std::vector<Operation<F>*> operations_;
  };
 
 
  // *********************************************** PiecewiseFormulaCombine **
 
  template<typename F, typename G, typename H = G,
     typename I = multiplies<G,H,F> >
  class PiecewiseFormulaCombine : public PiecewiseFormulaBase<F> {
  public:
    PiecewiseFormulaCombine(const PiecewiseFormulaBase<G>& formula1,
                            const PiecewiseFormulaBase<H>& formula2);
    virtual ~PiecewiseFormulaCombine() {}
    virtual PiecewiseFormulaCombine<F,G,H,I>* clone() const;
    virtual F operator()(const Connector& cell, const Real p,
       const Real t = 0.0) const;
    virtual F operator()(const Connector& cell, const Real2d& p,
       const Real t = 0.0) const;
    virtual F operator()(const Connector& cell, const Real3d& p,
       const Real t = 0.0) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    std::unique_ptr<const PiecewiseFormulaBase<G> > formula1_;
    std::unique_ptr<const PiecewiseFormulaBase<H> > formula2_;
    const I fun_;
  };
 
  // ************************************************ CurvatureElementFormula **
 
  class CurvatureElementFormula : public ElementFormula<Real> {
  public:
    CurvatureElementFormula(uint n = 0) : n_(n) {}
    virtual Real operator() (const ElementWithCell<Real>& elm, const Real p,
                             const Real t = 0.0) const;
    virtual Real operator() (const ElementWithCell<Real>& elm, const Real2d& p,
                             const Real t = 0.0) const;
    virtual Real operator() (const ElementWithCell<Real>& elm, const Real3d& p,
                             const Real t = 0.0) const;
    virtual CurvatureElementFormula* clone() const {
      return new CurvatureElementFormula();
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    uint n_;
  };
 
 
} // namespace concepts
 
#endif // geomFormula_hh