spacePreBuilder.hh

/* space pre builder classes for 2D spaces
*/
 
#ifndef hpSpacePreBuilder2D_hh
#define hpSpacePreBuilder2D_hh
 
#include <set>
#include <queue>
#include "basics/outputOperator.hh"
#include "geometry/cell2D.hh"
#include "space/spacePreBuilder.hh"
#include "space/space.hh"
#include "space/hpMethod.hh"          // for AdaptiveControlP
#include "toolbox/stiffArray.hh"
#include "toolbox/multiArray.hh"
#include "toolbox/hashMap.hh"
#include "toolbox/sequence.hh"
#include "toolbox/set.hh"
#include "hp2D/edgeMinMax.hh"
 
// If flag is set to 1, edgelist_ and vertexlist_ collect keys, not
// local indices. That is easier to debug, but lightly slower.
#define KeyLists 1
 
namespace hp2D {
 
  // *********************************************************** NotValidDof **
 
  //   class NotValidDof : public concepts::ExceptionBase {};
 
  // ******************************************************* SpacePreBuilder **
 
  class SpacePreBuilder : public concepts::SpacePreBuilder {
  public:
    SpacePreBuilder(concepts::Mesh2& msh) : msh_(msh) {}
    virtual ~SpacePreBuilder() {}
 
    virtual concepts::Mesh2& mesh() const { return msh_; }
    virtual bool vtxDof(const concepts::Connector2& cntr, uint i) const = 0;
    virtual ushort edgeDof(const concepts::Connector2& cntr, uint i) const = 0;
    virtual const ushort* innerDof(const concepts::Connector2& cntr) const = 0;
    virtual const ushort* pMax(const concepts::Connector2& cntr,
             ushort* p = 0) const = 0;
    virtual ushort edgeP(const concepts::Connector2& cntr, uint i) const = 0;
 
    //not necessary declaration, but usefull for pyplusplus
    virtual concepts::IndexRange&
      setIndex(uint& firstIdx, uint noIdx, uint dim,
               const concepts::Connector& cntr, uint i = 0,
         uint spcNo = 0) = 0;
 
    //not necessary declaration, but usefull for pyplusplus
    virtual concepts::Set<concepts::IndexRange> indices
      (uint dim, const concepts::Connector& cntr, uint spcNo = 0) const = 0;
 
  protected:
    concepts::Mesh2& msh_;
  };
  
  typedef concepts::Adaptivity<concepts::Connector,
                               concepts::AdaptiveAdjustP<2> > Adaptivity;
 
  // *********************************************************** DegreeDim **
 
  template<typename F>
  struct DegreeDim {};
 
  template<>
  struct DegreeDim<concepts::Quad>         { static const ushort dim = 2; };
 
  template<>
  struct DegreeDim<concepts::InfiniteQuad> { static const ushort dim = 2; };
 
  template<>
  struct DegreeDim<concepts::Triangle>     { static const ushort dim = 1; };
  
 
  // **************************************************************** hpFull **
 
  class hpFull : public SpacePreBuilder, public Adaptivity {
  public:
    hpFull(concepts::Mesh2& msh, uint l = 0, uint p = 1, bool build = false);
    hpFull(const hpFull& preBuild);
    virtual ~hpFull() {}
 
    void setPolyEdgeRule(const PolyEdgeRule& rule);
 
    virtual void adjust(const concepts::Connector& cell,
      const concepts::AdaptiveAdjustP<2>& a);
    void adjust(const uint nrCell, const concepts::AdaptiveAdjustP<2>& a);
    void adjust(const concepts::AdaptiveAdjustP<2>& a);
 
    uint available() const { return rebuild_ ? 0 : build_; }
    virtual void rebuildMesh();
    virtual void rebuildDof();
    virtual bool vtxDof(const concepts::Connector2& cntr, uint i) const;
    virtual ushort edgeDof(const concepts::Connector2& cntr, uint i) const;
    virtual const ushort* innerDof(const concepts::Connector2& cntr) const;
    virtual const ushort* pMax(const concepts::Connector2& cntr,
             ushort* p = 0) const;
    virtual ushort edgeP(const concepts::Connector2& cntr, uint i) const;
    const concepts::IndexRange& index(uint dim,
              const concepts::Connector2& cntr,
              uint i = 0, uint spcNo = 0) const;
    virtual concepts::IndexRange&
      setIndex(uint& firstIdx, uint noIdx, uint dim,
         const concepts::Connector& cntr, uint i = 0,
         uint spcNo = 0);
 
    virtual concepts::Set<concepts::IndexRange> indices
      (uint dim, const concepts::Connector& cntr, uint spcNo = 0) const;
 
    virtual void setP(uint p, const concepts::Set<uint>& donts);
 
    void clearAllIndices(uint spcNo = 0);
 
    uint ncells() { return cells_.size(); }
 
    bool isHangEdge(const concepts::Connector1& cntr) const;
    void storeMatlab(const std::string name);
 
