hp2Dedge/quad.hh

Go to the documentation of this file.

/* quadrilateral element for 2D hp edgeFEM
 */

#ifndef quadEdge_hh
#define quadEdge_hh

#include "basics/typedefs.hh"
#include "basics/vectorsMatrices.hh"
#include "integration/quadrature.hh"
#include "hp2D/quad.hh"
#include "integration/karniadakis.hh"

namespace hp2Dedge {

  using concepts::Real;

  // ***************************************************** KarniadakisDeriv2 **

  class KarniadakisDeriv2 : public concepts::ShapeFunction1D<Real> {
  public:
    KarniadakisDeriv2(const int P, const Real* xP, const int NxP,
          const bool cache = true);

    KarniadakisDeriv2(const KarniadakisDeriv2& arg);

    ~KarniadakisDeriv2();
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    concepts::Karniadakis<1,1>* tmp_;
  };

  // ***************************************************** QuadEdgeFunctions **

  class QuadEdgeFunctions {
  public:
    QuadEdgeFunctions(ushort p, const concepts::QuadratureRule *intX,
          const concepts::QuadratureRule *intY);
    QuadEdgeFunctions(ushort* p, const concepts::QuadratureRule *intX,
          const concepts::QuadratureRule *intY);
    virtual ~QuadEdgeFunctions();

    inline const ushort* p() const { return p_; }

    /* Returns the tangential shape functions in x direction
       that are multiples of Legendre Polynoms of order 0..p[0]
    */
    inline const KarniadakisDeriv2* shpfctX_t() const {
      return shpfctX_t_.get(); }
    /* Returns the tangential shape functions in y direction
       that are multiples of Legendre Polynoms of order 0..p[1]
    */
    inline const KarniadakisDeriv2* shpfctY_t() const {
      return shpfctY_t_.get(); }
    /* Returns the normal shape functions in x direction
       that are multiples of Integrated Legendre Polynoms of order 1..p[0]
       order 0 is not included, but not used
     */
    inline const concepts::Karniadakis<1,0>* shpfctX_n() const {
      return shpfctX_n_.get(); }
    /* Returns the shape functions in y direction
       that are multiples of Integrated Legendre Polynoms of order 1..p[1]
       order 0 is not included, but not used
     */
    inline const concepts::Karniadakis<1,0>* shpfctY_n() const {
      return shpfctY_n_.get(); }

    /* Returns the derivatives of the normal shape functions in x direction
       that are multiples of Legendre Polynoms of order 0..p[0]
     */
    inline const concepts::Karniadakis<1,1>* shpfctDX_n() const {
      return shpfctDX_n_.get(); }
    /* Returns the derivatives of the  normal shape functions in y direction
       that are multiples of Legendre Polynoms of order 0..p[1]       
     */
    inline const concepts::Karniadakis<1,1>* shpfctDY_n() const {
      return shpfctDY_n_.get(); }
  protected:
    void computeShapefunctions_(const concepts::QuadratureRule *intX,
        const concepts::QuadratureRule *intY);
  private:
    ushort p_[2];
    std::auto_ptr<KarniadakisDeriv2> shpfctX_t_, shpfctY_t_;
    std::auto_ptr<concepts::Karniadakis<1,0> > shpfctX_n_, shpfctY_n_;
    std::auto_ptr<concepts::Karniadakis<1,1> > shpfctDX_t_, shpfctDY_t_,
      shpfctDX_n_, shpfctDY_n_;
  };

  // ****************************************************************** Quad **

  template<class F = Real>
  class Quad : public hp2D::BaseQuad<F>, public QuadEdgeFunctions {
  public:
    Quad(concepts::Quad2d& cell, ushort* p, 
   concepts::TColumn<F>* T0, concepts::TColumn<F>* T1);

    virtual const concepts::ElementGraphics<F>* graphics() const;
    void recomputeShapefunctions();
    void recomputeShapefunctions(const uint nq[2]);
    virtual ~Quad();
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    static std::auto_ptr<concepts::ElementGraphics<F> > graphics_;
  };

  // ****************************************************************** Edge **

  template<class F = Real>
  class Edge : public concepts::Element<F> {
  public:
    Edge(const Quad<F>& elm, const ushort k);

    virtual ~Edge();

    inline concepts::Real2d vertex(uint i) const {
      conceptsAssert(i < 2, concepts::Assertion());
      return i ? chi(1.0) : chi(0.0);
    }

    virtual const concepts::Edge& support() const { 
      return *elm_.cell().connector().edge(k_); }

    inline const ushort edge() const { return k_;}

    inline const Quad<F>& elm() const { return elm_;}

    inline const ushort direction() const { return l_;}

    inline const Real sign() const { return k_%3==0 ? -1.0 : 1.0;}

    /* Returns the shape functions
       that are multiples of Legendre Polynoms of order 0..p
    */
    inline const KarniadakisDeriv2* shpfct() const {
      return shpfct_; }
    inline const concepts::QuadratureRule* integration() const {
      return int_; }

    inline concepts::Real2d localCoords(const Real t) const {
      concepts::Real2d x;
      x[l_] = x_;
      x[(l_+1)%2] = t;
      return x;
    }

    inline concepts::Real2d chi(const Real t) const {
      concepts::Real2d x = localCoords(t);
      return elm_.chi(x[0], x[1]);
    }

    inline concepts::MapReal2d jacobian(const Real t) const {
      concepts::Real2d x = localCoords(t);
      return elm_.jacobian(x[0],x[1]);
    }

    inline concepts::MapReal2d jacobianInverse(const Real t) const {
      return jacobian(t).inverse(); }

    inline Real jacobianDeterminant(const Real t) const {
      return jacobian(t).determinant(); }

    inline Real diffElement(const Real t) const {
      // factor comes due to integration over [-1,1], but differential element
      // was defined for reference element [0,1]
      return elm_.cell().lineElement(t, k_) * 0.5;
    }

    virtual const concepts::TMatrixBase<F>& T() const { return *T_;};
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const Quad<F>& elm_;
    const ushort k_;
    const ushort l_;
    concepts::Real x_;
    const KarniadakisDeriv2* shpfct_;
    const concepts::QuadratureRule *int_;
    const concepts::TMatrixBase<F>* T_;
  };

} // namespace hp2Dedge

#endif // quadEdge_hh