geometry/formula.hh

Go to the documentation of this file.

/* Formulas 
   Formulas defined on the cells
 */

#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 << "PiecewiseFormula0(" << *formula_ << ")"; }
  private:
    std::auto_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_(0) {}
      FormulaFactor(const FormulaFactor& frm) 
        : formula_(0), factor_(frm.factor_) { 
        formula_.reset(frm.formula_.get());
      }
      std::auto_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 << "PiecewiseElementFormula( " << 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::auto_ptr<const PiecewiseFormulaBase<G> > formula1_;
    std::auto_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