// 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 __serendipity_polynomial_h__
#define __serendipity_polynomial_h__

#include "EnumElemType.hpp"
#include "Eigen/Dense"

struct SerendipityPoly 
{
  static void computeBasisIndex(Eigen::Matrix<Eigen::Matrix<unsigned int,1,1>,Eigen::Dynamic,1>&  basisIndex, const ElemType itemType, const unsigned int order)
  {
    basisIndex.resize(order+1);
    for (unsigned int i = 0; i < (order + 1); i++)
    {
      basisIndex(i)(0) = i; 
    }
    return;
  }

  static void computeBasisIndex(Eigen::Matrix<Eigen::Matrix<unsigned int,2,1>,Eigen::Dynamic,1>&  basisIndex, const ElemType itemType, const unsigned int order)
  {
    switch(itemType)
    {
      case TRI3:
      case TRI6:
      {
        basisIndex.resize((order+1)*(order+2)/2);
        unsigned int v = 0;
        for (unsigned int i = 0; i <= order; i++)
        {
          for (unsigned int j = 0; j <= i; j++)
          {
            basisIndex(v)(0)= i - j;
            basisIndex(v)(1)= j;
            v++;
          }
        }
        return;
      }
      case QUAD4:
      case QUAD8:
      case QUAD9:
      {
	if (order == 1)
          basisIndex.resize(4);
	else
          basisIndex.resize(((order+1)*(order+2)/2)+2);
        unsigned int v = 0;
        for (unsigned int o = 0; o <= order; o++)
          for (unsigned int i = 0; i <= o; i++)
            for (unsigned int j = 0; j <= o; j++)
              if (i + j == o || (i == 1 && j == order) || (j == 1 && i == order))
              {
                basisIndex(v)(0)= i;
                basisIndex(v)(1)= j;
                v++;
              }
        return;
      }
      default:
      {
        std::cerr<< "SerendipityPoly: Unsupported Itemtype "<<itemType<<std::endl;
        std::logic_error("halt");
        return;
      }
    }
  }
  
  static void computeBasisIndex(Eigen::Matrix<Eigen::Matrix<unsigned int,3,1>,Eigen::Dynamic,1>&  basisIndex, const ElemType itemType, const unsigned int order)
  {
    switch(itemType)
    {
      case TET4:
      case TET10:
      {
        basisIndex.resize((order+1)*(order+2)*(order+3)/6);
        unsigned int v = 0;
        for (unsigned int i = 0; i <= order; i++)
        {
          for (unsigned int j = 0; j <= i; j++)
          {
            for (unsigned int k = 0; k <= j; k++)
            {
              basisIndex(v)(0)= i - j;
              basisIndex(v)(1)= j - k;
              basisIndex(v)(2)= k;
              v++;
            }
          }
        }
        return;
      }
      case HEX8:
      case HEX20:
      case HEX27:
      {
        unsigned int v = 0;
        for (unsigned int o = 0; o <= order; o++)
          for (unsigned int i = 0; i <= o; i++)
            for (unsigned int j = 0; j <= o; j++)
              for (unsigned int k = 0; k <= o; k++)
              {
                int deg2=0;
                if (i/2 >= 1)
                 deg2 += i;
                if (j/2 >= 1)
                 deg2 += j;
                if (k/2 >= 1)
                 deg2 += k;
                if (deg2<=order && (i == o || j == o || k == o))
                {
		  basisIndex.conservativeResize(v+1);
                  basisIndex(v)(0)= i;
                  basisIndex(v)(1)= j;
                  basisIndex(v)(2)= k;
                  v++;
                }
              }
        return;
      }
      default:
      {
        std::cerr<< "SerendipityPoly: Unsupported Itemtype "<<itemType<<std::endl;
        std::logic_error("halt");
        return;
      }
    }
  }

  static unsigned int modesNumber(const unsigned int Dim, const unsigned int order, const ElemType elemType)
  {
    switch (elemType)
    {
      case EDGE2:
      case EDGE3:
      {
        return (order+1);
      }
      case TRI3:
      case TRI6:
      {
        return ((order+1)*(order+2)/2);
      }
      case QUAD8:
      case QUAD9:
      {
        return ((order+1)*(order+2)/2)+2;
      }
      case TET4:
      case TET10:
      {
        return ((order+1)*(order+2)*(order+3)/6);
      }
      case HEX8:
      case HEX20:
      case HEX27:
      {
        return (order+1)*(order+1)*(order+1);
      }
      default:
      {
        std::cerr<< "SerendipityPoly: Unsupported itemType "<<elemType<<std::endl;
        std::logic_error("halt");
	return 0;
      }
    }
  }
};
 
#endif
