quad.hh

Go to the documentation of this file.

 
#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::QuadNd& cell);
 
    inline concepts::Real2d chi(const Real x, const Real y) const {
      return cell_.elemMap( concepts::Real2d(x,y) ); 
    }
 
    inline concepts::MapReal2d jacobian(const Real x, const Real y) const 
    {
      concepts::Quad2d* cell = dynamic_cast<concepts::Quad2d*>(&cell_);
      conceptsAssert(cell, concepts::Assertion());
      concepts::MapReal2d jac = cell->jacobian(x, y);
 
      conceptsAssert(intRule_.get(), concepts::Assertion());
      if(intRule_.get()->domain())
        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(); 
    }
 
    // Computes the second partial derivatives of the inverse of the Mapping.
    inline concepts::MapReal2d inverseLaplace(const Real x, const Real y) const
    {
     // TODO: move up the inverseLaplace to QuadNd and then
     // remove this assert and dynamic cast again
     const concepts::Quad2d *cell = dynamic_cast<concepts::Quad2d*>(&cell_);
     conceptsAssert(cell, concepts::Assertion());
 
     concepts::MapReal2d invLap = cell->inverseLaplace(x,y);
 
     conceptsAssert(intRule_.get(), concepts::Assertion());
     if(intRule_.get()->domain())
       invLap *= 2;
 
     return invLap;
    }
 
 
