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

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

unsigned int 
JacobiGaussFiniteElementQuadrature::
nPoints(const unsigned int dim, const unsigned int q_Order, const ElemType itemType)
{
  return pow(q_Order/2 + 1, dim);
}

void 
JacobiGaussFiniteElementQuadrature::
initQRule(Quadrature& quadrature)
{
  ElemType& q_itemType = quadrature.getItemType();
  unsigned int& q_Order = quadrature.getOrder();
  PointVector& q_Qpoint = quadrature.getQpoint();
  RealVector& q_Weight = quadrature.getWeight();
  switch (q_itemType)
  {
    case TET4:
    case TET10:
    { 
      tetQuadrature(q_Order/2 + 1, q_Qpoint, q_Weight);
      return;
    } 
    case PRISM6:
    case PRISM15:
    { 
      prismQuadrature(q_Order/2 + 1, q_Qpoint, q_Weight);
      return;
    } 
    case PYRAMID5:
    case PYRAMID13:
    { 
      pyramidQuadrature(q_Order/2 + 1, q_Qpoint, q_Weight);
      return;
    } 
    case TRI3:
    case TRI6:
    {
      triQuadrature(q_Order/2 + 1, q_Qpoint, q_Weight);
      return;
    }
    case QUAD4:
    case QUAD8:
    case QUAD9:
    {
      quadQuadrature(q_Order/2 + 1, q_Qpoint, q_Weight);
      return;
    }
    case HEX8:
    case HEX20:
    case HEX27:
    {
      hexQuadrature(q_Order/2 + 1, q_Qpoint, q_Weight);
      return;
    }
    case EDGE2:
    case EDGE3:
    {
      lineQuadrature(q_Order/2 + 1, q_Qpoint, q_Weight);
      return;
    }
    case NODEELEM:
    {
      vertexQuadrature(q_Order/2 + 1, q_Qpoint, q_Weight);
      return;
    }
    default:
    {
      std::cerr<< "JacobiGauss Quadrature: Unsupported Itemtype "<<q_itemType<<std::endl;
      std::logic_error("halt");
      return;
    }
  }
}

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

