function/formula.hh

Go to the documentation of this file.

/* Formulas 
   Formulas defined on FE solution
 */
#ifndef funFormula_hh
#define funFormula_hh

#include <map>
#include "basics/debug.hh"
#include "basics/outputMatlab.hh"
#include "space/formula.hh"
#include "space/function.hh"
#include "function/vector.hh"

#define PWFormVecBaseConstr_D 0

namespace concepts {

  // *************************************************** PiecewiseVectorBase **

  template<class F, class G, class H, class I>
  class PiecewiseFormulaVectorBase : public PiecewiseFormulaBase<H> {
  public:
    PiecewiseFormulaVectorBase(const Space<G>& spc, const Vector<F>& coeff,
                               const ElementFunction<I,G>& fun);
  protected:
    const Space<G>& spc_;
    Array<F> coeff_;
    const ElementFunction<I,G>* fun_;
    std::map<uint, const ElementWithCell<G>*> elm_;
    const ElementWithCell<G>* element_(const Connector& cell) const;
  };

  template<class F, class G, class H, class I>
  PiecewiseFormulaVectorBase<F,G,H,I>::PiecewiseFormulaVectorBase
  (const Space<G>& spc, const Vector<F>& coeff,
   const ElementFunction<I,G>& fun) 
    : spc_(spc), coeff_(spc_->dim()), fun_(&fun)
  {
    memorycpy((F*)coeff_, (const F*)coeff, spc_->dim());
    DEBUGL(PWFormVecBaseConstr_D, "Space = " << *spc_);
    // scanner over the elements
    Scan<Element<G> >* sc = spc_->scan();
    while(*sc) {
      const Element<G>& e = (*sc)++;
      DEBUGL(PWFormVecBaseConstr_D, "Element = " << e);
      const ElementWithCell<G>* elm = 
        dynamic_cast<const ElementWithCell<G>*>(&e);
      conceptsAssert(elm, concepts::Assertion());
      elm_[elm->cell().connector().key()] = elm;
    }
    DEBUGL(PWFormVecBaseConstr_D, "Mapping Cell -> Element = " <<
     (OutputMatlab<std::map<uint, const ElementWithCell<G>*> >(elm_)));
  }

  template<class F, class G, class H, class I>
  const ElementWithCell<G>* 
  PiecewiseFormulaVectorBase<F,G,H,I>::element_(const Connector& cell) const {
    // look for element to which the cell belongs to
    typename std::map<uint, const ElementWithCell<G>*>::
      const_iterator i = elm_.find(cell.key());
    // if such a cell exists return the element
    if (i != elm_.end()) return i->second;
    // look at the children cell and call the method recursivly
    const Connector* child = 0;
    for(uint j = 0; child = cell.child(j); ++j) {
      const ElementWithCell<G>* elm = element_(*child);
      // if a children cell belongs to a element return the element
      if (elm) return elm;
    }
    // no element found
    return 0;
  }

  // ************************************************** PiecewiseVector<dim> **

  template<uint dim, class F, class G, class H>
  class PiecewiseFormulaVector : 
    public PiecewiseFormulaVectorBase<F,G,Point<H,dim>,H> {
    public:
    PiecewiseFormulaVector(const Space<G>& spc, const Vector<F>& coeff,
                           const ElementFunction<H,G>& fun);
    virtual PiecewiseFormulaVector<dim,F,G,H>* clone() const {
      Array<F> c((const F*)this->coeff_, this->spc_.dim());
      const Vector<F> coeff(*this->spc_, (F*)c);
      return new PiecewiseFormulaVector<dim,F,G,H>(coeff, *this->fun_);
    }
    virtual Point<H,dim> operator()(const Connector& cell, const Real p,
            const Real t = 0.0) const;
    virtual Point<H,dim> operator()(const Connector& cell, const Real2d& p,
            const Real t = 0.0) const;
    virtual Point<H,dim> operator()(const Connector& cell, const Real3d& p,
            const Real t = 0.0) const;
  };

  template<uint dim, class F, class G, class H>
  PiecewiseFormulaVector<dim,F,G,H>::PiecewiseFormulaVector
  (const Space<G>& spc, const Vector<F>& coeff,
   const ElementFunction<H,G>& fun) :
    PiecewiseFormulaVectorBase<F,G,Point<H,dim>,H>(spc, coeff, fun) {}

  template<uint dim, class F, class G, class H>
  Point<H,dim> PiecewiseFormulaVector<dim,F,G,H>::operator()
    (const Connector& cell, const Real p, const Real t) const
  {
    const ElementWithCell<G>* elm = this->element_(cell);
    Point<H,dim> res;
    if (elm) {
      // get's function value
      Array<H> val;
      (*this->fun_)(*elm, this->coeff_, val, p, t);
      conceptsAssert(val.size() == dim, Assertion());
      // copy it to a point
      memorycpy((H*)res, (const H*)val, dim);
    }
    return res;
  }

