// SpaFEDTe, a Template based C++ library for creating 
// Discontinuous Finite Element Spaces,
// Copyright (C) 2012 Lorenzo Alessio Botti

/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Lesser General Public */ 
/* License as published by the Free Software Foundation either */ 
/* version 3.0 of the License, or (at your option) any later version. */

/* This software is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU */
/* Lesser General Public License for more details. */

/* You should have received a copy of the GNU Lesser General Public */
/* License along with this software; if not, a copy of the full */
/* GNU Lesser General Public License can be found at */
/* http://www.gnu.org/licenses/ */

#ifndef __jacobi_gauss_finite_element_quadrature_h__
#define __jacobi_gauss_finite_element_quadrature_h__

#include "FiniteElementQuadrature.hpp"
#include <vector>
#include "JacobiPoly.hpp"
#include "CollapsedMapping.hpp"
#include "Gamma.hpp"
#include "Tuple.hpp"
#include "EnumElemType.hpp"
#include "Quadrature.hpp"
#include <stdexcept>

class JacobiGaussFiniteElementQuadrature 
{
public:
  typedef Quadrature::PointVector PointVector;
  typedef Quadrature::RealVector RealVector;

  JacobiGaussFiniteElementQuadrature() 
    {
      ;
    }
  
  ~JacobiGaussFiniteElementQuadrature() {}
  
public:
  static void initQRule(Quadrature& quadrature);
  
  static unsigned int nPoints(const unsigned int dim, const unsigned int q_Order, const ElemType itemType);

private:
  static void tetQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight);
  
  static void prismQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight);
  
  static void pyramidQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight);

  static void hexQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight);

  static void triQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight);

  static void quadQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight);
  
  static void lineQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight);
  
  static void vertexQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight);
 
  static void tensorProduct2D(const RealVector& qPoints, const RealVector& qWeights, PointVector& q_Qpoint, RealVector& q_Weight);
  
  static void tensorProduct2D(const RealVector& qPoints1, const RealVector& qPoints2,
                              const RealVector& qWeights1, const RealVector& qWeights2, PointVector& q_Qpoint, RealVector& q_Weight);

  static void tensorProduct3D(const RealVector& qPoints, const RealVector& qWeights, PointVector& q_Qpoint, RealVector& q_Weight);
  
  static void tensorProduct3D(const RealVector& qPoints1, const RealVector& qPoints2, const RealVector& qPoints3,
                              const RealVector& qWeights1, const RealVector& qWeights2, const RealVector& qWeights3, 
                              PointVector& q_Qpoint, RealVector& q_Weight);

  static void jacobiGaussQweights(double a, double b, unsigned int n, RealVector& qWeights, const RealVector& qPoints);

  static void jacobiGaussQweights(double a, unsigned int n, RealVector& qWeights, const RealVector& qPoints);

};

#endif 
