hp2D/quad.hh

Go to the documentation of this file.

/* quadrilateral element for 2D hp FEM
*/

#ifndef quad_hh
#define quad_hh

#include <iostream>
#include <memory>
#include "basics/typedefs.hh"
#include "basics/outputOperator.hh"
#include "basics/vectorsMatrices.hh"
#include "geometry/cell2D.hh"
#include "geometry/integral.hh"
#include "integration/quadRule.hh"
#include "space/hpMethod.hh"
#include "element.hh"
#include "integration/karniadakis.hh"
#include "integration/laguerre.hh"

namespace hp2D {

  using concepts::Real;

  // ******************************************************** IntegrableQuad **

  class IntegrableQuad : public concepts::IntegrationCell {
  public:
    IntegrableQuad(concepts::Quad2d& cell);

    inline const concepts::QuadratureRule* integrationX() const 
    { return intX_.get(); }

    inline const concepts::QuadratureRule* integrationY() const 
    { return intY_.get(); }

    static concepts::QuadRuleFactory& rule() { return rule_; }

    inline concepts::Real2d chi(const Real x, const Real y) const {
      return cell_.chi(x, y); 
    }

    inline concepts::MapReal2d jacobian(const Real x, const Real y) const 
    {
      concepts::MapReal2d jac = cell_.jacobian(x, y);
      jac *= 0.5;         // factor comes due to integration over [-1,1]^2
      return jac; 
    }

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

    inline Real jacobianDeterminant(const Real x, const Real y) const 
    {
      return jacobian(x, y).determinant(); 
    }

    void setStrategy(const concepts::Quad2dSubdivision* strategy = 0)
      throw(concepts::StrategyChange) 
    { 
      cell_.setStrategy(strategy); 
    }

    const concepts::Quad2dSubdivision* getStrategy() const
    { 
      return cell_.getStrategy(); 
    }

    virtual bool quadraturePoint(uint i, intPoint& p,
                                 intFormType form = ZERO,
                                 bool localCoord = false) const;
  protected:
    std::auto_ptr<concepts::QuadratureRule> intX_;
    std::auto_ptr<concepts::QuadratureRule> intY_;

    concepts::Quad2d& cell_;

    // Factory for the integration rule
    static concepts::QuadRuleFactory rule_;
  };

  // ************************************************************** BaseQuad **

  template<class F = Real>
  class BaseQuad : public Element<F>, public IntegrableQuad {
  public:
    BaseQuad(concepts::Quad2d& cell, const ushort* p, 
             concepts::TColumn<F>* T0, concepts::TColumn<F>* T1);
    virtual ~BaseQuad();

    virtual const concepts::Quad& support() const { return cell_.connector(); }
    virtual concepts::Real3d vertex(uint i) const { return cell_.vertex(i); }
    virtual const concepts::Quad2d& cell() const  { return cell_; }

    virtual const concepts::ElementGraphics<F>* graphics() const = 0;

  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };

  // **************************************************** QuadShapeFunctions **

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

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

      inline const concepts::Karniadakis<1,0>* shpfctX() const {
        return shpfctX_.get(); 
      }
      inline const concepts::Karniadakis<1,0>* shpfctY() const {
        return shpfctY_.get(); }

