function/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"


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);
  protected:
    const Space<G>& spc_;
    const Vector<F>& v_;
    std::auto_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:
    __gnu_cxx::hash_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)
  {
    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));
    }
  }

  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 __gnu_cxx::hash_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());
    elmOur->T().extract(v_, coeff_);
    (*f_)(*elmOur, coeff_, val, p, t);
  }

  // ********************************************* 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 << "ElementFormulaVector<" << dim << ">(" << *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