elementFormula.hh

Go to the documentation of this file.

 
#ifndef elementFormVector_hh
#define elementFormVector_hh
 
#include "toolbox/hashMap.hh"
#include "formula/exceptions.hh"
#include "geometry/cell.hh"
#include "space/element.hh"
#include "space/formula.hh"
#include "function/vector.hh"
 
#include "basics/debug.hh"
 
#define ElmFormVectorBaseConstr_D 0
 
 
namespace concepts {
 
  template<typename F, typename G>
  class ElementFunction;     // declared in space/function.hh
 
  // **************************************************** ElementFormulaBase **
 
  template<typename F, class G, class H, class I>
  class ElementFormulaVectorBase : public ElementFormula<H,G> {
  public:
    ElementFormulaVectorBase(const Space<G>& spc, const Vector<F>& v,
                             const ElementFunction<I,G>& f);
    virtual ElementFormulaVectorBase<F,G,H,I>* clone() const = 0;
  protected:
    const Space<G>& spc_;
    const Vector<F> v_;
    std::unique_ptr<const ElementFunction<I,G> > f_;
 
    mutable const Element<G> *lastElm_, *lastElmOur_;
 
    mutable Array<F> coeff_;
 
    template<typename J, typename P>
    void compute_(const Element<G>& elm, Array<J>& val, const P& p,
                  const Real t = 0.0) const;
  private:
    std::unordered_map<uint, const Element<G>*> elementPointer_;
 
    uint getKey_(const Element<G>& elm) const;
    void getElement_(const Element<G>& elm, const Element<G>*& elmOur) const;
  };
 
  template<class F, class G, class H, class I>
  ElementFormulaVectorBase<F,G,H,I>::ElementFormulaVectorBase
  (const Space<G>& spc, const Vector<F>& v, const ElementFunction<I,G>& f) 
    : spc_(spc), v_(v), f_(f.clone()), lastElm_(0), lastElmOur_(0), coeff_(0)
  {
    DEBUGL(ElmFormVectorBaseConstr_D, "start");
    Scan<Element<G> >* sc = spc_.scan();
    while (!sc->eos()) {
      const Element<G>& elm = (*sc)++;
      uint key = getKey_(elm);
      elementPointer_.insert
        (std::pair<uint, const Element<G>*>(key, &elm));
      DEBUGL(ElmFormVectorBaseConstr_D, 
             "Inserted key = " << key << " -> elm = " << elm);
    }
    DEBUGL(ElmFormVectorBaseConstr_D, "done.");
  }
 
  template<class F, class G, class H, class I>
  uint ElementFormulaVectorBase<F,G,H,I>::getKey_(const Element<G>& elm) const{
    const ElementWithCell<G>* elmC =
      dynamic_cast<const ElementWithCell<G>*>(&elm);
    conceptsAssert(elmC, Assertion());
    const Cell& cell = elmC->cell();
    uint key = cell.connector().key().key();
    return key;
  }
 
  template<class F, class G, class H, class I>
  void ElementFormulaVectorBase<F,G,H,I>::getElement_
  (const Element<G>& elm, const Element<G>*& elmOur) const {
    uint key = getKey_(elm);
    typename std::unordered_map<uint, const Element<G>*>::const_iterator i =
      elementPointer_.find(key);
    if (i == elementPointer_.end())
      throw(ElementNotInDomainOfFormula(elm)); // this
    
    lastElmOur_ = elmOur = i->second;
    lastElm_    = &elm;
  }
 
  template<class F, class G, class H, class I>
  template<typename J, typename P>
  void ElementFormulaVectorBase<F,G,H,I>::compute_
  (const Element<G>& elm, Array<J>& val, const P& p, const Real t) const {
    const Element<G>* elmOur = 0;
    getElement_(elm, elmOur);
    conceptsAssert(elmOur, Assertion());
    DEBUGL(0, "elmOur = " << *elmOur);
    DEBUGL(0, "v_ = " << v_);
    elmOur->T().extract(v_, coeff_);
    DEBUGL(0, "coeff_ = " << coeff_);
    (*f_)(*elmOur, coeff_, val, p, t);
    DEBUGL(0, "val = " << val);
  }
 
  // ********************************************* ElementFormulaVector<dim> **
 