//    /// Computes the determinant of the Jacobian
//    inline Real jacobianDeterminant(const Real x, const Real y) const 
//    {
//      return jacobian(x, y).determinant(); 
//    }
    
    inline Real gramDeterminantRoot(const Real x, const Real y) const
    {
      Real gramian = cell_.gramDeterminantRoot(x,y);
 
      conceptsAssert(intRule_.get(), concepts::Assertion());
      //if integration is not on [0,1]^2, we need to scale
      if(intRule_.get()->domain())
        gramian *= 0.25;   // factor due to the integration over [-1,1]^2
      return gramian;    
    }
 
    void setStrategy(const concepts::Quad2dSubdivision* strategy = 0)
      throw(concepts::StrategyChange) 
    { 
      // TODO: move up the subdivision strategy to QuadNd and then
      // remove this assert and dynamic cast again
      concepts::Quad2d *cell = dynamic_cast<concepts::Quad2d*>(&cell_);
      conceptsAssert(cell, concepts::Assertion());
      cell->setStrategy(strategy); 
    }
 
    const concepts::Quad2dSubdivision* getStrategy() const
    { 
      // TODO: move up the subdivision strategy to QuadNd and then
      // remove this assert and dynamic cast again
      const concepts::Quad2d *cell = dynamic_cast<concepts::Quad2d*>(&cell_);
      conceptsAssert(cell, concepts::Assertion());
      return cell->getStrategy(); 
    }
 
    virtual bool quadraturePoint(uint i, intPoint& p, intFormType form = ZERO,
                                 bool localCoord = false) const;
 
 
    // Returns the integration rule
    const concepts::QuadratureRule2d* const intRule() const { return intRule_.get(); }
 
    static std::unique_ptr<concepts::QuadRuleFactoryTensor2d>& factory() { return factory_; }
 
    static concepts::QuadRuleFactoryTensor2d* factory_rp() { return factory_.get(); }
 
  protected:
    //std::unique_ptr<concepts::QuadratureRule1d> intX_;
   // std::unique_ptr<concepts::QuadratureRule1d> intY_;
    // Factory for the integration rule
    //static concepts::QuadRuleFactory rule_;
 
  // integration rule
    std::unique_ptr<concepts::QuadratureRule2d> intRule_;
 
    concepts::QuadNd& cell_;
 
  
    // factory for the integration rule, i.e tensored structural
    static std::unique_ptr<concepts::QuadRuleFactoryTensor2d> factory_;
  };
 
 
  void setQuadrature(enum concepts::intRule rule, uint noP);
 
  // ************************************************************** BaseQuad **
 
  template<class F = Real>
  class BaseQuad : public Element<F>, public IntegrableQuad {
  public:
    BaseQuad(concepts::QuadNd& 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::QuadNd& 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::QuadratureRule2d *intRule);
      QuadShapeFunctions(const ushort* p, const concepts::QuadratureRule2d *intRule);
      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::QuadratureRule1d *intX,
  //                                const concepts::QuadratureRule1d *intY);
 
     void computeShapefunctions_(const concepts::QuadratureRule2d *intRule);
 
    private:
      ushort p_[2];
      std::unique_ptr<concepts::Karniadakis<1,0> > shpfctX_, shpfctY_;
      std::unique_ptr<concepts::Karniadakis<1,1> > shpfctDX_, shpfctDY_;
  };
 
  // ****************************************************************** Quad **
 
  template<class F = Real>
  class Quad : public BaseQuad<F>, public QuadShapeFunctions {
  public:
    Quad(concepts::QuadNd& 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::unique_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::QuadratureRule1d* integrationX() const
    { return intX_.get(); }
    
    inline const concepts::QuadratureRule1d* 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::unique_ptr<concepts::QuadratureRule1d> intX_;
    std::unique_ptr<concepts::QuadratureRule1d> 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::unique_ptr<concepts::Karniadakis<1,0> > shpfctX_;
    std::unique_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::unique_ptr<concepts::ElementGraphics<Real> > graphics_;
 
    std::unique_ptr<concepts::LaguerreBasis<0> > shpfctY_;
    std::unique_ptr<concepts::LaguerreBasis<1> > shpfctDY_;
 
    void computeShapefunctionsY_();
  };
 
  // ****************************************************** ArrayQuadWeights **
 
  class ArrayQuadWeights : public concepts::Array<Real> {
  public:
    ArrayQuadWeights(const Quad<>& quad);
  };
 
  // ***************************************************** KarniadakisDeriv2 **
 
  class KarniadakisDeriv2 : public concepts::ShapeFunction1D<Real> {
  public:
    KarniadakisDeriv2(const int P, const Real* xP, const int NxP,
                      const bool cache = true);
    
    // commented, because there is not an assignement operator= . See Stroustrup "The Meaning of Copy" and "Explicit Defaults" chapters.
    //KarniadakisDeriv2(const KarniadakisDeriv2& arg);
    
    ~KarniadakisDeriv2();
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    concepts::Karniadakis<1,1>* tmp_;
  };
 
  namespace l2 {
 
    // **************************************************** QuadShapeFunctions **
 
    class QuadShapeFunctions {
    public:
      QuadShapeFunctions(const ushort p, const concepts::QuadratureRule2d *intRule);
      QuadShapeFunctions(const ushort* p, const concepts::QuadratureRule2d *intRule);
      virtual ~QuadShapeFunctions();
 
      inline const ushort* p() const { return p_; }
 
      inline const KarniadakisDeriv2* shpfctX() const {
        return shpfctX_.get(); 
      }
      inline const KarniadakisDeriv2* shpfctY() const {
        return shpfctY_.get(); }
    protected:
      void computeShapefunctions_(const concepts::QuadratureRule2d *intRule);
    private:
      ushort p_[2];
      std::unique_ptr<KarniadakisDeriv2 > shpfctX_, shpfctY_;
    };
 
    // ****************************************************************** Quad **
 
    template<class F = Real>
    class Quad : public BaseQuad<F>, public QuadShapeFunctions {
    public:
      Quad(concepts::QuadNd& 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]);
    protected:
      virtual std::ostream& info(std::ostream& os) const;
    private:
      static std::unique_ptr<concepts::ElementGraphics<F> > graphics_;
    };
 
  } // namespace hp2D::l2
 
} // namespace hp2D
 
#endif // quad_hh