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

#include "EnumElemType.hpp"
#include "EntityTypesMesh.hpp"
#include "Tuple.hpp"
#include <stdexcept>
#include <cassert>
#include <iostream>
#include "Quadrature.hpp"

struct ReferenceFiniteElement
{
  typedef Quadrature::PointVector PointVector;

  static ElemType elemTypeFromEntityType (const SpaMesh::ElemTypesMesh& itemType, const int num_nodes)
  {
    switch(itemType)
    {
      case SpaMesh::node:
      {
         switch(num_nodes)
         {
           case 1:
           {
              return NODEELEM;
	   }
           default:
           {
             std::cerr<<"Error in ReferenceFiniteElement: vertex with "<<num_nodes<<" nodes not implemented"<<std::endl;
             std::logic_error("halt"); 
             return INVALID_ELEM;
           }
	 }
      }
      case SpaMesh::line:
      {
         switch(num_nodes)
         {
           case 2:
           {
              return EDGE2;
	   }
           case 3:
           {
              return EDGE3;
	   }
           default:
           {
             std::cerr<<"Error in ReferenceFiniteElement: edge with "<<num_nodes<<" nodes not implemented"<<std::endl;
             std::logic_error("halt"); 
             return INVALID_ELEM;
           }
	 }
      }
      case SpaMesh::tri:
      {
         switch(num_nodes)
         {
           case 3:
           {
              return TRI3;
	   }
           case 6:
           {
              return TRI6;
	   }
           default:
           {
             std::cerr<<"Error in ReferenceFiniteElement: tri with "<<num_nodes<<" nodes not implemented"<<std::endl;
             std::logic_error("halt"); 
             return INVALID_ELEM;
           }
	 }
      }
      case SpaMesh::quad:
      {
         switch(num_nodes)
         {
           case 4:
           {
              return QUAD4;
	   }
           case 8:
           {
              return QUAD8;
	   }
           case 9:
           {
              return QUAD9;
	   }
           default:
           {
             std::cerr<<"Error in ReferenceFiniteElement: quad with "<<num_nodes<<" nodes not implemented"<<std::endl;
             std::logic_error("halt"); 
             return INVALID_ELEM;
           }
	 }
      }
      case SpaMesh::hex:
      {
         switch(num_nodes)
         {
           case 8:
           {
              return HEX8;
	   }
           case 20:
           {
              return HEX20;
	   }
           case 27:
           {
              return HEX27;
	   }
           default:
           {
             std::cerr<<"Error in ReferenceFiniteElement: quad with "<<num_nodes<<" nodes not implemented"<<std::endl;
             std::logic_error("halt"); 
             return INVALID_ELEM;
           }
	 }
      }
      case SpaMesh::tet:
      {
         switch(num_nodes)
         {
           case 4:
           {
              return TET4;
	   }
           case 10:
           {
              return TET10;
	   }
           default:
           {
             std::cerr<<"Error in ReferenceFiniteElement: quad with "<<num_nodes<<" nodes not implemented"<<std::endl;
             std::logic_error("halt"); 
             return INVALID_ELEM;
           }
	 }
      }
      case SpaMesh::prism:
      {
         switch(num_nodes)
         {
           case 6:
           {
              return PRISM6;
	   }
           case 15:
           {
              return PRISM15;
	   }
           default:
           {
             std::cerr<<"Error in ReferenceFiniteElement: quad with "<<num_nodes<<" nodes not implemented"<<std::endl;
             std::logic_error("halt"); 
             return INVALID_ELEM;
           }
	 }
      }
      case SpaMesh::pyram:
      {
         switch(num_nodes)
         {
           case 5:
           {
              return PYRAMID5;
	   }
           case 13:
           {
              return PYRAMID13;
	   }
           default:
           {
             std::cerr<<"Error in ReferenceFiniteElement: quad with "<<num_nodes<<" nodes not implemented"<<std::endl;
             std::logic_error("halt"); 
             return INVALID_ELEM;
           }
	 }
      }
      default:
      {
        std::cerr<<"Error in ReferenceFiniteElement: itemtype "<<itemType<<" not implemented"<<std::endl;
        std::logic_error("halt"); 
        return INVALID_ELEM;
      }
    }
  }
  
