// 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/ */

#include "JacobiGaussFiniteElementQuadrature.hpp"

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

void 
JacobiGaussFiniteElementQuadrature::
tetQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight)
{
  RealVector qPointsN1;
  RealVector qPointsN2;
  RealVector qPointsN3;
  RealVector qWeights1;
  RealVector qWeights2;
  RealVector qWeights3;
  JacobiPoly::jacobiGaussQpoints(0.0, 0.0, order, qPointsN1); 
  JacobiPoly::jacobiGaussQpoints(1.0, 0.0, order, qPointsN2); 
  JacobiPoly::jacobiGaussQpoints(2.0, 0.0, order, qPointsN3); 
  jacobiGaussQweights(0.0, order, qWeights1, qPointsN1);
  jacobiGaussQweights(1.0, order, qWeights2, qPointsN2);
  jacobiGaussQweights(2.0, order, qWeights3, qPointsN3);
  tensorProduct3D(qPointsN1, qPointsN2, qPointsN3, qWeights1, qWeights2, qWeights3, q_Qpoint, q_Weight);
  CollapsedMapping::mapTetNE(q_Qpoint);
}

void 
JacobiGaussFiniteElementQuadrature::
prismQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight)
{
  RealVector qPointsN1;
  RealVector qPointsN2;
  RealVector qWeights1;
  RealVector qWeights2;
  JacobiPoly::jacobiGaussQpoints(0.0, 0.0, order, qPointsN1); 
  JacobiPoly::jacobiGaussQpoints(1.0, 0.0, order, qPointsN2); 
  jacobiGaussQweights(0.0, order, qWeights1, qPointsN1);
  jacobiGaussQweights(1.0, order, qWeights2, qPointsN2);
  tensorProduct3D(qPointsN1, qPointsN2, qPointsN1, qWeights1, qWeights2, qWeights1, q_Qpoint, q_Weight);
  CollapsedMapping::mapPrismNE(q_Qpoint);
}

void 
JacobiGaussFiniteElementQuadrature::
pyramidQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight)
{
  RealVector qPointsN1;
  RealVector qPointsN2;
  RealVector qWeights1;
  RealVector qWeights2;
  JacobiPoly::jacobiGaussQpoints(0.0, 0.0, order, qPointsN1); 
  JacobiPoly::jacobiGaussQpoints(2.0, 0.0, order, qPointsN2); 
  jacobiGaussQweights(0.0, order, qWeights1, qPointsN1);
  jacobiGaussQweights(2.0, order, qWeights2, qPointsN2);
  tensorProduct3D(qPointsN1, qPointsN1, qPointsN2, qWeights1, qWeights1, qWeights2, q_Qpoint, q_Weight);
  CollapsedMapping::mapPyramidNE(q_Qpoint);
}

void 
JacobiGaussFiniteElementQuadrature::
hexQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight)
{
  RealVector qPointsE;
  RealVector qWeights;
  JacobiPoly::jacobiGaussQpoints(0.0, 0.0, order, qPointsE); 
  jacobiGaussQweights(0.0, order, qWeights, qPointsE);
  tensorProduct3D(qPointsE, qWeights, q_Qpoint, q_Weight);
}

void 
JacobiGaussFiniteElementQuadrature::
triQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight)
{
  RealVector qPointsN1;
  RealVector qPointsN2;
  RealVector qWeights1;
  RealVector qWeights2;
  JacobiPoly::jacobiGaussQpoints(0.0, 0.0, order, qPointsN1); 
  JacobiPoly::jacobiGaussQpoints(1.0, 0.0, order, qPointsN2); 
  jacobiGaussQweights(0.0, order, qWeights1, qPointsN1);
  jacobiGaussQweights(1.0, order, qWeights2, qPointsN2);
  tensorProduct2D(qPointsN1, qPointsN2, qWeights1, qWeights2, q_Qpoint, q_Weight);
  CollapsedMapping::mapTriNE(q_Qpoint);
}

void 
JacobiGaussFiniteElementQuadrature::
quadQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight)
{
  RealVector qPointsE;
  RealVector qWeights;
  JacobiPoly::jacobiGaussQpoints(0.0, 0.0, order, qPointsE); 
  jacobiGaussQweights(0.0, order, qWeights, qPointsE);
  tensorProduct2D(qPointsE, qWeights, q_Qpoint, q_Weight);
}

void 
JacobiGaussFiniteElementQuadrature::
lineQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight)
{
  RealVector qPointsE;
  RealVector qWeights;
  JacobiPoly::jacobiGaussQpoints(0.0, 0.0, order, qPointsE); 
  jacobiGaussQweights(0.0, order, qWeights, qPointsE);
  //JacobiPoly::legendreGaussQpointsWeights(order, qPointsE, qWeights); 
  const unsigned int nPoints = qPointsE.size();
  q_Qpoint.resize(nPoints,Eigen::NoChange);
  q_Weight.resize(nPoints);
  unsigned int qp = 0;
  for (unsigned int i = 0; i < nPoints; i++)
  {
    q_Qpoint(qp,0) = qPointsE(i);
    q_Qpoint(qp,1) = 0.;
    q_Qpoint(qp,2) = 0.;
    q_Weight(qp) = qWeights(i);
    qp++;
  }
}

void
JacobiGaussFiniteElementQuadrature::
vertexQuadrature(unsigned int order, PointVector& q_Qpoint, RealVector& q_Weight)
{
  RealVector qPointsE;
  RealVector qWeights;
  q_Qpoint.resize(1,Eigen::NoChange);
  q_Weight.resize(1);
  q_Qpoint(0,0) = 0.;
  q_Qpoint(0,1) = 0.;
  q_Qpoint(0,2) = 0.;
  q_Weight(0) = 1.;
}

