/*
Copyright (c) 2015, Michael Kazhdan and Fabian Prada
All rights reserved.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of
conditions and the following disclaimer. Redistributions in binary form must reproduce
the above copyright notice, this list of conditions and the following disclaimer
in the documentation and/or other materials provided with the distribution. 

Neither the name of the Johns Hopkins University nor the names of its contributors
may be used to endorse or promote products derived from this software without specific
prior written permission. 

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES 
OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
TO, PROCUREMENT OF SUBSTITUTE  GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
*/

#ifndef FEM_INCLUDED
#define FEM_INCLUDED

#define NEW_FEM
#define NEW_FEM_CODE

#undef INCLUDE_EXTRA_FEM

#include <string.h>
#include <algorithm>
#include "SparseMatrix.h"
#include "Miscellany.h"
#include "Geometry.h"
#include "Array.h"


/************ Notes ****************
-.The evaluation w.r.t. to some bases can be metric dependent:
-- BASIS_1_CONFORMING requires the metric to define a rotation.
-- BASIS_2_WHITNEY requires the metric to transform the integral into a value.
. Using the DYNAMIC_MESH flag, we do not copy over the triangles, so edge-flips are performed on the input array.
. Althought we are calling this "FEM", the definition of the stiffness matrix is done using "DEC". (Probably should call it "DDG".)
. Assumes the mesh is oriented
. The values of the vector fields live in the cotangent space. (That is given a vector field v and a tangent vector w, the inner-product is simple < v , w > not g(v,w).)
. Throughout, "tensor" refers to the inner-product on the tangent space.
. There is some probably (parallelization) resulting in inconsistent runs (and undefined error values).
. Add an explicit method for computing the mass-matrix inverse.
. Modify the code to copy the triangles (?). This would enable computing edge-flips.
. Note that the method will fail if we try computing the 2-form Laplacian where the 1-form basis is BASIS_1_TRIANGLE_CONSTANT because we haven't defined the diagonal correctly.
. Make the Point2DWrapper its own class and allow it to have coefficients that are themselves vectors.
. Merge GetMassMatrix and GetLumpedMassMatrix
. Check that flow/exponentiation uses tangent vectors where appropriate (and overload the coordinate transform function to support processing of tangent vectors)
 ***********************************/

namespace FEM
{	
	//////////////
	// Elements //
	//////////////
	enum
	{
		ELEMENT_VERTEX ,
		ELEMENT_EDGE ,
		ELEMENT_TRIANGLE ,
		ELEMENT_COUNT
	};
	static const char* ElementNames[] = { "vertex" , "edge" , "triangle" };
	template< unsigned int Type > struct ElementInfo
	{
		static const unsigned int ElementsPerTriangle;
	};
	static void TestElementType( unsigned int ElementType , const char* header , bool forceFailure )
	{
		if( ElementType>=ELEMENT_COUNT ) Miscellany::ErrorOut( "%s: Unrecognized element type: 0 <= %d < %d" , header , ElementType , ELEMENT_COUNT );
		else if( forceFailure ) Miscellany::ErrorOut( "%s: Element type unsupported: %d (%s)\n" , header , ElementType , ElementNames[ElementType] );
	}

	//////////////////////////////
	// Types of Basis Functions //
	//////////////////////////////
	enum
	{
		BASIS_0_WHITNEY ,				// 0-form represented as the linear combination of hat functions (1 value per vertices)
		BASIS_1_CONFORMING ,			// 1-form represented as the linear combination of gradients and 90-degree rotated gradients of hat functions (2 values per vertex)
		BASIS_1_WHITNEY ,				// 1-form represented as the linear combination of the Whitney 1-form basis functions (1 value per edge)
		BASIS_1_TRIANGLE_CONSTANT ,		// 1-form represented as piecewise (per triangle) constant cotangent vectors (2 values per triangle)
		BASIS_2_WHITNEY ,				// 2-form represented as the linear combination of the Whitney  2-form basis functions (1 value per triangle)
		BASIS_2_VERTEX_CONSTANT ,		// 2-form represented as piecewise (around vertex) constant density fields (1 value per vertex)
		BASIS_COUNT
	};