  template<uint dim, class F = Real, class G = F,
           class H = typename Realtype<F>::type>
  class ElementFormulaVector : 
    public ElementFormulaVectorBase<F,H,Point<G,dim>,G> {
  public:
    ElementFormulaVector(const Space<H>& spc, const Vector<F>& v,
                         const ElementFunction<G,H>& f);
    virtual Point<G,dim> operator() (const ElementWithCell<H>& elm,
                                     const Real p, const Real t = 0.0) const;
    virtual Point<G,dim> operator() (const ElementWithCell<H>& elm,
                                     const Real2d& p, const Real t = 0.0) const;
    virtual Point<G,dim> operator() (const ElementWithCell<H>& elm,
                                     const Real3d& p, const Real t = 0.0) const;
 
    virtual ElementFormulaVector<dim,F,G,H>* clone() const {
      return new ElementFormulaVector(this->spc_, this->v_, *this->f_);
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };
 
  template<uint dim, class F, class G, class H>
  ElementFormulaVector<dim,F,G,H>::ElementFormulaVector
  (const Space<H>& spc, const Vector<F>& v, const ElementFunction<G,H>& f) : 
    ElementFormulaVectorBase<F,H,Point<G,dim>,G>(spc, v, f)
  {
    conceptsAssert3(this->f_->n() == dim, Assertion(), 
                    "not matching number of vector components of element function");
  }
 
  template<uint dim, class F, class G, class H>
  Point<G,dim> ElementFormulaVector<dim,F,G,H>::operator()
    (const ElementWithCell<H>& elm, const Real p, const Real t) const {
    // get's function value
    Array<G> val;
    this->compute_(elm, val, p, t);
    conceptsAssert(val.size() == dim, Assertion());
    // copy it to a point
    Point<G,dim> res;
    memorycpy((G*)res, (const G*)val, dim);
    return res;
  }
 
  template<uint dim, class F, class G, class H>
  Point<G,dim> ElementFormulaVector<dim,F,G,H>::operator()
    (const ElementWithCell<H>& elm, const Real2d& p, const Real t) const {
    // get's function value
    Array<G> val;
    this->compute_(elm, val, p, t);
    conceptsAssert(val.size() == dim, Assertion());
    // copy it to a point
    Point<G,dim> res;
    memorycpy((G*)res, (const G*)val, dim);
    return res;
  }
 
 
  template<uint dim, class F, class G, class H>
  Point<G,dim> ElementFormulaVector<dim,F,G,H>::operator()
    (const ElementWithCell<H>& elm, const Real3d& p, const Real t) const {
    // get's function value
    Array<G> val;
    this->compute_(elm, val, p, t);
    conceptsAssert(val.size() == dim, Assertion());
    // copy it to a point
    Point<G,dim> res;
    memorycpy((G*)res, (const G*)val, dim);
    return res;
  }
 
  template<uint dim, class F, class G, class H>
  std::ostream& ElementFormulaVector<dim,F,G,H>::info(std::ostream& os) const {
    return os << concepts::typeOf(*this) << "(" << *this->f_ << ", "
              << this->v_ << ")";
  }
 
  // *********************************************** ElementFormulaVector<1> **
 
  template<class F, class G, class H>
  class ElementFormulaVector<1,F,G,H> : 
    public ElementFormulaVectorBase<F,H,G,G> {
  public:
    ElementFormulaVector(const Space<H>& spc, const Vector<F>& v,
                         const ElementFunction<G,H>& f);
    virtual G operator() (const ElementWithCell<H>& elm, const Real p,
                          const Real t = 0.0) const;
    virtual G operator() (const ElementWithCell<H>& elm, const Real2d& p,
                          const Real t = 0.0) const;
    virtual G operator() (const ElementWithCell<H>& elm, const Real3d& p,
                          const Real t = 0.0) const;
 
    virtual ElementFormulaVector<1,F,G,H>* clone() const {
      return new ElementFormulaVector(this->spc_, this->v_, *this->f_); 
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };
 
  template<class F, class G, class H>
  ElementFormulaVector<1,F,G,H>::ElementFormulaVector
  (const Space<H>& spc, const Vector<F>& v, const ElementFunction<G,H>& f) : 
    ElementFormulaVectorBase<F,H,G,G>(spc, v, f) {
    conceptsAssert3(this->f_->n() == 1, Assertion(),
                    "only scalar element functions");
  }
 
  template<class F, class G, class H>
  G ElementFormulaVector<1,F,G,H>::operator()
    (const ElementWithCell<H>& elm, const Real p, const Real t) const {
    Array<G> res;
    this->compute_(elm, res, p, t);
    return res[0];
  }
 
  template<class F, class G, class H>
  G ElementFormulaVector<1,F,G,H>::operator()
    (const ElementWithCell<H>& elm, const Real2d& p, const Real t) const {
    Array<G> res;
    this->compute_(elm, res, p, t);
    return res[0];
  }
 
 
  template<class F, class G, class H>
  G ElementFormulaVector<1,F,G,H>::operator()
    (const ElementWithCell<H>& elm, const Real3d& p, const Real t) const {
    Array<G> res;
    this->compute_(elm, res, p, t);
    return res[0];
  }
 
  template<class F, class G, class H>
  std::ostream& ElementFormulaVector<1,F,G,H>::info(std::ostream& os) const {
    
    return os << "ElementFormulaVector<1>(" << *this->f_ << ", " << this->v_
        << ")";
  }
 
} // namespace concepts
 
#endif // elementFormVector_hh