void
JacobiGaussFiniteElementQuadrature::
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)
{
  assert( qPoints1.size() == qPoints2.size() &&  qPoints2.size() == qPoints3.size());//, 
  const unsigned int nPoints = qPoints1.size();
  q_Qpoint.resize(nPoints*nPoints*nPoints, Eigen::NoChange);
  q_Weight.resize(nPoints*nPoints*nPoints);
  unsigned int qp = 0;
  for (unsigned int k = 0; k < nPoints; k++)
  {
    for (unsigned int i = 0; i < nPoints; i++)
    {
      for (unsigned int j = 0; j < nPoints; j++)
      {
         q_Qpoint(qp,0) = qPoints1(k);
         q_Qpoint(qp,1) = qPoints2(i);
         q_Qpoint(qp,2) = qPoints3(j);
         q_Weight(qp) = qWeights1(k) * qWeights2(i) * qWeights3(j);
         qp++;
      }
    }
  }
}

void
JacobiGaussFiniteElementQuadrature::
tensorProduct3D(const RealVector& qPoints, const RealVector& qWeights, PointVector& q_Qpoint, RealVector& q_Weight)
{
  const unsigned int nPoints = qPoints.size();
  q_Qpoint.resize(nPoints*nPoints*nPoints, Eigen::NoChange);
  q_Weight.resize(nPoints*nPoints*nPoints);
  unsigned int qp = 0;
  for (unsigned int k = 0; k < nPoints; k++)
    for (unsigned int i = 0; i < nPoints; i++)
      for (unsigned int j = 0; j < nPoints; j++)
      {
         q_Qpoint(qp,0) = qPoints(k);
         q_Qpoint(qp,1) = qPoints(i);
         q_Qpoint(qp,2) = qPoints(j);
         q_Weight(qp) = qWeights(k) * qWeights(i) * qWeights(j);
         qp++;
      }
}

void
JacobiGaussFiniteElementQuadrature::
tensorProduct2D(const RealVector& qPoints1, const RealVector& qPoints2,
                const RealVector& qWeights1, const RealVector& qWeights2, PointVector& q_Qpoint, RealVector& q_Weight)
{
  assert( qPoints1.size() == qPoints2.size());//, ("Error in JacobiGauss Quadrature: vector size doesn't match in 2D tri tensor product "));
  const unsigned int nPoints = qPoints1.size();
  q_Qpoint.resize(nPoints*nPoints, Eigen::NoChange);
  q_Weight.resize(nPoints*nPoints);
  unsigned int qp = 0;
  for (unsigned int i = 0; i < nPoints; i++)
    for (unsigned int j = 0; j < nPoints; j++)
    {
       q_Qpoint(qp,0) = qPoints1(i);
       q_Qpoint(qp,1) = qPoints2(j);
       q_Qpoint(qp,2) = 0.;
       q_Weight(qp) = qWeights1(i) * qWeights2(j);
       qp++;
    }
}

void
JacobiGaussFiniteElementQuadrature::
tensorProduct2D(const RealVector& qPoints, const RealVector& qWeights, PointVector& q_Qpoint, RealVector& q_Weight)
{
  const unsigned int nPoints = qPoints.size();
  q_Qpoint.resize(nPoints*nPoints, Eigen::NoChange);
  q_Weight.resize(nPoints*nPoints);
  unsigned int qp = 0;
  for (unsigned int i = 0; i < nPoints; i++)
  {
    for (unsigned int j = 0; j < nPoints; j++)
    {
       q_Qpoint(qp,0) = qPoints(i);
       q_Qpoint(qp,1) = qPoints(j);
       q_Qpoint(qp,2) = 0.;
       q_Weight(qp) = qWeights(i) * qWeights(j);
       qp++;
    }
  }
}

void 
JacobiGaussFiniteElementQuadrature::
jacobiGaussQweights(double a , double b, unsigned int n, RealVector& qWeights, const RealVector& qPoints)
{
  // this function for gauss lobatto quadrature weights is not used
  // it should be modified because it is not accurate for large n
  qWeights.resize(n);
  double scale = 1.0;
  for (unsigned int i = 0; i < a; i++)
  { 
    scale *= 0.5; 
  }
  double a1 = pow(2,a+b+1);
  double a2 = Gamma::gammaF(a + n + 1);
  double a3 = Gamma::gammaF(b + n + 1);
  double a4 = Gamma::gammaF(a + b + n + 1);
  double a5 = Gamma::fact(n);
  double a6 = a1 * a2 * a3 / a4 / a5;
  

  for (unsigned int i = 0; i<n; i++)
  {
      double jp = JacobiPoly::evalJacobiDeriv(a,b,n,qPoints(i));
      qWeights(i) = ((a6 * scale) / ( 1.0  - pow(qPoints(i),2.0)) / pow(jp,2.0));
  }
}

void 
JacobiGaussFiniteElementQuadrature::
jacobiGaussQweights(double a, unsigned int n, RealVector& qWeights, const RealVector& qPoints)
{
  qWeights.resize(n);
  for (unsigned int i = 0; i<n; i++)
  {
    double jp = JacobiPoly::evalJacobiDeriv(a,0.0,n,qPoints(i));
    qWeights(i) = 2. / (( 1.0  - qPoints(i)*qPoints(i)) * jp*jp);
  }
}