	static const char* BasisNames[] = { "scalar whitney" , "vector conforming" , "vector whitney" , "vector triangle constant" , "density whitney" , "density vertex constant" };
	template< unsigned int Type > struct BasisInfo
	{
		static const unsigned int CoefficientsPerElement;
		static const unsigned int Dimension;
		static const bool Lumpable;
		static const unsigned int ElementType;
		static const unsigned int Coefficients;
		static const bool Singular;
	};
#include "FEM.BasisInfo.h"
	template< class Real , unsigned int Type > struct BasisInfoSystem
	{
		typedef ::SquareMatrix<          Real , BasisInfo< Type >::Coefficients >     Matrix;
		typedef ::SquareMatrix< unsigned char , BasisInfo< Type >::Coefficients > MaskMatrix;
		typedef ::Point       <          Real , BasisInfo< Type >::Coefficients >      Point;
	};
	template< class Real , unsigned int InType , unsigned int OutType > struct BasisInfoSystem2
	{
		typedef ::Matrix< unsigned char , BasisInfo< InType >::Coefficients , BasisInfo< OutType >::Coefficients > MaskMatrix;
		typedef ::Matrix<          Real , BasisInfo< InType >::Coefficients , BasisInfo< OutType >::Coefficients >     Matrix;
	};

	static void TestBasisType( unsigned int BasisType , const char* header , bool forceFailure )
	{
		if( BasisType>=BASIS_COUNT ) Miscellany::ErrorOut( "%s: Unrecognized basis type: 0 <= %d < %d" , header , BasisType , BASIS_COUNT );
		else if( forceFailure ) Miscellany::ErrorOut( "%s: Basis type unsupported: %d (%s)" , header , BasisType , BasisNames[BasisType] );
	}

	///////////////////////////////
	// Tangent/Cotangent Vectors //
	///////////////////////////////
	enum
	{
		VECTOR_PRIMAL ,
		VECTOR_DUAL
	};
	// Defining Tangent/Cotangent vector classes
	template< class Real , int ID >
	struct Point2DWrapper : public Point2D< Real >
	{
		Point2DWrapper( void ) : Point2D< Real >() {;}
		Point2DWrapper( Real v1 , Real v2 ) : Point2D< Real >( v1 , v2 ) {;}
		Point2DWrapper( const Point2D< Real >& v ) : Point2D< Real >( v ) {;}
		Point2DWrapper  operator +  ( const Point2DWrapper& v ) const { return Point2D< Real >::operator + ( v ); }
		Point2DWrapper  operator -  ( const Point2DWrapper& v ) const { return Point2D< Real >::operator - ( v ); }
		Point2DWrapper& operator += ( const Point2DWrapper& v )       { Point2D< Real >::operator += ( v ) ; return *this; }
		Point2DWrapper& operator -= ( const Point2DWrapper& v )       { Point2D< Real >::operator -= ( v ) ; return *this; }
		Point2DWrapper  operator *  ( Real s ) const { return Point2D< Real >::operator * ( s ); }
		Point2DWrapper  operator /  ( Real s ) const { return Point2D< Real >::operator / ( s ); }
		Point2DWrapper& operator *= ( Real s )       { Point2D< Real >::operator *= ( s ) ; return *this; }
		Point2DWrapper& operator /= ( Real s )       { Point2D< Real >::operator /= ( s ) ; return *this; }
		Point2DWrapper  operator -  ( void ) const { return Point2D< Real >::operator - ( ); }
	};
	template< class Real > using   TangentVector = Point2DWrapper< Real , VECTOR_PRIMAL >;
	template< class Real > using CotangentVector = Point2DWrapper< Real , VECTOR_DUAL   >;

	// Operators for moving between the primal and dual spaces
	template< class Real >   TangentVector< Real > Dual( const SquareMatrix< Real , 2 >& tensor , CotangentVector< Real > v ){ return tensor.inverse() * v; }
	template< class Real > CotangentVector< Real > Dual( const SquareMatrix< Real , 2 >& tensor ,   TangentVector< Real > v ){ return tensor * v; }