      inline const concepts::Karniadakis<1,1>* shpfctDX() const {
        return shpfctDX_.get(); }
      inline const concepts::Karniadakis<1,1>* shpfctDY() const {
        return shpfctDY_.get(); }
    protected:
      void computeShapefunctions_(const concepts::QuadratureRule *intX,
                                  const concepts::QuadratureRule *intY);
    private:
      ushort p_[2];
      std::auto_ptr<concepts::Karniadakis<1,0> > shpfctX_, shpfctY_;
      std::auto_ptr<concepts::Karniadakis<1,1> > shpfctDX_, shpfctDY_;
  };

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

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

    virtual ~Quad();

    virtual const concepts::ElementGraphics<F>* graphics() const;
    void recomputeShapefunctions();
    void recomputeShapefunctions(const uint nq[2]);

    void edgeP(const uint i, const uint& p) {
      edges_[i] = &p;
    }
    const uint edgeP(const uint i) const { 
      conceptsAssert(i < 4, concepts::Assertion());
      conceptsAssert(edges_[i], concepts::Assertion());
      return *edges_[i];
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    static std::auto_ptr<concepts::ElementGraphics<F> > graphics_;
    const uint* edges_[4];
  };


  // ********************************************************** InfiniteQuad **

  class InfiniteQuad : public Element<Real> {
  public:
    InfiniteQuad(concepts::InfiniteQuad2d& cell, const ushort* p, 
                 concepts::TColumn<Real>* T0, concepts::TColumn<Real>* T1);

    virtual ~InfiniteQuad();

    const concepts::Connector2& support()  const { return cell_.connector(); }
    concepts::Real3d vertex(uint i)        const { return cell_.vertex(i); }
    const concepts::InfiniteQuad2d& cell() const { return cell_; }

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

    inline const concepts::QuadratureRule* integrationX() const 
    { return intX_.get(); }
    
    inline const concepts::QuadratureRule* integrationY() const 
    { return intY_.get(); }

    static void setIntegrationRuleX(enum concepts::intRule rule, bool constant,
                                    uint points);

    static void setIntegrationRuleY(bool constant, uint points,
                                    const Real border);

    static void resetIntegrationRule();

    static std::string integrationRule();

    inline const concepts::Karniadakis<1,0>* shpfctX() const {
      return shpfctX_.get(); 
    }

    inline const concepts::Karniadakis<1,1>* shpfctDX() const {
      return shpfctDX_.get(); }

    virtual const concepts::LaguerreBasis<0>* shpfctY() const = 0;

    virtual const concepts::LaguerreBasis<1>* shpfctDY() const = 0;

    virtual const concepts::ElementGraphics<Real>* graphics() const = 0;

    static Real& borderY() { return borderY_; }
  protected:
    void computeIntegrationRuleFromP_(const ushort* p);

    void computeIntegrationRule_(const uint nq[2]);

    void computeShapefunctionsX_();

    virtual std::ostream& info(std::ostream& os) const;

    ushort p_[2];

    std::auto_ptr<concepts::QuadratureRule> intX_;
    std::auto_ptr<concepts::QuadratureRule> intY_;
  private:
    concepts::InfiniteQuad2d& cell_;

    static enum concepts::intRule integrationTypeX_;
    static enum concepts::intRule integrationTypeY_;
    static uint constNumerOfPointsX_, constNumerOfPointsY_;
    static uint addNumberOfPointsX_, addNumberOfPointsY_;
    static bool useConstantNumberOfPointsX_, useConstantNumberOfPointsY_;
    static Real borderY_;

    std::auto_ptr<concepts::Karniadakis<1,0> > shpfctX_;
    std::auto_ptr<concepts::Karniadakis<1,1> > shpfctDX_;
  };

  // ************************************************** InfiniteLaguerreQuad **

  class InfiniteLaguerreQuad : public InfiniteQuad {
  public:
    InfiniteLaguerreQuad(concepts::InfiniteQuad2d& cell, const ushort* p, 
                         concepts::TColumn<Real>* T0,
                         concepts::TColumn<Real>* T1);
    virtual ~InfiniteLaguerreQuad();

    inline const concepts::LaguerreBasis<0>* shpfctY() const {
      return shpfctY_.get(); }

    inline const concepts::LaguerreBasis<1>* shpfctDY() const {
      return shpfctDY_.get(); }

    void recomputeShapefunctions();

    virtual const concepts::ElementGraphics<Real>* graphics() const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    static std::auto_ptr<concepts::ElementGraphics<Real> > graphics_;

    std::auto_ptr<concepts::LaguerreBasis<0> > shpfctY_;
    std::auto_ptr<concepts::LaguerreBasis<1> > shpfctDY_;

    void computeShapefunctionsY_();
  };


} // namespace hp2D

#endif // quad_hh