formulafrmEformula.hh

Go to the documentation of this file.

 
#ifndef formulatfrmEformula_hh
#define formulatfrmEformula_hh
 
#include "geometry.hh"
#include "formula.hh"
#include "toolbox/sharedPointer.hh"
 
#include <algorithm>
 
#if __GNUC__ == 2
#  include <float.h>
#  define EPS DBL_EPSILON
#else
#  include <limits>
#  define EPS std::numeric_limits<double>::epsilon()
#endif
 
namespace concepts{
 
  // ************************************************************* JacobianCell *
 
  template<uint dimC> class JacobianCell;
 
  template<> class JacobianCell<1>{
  public:
    typedef concepts::Edge1d cell;
  };
 
  template<> class JacobianCell<2>{
  public:
    typedef concepts::Quad2d cell;
  };
 
  template<> class JacobianCell<3>{
  public:
    typedef concepts::Hexahedron3d cell;
  };
 
 
  // ********************************************************** TensorVertexMap *
 
  template<uint dimC> class TensorVertexMap;
 
  template<> class TensorVertexMap<1>{
  public:
    typedef concepts::MapEdge1d vtxmap;
  };
 
  template<> class TensorVertexMap<2>{
  public:
    typedef concepts::VertexQuad2d vtxmap;
  };
 
  template<> class TensorVertexMap<3>{
  public:
    typedef concepts::MapHexahedron3d vtxmap;
  };
 
 
 
  template<uint dimC>
  struct CellStripeElement : public OutputOperator{
 
    CellStripeElement(const typename JacobianCell<dimC>::cell& cell)
      : cell(cell), key(cell.connector().key()){};
 
    virtual bool operator<(const CellStripeElement<dimC>& sb) const{
      return this->key < sb.key;
    }
 
    virtual bool operator==(const CellStripeElement<dimC>& sb) const{
      return key == sb.key;
    }
 
    const typename JacobianCell<dimC>::cell& cell;
    const uint key;
 
 
  protected:
    virtual std::ostream& info(std::ostream& os) const{
      return os<< concepts::typeOf(*this)<<"[ "<< key << " ] ";
    }
  };
 
 
 
  template<uint dimC>
  class CellBox : public OutputOperator{
  public:
    CellBox(const typename JacobianCell<dimC>::cell& cell):
      cell(cell), isDummy(false), key(cell.connector().key().key()){
 
      for(uint j = 0; j < dimC; ++j)
        minX[j]  = maxX[j] = cell.vertex(0)[j];
 
      for(uint i = 1; i < pow(2,dimC); ++i){
 
        for(uint j = 0; j < dimC; ++j){
          maxX[j] = std::max(maxX[j], cell.vertex(i)[j] );
          minX[j] = std::min(minX[j], cell.vertex(i)[j] );
        }
      }//for i
    }
 
 
    //  //dummy constructor
    CellBox(const typename JacobianCell<dimC>::cell& cell, Real xMax):
      cell(cell), isDummy(true), key(cell.connector().key().key()){
 
      minX[0] = - std::numeric_limits<double>::max();
      maxX[0] = xMax;
 
      for(uint i = 1 ; i < dimC; ++i){
        minX[i] = maxX[i] = - std::numeric_limits<double>::max();
      }
 
    }
 
    virtual bool operator<(const CellBox<dimC>& box) const;
 
 
    virtual bool operator==(const CellBox<dimC>& box) const{
      //safety hold ask for isDummy
      return this->key == box.key && isDummy == box.isDummy;
    }
 
    Real minX[dimC];
    Real maxX[dimC];
 
    const typename JacobianCell<dimC>::cell& cell;
    bool isDummy;
    uint key;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const{
      os<< concepts::typeOf(*this)<<"[key = "<< key<<"] << box frame = ";
      for(uint i = 0 ; i < dimC-1; ++i)
        os << "["<<minX[i] << " , "<< maxX[i] << "] x ";
      os << "["<<minX[dimC-1] << " , "<< maxX[dimC-1] << "]";
      return os;
    }
  };
 
 
  template<uint dimC>
  struct CCell_F : public OutputOperator{
 
    CCell_F(const typename JacobianCell<dimC>::cell& cell, const uint cardF=2):
      cell(cell), key(cell.connector().key()), cardF(cardF){
      conceptsAssert(cardF==2, concepts::Assertion());
      computeF();
    }
 
 
    bool operator==(CCell_F cCF) const{
      return this->key == cCF.key;
    }
 
    bool operator<(CCell_F cCF) const{
      return this->key < cCF.key;
    }
 
 
    const typename JacobianCell<dimC>::cell& cell;
    const uint key;
    const uint cardF;
 