  static ElemType elemFirstOrderTypeFromEntityType (const SpaMesh::ElemTypesMesh& itemType)
  {
    switch(itemType)
    {
      case SpaMesh::node:
      {
        return NODEELEM;
      }
      case SpaMesh::line:
      {
        return EDGE2;
      }
      case SpaMesh::tri:
      {
        return TRI3;
      }
      case SpaMesh::quad:
      {
        return QUAD4;
      }
      case SpaMesh::hex:
      {
        return HEX8;
      }
      case SpaMesh::tet:
      {
        return TET4;
      }
      case SpaMesh::prism:
      {
        return PRISM6;
      }
      case SpaMesh::pyram:
      {
        return PYRAMID5;
      }
      default:
      {
        std::cerr<<"Error in ReferenceFiniteElement: itemtype "<<itemType<<" not implemented"<<std::endl;
        std::logic_error("halt"); 
        return INVALID_ELEM;
      }
    }
  }
  
  static unsigned int getGeoOrder(const ElemType itemType)
  {
    switch(itemType)
    {
      case EDGE2:
      case TRI3:
      case QUAD4:
      case TET4:
      case PYRAMID5:
      case PRISM6:
      case HEX8:
      case AGGLOELEM:
      {
	return 1;
      }
      case EDGE3:
      case PRISM15:
      case PYRAMID13:
      case TRI6:
      case QUAD8:
      case QUAD9:
      case TET10:
      case HEX20:
      case HEX27:
      {
        return 2;
      }
      default:
      {
        std::cerr<<"Error in ReferenceFiniteElement: itemtype "<<itemType<<" not implemented"<<std::endl;
        std::logic_error("halt"); 
        return 0;
      }
    }
  }