  template<uint dim, class F, class G, class H>
  Point<H,dim> PiecewiseFormulaVector<dim,F,G,H>::operator()
    (const Connector& cell, const Real2d& p, const Real t) const
  {
    const ElementWithCell<G>* elm = this->element_(cell);
    Point<H,dim> res;
    if (elm) {
      // get's function value
      Array<H> val;
      (*this->fun_)(*elm, this->coeff_, val, p, t);
      conceptsAssert(val.size() == dim, Assertion());
      // copy it to a point
      memorycpy((H*)res, (const H*)val, dim);
    }
    return res;
  }

  template<uint dim, class F, class G, class H>
  Point<H,dim> PiecewiseFormulaVector<dim,F,G,H>::operator()
    (const Connector& cell, const Real3d& p, const Real t) const
  {
    const ElementWithCell<G>* elm = this->element_(cell);
    Point<H,dim> res;
    if (elm) {
      // get's function value
      Array<H> val;
      (*this->fun_)(*elm, this->coeff_, val, p, t);
      conceptsAssert(val.size() == dim, Assertion());
      // copy it to a point
      memorycpy((H*)res, (const H*)val, dim);
    }
    return res;
  }

  // ********************************************* PiecewiseFormulaVector<1> **

  template<class F, class G, class H>
  class PiecewiseFormulaVector<1,F,G,H> : 
    public PiecewiseFormulaVectorBase<F,G,H,H> {
    public:
    PiecewiseFormulaVector(const Space<G>& spc, const Vector<F>& coeff,
                           const ElementFunction<H,G>& fun);
    virtual PiecewiseFormulaVector<1,F,G,H>* clone() const {
      Array<F> c((const F*)this->coeff_, this->spc_.dim());
      const Vector<F> coeff(*this->spc_, (F*)c);
      return new PiecewiseFormulaVector<1,F,G,H>(coeff, *this->fun_);
    }
    virtual H operator()(const Connector& cell, const Real p,
       const Real t = 0.0) const;
    virtual H operator()(const Connector& cell, const Real2d& p,
       const Real t = 0.0) const;
    virtual H operator()(const Connector& cell, const Real3d& p,
       const Real t = 0.0) const;
  };

  template<class F, class G, class H>
  PiecewiseFormulaVector<1,F,G,H>::PiecewiseFormulaVector
  (const Space<G>& spc, const Vector<F>& coeff,
   const ElementFunction<H,G>& fun) :
    PiecewiseFormulaVectorBase<F,G,H,H>(spc, coeff, fun) {}

  template<class F, class G, class H>
  H PiecewiseFormulaVector<1,F,G,H>::operator()
    (const Connector& cell, const Real p, const Real t) const
  {
    const ElementWithCell<G>* elm = this->element_(cell);
    if (elm) {
      // get's function value
      Array<H> val;
      (*this->fun_)(*elm, this->coeff_, val, p, t);
      conceptsAssert(val.size() == 1, Assertion());
      return val[0];
    }
    return H(0);
  }

  template<class F, class G, class H>
  H PiecewiseFormulaVector<1,F,G,H>::operator()
    (const Connector& cell, const Real2d& p, const Real t) const
  {
    const ElementWithCell<G>* elm = this->element_(cell);
    DEBUGL(0, "Cell = " << cell);
    if (elm) {
      DEBUGL(0, "Element = " << *elm);
      // get's function value
      Array<H> val;
      (*this->fun_)(*elm, this->coeff_, val, p, t);
      conceptsAssert(val.size() == 1, Assertion());
      DEBUGL(0, "p = " << p << ", Value = " << val[0]);
      return val[0];
    }
    return H(0);
  }

  template<class F, class G, class H>
  H PiecewiseFormulaVector<1,F,G,H>::operator()
    (const Connector& cell, const Real3d& p, const Real t) const
  {
    DEBUGL(0, "Cell = " << cell);
    const ElementWithCell<G>* elm = this->element_(cell);
    if (elm) {
      DEBUGL(0, "Element = " << *elm);
      // get's function value
      Array<H> val;
      (*this->fun_)(*elm, this->coeff_, val, p, t);
      conceptsAssert(val.size() == 1, Assertion());
      DEBUGL(0, "p = " << p << ", Value = " << val[0]);
      return val[0];
    }
    return H(0);
  }

} // namespace concepts

#endif // funFormula_hh