    //Finite set of physical points
    Sequence<Point<Real,dimC> > F;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const{
      return os << concepts::typeOf(*this)<<"[ Cellkey = "<< key << " ], F = " << F;
    }
    void computeF();
  };
 
 
  template<uint dimC>
  struct CCell_F_dist : public OutputOperator{
    CCell_F_dist(const CCell_F<dimC>& ccF, const Point<Real,dimC>& P):
      ccF(ccF), dist_2((ccF.F[0]-P).l2_2()) {
 
      for(uint i = 1; i < ccF.F.size(); ++i){
        dist_2 = std::min(dist_2, (ccF.F[i]-P).l2_2());
      }
    }
    bool operator==(CCell_F_dist ccF_d) const{
      return this->ccF.key == ccF_d.ccF.key;
    }
 
    bool operator<(CCell_F_dist ccF_d) const{
      return this->dist_2 < ccF_d.dist_2;
    }
    const CCell_F<dimC>& ccF;
    //squared distance of P and F
    Real dist_2;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const{
      return os << concepts::typeOf(*this)<<"[] ";
    }
  };
 
 
 
  template<uint dimC,class G> class MultiVector;
 
 
  // ********************************* specific MultiVector classes **
 
  template<> class MultiVector<0,const ElementWithCell<Real>* >{
  public:
    typedef const ElementWithCell<Real>* map;
  };
  template<> class MultiVector<0,const ElementWithCell<Cmplx>* >{
  public:
    typedef const ElementWithCell<Cmplx>* map;
  };
 
 
 
  template<> class MultiVector<0,concepts::Set<CellBox<1>>>{
  public:
    typedef std::pair<concepts::Set<CellBox<1>>, Real>  map;
  };
  template<> class MultiVector<0,concepts::Set<CellBox<2>>>{
  public:
    typedef std::pair<concepts::Set<CellBox<2>>, Real>  map;
  };
  template<> class MultiVector<0,concepts::Set<CellBox<3>>>{
  public:
    typedef std::pair<concepts::Set<CellBox<3>>, Real>  map;
  };
 
 
  template<uint dimC,class G> class MultiVector{
  public:
    typedef  std::vector<typename MultiVector<dimC-1,G>::map > map;
  };
 
 
  // **************************************************** CellMap **
 
  template<uint dimC, class H>
  struct CellMap: public OutputOperator{
 
    CellMap();
 
    CellMap(const CellMap<dimC,H>& map);
 
    void resize(const uint (&maxlevel)[dimC]);
 
    //add information at given index
    void add(const ElementWithCell<H>*  elm, const uint (&idx)[dimC]);
 
    //@ pre idx is accessable
    const ElementWithCell<H>* get(const uint (&idx)[dimC]);
 
    uint sizes[dimC];
 
    //a better range safe multiarray, since 2^64 cannot be holded in a ulonglong!
    typename MultiVector<dimC, const ElementWithCell<H>* >::map map;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const{
      return os << concepts::typeOf(*this)<<"["<<dimC<<"]";
    }
  };
 
 
 
 
  // **************************************************** FFEF_Error **
 
  class FFEF_Error: public ExceptionBase {
  public:
    FFEF_Error(const std::string& file, const unsigned int line,
        const std::string& function, const std::string& errorMessage) throw ();
 
 
    virtual const char* what() const throw();
 
      virtual ~FFEF_Error() throw () {};
  protected:
      virtual std::ostream& info(std::ostream& os) const throw ();
  private:
      std::string errorMessage_;
      mutable std::string outputMessage_ ;
  };
 
 
 
  // ************************************************* FormulaFromElementFormula *
 
  template<uint dim, class F, typename G = typename Realtype<F>::type>
  class FormulaFromElementFormula: public concepts::Formula<F>{
 
  private:
    //this is the maximal refinement level allowed here, maximal practical level should be around < 32.
    static uint max_lvl;
 
  public:
 
    enum PNF {Exc = 0, NaN, Zero};
 
 
    //template<class F, class H>
    FormulaFromElementFormula(const concepts::ElementFormulaContainer<F,G> frm,
                   const concepts::SpaceOnCoarseCells<dim, G>& sSpc,
                   const Set<uint>* s_attrb = 0,
                   const Set<uint>* d_attrb = 0,
                   const Point<Real, dim> scale = Point<Real, dim>(1),
                   const Point<Real, dim> shift = Point<Real, dim>(0),
                   const uint N = 10,
                   const uint cardF_2 = 2,
                   const Real tol = EPS,
                   const Real interior = 0.00000001,
                   const uint maxDoI = 3
     );
 
    virtual ~FormulaFromElementFormula(){};
 
 
    virtual F operator() (const Real p,    const Real t = 0.0) const;
    virtual F operator() (const Real2d& p, const Real t = 0.0) const;
    virtual F operator() (const Real3d& p, const Real t = 0.0) const;
 