	// Operators for measuring stuff
	template< class Real > Real Dot( const SquareMatrix< Real , 2 >& tensor ,   TangentVector< Real > v1 ,   TangentVector< Real > v2 ){ return Point2D< Real >::Dot( v1 , tensor*v2 ); }
	template< class Real > Real Dot( const SquareMatrix< Real , 2 >& tensor , CotangentVector< Real > v1 , CotangentVector< Real > v2 ){ return Point2D< Real >::Dot( v1 , tensor.inverse()*v2 ); }
	template< class Real > Real SquareLength( const SquareMatrix< Real , 2 >& tensor ,   TangentVector< Real > v ){ return Dot( tensor , v , v ); }
	template< class Real > Real SquareLength( const SquareMatrix< Real , 2 >& tensor , CotangentVector< Real > v ){ return Dot( tensor , v , v ); }
	template< class Real > Real Angle( const SquareMatrix< Real , 2 >& tensor ,   TangentVector< Real > v1 ,   TangentVector< Real > v2 );
	template< class Real > Real Angle( const SquareMatrix< Real , 2 >& tensor , CotangentVector< Real > v1 , CotangentVector< Real > v2 );
	template< class Real >   TangentVector< Real > Rotate90( const SquareMatrix< Real , 2 >& tensor ,   TangentVector< Real > v );
	template< class Real > CotangentVector< Real > Rotate90( const SquareMatrix< Real , 2 >& tensor , CotangentVector< Real > v );
	template< class Real > SquareMatrix< Real , 2 >   J( const SquareMatrix< Real , 2 >& tensor );
	template< class Real > SquareMatrix< Real , 2 > CoJ( const SquareMatrix< Real , 2 >& tensor );
	template< class Real > Real Area( const SquareMatrix< Real , 2 >& tensor , const Point2D< Real > triangle[3] );
	template< class Real > Real Area( const SquareMatrix< Real , 2 >& tensor , Point2D< Real > v1 , Point2D< Real > v2 , Point2D< Real > v3 );

	// Return the conformal/authalic metric that is closest to the given metric
	template< class Real > SquareMatrix< Real , 2 > MakeConformal( const SquareMatrix< Real , 2 >& sourceTensor , const SquareMatrix< Real , 2 >& targetTensor );
	template< class Real > SquareMatrix< Real , 2 > MakeAuthalic ( const SquareMatrix< Real , 2 >& sourceTensor , const SquareMatrix< Real , 2 >& targetGensor );

	// Compute a symmetric square root of the metric tensor
	template< class Real > SquareMatrix< Real , 2 > TensorRoot( const SquareMatrix< Real , 2 >& tensor );

	///////////////////////////////////////////
	// Right Triangle (with a metric tensor) //
	///////////////////////////////////////////
	template< class Real >
	struct RightTriangle
	{
		// (Some of) the different centers one can define for a triangle
		enum
		{
			CENTER_BARYCENTRIC ,
			CENTER_CIRCUMCENTRIC ,
			CENTER_INCENTRIC ,
			CENTER_ISOGONIC ,
			CENTER_COUNT
		};
		static const char* CenterNames[];

		// Coordinates of the corners/edge-centers/edge-directions
		// The vertices are ordered { (0,0) , (1,0) , (0,1) }.
		// The edges are ordered so that the i-th edge is opposite the i-th vertex (w/ ccw orientation)
		static const Point2D< Real > Corners[];
		static const Point2D< Real > EdgeMidpoints[];
		static const CotangentVector< Real > CornerGradients[];
		static const   TangentVector< Real > EdgeDirections[];

		// Evaluate the scalar/vector/density fields represented by the coefficients.
		// Note that the tensor is required for the conforming basis to define the 90-degree rotation and for the 2-forms to remove the volume
		template< unsigned int BasisType , class V > static V                    EvaluateScalarField        (                                          ConstPointer( V ) coefficients , const Point2D< Real >& position );
		template< unsigned int BasisType , class V > static CotangentVector< V > EvaluateScalarFieldGradient(                                          ConstPointer( V ) coefficients , const Point2D< Real >& position );
		template< unsigned int BasisType , class V > static CotangentVector< V > EvaluateCovectorField      ( const SquareMatrix< Real , 2 >& tensor , ConstPointer( V ) coefficients , const Point2D< Real >& position );
		template< unsigned int BasisType , class V > static V                    EvaluateDensityField       ( const SquareMatrix< Real , 2 >& tensor , ConstPointer( V ) coefficients , const Point2D< Real >& position );

