hp2D/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 hp2D {

  using concepts::Real;

  // forward declaration
  template<typename F>
  class Quad;

  // ************************************************************ TraceSpace **

  class TraceSpace : public concepts::SpaceOnCells<Real> {
  public:
    typedef concepts::Scan<hp1D::Element<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);

    TraceSpace(concepts::SpaceOnCells<Real>& spc,
               const concepts::CellConditions* cc = 0,
               enum traceTypes type = FIRST);
    TraceSpace(const concepts::SpaceOnCells<Real>& spc,
               const concepts::CellConditions* cc = 0,
               enum traceTypes type = FIRST);

    virtual ~TraceSpace();

    virtual uint dim() const { return dim_; }

    virtual uint nelm() const { return nelm_; }

    virtual Scan* scan() const {
      return new concepts::PListScan<hp1D::Element<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::Element<Real>*, 1>* elm_;

    concepts::HashMap<concepts::Sequence<UnderlyingElement> > uelm_;

    concepts::HashMap<hp1D::Element<Real>* > edges_;

    std::auto_ptr<QuadEdgeBase> quadEdge_;

    bool warn_edgeBuildMissElem_;

    void constructType_(enum traceTypes type);

    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