    virtual F operator() (const Connector& cntr,
                          const Real p,
                          const Real t = 0.0) const;
 
    virtual F operator() (const Connector& cntr,
                          const Real2d& p,
                          const Real t = 0.0) const;
 
    virtual F operator() (const Connector& cntr,
                          const Real3d& p,
                          const Real t = 0.0) const;
 
 
    void allowPeriodicity(const Point<Real,dim>& O, const Point<Real, dim>& eL);
 
 
    inline void setPNF(PNF pnf) {pnf_ = pnf;}
 
 
    virtual FormulaFromElementFormula<dim, F,G>* clone() const;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const {
      return os << "u_ = " << this->u_ << ", pnf_ = " << this->pnf_ ;
    }
 
 
  private:
 
    const concepts::ElementFormulaContainer<F,G> u_;
 
    const concepts::RCP<const concepts::Set<uint> > s_attrb_;
    const concepts::RCP<const concepts::Set<uint> > d_attrb_;
 
    const Point<Real, dim> scale_;
    const Point<Real, dim> shift_;
 
    const Real tol_;
    const uint N_;
    const uint cardF_2_;
    const Real interior_;
 
    //heuristic approach due to integration or graphical application
    //often next local coordinate is in the same element as the last one
    mutable const ElementWithCell<G>* lastElm_;
 
    mutable const typename JacobianCell<dim>::cell* lastCell_;
 
    mutable const typename JacobianCell<dim>::cell* lastCoarseCell_;
 
    mutable bool wasCurved_;
 
    mutable bool wasCurved_c_;
 
 
    mutable Z2 cCellType_;
 
    bool allowPer_;
 
    Point<Real, dim> O_;
    Point<Real, dim> eL_;
 
    const uint maxDoI_;
 
    //Flag for Point not found behaviour of underlying search routine.
    PNF pnf_;
 
 
 
    Set<const typename JacobianCell<dim>::cell* > coarsestCell_;
 
    Set<const typename JacobianCell<dim>::cell* > irrCCells_;
 
 
 
 
    //maps I in IN^dim to the coressponding child cells key
    //it holds 0<= I_i maxLevel_i
    // where maxLevel is the highest number of levels in diretion x_i
 
    //example :
    //
    //  coarsest cell        href4
    //   1_________         ________                     _________
    //    |       |         |   |   |  maxlevel_0 = 1    |   |   |  maxlevel_0 = 2
    //    |       |         |___|___|  maxlevel_1 = 1    |___|___|  maxlevel_1 = 1
    //    |       |         |   |   |                    | | |   |
    //    |_______|         |___|___|                    |_|_|___|
    //   0        1
 
    //  We devide the coarsest cell into sectors along coord i
    // 0 <=  l_i / (2^maxLevel_i) <= 1 ,  l_i = 0,..., 2^maxLevel_i
    //
    // Then  map(l_0, l_1, l_2) maps into a sector  surjectiv.
    //
    // Example :
    //      _____________
    //    1 |  e  |  d  |
    //      |_____|_____|
    //    0 | a| b| c   |
    //      |__|__|_____|
    //       0  1   2 3
    //
    // then for example map(0,0) = a, map(2,0) = map(3,0) = c
 
    //the map is from coarse elms to its finest children
    mutable concepts::HashMap<CellMap<dim,G> > mapToFinestElm_;
 
 
    mutable HashMap<HashMap<const concepts::ElementWithCell<G>* > > coarseToElm_;
 
 
    concepts::Set<CellBox<dim> > cBoxes_;
 
    concepts::Set<CCell_F<dim> > rCC_;
    Set<uint> mapCCells_;
 
 
    void build_directMap_(const typename JacobianCell<dim>::cell* cell,
        const concepts::HashMap<const concepts::ElementWithCell<G>* >& mapToElm);
 
 
    template<uint dim1 = dim,
          typename std::enable_if<dim1 == 1u, uint>::type = 0>
    void setLocalMap_(const typename JacobianCell<dim>::cell* coarseCell,
                     const HashMap<const ElementWithCell<G>* >& mapToElm,
                     CellMap<dim, G>& map,
                     const std::array<uint,dim> idxL,
                     const std::array<uint,dim> idxR) const;
 
    template<uint dim2 = dim,
          typename std::enable_if<dim2 == 2u, uint>::type = 0>
    void setLocalMap_(const typename JacobianCell<dim>::cell* coarseCell,
                     const HashMap<const ElementWithCell<G>* >& mapToElm,
                     CellMap<dim, G>& map,
                     const std::array<uint,dim> idxL,
                     const std::array<uint,dim> idxR) const ;
 
