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 hp2D
 
#endif // hp2dtraces_hh