traces.hh

Go to the documentation of this file.

 
#ifndef hp2dtraces_hh
#define hp2dtraces_hh
 
#include "toolbox/sequence.hh"
#include "formula/boundary.hh"
#include "formula/exceptions.hh"
#include "geometry/boundaryConditions.hh"
#include "geometry/cellConditions.hh"
#include "space/space.hh"
#include "hp1D/element.hh"
#include "hp2D/element.hh"
#include "hp2D/edge.hh"
 
namespace concepts {
 
  // forward declaration
  class EdgeNormalVectorRule;
}
 
namespace hp2D {
 
  using concepts::Real;
 
  // forward declaration
  template<typename F>
  class Quad;
 
 
  // ************************************************************ TraceSpace **
 
  class TraceSpace : public concepts::SpaceOnCells<Real> {
  public:
    typedef concepts::Scan<hp1D::BaseElement<Real> > Scan;
    typedef concepts::ElementAndFacette<hp2D::Element<Real> > UnderlyingElement;
 
    enum traceTypes { FIRST, MEAN, JUMP };
 
    TraceSpace(const concepts::SpaceOnCells<Real>& spc,
               const concepts::Set<uint> edgeAttr,
               enum traceTypes type = FIRST,
               const concepts::EdgeNormalVectorRule&
               normalVectorRule = concepts::EdgeNormalVectorRule());
 
    TraceSpace(concepts::SpaceOnCells<Real>& spc,
               const concepts::CellConditions* cc = 0,
               enum traceTypes type = FIRST,
               const concepts::EdgeNormalVectorRule
               normalVectorRule = concepts::EdgeNormalVectorRule());
    TraceSpace(const concepts::SpaceOnCells<Real>& spc,
               const concepts::CellConditions* cc = 0,
               enum traceTypes type = FIRST,
               const concepts::EdgeNormalVectorRule
               normalVectorRule = concepts::EdgeNormalVectorRule());
 
    virtual ~TraceSpace();
 
    virtual uint dim() const { return dim_; }
 
    virtual uint nelm() const { return nelm_; }
 
    virtual Scan* scan() const
    {
      return new concepts::PListScan<hp1D::BaseElement<Real> >(*elm_);
    }
 
    const concepts::Sequence<UnderlyingElement>
    uelm(const concepts::Edge edge) const
    {
      concepts::HashMap<concepts::Sequence<UnderlyingElement> >::
        const_iterator i = uelm_.find(edge.key());
      if (i == uelm_.end())
        return concepts::Sequence<UnderlyingElement>();
      return i->second;
    }
 
    const concepts::HashMap<concepts::Sequence<UnderlyingElement> > 
    uelm() const
    {
      return uelm_;
    }
 
    virtual void recomputeShapefunctions();
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const uint dim_;
 
    uint nelm_;
 
    concepts::Joiner<hp1D::BaseElement<Real>*, 1>* elm_;
 
    concepts::HashMap<concepts::Sequence<UnderlyingElement> > uelm_;
 
    concepts::HashMap<hp1D::Element<Real>* > edges_;
    
    std::unique_ptr<QuadEdgeBase> quadEdge_;
 
    bool warn_edgeBuildMissElem_;
 
    void constructType_(enum traceTypes type,
                        const concepts::EdgeNormalVectorRule& normalVectorRule);
 
    template<class F>
    bool build_(const Quad<Real>* elm, F condition);
 
    template<class F>
    void test_(const concepts::ElementWithCell<Real>& elm, F condition);
  };
 
} // namespace hp2D
 
namespace concepts {
 
  // *********************************************** DirichletElementFormula **
 
  template<typename F, typename G = typename Realtype<F>::type>
  class DirichletElementFormula : public ElementFormula<F,G> {
  public:
    DirichletElementFormula(const BoundaryConditions bc) : bc_(bc) {}
 