    template<uint dim3 = dim,
          typename std::enable_if<dim3 == 3u, uint>::type = 0>
    void setLocalMap_(const typename JacobianCell<dim>::cell* coarseCell,
                     const HashMap<const ElementWithCell<G>* >& mapToElm,
                     CellMap<dim, G>& map,
                     const std::array<uint,dim> idxL,
                     const std::array<uint,dim> idxR) const;
 
 
 
    FormulaFromElementFormula(const concepts::ElementFormulaContainer<F,G> frm,
                   const concepts::RCP<const concepts::Set<uint> > s_attrb,
                   const concepts::RCP<const concepts::Set<uint> > d_attrb,
                   const Point<Real, dim> scale,
                   const Point<Real, dim> shift,
                   const uint N,
                   const uint cardF_2 ,
                   const Real tol,
                   const Real interior,
                   const concepts::Set<CellBox<dim> > cBoxes,
                   const concepts::Set<CCell_F<dim> > rCC,
                   const bool allowPer,
                   const Point<Real, dim> O,
                   const Point<Real, dim> eL,
                   const Set<uint> mapCCells,
                   const Set<const typename JacobianCell<dim>::cell* > coarsestCell,
                   const Set<const typename JacobianCell<dim>::cell* > irrCCells,
                   const concepts::HashMap<CellMap<dim,G> > mapToFinestElm,
                   const uint maxDoI,
                   const PNF pnf,
                   const HashMap<HashMap<const concepts::ElementWithCell<G>* > > coarseToElm);
 
    F compute_(const Point<Real,dim>& p, const Real t) const;
 
    void buildBoxCandidates_(const Point<Real, dim> P, concepts::Set<CellStripeElement<dim> >& box) const;
 
    bool getOrigPoint_(//const concepts::ElementWithCell<G>& elm,
                       const typename JacobianCell<dim>::cell& cell,
                       const Point<Real, dim>  P, Point<Real, dim>& eta, const uint N,
                       const Real interior,
                       const Point<Real, dim>&  x0 = Real3d(0.49,0.51,0.49)) const;
 
 
    //
    // then applies original shift/scale to p to build P
 
    void shift_DBox_(const Point<Real,dim>& p, Point<Real,dim>& P) const;
 
 
 
    //a cell is called coarse iff its level of refinement is zero in each direction
    bool isCoarse_(const typename JacobianCell<dim>::cell& cell);
 
 
    void addElements_(const uint cKey,
                      const typename JacobianCell<dim>::cell& cell,
                      const HashMap<const concepts::ElementWithCell<G>* >& mapToElm,
                      HashMap<HashMap<const concepts::ElementWithCell<G>* > >& cCtE );
 
 
    //@pre arrays are initialized with zero entries
    //template<uint dim>
    void getExtremalLevel(const typename JacobianCell<dim>::cell* cell,
                          uint (&maxL)[dim], uint (&minL)[dim]);
 
 
    bool compute0_(const typename JacobianCell<dim>::cell& cCell, const Point<Real, dim>& P,
                   const Point<Real, dim>& cEta,const Real t, F& val) const;
 
 
    bool compute1_(const typename JacobianCell<dim>::cell& cCell, const Point<Real, dim>& P,
                   const Point<Real, dim>& cEta,const Real t, F& val) const;
 
    bool checkCCells_(const std::set<typename CellType<dim>::cell *>& allCell,
        typename std::set<typename CellType<dim>::cell *>::const_iterator iter);
 
 
    F handlePNF_() const;
 
 
 
 
 
  };
 
#ifdef __clang__
  void getChild(const concepts::Edge1d& cCell,const Point<Real,1u>& eta,
                const std::array<Real,std::size_t(1u)> x0,
                const std::array<Real,std::size_t(1u)> xN,
                const concepts::Edge1d*& cld);
  void getChild(const concepts::Quad2d& cCell,const Point<Real,2u>& eta,
                const std::array<Real,std::size_t(2u)> x0,
                const std::array<Real,std::size_t(2u)> xN,
                const concepts::Quad2d*& cld);
  void getChild(const concepts::Hexahedron3d& cCell,const Point<Real,3u>& eta,
                const std::array<Real,std::size_t(3u)> x0,
                const std::array<Real,std::size_t(3u)> xN,
                const concepts::Hexahedron3d*& cld);
 
#else
 
  template<uint dimC>
  void getChild(const typename JacobianCell<dimC>::cell& cCell,const Point<Real,dimC>& eta,
                const std::array<Real,std::size_t(dimC)> x0,
                const std::array<Real,std::size_t(dimC)> xN,
                const typename JacobianCell<dimC>::cell*& cld);
#endif
 
 
 
 
} //namespace concepts
 
#endif // formulatfrmEformula.hh