		// Compute the (possibly lumped/weighted) mass matrix
#ifdef NEW_FEM_CODE
		template< unsigned int BasisType > static typename BasisInfoSystem< Real , BasisType >::Matrix        GetMassMatrix( const SquareMatrix< Real , 2 >& tensor );
		template< unsigned int BasisType > static typename BasisInfoSystem< Real , BasisType >::Matrix        GetMassMatrix( const SquareMatrix< Real , 2 >& tensor , const SquareMatrix< Real , 2 >& newTensor );
		template< unsigned int BasisType > static typename BasisInfoSystem< Real , BasisType >::Point GetDiagonalMassMatrix( const SquareMatrix< Real , 2 >& tensor );

		// Compute the differential operator
		template< unsigned int InBasisType , unsigned int OutBasisType > static typename BasisInfoSystem2< Real , InBasisType , OutBasisType >::Matrix GetDMatrix( const SquareMatrix< Real , 2 >& tensor );

		// Compute the masks that indicate which entries could be non-zero
		template< unsigned int BasisType > static typename BasisInfoSystem< Real , BasisType >::MaskMatrix GetMassMask( bool useTensor );
		template< unsigned int InBasisType , unsigned int OutBasisType > static typename BasisInfoSystem2< Real , InBasisType , OutBasisType >::MaskMatrix GetDMask( bool& redundant );
#else // !NEW_FEM_CODE
		template< unsigned int BasisType > static SquareMatrix< Real , BasisInfo< BasisType >::Coefficients >  GetMassMatrix( const SquareMatrix< Real , 2 >& tensor );
		template< unsigned int BasisType > static SquareMatrix< Real , BasisInfo< BasisType >::Coefficients >  GetMassMatrix( const SquareMatrix< Real , 2 >& tensor , const SquareMatrix< Real , 2 >& newTensor );
		template< unsigned int BasisType > static Point< Real , BasisInfo< BasisType >::Coefficients > GetDiagonalMassMatrix( const SquareMatrix< Real , 2 >& tensor );

		// Compute the differential operator
		template< unsigned int InBasisType , unsigned int OutBasisType > static Matrix< Real , BasisInfo< InBasisType >::Coefficients , BasisInfo< OutBasisType >::Coefficients > GetDMatrix( const SquareMatrix< Real , 2 >& tensor );

		// Compute the masks that indicate which entries could be non-zero
		template< unsigned int BasisType > static SquareMatrix< unsigned char , BasisInfo< BasisType >::Coefficients > GetMassMask( bool useTensor );
		template< unsigned int InBasisType , unsigned int OutBasisType > static Matrix< unsigned char , BasisInfo< InBasisType >::Coefficients , BasisInfo< OutBasisType >::Coefficients > GetDMask( bool& redundant );
#endif // NEW_FEM_CODE

		// Compute the integrals
		template< unsigned int BasisType > static Real Integrate( const SquareMatrix< Real , 2 >& tensor , ConstPointer( Real ) linear );
#ifdef NEW_FEM_CODE
		template< unsigned int BasisType > static typename BasisInfoSystem< Real , BasisType >::Point IntegrationDual( const SquareMatrix< Real , 2 >& tensor , ConstPointer( Real ) linear );
		template< unsigned int BasisType > static typename BasisInfoSystem< Real , BasisType >::Point IntegrationDual( const SquareMatrix< Real , 2 >& tensor , ConstPointer( CotangentVector< Real > ) linear );
		template< unsigned int BasisType > static typename BasisInfoSystem< Real , BasisType >::Point IntegrationDual( const SquareMatrix< Real , 2 >& tensor , Point2D< Real > p , Point2D< Real > q );
#else // !NEW_FEM_CODE
		template< unsigned int BasisType > static Point< Real , BasisInfo< BasisType >::Coefficients > IntegrationDual( const SquareMatrix< Real , 2 >& tensor , ConstPointer( Real ) linear );
		template< unsigned int BasisType > static Point< Real , BasisInfo< BasisType >::Coefficients > IntegrationDual( const SquareMatrix< Real , 2 >& tensor , ConstPointer( CotangentVector< Real > ) linear );
		template< unsigned int BasisType > static Point< Real , BasisInfo< BasisType >::Coefficients > IntegrationDual( const SquareMatrix< Real , 2 >& tensor , Point2D< Real > p , Point2D< Real > q );
#endif // NEW_FEM_CODE

