// 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 "EquiSpacedFiniteElementNodes.hpp"

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

void 
EquiSpacedFiniteElementNodes::
tetQuadrature(unsigned int order, PointVector& q_Qpoint)
{
  RealVector qPointsN;
  equispacedNodes(order, qPointsN); 
  tensorProduct3D(qPointsN, q_Qpoint);
  CollapsedMapping::mapTetNE(q_Qpoint);
}

void 
EquiSpacedFiniteElementNodes::
prismQuadrature(unsigned int order, PointVector& q_Qpoint)
{
  RealVector qPointsN;
  equispacedNodes(order, qPointsN); 
  tensorProduct3D(qPointsN, q_Qpoint);
  CollapsedMapping::mapPrismNE(q_Qpoint);
}

void 
EquiSpacedFiniteElementNodes::
pyramidQuadrature(unsigned int order, PointVector& q_Qpoint)
{
  RealVector qPointsN;
  equispacedNodes(order, qPointsN); 
  tensorProduct3D(qPointsN, q_Qpoint);
  CollapsedMapping::mapPyramidNE(q_Qpoint);
}

void 
EquiSpacedFiniteElementNodes::
hexQuadrature(unsigned int order, PointVector& q_Qpoint)
{
  RealVector qPointsE;
  equispacedNodes(order, qPointsE); 
  tensorProduct3D(qPointsE, q_Qpoint);
}

void 
EquiSpacedFiniteElementNodes::
triQuadrature(unsigned int order, PointVector& q_Qpoint)
{
  q_Qpoint.resize(((order)*(order+1))/2,Eigen::NoChange);
  if (order==1)
  {
    q_Qpoint(0,0) = -0.25;
    q_Qpoint(0,1) = -0.25;
  }
  else
  {
    double step = 2./((double)order);
    q_Qpoint(0,0) = step/4. - 1.;
    q_Qpoint(0,1) = step/4. - 1.;
    step += step/(4.*((double)order-1.));
    unsigned int pt = 0;
    for (unsigned int i = 0; i != order; i++)
      for (int j = 0; j != order-i; j++)
      {
        q_Qpoint(pt,0) = q_Qpoint(0,0)+i*step;
        q_Qpoint(pt,1) = q_Qpoint(0,1)+j*step;
        pt++;
      }
  }
}

void 
EquiSpacedFiniteElementNodes::
quadQuadrature(unsigned int order, PointVector& q_Qpoint)
{
  RealVector qPointsE;
  equispacedNodes(order, qPointsE); 
  tensorProduct2D(qPointsE, q_Qpoint);
}

void 
EquiSpacedFiniteElementNodes::
lineQuadrature(unsigned int order, PointVector& q_Qpoint)
{
  RealVector qPointsE;
  equispacedNodes(order, qPointsE); 
  q_Qpoint = Quadrature::PointVector::Zero(qPointsE.size(),3);
  q_Qpoint.col(0) = qPointsE;
}

void
EquiSpacedFiniteElementNodes::
tensorProduct3D(const RealVector& qPoints1, const RealVector& qPoints2, const RealVector& qPoints3, PointVector& q_Qpoint)
{
  assert( qPoints1.size() == qPoints2.size() &&  qPoints2.size() == qPoints3.size()); 
  const unsigned int nPoints = qPoints1.size();
  q_Qpoint.resize(nPoints*nPoints*nPoints, Eigen::NoChange);
  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);
         qp++;
      }
    }
  }
}

void
EquiSpacedFiniteElementNodes::
tensorProduct3D(const RealVector& qPoints, PointVector& q_Qpoint)
{
  const unsigned int nPoints = qPoints.size();
  q_Qpoint.resize(nPoints*nPoints*nPoints, Eigen::NoChange);
  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);
         qp++;
      }
}

void
EquiSpacedFiniteElementNodes::
tensorProduct2D(const RealVector& qPoints1, const RealVector& qPoints2, PointVector& q_Qpoint)
{
  assert( qPoints1.size() == qPoints2.size());
  const unsigned int nPoints = qPoints1.size();
  q_Qpoint.resize(nPoints*nPoints, Eigen::NoChange);
  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.;
       qp++;
    }
}

void
EquiSpacedFiniteElementNodes::
tensorProduct2D(const RealVector& qPoints, PointVector& q_Qpoint)
{
  const unsigned int nPoints = qPoints.size();
  q_Qpoint.resize(nPoints*nPoints, Eigen::NoChange);
  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.;
       qp++;
    }
  }
}

void 
EquiSpacedFiniteElementNodes::
equispacedNodes(const unsigned int n, RealVector& qPoints)
{
  qPoints.resize(n);
  if (n==1)
   qPoints(0) = 0.;
  else
  {
    double step = 2./(double)n;
    qPoints(0) = step/4. - 1.;
    step += step/(2.*((double)n-1.));
    for (unsigned int i = 1; i != n; i++)
      qPoints(i) = qPoints(i-1)+step;
  }
}