    virtual F operator() (const ElementWithCell<G>& elm, const Real p,
        const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<G>& elm, const Real2d& p,
        const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<G>& elm, const Real3d& p,
        const Real t = 0.0) const;
    virtual DirichletElementFormula<F,G>* clone() const {
      return new DirichletElementFormula<F,G>(bc_);
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const BoundaryConditions bc_;
 
    template<typename P>
    F compute_(const Element<G>& elm, const P& p, const Real t) const {
      const ElementWithCell<G>* e = 
        dynamic_cast<const ElementWithCell<G>*>(&elm);
      conceptsAssert(e, Assertion());
      const Cell& cell = e->cell();
      const Boundary b(bc_(cell.connector().attrib()));
      if (b.type() != Boundary::DIRICHLET) 
        throw(ElementNotInDomainOfFormula(elm, *this));
 
      F val = 0.0;
      if (!compute_(dynamic_cast<const Edge2d*>(&cell), b, p, t, val))
        throw conceptsException(concepts::MissingFeature
        ("element not supported"));
      
      return val;
    }
 
    bool compute_(const Edge2d* cell, const Boundary b,
      const Real p, const Real t, F& val) const {
      if (cell) {
        val = b(cell->chi(p), t);
        return true;
      }
      return false;
    }
 
    template<typename P>
    bool compute_(const Edge2d* cell, const Boundary b,
                  const P& p, const Real t, F& val) const {
      throw conceptsException(concepts::MissingFeature("p should be a real"));
      return false;
    }
  };
 
} // namespace concepts
 
 
namespace hp2D {
 
  namespace l2 {
 
    using concepts::Real;
 
    // forward declaration
    template<typename F>
    class Quad;
 
    // ************************************************************ TraceSpace **
 
    class TraceSpace : public concepts::SpaceOnCells<Real> {
    public:
      typedef concepts::Scan<hp1D::BaseElement<Real> > Scan;
      typedef concepts::ElementAndFacette<hp2D::Element<Real> > UnderlyingElement;
 
      enum traceTypes { FIRST, MEAN, JUMP };
 
      TraceSpace(const concepts::SpaceOnCells<Real>& spc,
                 const concepts::Set<uint> edgeAttr,
                 enum traceTypes type = FIRST,
                 const concepts::EdgeNormalVectorRule&
                 normalVectorRule = concepts::EdgeNormalVectorRule());
 
      TraceSpace(concepts::SpaceOnCells<Real>& spc,
                 const concepts::CellConditions* cc = 0,
                 enum traceTypes type = FIRST,
                 const concepts::EdgeNormalVectorRule
                 normalVectorRule = concepts::EdgeNormalVectorRule());
      TraceSpace(const concepts::SpaceOnCells<Real>& spc,
                 const concepts::CellConditions* cc = 0,
                 enum traceTypes type = FIRST,
                 const concepts::EdgeNormalVectorRule
                 normalVectorRule = concepts::EdgeNormalVectorRule());
 
      virtual ~TraceSpace();
 
      virtual uint dim() const { return dim_; }
 
      virtual uint nelm() const { return nelm_; }
 
      virtual Scan* scan() const
      {
        return new concepts::PListScan<hp1D::BaseElement<Real> >(*elm_);
      }
 
      const concepts::Sequence<UnderlyingElement>
      uelm(const concepts::Edge edge) const
      {
        concepts::HashMap<concepts::Sequence<UnderlyingElement> >::
          const_iterator i = uelm_.find(edge.key());
        if (i == uelm_.end())
          return concepts::Sequence<UnderlyingElement>();
        return i->second;
      }
 
      const concepts::HashMap<concepts::Sequence<UnderlyingElement> >
      uelm() const
      {
        return uelm_;
      }
 
      virtual void recomputeShapefunctions();
    protected:
      virtual std::ostream& info(std::ostream& os) const;
    private:
      const uint dim_;
 
      uint nelm_;
 
      concepts::Joiner<hp1D::BaseElement<Real>*, 1>* elm_;
 
      concepts::HashMap<concepts::Sequence<UnderlyingElement> > uelm_;
 
      concepts::HashMap<hp1D::LegendreElement<Real>* > edges_;
 
      std::unique_ptr<QuadEdgeBase> quadEdge_;
 
      bool warn_edgeBuildMissElem_;
 
      void constructType_(enum traceTypes type,
                          const concepts::EdgeNormalVectorRule& normalVectorRule);
 
      template<class F>
      bool build_(const Quad<Real>* elm, F condition);
 
      template<class F>
      void test_(const concepts::ElementWithCell<Real>& elm, F condition);
    };
 
  } // namespace l2
 
 
} // namespace hp2D
 
#endif // hp2dtraces_hh