		static Point2D< Real > Center( const SquareMatrix< Real , 2 >& tensor , int centerType );
		static Point3D< Real > CenterAreas( const SquareMatrix< Real , 2 >& tensor , int centerType );

		static Point3D< Real > SubTriangleAreas( const SquareMatrix< Real , 2 >& tensor , Point2D< Real > center );

		// Reflect a point across the specified edge
		static Point2D< Real > EdgeReflect( const SquareMatrix< Real , 2 >& tensor , int e , Point2D< Real > p );
	};

	// Position on a mesh
	template< class Real >
	struct SamplePoint
	{
		int tIdx;
		Point2D< Real > p;
		SamplePoint( void ){ ; }
		SamplePoint( int tIdx , const Point2D< Real > p ){ this->tIdx = tIdx , this->p = p; }
	};
	// Position and tangent direction on a mesh
	template< class Real >
	struct HermiteSamplePoint : public SamplePoint< Real >
	{
		using SamplePoint< Real >::tIdx;
		using SamplePoint< Real >::p;
		TangentVector< Real > v;
		HermiteSamplePoint( void ){ ; }
		HermiteSamplePoint( int tIdx , const Point2D< Real >& p , const TangentVector< Real > v=TangentVector< Real >() ){ this->tIdx = tIdx , this->p = p , this->v = v; }
		HermiteSamplePoint( const SamplePoint< Real >& p , const TangentVector< Real > v=TangentVector< Real >() ){ tIdx = p.tIdx , this->p = p.p , this->v = v; }
	};
	// Transformations taking points/directions in one coordinate frame to points/directions in the other
	template< class Real >
	struct CoordinateXForm : public Group< CoordinateXForm< Real > >
	{
		SquareMatrix< Real , 2 > linear;
		Point2D< Real > constant;
		CoordinateXForm( void ){ constant = Point2D< Real >() , linear = SquareMatrix< Real , 2 >::Identity(); }
		// (A,s) * p = A*p + s
		Point2D< Real > operator() ( const Point2D< Real >& p ) const { return linear*p + constant; }
		Point2D< Real > operator * ( const Point2D< Real >& p ) const { return linear*p + constant; }

		// (A,s) * (B,t) * p = (A,s) * (B*p + t) = A*B*p + (A*t + s)
		void SetIdentity( void ){ constant = Point2D< Real >() , linear = SquareMatrix< Real , 2 >::Identity(); }
		void Multiply( const CoordinateXForm& xForm ){ constant += linear * xForm.constant , linear *= xForm.linear; }
		void Invert( void ){ linear = linear.inverse() , constant = - linear * constant; }
	};

	// An abstract class that can be sampled at a point on the mesh and returns a tangent vector
	template< class Real >
	struct TangentVectorField{ virtual TangentVector< Real > operator() ( const SamplePoint< Real >& p ) const = 0; };
	// A realization of the abstract class by a vector field that is constant per triangle.
	template< class Real >
	struct SampledTangentVectorField : public TangentVectorField< Real >
	{
		SampledTangentVectorField( ConstPointer( TangentVector< Real > ) samples ) : _samples( samples ){ ; }
		TangentVector< Real > operator()( const SamplePoint< Real >& p ) const { return _samples[p.tIdx]; }
	protected:
		ConstPointer( TangentVector< Real > ) _samples;
	};