  static void getVerticesN(const ElemType itemType, PointVector& vertices)
  {
    switch(itemType)
    {
      case EDGE2:
      {
         vertices.resize(2,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = 0.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.; 
         vertices(1,1) = 0.; 
         vertices(1,2) = 0.; 
         return;
      }
      case EDGE3:
      {
         vertices.resize(3,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = 0.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.; 
         vertices(1,1) = 0.; 
         vertices(1,2) = 0.; 
         vertices(2,0) = 0.;
         vertices(2,1) = 0.;
         vertices(2,2) = 0.;
         return;
      }
      case TRI3:
      {
         vertices.resize(3,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         return;
      } 
      case PRISM6:
      {
         vertices.resize(6,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = -1.;
         vertices(3,2) = 1.;
         vertices(4,0) = 1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = -1.;
         vertices(5,1) = 1.;
         vertices(5,2) = 1.;
         return;
      } 
      case PRISM15:
      {
         vertices.resize(15,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = -1.;
         vertices(3,2) = 1.;
         vertices(4,0) = 1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = -1.;
         vertices(5,1) = 1.;
         vertices(5,2) = 1.;
         vertices(6,0) = 0.;
         vertices(6,1) = -1.;
         vertices(6,2) = -1.;
         vertices(7,0) = 1.;
         vertices(7,1) = 0.;
         vertices(7,2) = -1.;
         vertices(8,0) = -1.;
         vertices(8,1) = 0.;
         vertices(8,2) = -1.;
         vertices(9,0) = -1.;
         vertices(9,1) = -1.;
         vertices(9,2) = 0.;
         vertices(10,0) = 1.;
         vertices(10,1) = -1.;
         vertices(10,2) = 0.;
         vertices(11,0) = -1.;
         vertices(11,1) = 1.;
         vertices(11,2) = 0.;
         vertices(12,0) = 0.;
         vertices(12,1) = -1.;
         vertices(12,2) = 1.;
         vertices(13,0) = 1.;
         vertices(13,1) = 0.;
         vertices(13,2) = 1.;
         vertices(14,0) = -1.;
         vertices(14,1) = 0.;
         vertices(14,2) = 1.;
         return;
      } 
      case PYRAMID5:
      {
         vertices.resize(5,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         return;
      } 
      case PYRAMID13:
      {
         vertices.resize(13,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = 0.;
         vertices(5,1) = -1.;
         vertices(5,2) = -1.;
         vertices(6,0) = 1.;
         vertices(6,1) = 0.;
         vertices(6,2) = -1.;
         vertices(7,0) = 0.;
         vertices(7,1) = 1.;
         vertices(7,2) = -1.;
         vertices(8,0) = -1.;
         vertices(8,1) = 0.;
         vertices(8,2) = -1.;
         vertices(9,0) = -1.;
         vertices(9,1) = -1.;
         vertices(9,2) = 0.;
         vertices(10,0) = 1.;
         vertices(10,1) = -1.;
         vertices(10,2) = 0.;
         vertices(11,0) = 1.;
         vertices(11,1) = 1.;
         vertices(11,2) = 0.;
         vertices(12,0) = -1.;
         vertices(12,1) = 1.;
         vertices(12,2) = 0.;
         return;
      } 
      case TRI6:
      {
         vertices.resize(6,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         vertices(3,0) = 0.;
         vertices(3,1) = -1.;
         vertices(3,2) = 0.;
         vertices(4,0) = 1.;
         vertices(4,1) = 0.;
         vertices(4,2) = 0.;
         vertices(5,0) = -1.;
         vertices(5,1) = 0.;
         vertices(5,2) = 0.;
         return;
      } 
      case QUAD4:
      {
         vertices.resize(4,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = 0.;
         return;
      }
      case QUAD8:
      {
         vertices.resize(8,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = 0.;
         vertices(4,0) = 0.;
         vertices(4,1) = -1.;
         vertices(4,2) = 0.;
         vertices(5,0) = 1.;
         vertices(5,1) = 0.;
         vertices(5,2) = 0.;
         vertices(6,0) = 0.;
         vertices(6,1) = 1.;
         vertices(6,2) = 0.;
         vertices(7,0) = -1.;
         vertices(7,1) = 0.;
         vertices(7,2) = 0.;
         return;
      }
      case QUAD9:
      {
         vertices.resize(9,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = 0.;
         vertices(4,0) = 0.;
         vertices(4,1) = -1.;
         vertices(4,2) = 0.;
         vertices(5,0) = 1.;
         vertices(5,1) = 0.;
         vertices(5,2) = 0.;
         vertices(6,0) = 0.;
         vertices(6,1) = 1.;
         vertices(6,2) = 0.;
         vertices(7,0) = -1.;
         vertices(7,1) = 0.;
         vertices(7,2) = 0.;
         vertices(8,0) = 0.;
         vertices(8,1) = 0.;
         vertices(8,2) = 0.;
         return;
      }
      case TET4:
      {
         vertices.resize(4,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = -1.;
         vertices(3,2) = 1.;
         return;
      } 
      case TET10:
      {
         vertices.resize(10,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = -1.;
         vertices(3,2) = 1.;
         vertices(4,0) = 0.;
         vertices(4,1) = -1.;
         vertices(4,2) = -1.;
         vertices(5,0) = 1.;
         vertices(5,1) = 0.;
         vertices(5,2) = -1.;
         vertices(6,0) = -1.;
         vertices(6,1) = 0.;
         vertices(6,2) = -1.;
         vertices(7,0) = -1.;
         vertices(7,1) = -1.;
         vertices(7,2) = 0.;
         vertices(8,0) = 1.;
         vertices(8,1) = -1.;
         vertices(8,2) = 0.;
         vertices(9,0) = 1.;
         vertices(9,1) = 1.;
         vertices(9,2) = 0.;
         return;
      } 
      case HEX8:
      {
         vertices.resize(8,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = 1.;
         vertices(5,1) = -1.;
         vertices(5,2) = 1.;
         vertices(6,0) = 1.;
         vertices(6,1) = 1.;
         vertices(6,2) = 1.;
         vertices(7,0) = -1.;
         vertices(7,1) = 1.;
         vertices(7,2) = 1.;
         return;
      }
      case HEX20:
      {
         vertices.resize(20,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = 1.;
         vertices(5,1) = -1.;
         vertices(5,2) = 1.;
         vertices(6,0) = 1.;
         vertices(6,1) = 1.;
         vertices(6,2) = 1.;
         vertices(7,0) = -1.;
         vertices(7,1) = 1.;
         vertices(7,2) = 1.;
         vertices(8,0) = 0.;
         vertices(8,1) = -1.;
         vertices(8,2) = -1.;
         vertices(9,0) = 1.;
         vertices(9,1) = 0.;
         vertices(9,2) = -1.;
         vertices(10,0) = 0.;
         vertices(10,1) = 1.;
         vertices(10,2) = -1.;
         vertices(11,0) = -1.;
         vertices(11,1) = 0.;
         vertices(11,2) = -1.;
         vertices(12,0) = -1.;
         vertices(12,1) = -1.;
         vertices(12,2) = 0.;
         vertices(13,0) = 1.;
         vertices(13,1) = -1.;
         vertices(13,2) = 0.;
         vertices(14,0) = 1.;
         vertices(14,1) = 1.;
         vertices(14,2) = 0.;
         vertices(15,0) = -1.;
         vertices(15,1) = 1.;
         vertices(15,2) = 0.;
         vertices(16,0) = 0.;
         vertices(16,1) = -1.;
         vertices(16,2) = 1.;
         vertices(17,0) = 1.;
         vertices(17,1) = 0.;
         vertices(17,2) = 1.;
         vertices(18,0) = 0.;
         vertices(18,1) = 1.;
         vertices(18,2) = 1.;
         vertices(19,0) = -1.;
         vertices(19,1) = 0.;
         vertices(19,2) = 1.;
         return;
      }
      case HEX27:
      {
         vertices.resize(27,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = 1.;
         vertices(5,1) = -1.;
         vertices(5,2) = 1.;
         vertices(6,0) = 1.;
         vertices(6,1) = 1.;
         vertices(6,2) = 1.;
         vertices(7,0) = -1.;
         vertices(7,1) = 1.;
         vertices(7,2) = 1.;
         vertices(8,0) = 0.;
         vertices(8,1) = -1.;
         vertices(8,2) = -1.;
         vertices(9,0) = 1.;
         vertices(9,1) = 0.;
         vertices(9,2) = -1.;
         vertices(10,0) = 0.;
         vertices(10,1) = 1.;
         vertices(10,2) = -1.;
         vertices(11,0) = -1.;
         vertices(11,1) = 0.;
         vertices(11,2) = -1.;
         vertices(12,0) = -1.;
         vertices(12,1) = -1.;
         vertices(12,2) = 0.;
         vertices(13,0) = 1.;
         vertices(13,1) = -1.;
         vertices(13,2) = 0.;
         vertices(14,0) = 1.;
         vertices(14,1) = 1.;
         vertices(14,2) = 0.;
         vertices(15,0) = -1.;
         vertices(15,1) = 1.;
         vertices(15,2) = 0.;
         vertices(16,0) = 0.;
         vertices(16,1) = -1.;
         vertices(16,2) = 1.;
         vertices(17,0) = 1.;
         vertices(17,1) = 0.;
         vertices(17,2) = 1.;
         vertices(18,0) = 0.;
         vertices(18,1) = 1.;
         vertices(18,2) = 1.;
         vertices(19,0) = -1.;
         vertices(19,1) = 0.;
         vertices(19,2) = 1.;
         vertices(20,0) = 0.;
         vertices(20,1) = -1.;
         vertices(20,2) = 0.;
         vertices(21,0) = 1.;
         vertices(21,1) = 0.;
         vertices(21,2) = 0.;
         vertices(22,0) = 0.;
         vertices(22,1) = 1.;
         vertices(22,2) = 0.;
         vertices(23,0) = -1.;
         vertices(23,1) = 0.;
         vertices(23,2) = 0.;
         vertices(24,0) = 0.;
         vertices(24,1) = 0.;
         vertices(24,2) = -1.;
         vertices(25,0) = 0.;
         vertices(25,1) = 0.;
         vertices(25,2) = 1.;
         vertices(26,0) = 0.;
         vertices(26,1) = 0.;
         vertices(26,2) = 0.;
         return;
      }
      default:
      {
        std::cerr<<"Error in ReferenceFiniteElement: itemtype "<<itemType<<" not implemented"<<std::endl;
        std::logic_error("halt"); 
        return;
      }
    }
  }
  
  static void getVerticesOnlyE(const ElemType itemType, PointVector& vertices)
  {
    switch(itemType)
    {
      case EDGE2:
      case EDGE3:
      {
         vertices.resize(2,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = 0.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.; 
         vertices(1,1) = 0.; 
         vertices(1,2) = 0.; 
         return;
      }
      case TRI3:
      case TRI6:
      {
         vertices.resize(3,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         return;
      } 
      case QUAD4:
      case QUAD8:
      case QUAD9:
      {
         vertices.resize(4,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = 0.;
         return;
      }
      case TET4:
      case TET10:
      {
         vertices.resize(4,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = -1.;
         vertices(3,2) = 1.;
         return;
      } 
      case HEX8:
      case HEX20:
      case HEX27:
      {
         vertices.resize(8,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = 1.;
         vertices(5,1) = -1.;
         vertices(5,2) = 1.;
         vertices(6,0) = 1.;
         vertices(6,1) = 1.;
         vertices(6,2) = 1.;
         vertices(7,0) = -1.;
         vertices(7,1) = 1.;
         vertices(7,2) = 1.;
         return;
      }
      case PRISM6:
      case PRISM15:
      {
         vertices.resize(6,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = -1.;
         vertices(3,2) = 1.;
         vertices(4,0) = 1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = -1.;
         vertices(5,1) = 1.;
         vertices(5,2) = 1.;
         return;
      }
      case PYRAMID5:
      case PYRAMID13:
      {
         vertices.resize(5,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         return;
      }
      default:
      {
        std::cerr<<"Error in ReferenceFiniteElement: itemtype "<<itemType<<" not implemented"<<std::endl;
        std::logic_error("halt"); 
        return;
      }
    }
  }
  
  static void getVerticesE(const ElemType itemType, PointVector& vertices)
  {
    switch(itemType)
    {
      case EDGE2:
      {
         vertices.resize(2,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = 0.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.; 
         vertices(1,1) = 0.; 
         vertices(1,2) = 0.; 
         return;
      }
      case EDGE3:
      {
         vertices.resize(3,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = 0.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.; 
         vertices(1,1) = 0.; 
         vertices(1,2) = 0.; 
         vertices(2,0) = 0.;
         vertices(2,1) = 0.;
         vertices(2,2) = 0.;
         return;
      }
      case TRI3:
      {
         vertices.resize(3,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         return;
      } 
      case TRI6:
      {
         vertices.resize(6,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         vertices(3,0) = 0.;
         vertices(3,1) = -1.;
         vertices(3,2) = 0.;
         vertices(4,0) = 0.;
         vertices(4,1) = 0.;
         vertices(4,2) = 0.;
         vertices(5,0) = -1.;
         vertices(5,1) = 0.;
         vertices(5,2) = 0.;
         return;
      } 
      case QUAD4:
      {
         vertices.resize(4,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = 0.;
         return;
      }
      case QUAD8:
      {
         vertices.resize(8,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = 0.;
         vertices(4,0) = 0.;
         vertices(4,1) = -1.;
         vertices(4,2) = 0.;
         vertices(5,0) = 1.;
         vertices(5,1) = 0.;
         vertices(5,2) = 0.;
         vertices(6,0) = 0.;
         vertices(6,1) = 1.;
         vertices(6,2) = 0.;
         vertices(7,0) = -1.;
         vertices(7,1) = 0.;
         vertices(7,2) = 0.;
         return;
      }
      case QUAD9:
      {
         vertices.resize(9,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = 0.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = 0.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = 0.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = 0.;
         vertices(4,0) = 0.;
         vertices(4,1) = -1.;
         vertices(4,2) = 0.;
         vertices(5,0) = 1.;
         vertices(5,1) = 0.;
         vertices(5,2) = 0.;
         vertices(6,0) = 0.;
         vertices(6,1) = 1.;
         vertices(6,2) = 0.;
         vertices(7,0) = -1.;
         vertices(7,1) = 0.;
         vertices(7,2) = 0.;
         vertices(8,0) = 0.;
         vertices(8,1) = 0.;
         vertices(8,2) = 0.;
         return;
      }
      case TET4:
      {
         vertices.resize(4,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = -1.;
         vertices(3,2) = 1.;
         return;
      } 
      case PRISM6:
      {
         vertices.resize(6,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = -1.;
         vertices(3,2) = 1.;
         vertices(4,0) = 1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = -1.;
         vertices(5,1) = 1.;
         vertices(5,2) = 1.;
         return;
      } 
      case PYRAMID5:
      {
         vertices.resize(5,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         return;
      } 
      case PRISM15:
      {
         vertices.resize(15,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = -1.;
         vertices(3,2) = 1.;
         vertices(4,0) = 1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = -1.;
         vertices(5,1) = 1.;
         vertices(5,2) = 1.;
         vertices(6,0) = 0.;
         vertices(6,1) = -1.;
         vertices(6,2) = -1.;
         vertices(7,0) = 0.;
         vertices(7,1) = 0.;
         vertices(7,2) = -1.;
         vertices(8,0) = -1.;
         vertices(8,1) = 0.;
         vertices(8,2) = -1.;
         vertices(9,0) = -1.;
         vertices(9,1) = -1.;
         vertices(9,2) = 0.;
         vertices(10,0) = 1.;
         vertices(10,1) = -1.;
         vertices(10,2) = 0.;
         vertices(11,0) = -1.;
         vertices(11,1) = 1.;
         vertices(11,2) = 0.;
         vertices(12,0) = 0.;
         vertices(12,1) = -1.;
         vertices(12,2) = 1.;
         vertices(13,0) = 0.;
         vertices(13,1) = 0.;
         vertices(13,2) = 1.;
         vertices(14,0) = -1.;
         vertices(14,1) = 0.;
         vertices(14,2) = 1.;
         return;
      } 
      case PYRAMID13:
      {
         vertices.resize(13,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = 0.;
         vertices(5,1) = -1.;
         vertices(5,2) = -1.;
         vertices(6,0) = 1.;
         vertices(6,1) = 0.;
         vertices(6,2) = -1.;
         vertices(7,0) = 0.;
         vertices(7,1) = 1.;
         vertices(7,2) = -1.;
         vertices(8,0) = -1.;
         vertices(8,1) = 0.;
         vertices(8,2) = -1.;
         vertices(9,0) = -1.;
         vertices(9,1) = -1.;
         vertices(9,2) = 0.;
         vertices(10,0) = 0.;
         vertices(10,1) = -1.;
         vertices(10,2) = 0.;
         vertices(11,0) = 0.;
         vertices(11,1) = 0.;
         vertices(11,2) = 0.;
         vertices(12,0) = -1.;
         vertices(12,1) = 0.;
         vertices(12,2) = 0.;
         return;
      } 
      case TET10:
      {
         vertices.resize(10,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = -1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = -1.;
         vertices(3,2) = 1.;
         vertices(4,0) = 0.;
         vertices(4,1) = -1.;
         vertices(4,2) = -1.;
         vertices(5,0) = 0.;
         vertices(5,1) = 0.;
         vertices(5,2) = -1.;
         vertices(6,0) = -1.;
         vertices(6,1) = 0.;
         vertices(6,2) = -1.;
         vertices(7,0) = -1.;
         vertices(7,1) = -1.;
         vertices(7,2) = 0.;
         vertices(8,0) = 0.;
         vertices(8,1) = -1.;
         vertices(8,2) = 0.;
         vertices(9,0) = -1.;
         vertices(9,1) = 0.;
         vertices(9,2) = 0.;
         return;
      } 
      case HEX8:
      {
         vertices.resize(8,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = 1.;
         vertices(5,1) = -1.;
         vertices(5,2) = 1.;
         vertices(6,0) = 1.;
         vertices(6,1) = 1.;
         vertices(6,2) = 1.;
         vertices(7,0) = -1.;
         vertices(7,1) = 1.;
         vertices(7,2) = 1.;
         return;
      }
      case HEX20:
      {
         vertices.resize(20,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = 1.;
         vertices(5,1) = -1.;
         vertices(5,2) = 1.;
         vertices(6,0) = 1.;
         vertices(6,1) = 1.;
         vertices(6,2) = 1.;
         vertices(7,0) = -1.;
         vertices(7,1) = 1.;
         vertices(7,2) = 1.;
         vertices(8,0) = 0.;
         vertices(8,1) = -1.;
         vertices(8,2) = -1.;
         vertices(9,0) = 1.;
         vertices(9,1) = 0.;
         vertices(9,2) = -1.;
         vertices(10,0) = 0.;
         vertices(10,1) = 1.;
         vertices(10,2) = -1.;
         vertices(11,0) = -1.;
         vertices(11,1) = 0.;
         vertices(11,2) = -1.;
         vertices(12,0) = -1.;
         vertices(12,1) = -1.;
         vertices(12,2) = 0.;
         vertices(13,0) = 1.;
         vertices(13,1) = -1.;
         vertices(13,2) = 0.;
         vertices(14,0) = 1.;
         vertices(14,1) = 1.;
         vertices(14,2) = 0.;
         vertices(15,0) = -1.;
         vertices(15,1) = 1.;
         vertices(15,2) = 0.;
         vertices(16,0) = 0.;
         vertices(16,1) = -1.;
         vertices(16,2) = 1.;
         vertices(17,0) = 1.;
         vertices(17,1) = 0.;
         vertices(17,2) = 1.;
         vertices(18,0) = 0.;
         vertices(18,1) = 1.;
         vertices(18,2) = 1.;
         vertices(19,0) = -1.;
         vertices(19,1) = 0.;
         vertices(19,2) = 1.;
         return;
      }
      case HEX27:
      {
         vertices.resize(27,Eigen::NoChange);
         vertices(0,0) = -1.;
         vertices(0,1) = -1.;
         vertices(0,2) = -1.;
         vertices(1,0) = 1.;
         vertices(1,1) = -1.;
         vertices(1,2) = -1.;
         vertices(2,0) = 1.;
         vertices(2,1) = 1.;
         vertices(2,2) = -1.;
         vertices(3,0) = -1.;
         vertices(3,1) = 1.;
         vertices(3,2) = -1.;
         vertices(4,0) = -1.;
         vertices(4,1) = -1.;
         vertices(4,2) = 1.;
         vertices(5,0) = 1.;
         vertices(5,1) = -1.;
         vertices(5,2) = 1.;
         vertices(6,0) = 1.;
         vertices(6,1) = 1.;
         vertices(6,2) = 1.;
         vertices(7,0) = -1.;
         vertices(7,1) = 1.;
         vertices(7,2) = 1.;
         vertices(8,0) = 0.;
         vertices(8,1) = -1.;
         vertices(8,2) = -1.;
         vertices(9,0) = 1.;
         vertices(9,1) = 0.;
         vertices(9,2) = -1.;
         vertices(10,0) = 0.;
         vertices(10,1) = 1.;
         vertices(10,2) = -1.;
         vertices(11,0) = -1.;
         vertices(11,1) = 0.;
         vertices(11,2) = -1.;
         vertices(12,0) = -1.;
         vertices(12,1) = -1.;
         vertices(12,2) = 0.;
         vertices(13,0) = 1.;
         vertices(13,1) = -1.;
         vertices(13,2) = 0.;
         vertices(14,0) = 1.;
         vertices(14,1) = 1.;
         vertices(14,2) = 0.;
         vertices(15,0) = -1.;
         vertices(15,1) = 1.;
         vertices(15,2) = 0.;
         vertices(16,0) = 0.;
         vertices(16,1) = -1.;
         vertices(16,2) = 1.;
         vertices(17,0) = 1.;
         vertices(17,1) = 0.;
         vertices(17,2) = 1.;
         vertices(18,0) = 0.;
         vertices(18,1) = 1.;
         vertices(18,2) = 1.;
         vertices(19,0) = -1.;
         vertices(19,1) = 0.;
         vertices(19,2) = 1.;
         vertices(20,0) = 0.;
         vertices(20,1) = -1.;
         vertices(20,2) = 0.;
         vertices(21,0) = 1.;
         vertices(21,1) = 0.;
         vertices(21,2) = 0.;
         vertices(22,0) = 0.;
         vertices(22,1) = 1.;
         vertices(22,2) = 0.;
         vertices(23,0) = -1.;
         vertices(23,1) = 0.;
         vertices(23,2) = 0.;
         vertices(24,0) = 0.;
         vertices(24,1) = 0.;
         vertices(24,2) = -1.;
         vertices(25,0) = 0.;
         vertices(25,1) = 0.;
         vertices(25,2) = 1.;
         vertices(26,0) = 0.;
         vertices(26,1) = 0.;
         vertices(26,2) = 0.;
         return;
      }
      default:
      {
        std::cerr<<"Error in ReferenceFiniteElement: itemtype "<<itemType<<" not implemented"<<std::endl;
        std::logic_error("halt"); 
        return;
      }
    }
  }
  

  static bool onReferenceElement(const Point3D& point, const ElemType itemType, double epsilon = 1.e-6)
  {
    switch(itemType)
    {
      case EDGE2:
      case EDGE3:
      {
         if (point(0) >= -1.-epsilon && point(0) <= 1.+epsilon)
           return true;
         return false;
      }
      case TRI3:
      case TRI6:
      {
         if (point(0) >= -1.-epsilon && 
             point(1) >= -1.-epsilon && 
             point(0) + point(1) <= 0.+epsilon)
           return true;
         return false;
      } 
      case QUAD4:
      case QUAD8:
      case QUAD9:
      {
         if (point(0) >= -1.-epsilon && point(0) <= 1.+epsilon &&
             point(1) >= -1.-epsilon && point(1) <= 1.+epsilon)
           return true;
         return false;
      }
      case TET4:
      case TET10:
      {
         if (point(0) >= -1.-epsilon && 
             point(1) >= -1.-epsilon && 
             point(2) >= -1.-epsilon && 
             point(0) + point(1) + point(2) <= -1.+epsilon)
           return true;
         return false;
      } 
      case HEX8:
      case HEX20:
      case HEX27:
      {
         if (point(0) >= -1.-epsilon && point(0) <= 1.+epsilon && 
             point(1) >= -1.-epsilon && point(1) <= 1.+epsilon &&
             point(2) >= -1.-epsilon && point(2) <= 1.+epsilon)
           return true;
         return false;
      }
      case PRISM6:
      case PRISM15:
      {
         if (point(0) >= -1.-epsilon && 
             point(1) >= -1.-epsilon && 
             point(2) >= -1.-epsilon && point(2) <= 1.+epsilon &&
             point(0) + point(1) <= 0.+epsilon)
           return true;
         return false;
      } 
      case PYRAMID5:
      case PYRAMID13:
      {
         if (point(0) >= -1.-epsilon && point(0) <= 1.+epsilon && 
             point(1) >= -1.-epsilon && point(1) <= 1.+epsilon &&
             point(2) >= -1.-epsilon && 
             point(0) + point(2) <= 0.+epsilon &&
             point(1) + point(2) <= 0.+epsilon)
           return true;
         return false;
      } 
      default:
      {
        std::cerr<<"Error in ReferenceFiniteElement: itemtype "<<itemType<<" not implemented"<<std::endl;
        std::logic_error("halt"); 
        return false;
      }
    }
  }  
}; 
#endif