    concepts::Set<concepts::Cell2*> cells() const { return cells_; }
    std::set<concepts::Cell2*> allCells() const { return allCells_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    bool rebuild_;
    uint build_;
    std::unique_ptr<const PolyEdgeRule> polyEdgeRule_;
 
    concepts::HashMap<concepts::AdaptiveAdjustP<2> > adj_;
 
    // mapping from space number to an index range
    typedef std::map<uint, concepts::IndexRange> IndexType;
    // index range per space number for all 'member of space' vertices dof
    typedef concepts::AdaptiveControl<IndexType> CtrlType_0;
    // index range per space number for all 'member of space' edge dof
    typedef concepts::AdaptiveControlP<1, IndexType> CtrlType_1;
    // index range per space number for all inner dof
    typedef concepts::AdaptiveControlP<2, IndexType> CtrlType_2;
 
    concepts::HashMap<CtrlType_0> ctrl0_;
    concepts::HashMap<CtrlType_1> ctrl1_;
    concepts::HashMap<CtrlType_2> ctrl2_;
 
    concepts::HashMap<ushort> polyDep_;
    concepts::HashMap<ushort> polyHang_;
    mutable concepts::HashMap<concepts::StiffArray<2,ushort> > pmax_;
    mutable concepts::HashMap<concepts::StiffArray<2,ushort> > pminInner_;
 
    concepts::Set<concepts::Cell2*> cells_;
    std::set<concepts::Cell2*> allCells_;
    concepts::HashMap<std::set<uint> > vertexList_;
    concepts::HashMap<std::set<uint> > edgeList_;
    concepts::Set<uint> hangEdges_;
 
    bool isMember_(const concepts::Connector2& cntr) const;
    bool isMember_(const concepts::Connector1& cntr) const;
 
    const ushort* pMinInner_(const concepts::Quad& quad);
 
    void refine_(concepts::Quad2d& cell, concepts::Level<2> L, short* P,
                 concepts::HashMap<CtrlType_2>::iterator& i2);
    void refine_(concepts::InfiniteQuad2d& cell, concepts::Level<2> L,
                 short* P, concepts::HashMap<CtrlType_2>::iterator& i2);
    void refine_(concepts::Cell2& cell, ushort* P,
                 concepts::HashMap<CtrlType_2>::iterator& i2);
 
    void memberEdges_(const concepts::Connector2& cntr,
          std::set<const concepts::Connector1*>& edges);
 
    ushort takeOutchildrenEdges_
    (const concepts::Connector1& edge, const concepts::Connector1& cntr1,
     std::set<const concepts::Connector1*>& edges,
     std::set<const concepts::Connector0*>& vertices);
 
    template<class F>
    bool getEdgeSupport_(F cntr);
 
    class Cause : public std::pair<const concepts::Connector2*,
                                   const concepts::Connector1*>,
                  public concepts::OutputOperator
    {
    public:
      Cause(const concepts::Connector2* cell,
            const concepts::Connector1* edge) : 
        std::pair<const concepts::Connector2*,
                  const concepts::Connector1*>(cell, edge) {}
    protected:
      virtual std::ostream& info(std::ostream& os) const {
        return os << concepts::typeOf(*this)<<"(causing cell = " << *this->first
                  << ", appropiate edge = " << *this->second << ")";
      }
    };
 
    template<class F>
    bool getPretendedVertexSupport_
    (F cntr,
     concepts::HashMap<std::queue<Cause> >& nonLocal);
 
    void getRootCells_
    (const concepts::Connector2& cntr, const std::set<uint>& vertices,
     std::map<uint, concepts::Sequence<const concepts::Connector2*> >& rootCells);
 
    void getNonLocalVertexSupport_
    (concepts::HashMap<std::queue<Cause> >& nonLocal,
     const std::map<uint, concepts::Sequence<const concepts::Connector2*> >&
     rootCells);
 
    void getEdgeMappings_
    (uint key, uint cKey,
     const concepts::Sequence<const concepts::Connector2*>& root,
     concepts::MultiArray<2, const concepts::Connector2*>& edgesToCell, 
     std::map<uint,
       concepts::Sequence<const concepts::Connector2*> >& edgeToCell,
     std::map<uint, const concepts::Connector1*>& fatherEdges,
     std::map<uint, uint>& mosEdges) const;
 
    bool eliminateNonCauses_(std::queue<Cause>& causes);
 
    void dependentEdges_(const concepts::Quad& quad);
    void polyDependentEdges_(const concepts::Quad& quad,
           const concepts::Edge& edge,
           const concepts::Edge& oppEdge, const uint p);
 
    ushort* pMax_(const concepts::Quad& cntr) const;
    ushort* pMax_(const concepts::InfiniteQuad& cntr) const;
 
    void setPmax_(const concepts::Connector2& cntr, ushort* Pmax);
 
    void polyHangingEdges_(const concepts::Connector1& edge, const uint& p);
 
    template<class F>
    void findControl_(F*& ctrl, const concepts::Connector2& cntr, uint i);
 
          IndexType& findIndex_(uint dim, const concepts::Connector2& cntr, 
                                uint i);
    const IndexType& findIndex_(uint dim, const concepts::Connector2& cntr,
                                uint i) const;
  };
 
} // namespace hp2D
 
#endif // hpSpacePreBuilder2D_hh