	// A structure that maps from half-edges to edges and vice-versa
	// -- Given an edge index, e, the value returned by operator[e] is an array of integers corresponding to the two half-edges incident on the edge.
	//    (The orientation is given by the first and the second is -1 if the edge is on the boundary).
	// -- Given a half edge index, he, the value of edge( he ) is the associated edge.
	//    (To find if the half-edge is aligned with the orientation of the edge, a reference to a boolean can be passed in, or the isAligned method can be called.)
	// -- Given a half edge index, he, the value of opposite( he ) is the opposite half edge index.
	//    (Or -1 if the he is on the boundary.)
	struct EdgeMap
	{
		EdgeMap( ConstPointer( TriangleIndex ) triangles , size_t tCount );
		~EdgeMap( void ){ FreePointer( _e2he ) ; FreePointer( _he2e ); }
		size_t size( void ) const { return _eCount; }
		// Given an edge index, this returns a pointer to the two half-edges adjacent on that edge. (The first is positively aligned with the edge.)
		const int* operator[]( int eIdx ) const { return _e2he + eIdx*2; }
		// Given a half-edge, this returns the index of the associated edge and indicates whether or not the alignment is consistent
		int edge( int he , bool& aligned ) const { if( _he2e[he]<0 ){ aligned = false ; return -_he2e[he]-1; } else { aligned = true ; return _he2e[he]-1; } }
		int edge( int he ) const { return std::max< int >( _he2e[he] , -_he2e[he] )-1; }
		bool isAligned( int he ) const { return _he2e[he]>0; }
		int opposite( int he ) const { int e = edge(he) ; return he==_e2he[ 2*e+0 ] ? _e2he[ 2*e+1 ] : _e2he[ 2*e+0 ]; }
	protected:
		size_t _eCount;
		Pointer( int ) _e2he;
		Pointer( int ) _he2e;
		template< typename Real > friend struct RiemannianMesh;
	};
	// This structure represents a Riemmanian mesh, with the triangles giving the connectivity and the square (symmetric) matrices giving the metric
	template< class Real >
	struct RiemannianMesh
	{
	protected:
		Pointer( TriangleIndex ) _triangles;
		Pointer( SquareMatrix< Real , 2 > ) _g;
		size_t _tCount , _vCount;
		EdgeMap _edgeMap;
	public:
		ConstPointer( TriangleIndex ) triangles( void     ) const { return _triangles   ; }
		const        TriangleIndex&   triangles( size_t t ) const { return _triangles[t]; }
		size_t tCount( void ) const { return _tCount; }
		size_t vCount( void ) const { return _vCount; }
		size_t eCount( void ) const { return _edgeMap.size(); }

		void edgeVertices( int e , int& v1 , int& v2 ) const;
		bool oppositeEdgeVertices( int e , int& v1 , int& v2 ) const;

		SquareMatrix< Real , 2 >& g( size_t t ){ return _g[t]; }
		const SquareMatrix< Real , 2 >& g( size_t t ) const { return _g[t]; }
		int opposite( int he ) const { return _edgeMap.opposite( he ); }
		RiemannianMesh( Pointer( TriangleIndex ) t , size_t tC );
		~RiemannianMesh( void );

		template< class Vertex > void setMetricFromEmbedding( ConstPointer( Vertex ) vertices );
		void setMetricFromEdgeLengths( ConstPointer( Real ) edgeLengths );
		void setMetricFromSquareEdgeLengths( ConstPointer( Real ) squareEdgeLengths );
		void makeUnitArea( void ){ makeArea( (Real)1. ); }
		void makeArea( Real area );
		Real area( void ) const;
		Real area( int tIdx ) const;
		Real squareEdgeLength( int heIdx ) const;

		// Get the associated index, given the triangle and the index of the element within the triangle
		template< unsigned int BasisType > size_t dimension( void ) const;
		template< unsigned int BasisType > int index( int t , int idx ) const;
		template< unsigned int BasisType > int index( int t , int idx , bool& isAligned ) const;

		template< unsigned int BasisType , class V >                  V   evaluateScalarField        ( ConstPointer( V ) coefficients , const SamplePoint< Real >& p ) const;
		template< unsigned int BasisType , class V > CotangentVector< V > evaluateScalarFieldGradient( ConstPointer( V ) coefficients , const SamplePoint< Real >& p ) const;
		template< unsigned int BasisType , class V > CotangentVector< V > evaluateCovectorField      ( ConstPointer( V ) coefficients , const SamplePoint< Real >& p ) const;

		/////////////////////////
		// Geometric Operators //
		/////////////////////////
		template< unsigned int BasisType > SparseMatrix< Real , int > massMatrix( bool lump=false , ConstPointer( SquareMatrix< Real , 2 > ) newTensors = NullPointer< SquareMatrix< Real , 2 > >() ) const;
		template< unsigned int InBasisType , unsigned int OutBasisType > SparseMatrix< Real , int > dMatrix( void ) const;
		template< unsigned int BasisType , unsigned int PreBasisType , unsigned int PostBasisType > SparseMatrix< Real , int > stiffnessMatrix( ConstPointer( SquareMatrix< Real , 2 > ) newTensors = NullPointer< SquareMatrix< Real , 2 > >() ) const;
		template< unsigned int BasisType > SparseMatrix< Real , int > stiffnessMatrix( void ) const;

		// Integrate the piecewise linear function over the mesh
		Real getIntegral( ConstPointer( Real ) coefficients ) const;
		Real getDotProduct( ConstPointer( Real ) c1 , ConstPointer( Real ) c2 , bool lump ) const;

		CoordinateXForm< Real >  exp( ConstPointer( CoordinateXForm< Real > ) xForms , HermiteSamplePoint< Real >& p , Real eps=(Real)0 , bool noWarning=false ) const;
		CoordinateXForm< Real > flow( ConstPointer( CoordinateXForm< Real > ) xForms , const TangentVectorField< Real >& vf , Real flowTime , SamplePoint< Real >& p , Real minStepSize , Real eps=(Real)0 , std::vector< SamplePoint< Real > >* path=NULL , bool noWarning=false ) const;

		Pointer( CoordinateXForm< Real > ) getCoordinateXForms( void ) const;
		void                               setCoordinateXForms( Pointer( CoordinateXForm< Real > ) xForms ) const;
		CoordinateXForm< Real > getVertexCoordinateXForm( ConstPointer( CoordinateXForm< Real > ) xForms , int t , int v ) const;
		CoordinateXForm< Real > xForm( int he ) const;
		std::vector< int > getVertexCorners( int t , int v ) const;
		Real getVertexConeAngle( int t , int v ) const;

		bool edgeFlip( int e , Real eps=0 );
		bool isVoronoiEdge( int e , Real eps=0 ) const;

		std::vector< SamplePoint< Real > > randomSamples( unsigned int count ) const;
	};
	template< class Real , unsigned int BasisType >
	struct TangentVectorFieldWrapper : public TangentVectorField< Real >
	{
		TangentVectorFieldWrapper( const RiemannianMesh< Real >* mesh , ConstPointer( Real ) coefficients , bool precomputeInverses=false );
		~TangentVectorFieldWrapper( void );
		TangentVector< Real > operator() ( const SamplePoint< Real >& p ) const;
	protected:
		const RiemannianMesh< Real >* _mesh;
		Pointer( SquareMatrix< Real , 2 > ) _gInverse;
		ConstPointer( Real ) _coefficients;
	};
}
template< class Real > const char* FEM::RightTriangle< Real >::CenterNames[] = { "barycentric" , "circumcentric" , "incentric" , "isogonic" };

template< class Real > const Point2D< Real > FEM::RightTriangle< Real >::Corners        [] = { Point2D< Real >(0,0) , Point2D< Real >(1,0) , Point2D< Real >(0,1) };
template< class Real > const Point2D< Real > FEM::RightTriangle< Real >::EdgeMidpoints  [] = { Point2D< Real >((Real)0.5,(Real)0.5) , Point2D< Real >(0,(Real)0.5) , Point2D< Real >((Real)0.5,0) };
template< class Real > const FEM::CotangentVector< Real > FEM::RightTriangle< Real >::CornerGradients[] = { FEM::CotangentVector< Real >( Point2D< Real >(-1,-1) ) , FEM::CotangentVector< Real >( Point2D< Real >( 1, 0) ) , FEM::CotangentVector< Real >( Point2D< Real >( 0,1) ) };
template< class Real > const FEM::  TangentVector< Real > FEM::RightTriangle< Real >::EdgeDirections [] = { FEM::  TangentVector< Real >( Point2D< Real >(-1, 1) ) , FEM::  TangentVector< Real >( Point2D< Real >( 0,-1) ) , FEM::  TangentVector< Real >( Point2D< Real >( 1,0) ) };

template<> const unsigned int FEM::ElementInfo< FEM::ELEMENT_VERTEX   >::ElementsPerTriangle = 3;
template<> const unsigned int FEM::ElementInfo< FEM::ELEMENT_EDGE     >::ElementsPerTriangle = 3;
template<> const unsigned int FEM::ElementInfo< FEM::ELEMENT_TRIANGLE >::ElementsPerTriangle = 1;

#include "FEM.inl"

#ifdef INCLUDE_EXTRA_FEM
#include "FEM_Extra.h"
#endif // INCLUDE_EXTRA_FEM

#endif // FEM_INCLUDED