#include "FunctionBasis/PPolynomial.h"

//////////////////////
// FiniteElements2D //
//////////////////////
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::TensorMatrix(int inDim1,int outDim1,const SparseMatrix<Real>& m1,const SparseMatrix<Real>& m2,SparseMatrix<Real>& m)
{
	int stride;
	if(m1.rowMajor!=m2.rowMajor)	return  false;
	if(m1.rowMajor)	stride=inDim1;
	else			stride=outDim1;
	m.Resize(m1.groups*m2.groups,m1.rowMajor);
	for(int i=0;i<m1.groups;i++)
		for(int j=0;j<m2.groups;j++)
		{
			int idx=i+j*m1.groups;
			m.SetGroupSize(idx,m1.groupSizes[i]*m2.groupSizes[j]);
			for(int ii=0;ii<m1.groupSizes[i];ii++)
				for(int jj=0;jj<m2.groupSizes[j];jj++)
				{
					int iidx=ii+jj*m1.groupSizes[i];
					m.m_ppElements[idx][iidx].N=m1.m_ppElements[i][ii].N+m2.m_ppElements[j][jj].N*stride;
					m.m_ppElements[idx][iidx].Value=m1.m_ppElements[i][ii].Value*m2.m_ppElements[j][jj].Value;
				}
		}

	return true;
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
template<class Data>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::TensorMatrixMultiply(int inDim1,int outDim1,int inDim2,int outDim2,
																			  const SparseMatrix<Real>& m1,const SparseMatrix<Real>& m2,
																			  const Vector<Data>& in,Vector<Data>& out)
{
	int stride;
	if(m1.rowMajor!=m2.rowMajor || in.Dimensions()!=inDim1*inDim2)	return  false;
	if(m1.rowMajor)	stride= inDim1;
	else			stride=outDim1;
	out.Resize(outDim1*outDim2);
	for(int i=0;i<m1.groups;i++)
		for(int j=0;j<m2.groups;j++)
		{
			int idx=i+j*m1.groups;
			for(int ii=0;ii<m1.groupSizes[i];ii++)
			{
				if(m1.m_ppElements[i][ii].N<0)	continue;
				for(int jj=0;jj<m2.groupSizes[j];jj++)
				{
					if(m2.m_ppElements[j][jj].N<0)	continue;
					int N=m1.m_ppElements[i][ii].N+m2.m_ppElements[j][jj].N*stride;
					Real Value=m1.m_ppElements[i][ii].Value*m2.m_ppElements[j][jj].Value;
					if(m1.rowMajor)	out[idx]+=in[N]*Value;
					else			out[N]+=in[idx]*Value;
				}
			}
		}
	return true;
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
template<class Data>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::TensorMatrixMultiply(int inDim1,int outDim1,int inDim2,int outDim2,
																			  const SparseMatrix<Real>& m1,const SparseMatrix<Real>& m2,
																			  const Data* in,Data* out)
{
	int stride;
	if(m1.rowMajor!=m2.rowMajor)	return  false;
	if(m1.rowMajor)	stride= inDim1;
	else			stride=outDim1;
	for(int i=0;i<m1.groups;i++)
		for(int j=0;j<m2.groups;j++)
		{
			int idx=i+j*m1.groups;
			for(int ii=0;ii<m1.groupSizes[i];ii++)
			{
				if(m1.m_ppElements[i][ii].N<0)	continue;
				for(int jj=0;jj<m2.groupSizes[j];jj++)
				{
					if(m2.m_ppElements[j][jj].N<0)	continue;
					int N=m1.m_ppElements[i][ii].N+m2.m_ppElements[j][jj].N*stride;
					Real Value=m1.m_ppElements[i][ii].Value*m2.m_ppElements[j][jj].Value;
					if(m1.rowMajor)	out[idx]+=in[N]*Value;
					else			out[N]+=in[idx]*Value;
				}
			}
		}
	return true;
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::TensorMatrixMultiply(int inDim1,int outDim1,int inDim2,int outDim2,
																			  const SparseMatrix<Real>& m1,const SparseMatrix<Real>& m2,
																			  const Real* in,Real* out,int Channels)
{
	int stride;
	if(m1.rowMajor!=m2.rowMajor)	return  false;
	if(m1.rowMajor)	stride= inDim1;
	else			stride=outDim1;
	for(int c=0;c<Channels;c++)
		for(int i=0;i<m1.groups;i++)
			for(int j=0;j<m2.groups;j++)
			{
				int idx=i+j*m1.groups;
				for(int ii=0;ii<m1.groupSizes[i];ii++)
				{
					if(m1.m_ppElements[i][ii].N<0)	continue;
					for(int jj=0;jj<m2.groupSizes[j];jj++)
					{
						if(m2.m_ppElements[j][jj].N<0)	continue;
						Real Value=m1.m_ppElements[i][ii].Value*m2.m_ppElements[j][jj].Value;
#if 1
						int myIdx=i+j*m1.groups*Channels+c*m1.groups;
						int N=m1.m_ppElements[i][ii].N+m2.m_ppElements[j][jj].N*stride*Channels+c*stride;
						if(m1.rowMajor)	out[myIdx]+=in[N]*Value;
						else			out[N]+=in[myIdx]*Value;
#else
						int N=m1.m_ppElements[i][ii].N+m2.m_ppElements[j][jj].N*stride*Channels+c*stride;
						if(m1.rowMajor)	out[idx]+=in[N]*Value;
						else			out[N]+=in[idx]*Value;
#endif
					}
				}
			}
	return true;
}


template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::IsUpSamplable(int lowD1,int lowD2,int& highD1,int& highD2)
{
	return FiniteElements1D<Real,Type1,Degree1>::IsUpSamplable(lowD1,highD1) && FiniteElements1D<Real,Type2,Degree2>::IsUpSamplable(lowD2,highD2);
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::IsDownSamplable(int highD1,int highD2,int& lowD1,int& lowD2)
{
	return FiniteElements1D<Real,Type1,Degree1>::IsDownSamplable(highD1,lowD1) && FiniteElements1D<Real,Type2,Degree2>::IsDownSamplable(highD2,lowD2);
}

template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::DotProductMatrix(int dim1,int dim2,SparseMatrix<Real>& m)
{
	return DotProductMatrix(dim1,dim2,m,0,FiniteElements1D<Real,Type1,Degree1>::DomainSize(dim1),0,FiniteElements1D<Real,Type2,Degree2>::DomainSize(dim2));
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::DotProductMatrix(int dim1,int dim2,SparseMatrix<Real>& m,Real min1,Real max1,Real min2,Real max2)
{
	SparseMatrix<Real> dot1,d2Dot1,dot2,d2Dot2;
	if(	!FiniteElements1D<Real,Type1,Degree1>::DotProductMatrix(dim1,dot1,min1,max1) ||
		!FiniteElements1D<Real,Type2,Degree2>::DotProductMatrix(dim2,dot2,min2,max2) )
		return false;
	return TensorMatrix(dim1,dim1,dot1,dot2,m);
}


template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::LaplacianMatrix(int dim1,int dim2,SparseMatrix<Real>& m,bool weakForm,bool negate)
{
	return LaplacianMatrix(dim1,dim2,m,0,Real(FiniteElements1D<Real,Type1,Degree1>::DomainSize(dim1)),0,Real(FiniteElements1D<Real,Type2,Degree2>::DomainSize(dim2)),weakForm,negate);
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::LaplacianMatrix(int dim1,int dim2,SparseMatrix<Real>& m,Real min1,Real max1,Real min2,Real max2,bool weakForm,bool negate)
{
	m.Resize(dim1*dim2);
	for(int i=0;i<dim1*dim2;i++) m.SetGroupSize(i,(2*Degree1+1)*(2*Degree2+1));
	SparseMatrix<Real> dot1,d2Dot1,dot2,d2Dot2;
	FiniteElements1D<Real,Type1,Degree1>::DotProductMatrix(dim1,dot1,min1,max1);
	FiniteElements1D<Real,Type2,Degree2>::DotProductMatrix(dim2,dot2,min2,max2);
	FiniteElements1D<Real,Type1,Degree1>::LaplacianMatrix(dim1,d2Dot1,min1,max1,weakForm,negate);
	FiniteElements1D<Real,Type2,Degree2>::LaplacianMatrix(dim2,d2Dot2,min2,max2,weakForm,negate);

	for(int y=0;y<dim2;y++)
		for(int x=0;x<dim1;x++)
			for(int yy=0;yy<(2*Degree2+1);yy++)
				for(int xx=0;xx<(2*Degree1+1);xx++)
				{
					m.m_ppElements[x+y*dim1][xx+yy*(2*Degree1+1)].N=
						dot1.m_ppElements[x][xx].N+dot2.m_ppElements[y][yy].N*dim1;
					m.m_ppElements[x+y*dim1][xx+yy*(2*Degree1+1)].Value=
						  dot1.m_ppElements[x][xx].Value*d2Dot2.m_ppElements[y][yy].Value+
						d2Dot1.m_ppElements[x][xx].Value*  dot2.m_ppElements[y][yy].Value;
				}
	return true;
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::LaplacianMatrix(int dim1,int dim2,double iWeight,SparseMatrix<Real>& m,bool weakForm,bool negate)
{
	return LaplacianMatrix(dim1,dim2,iWeight,m,0,Real(FiniteElements1D<Real,Type1,Degree1>::DomainSize(dim1)),0,Real(FiniteElements1D<Real,Type2,Degree2>::DomainSize(dim2)),weakForm,negate);
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::LaplacianMatrix(int dim1,int dim2,double iWeight,SparseMatrix<Real>& m,Real min1,Real max1,Real min2,Real max2,bool weakForm,bool negate)
{
	m.Resize(dim1*dim2);
	for(int i=0;i<dim1*dim2;i++) m.SetGroupSize(i,(2*Degree1+1)*(2*Degree2+1));
	SparseMatrix<Real> dot1,d2Dot1,dot2,d2Dot2;
	FiniteElements1D<Real,Type1,Degree1>::DotProductMatrix(dim1,dot1,min1,max1);
	FiniteElements1D<Real,Type2,Degree2>::DotProductMatrix(dim2,dot2,min2,max2);
	FiniteElements1D<Real,Type1,Degree1>::LaplacianMatrix(dim1,d2Dot1,min1,max1,weakForm,negate);
	FiniteElements1D<Real,Type2,Degree2>::LaplacianMatrix(dim2,d2Dot2,min2,max2,weakForm,negate);

	for(int y=0;y<dim2;y++)
		for(int x=0;x<dim1;x++)
			for(int yy=0;yy<(2*Degree2+1);yy++)
				for(int xx=0;xx<(2*Degree1+1);xx++)
				{
					m.m_ppElements[x+y*dim1][xx+yy*(2*Degree1+1)].N=
						dot1.m_ppElements[x][xx].N+dot2.m_ppElements[y][yy].N*dim1;
					m.m_ppElements[x+y*dim1][xx+yy*(2*Degree1+1)].Value=
						  dot1.m_ppElements[x][xx].Value*d2Dot2.m_ppElements[y][yy].Value+
						d2Dot1.m_ppElements[x][xx].Value*  dot2.m_ppElements[y][yy].Value+
						(dot1.m_ppElements[x][xx].Value*dot2.m_ppElements[y][yy].Value)*iWeight;
				}
	return true;
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::DivergenceStencil(int dim1,int dim2,FullDivergenceStencil& s)
{
	typename FiniteElements1D< Real , Type1 , Degree1 >::template DotProduct<Type1,Degree1>::FullDotProductStencil dot1;
	typename FiniteElements1D< Real , Type2 , Degree2 >::template DotProduct<Type2,Degree2>::FullDotProductStencil dot2;
	typename FiniteElements1D< Real , DERIVATIVE(Type1) , Degree1-1 >::template DotProduct< Type1 , Degree1 >::FullDotProductStencil d1;
	typename FiniteElements1D< Real , DERIVATIVE(Type2) , Degree2-1 >::template DotProduct< Type2 , Degree2 >::FullDotProductStencil d2;

	FiniteElements1D< Real , Type1 , Degree1 >::template DotProduct< Type1 , Degree1 >::DotProductStencil( dim1 , dot1 , 0 , 0 );
	FiniteElements1D< Real , Type2 , Degree2 >::template DotProduct< Type2 , Degree2 >::DotProductStencil( dim2 , dot2 , 0 , 0 );
	FiniteElements1D< Real , DERIVATIVE(Type1) , Degree1-1 >::template DotProduct< Type1 , Degree1 >::DotProductStencil(FiniteElements1D< Real , DERIVATIVE(Type1) , Degree1-1 >::Dimension( FiniteElements1D< Real , Type1 , Degree1 >::DomainSize( dim1 ) ) , d1 , 0 , 1 );
	FiniteElements1D< Real , DERIVATIVE(Type2) , Degree2-1 >::template DotProduct< Type2 , Degree2 >::DotProductStencil(FiniteElements1D< Real , DERIVATIVE(Type2) , Degree2-1 >::Dimension( FiniteElements1D< Real , Type2 , Degree2 >::DomainSize( dim2 ) ) , d2 , 0 , 1 );

	for(int i=0;i<=2*Degree1;i++)
		for(int j=0;j<=2*Degree2;j++)
		{
			for( int k=0 ; k<2*Degree1 ; k++ ) for( int l=0 ; l<=2*Degree2 ; l++ )
				s.caseTable[i][j].values1[k][l]=d1.caseTable[i].values[k]*dot2.caseTable[j].values[l];
			for( int k=0 ; k<=2*Degree1 ; k++ ) for( int l=0 ; l<2*Degree2 ; l++ )
				s.caseTable[i][j].values2[k][l]=dot1.caseTable[i].values[k]*d2.caseTable[j].values[l];
		}
	return true;
}

template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::BiLaplacianStencil( int dim1 , int dim2 , FullMatrixStencil& s , bool weakForm , bool negate )
{
	typename FiniteElements1D< Real , Type1 , Degree1 >::template DotProduct< Type1 , Degree1 >::FullDotProductStencil dot1 , biLap1;
	typename FiniteElements1D< Real , Type2 , Degree2 >::template DotProduct< Type2 , Degree2 >::FullDotProductStencil dot2 , biLap2;

	FiniteElements1D< Real , Type1 , Degree1 >::template DotProduct< Type1 , Degree1 >::DotProductStencil( dim1 , dot1 , 0 , 0 );
	if( weakForm ) FiniteElements1D< Real , Type1 , Degree1 >::template DotProduct< Type1 , Degree1 >::DotProductStencil( dim1 , biLap1 , 2 , 2 , negate );
	else           FiniteElements1D< Real , Type1 , Degree1 >::template DotProduct< Type1 , Degree1 >::DotProductStencil( dim1 , biLap1 , 4 , 0 , negate );
	FiniteElements1D< Real , Type2 , Degree2 >::template DotProduct< Type2 , Degree2 >::DotProductStencil( dim2 , dot2 , 0 , 0 );
	if( weakForm ) FiniteElements1D< Real , Type2 , Degree2 >::template DotProduct< Type2 , Degree2 >::DotProductStencil( dim2 , biLap2 , 2 , 2 , negate );
	else           FiniteElements1D< Real , Type2 , Degree2 >::template DotProduct< Type2 , Degree2 >::DotProductStencil( dim2 , biLap2 , 4 , 0 , negate );

	for( int i=0 ; i<=2*Degree1 ; i++ )
		for( int j=0 ; j<=2*Degree2 ; j++ )
			for( int k=0 ; k<=2*Degree1 ; k++ )
				for( int l=0 ; l<=2*Degree2 ; l++ )
					s.caseTable[i][j].values[k][l] =
					dot1.caseTable[i].values[k]*biLap2.caseTable[j].values[l]+
					biLap1.caseTable[i].values[k]*dot2.caseTable[j].values[l];
	return true;
}

template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::LaplacianStencil(int dim1,int dim2,FullMatrixStencil& s,bool weakForm,bool negate)
{
	typename FiniteElements1D< Real , Type1 , Degree1 >::template DotProduct< Type1 , Degree1 >::FullDotProductStencil dot1,lap1;
	typename FiniteElements1D< Real , Type2 , Degree2 >::template DotProduct< Type2 , Degree2 >::FullDotProductStencil dot2,lap2;

	FiniteElements1D< Real , Type1 , Degree1 >::template DotProduct< Type1 , Degree1 >::DotProductStencil(dim1,dot1,0,0);
	if(weakForm)	FiniteElements1D< Real , Type1 , Degree1>::template DotProduct< Type1 , Degree1 >::DotProductStencil(dim1,lap1,1,1,negate);
	else			FiniteElements1D< Real , Type1 , Degree1>::template DotProduct< Type1 , Degree1 >::DotProductStencil(dim1,lap1,2,0,negate);
	FiniteElements1D<Real,Type2,Degree2>::template DotProduct<Type2,Degree2>::DotProductStencil(dim2,dot2,0,0);
	if(weakForm)	FiniteElements1D<Real,Type2,Degree2>::template DotProduct<Type2,Degree2>::DotProductStencil(dim2,lap2,1,1,negate);
	else			FiniteElements1D<Real,Type2,Degree2>::template DotProduct<Type2,Degree2>::DotProductStencil(dim2,lap2,2,0,negate);

	for(int i=0;i<=2*Degree1;i++)
		for(int j=0;j<=2*Degree2;j++)
			for(int k=0;k<=2*Degree1;k++)
				for(int l=0;l<=2*Degree2;l++)
					s.caseTable[i][j].values[k][l]=
					dot1.caseTable[i].values[k]*lap2.caseTable[j].values[l]+
					lap1.caseTable[i].values[k]*dot2.caseTable[j].values[l];
	return true;
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::DotProductStencil(int dim1,int dim2,FullMatrixStencil& s)
{
	typename FiniteElements1D< Real , Type1 , Degree1 >::template DotProduct< Type1 , Degree1 >::FullDotProductStencil dot1 , dot2;

	FiniteElements1D<Real,Type1,Degree1>::template DotProduct< Type1 , Degree1 >::DotProductStencil( dim1 , dot1 , 0 , 0 );
	FiniteElements1D<Real,Type2,Degree2>::template DotProduct< Type2 , Degree2 >::DotProductStencil( dim2 , dot2 , 0 , 0 );

	for(int i=0;i<=2*Degree1;i++)
		for(int j=0;j<=2*Degree2;j++)
			for(int k=0;k<=2*Degree1;k++)
				for(int l=0;l<=2*Degree2;l++)
					s.caseTable[i][j].values[k][l]=dot1.caseTable[i].values[k]*dot2.caseTable[j].values[l];
	return true;
}

template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::ProlongationStencil(int lowD1,int lowD2,FullProlongationStencil& s,int& highD1,int& highD2)
{
	typename FiniteElements1D<Real,Type1,Degree1>::FullProlongationStencil stencil1;
	typename FiniteElements1D<Real,Type2,Degree2>::FullProlongationStencil stencil2;

	FiniteElements1D<Real,Type1,Degree1>::ProlongationStencil(lowD1,stencil1,highD1);
	FiniteElements1D<Real,Type1,Degree2>::ProlongationStencil(lowD2,stencil2,highD2);

	for(int i=0;i<=2*Degree1;i++)
		for(int j=0;j<=2*Degree2;j++)
			for(int k=0;k<Degree1+2;k++)
				for(int l=0;l<Degree2+2;l++)
					s.caseTable[i][j].values[k][l]=stencil1.caseTable[i].values[k]*stencil2.caseTable[j].values[l];
	return true;
}
#if NEW_LAPLACIAN_CODE
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::RestrictionStencil(int highD1,int highD2,FullRestrictionStencil& s,int& lowD1,int& lowD2)
{
	typename FiniteElements1D<Real,Type1,Degree1>::FullRestrictionStencil stencil1;
	typename FiniteElements1D<Real,Type2,Degree2>::FullRestrictionStencil stencil2;

	FiniteElements1D<Real,Type1,Degree1>::RestrictionStencil(highD1,stencil1,lowD1);
	FiniteElements1D<Real,Type1,Degree2>::RestrictionStencil(highD2,stencil2,lowD2);

	for(int i=0;i<2*Degree1+2;i++)
		for(int j=0;j<2*Degree2+2;j++)
			for(int k=0;k<(Degree1+3)>>1;k++)
				for(int l=0;l<(Degree2+3)>>1;l++)
					s.caseTable[i][j].values[k][l]=stencil1.caseTable[i].values[k]*stencil2.caseTable[j].values[l];
	return true;
}
#endif // NEW_LAPLACIAN_CODE
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
void FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::StripDiagonal(const SparseMatrix<Real>& M,Vector<MatrixEntry<Real> >& D,bool clearDiagonal)
{
	int idx=((2*Degree1+1)*(2*Degree2+1))/2;
	D.Resize(M.groups);
	for(int i=0;i<M.groups;i++)
	{
		D(i)=M.m_ppElements[i][idx];
		if(clearDiagonal)	M.m_ppElements[i][idx].Value=0;
	}
}

template<class Real,int Type1,int Degree1,int Type2,int Degree2>
template<class Data>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::UpSample(Vector<Data>& in,Vector<Data>& out,int lowD1,int lowD2,int& highD1,int& highD2)
{
	SparseMatrix<Real> up1,up2;
	if(	!FiniteElements1D<Real,Type1,Degree1>::UpSampleMatrix(up1,lowD1,highD1) ||
		!FiniteElements1D<Real,Type2,Degree2>::UpSampleMatrix(up2,lowD2,highD2) )
		return false;
	return TensorMatrixMultiply(lowD1,highD1,lowD2,highD2,up1,up2,in,out);
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
template<class Data>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::DownSample(Vector<Data>& in,Vector<Data>& out,
																	int highD1,int highD2,int& lowD1,int& lowD2)
{
	SparseMatrix<Real> down1,down2;
	if(	!FiniteElements1D<Real,Type1,Degree1>::DownSampleMatrix(down1,highD1,lowD1) ||
		!FiniteElements1D<Real,Type2,Degree2>::DownSampleMatrix(down2,highD2,lowD2) )
		return false;
	return TensorMatrixMultiply(highD1,lowD1,highD2,lowD2,down1,down2,in,out);
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
template<class Data>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::UpSample(const Data* in,Data* out,int lowD1,int lowD2,int& highD1,int& highD2)
{
	SparseMatrix<Real> up1,up2;
	if(	!FiniteElements1D<Real,Type1,Degree1>::UpSampleMatrix(up1,lowD1,highD1) ||
		!FiniteElements1D<Real,Type2,Degree2>::UpSampleMatrix(up2,lowD2,highD2) )
		return false;
	return TensorMatrixMultiply(lowD1,highD1,lowD2,highD2,up1,up2,in,out);
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::UpSample(const Real* in,Real* out,int lowD1,int lowD2,int& highD1,int& highD2,int Channels)
{
	SparseMatrix<Real> up1,up2;
	if(	!FiniteElements1D<Real,Type1,Degree1>::UpSampleMatrix(up1,lowD1,highD1) ||
		!FiniteElements1D<Real,Type2,Degree2>::UpSampleMatrix(up2,lowD2,highD2) )
		return false;
	return TensorMatrixMultiply(lowD1,highD1,lowD2,highD2,up1,up2,in,out,Channels);
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
template<class Data>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::DownSample(const Data* in,Data* out,int highD1,int highD2,int& lowD1,int& lowD2)
{
	SparseMatrix<Real> down1,down2;
	if(	!FiniteElements1D<Real,Type1,Degree1>::DownSampleMatrix(down1,highD1,lowD1) ||
		!FiniteElements1D<Real,Type2,Degree2>::DownSampleMatrix(down2,highD2,lowD2) )
		return false;
	return TensorMatrixMultiply(highD1,lowD1,highD2,lowD2,down1,down2,in,out);
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::DownSample(const Real* in,Real* out,int highD1,int highD2,int& lowD1,int& lowD2,int Channels)
{
	SparseMatrix<Real> down1,down2;
	if(	!FiniteElements1D<Real,Type1,Degree1>::DownSampleMatrix(down1,highD1,lowD1) ||
		!FiniteElements1D<Real,Type2,Degree2>::DownSampleMatrix(down2,highD2,lowD2) )
		return false;
	return TensorMatrixMultiply(highD1,lowD1,highD2,lowD2,down1,down2,in,out,Channels);
}


template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::UpSampleMatrix(SparseMatrix<Real>& M,
																		int lowD1,int lowD2,int& highD1,int& highD2)
{
	SparseMatrix<Real> matrix1,matrix2;
	if(	!FiniteElements1D<Real,Type1,Degree1>::UpSampleMatrix(matrix1,lowD1,highD1) ||
		!FiniteElements1D<Real,Type2,Degree2>::UpSampleMatrix(matrix2,lowD2,highD2) )
		return false;
	return TensorMatrix(lowD1,highD1,matrix1,matrix2,M);
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::DownSampleMatrix(SparseMatrix<Real>& M,
																		  int highD1,int highD2,int& lowD1,int& lowD2)
{
	SparseMatrix<Real> matrix1,matrix2;
	if(	!FiniteElements1D<Real,Type1,Degree1>::DownSampleMatrix(matrix1,highD1,lowD1) ||
		!FiniteElements1D<Real,Type2,Degree2>::DownSampleMatrix(matrix2,highD2,lowD2) )
		return false;
	return TensorMatrix(highD1,lowD1,matrix1,matrix2,M);
}


template<class Real,int Type1,int Degree1,int Type2,int Degree2>
template<class Data>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::Gradient(const Vector<Data>& in,Vector<Data>& dX,Vector<Data>& dY,
																  int w,int h,int& ww,int& hh)
{
	SparseMatrix< Real > identity1,identity2,derivative1,derivative2;
	ww = FiniteElements1D< Real , Type1 , Degree1 >::DerivativeMatrix( w , derivative1 );
	hh = FiniteElements1D< Real , Type2 , Degree2 >::DerivativeMatrix( h , derivative2 );
	FiniteElements1D< Real , Type1 , Degree1 >::IdentityMatrix( w , identity1 , derivative1.rowMajor );
	FiniteElements1D< Real , Type2 , Degree2 >::IdentityMatrix( h , identity2 , derivative2.rowMajor );

	return
		FiniteElements2D< Real , Type1 , Degree1 >::TensorMatrixMultiply( w , ww , h ,h , derivative1 , identity2 , in , dX ) &&
		FiniteElements2D< Real , Type2 , Degree2 >::TensorMatrixMultiply( w , w , h , hh , identity1 , derivative2 , in , dY );
}
template<class Real,int Type1,int Degree1,int Type2,int Degree2>
template<class Data>
bool FiniteElements2D<Real,Type1,Degree1,Type2,Degree2>::Laplacian(const Vector<Data>& dX,const Vector<Data>& dY,Vector<Data>& out,
																   int w,int h,int ww,int hh)
{
#if 1
	FullDivergenceStencil stencil;
	DivergenceStencil(w,h,stencil);

	out.Resize(w*h);

	for(int i=0;i<h;i++)
	{
		int ii;
		if(i<Degree2)			ii=i;
		else if(i>=h-Degree2)	ii=2*Degree2+(i-(h-1));
		else					ii=Degree2;
		int dI=i+FiniteElements1D< Real , DERIVATIVE(Type2) , Degree2-1 >::template DotProduct< Type2 , Degree2 >::Helper::StartOffset();
		int I =i+FiniteElements1D< Real , Type2 , Degree2 >::template DotProduct< Type2 , Degree2 >::Helper::StartOffset();
		for(int j=0;j<w;j++)
		{
			Data temp=0;
			int jj;
			if(j<Degree1)			jj=j;
			else if(j>=w-Degree2)	jj=2*Degree1+(j-(w-1));
			else					jj=Degree1;
			int dJ=j+FiniteElements1D< Real , DERIVATIVE(Type1) , Degree1-1 >::template DotProduct< Type1 , Degree1 >::Helper::StartOffset();
			int J =j+FiniteElements1D< Real , Type1 , Degree1 >::template DotProduct< Type1 , Degree1 >::Helper::StartOffset();

			// Partial w.r.t minor index
			for(int di=0;di<2*Degree2;di++)
			{
				if(dI+di<0 || dI+di>=hh)	continue;
				for(int dj=0;dj<=2*Degree1;dj++)
				{
					if(J+dj<0 || J+dj>=w)	continue;
					temp+=dY[(dI+di)*w+(J+dj)]*stencil.caseTable[jj][ii].values2[dj][di];
				}
			}
			// Partial w.r.t major index
			for(int di=0;di<=2*Degree2;di++)
			{
				if(I+di<0 || I+di>=h)	continue;
				for(int dj=0;dj<2*Degree1;dj++)
				{
					if(dJ+dj<0 || dJ+dj>=ww)	continue;
					temp+=dX[(I+di)*ww+(dJ+dj)]*stencil.caseTable[jj][ii].values1[dj][di];
				}
			}
			out[i*w+j]=temp;
		}
	}
#else
	SparseMatrix<Real> dotX,dotY,dDotX,dDotY;
	Vector<Data> tempX,tempY;

	FiniteElements1D<Real,DERIVATIVE(Type1),Degree1-1>::DotProduct<Type1,Degree1>::DerivativeMatrix(ww,dDotX,0,1);
	FiniteElements1D<Real,DERIVATIVE(Type2),Degree2-1>::DotProduct<Type2,Degree2>::DerivativeMatrix(hh,dDotY,0,1);
	FiniteElements1D<Real,Type1,Degree1>::DotProductMatrix(w,dotX);
	FiniteElements1D<Real,Type2,Degree2>::DotProductMatrix(h,dotY);
	if(	!FiniteElements2D<Real,Type1,Degree1>::TensorMatrixMultiply(ww,w,h,h,dDotX,dotY,dX,tempX) ||
		!FiniteElements2D<Real,Type2,Degree2>::TensorMatrixMultiply(w,w,hh,h,dotX,dDotY,dY,tempY) )
		return false;
	out.Resize(tempX.Dimensions());
	for(size_t i=0;i<tempX.Dimensions();i++)	out[i]=tempX[i]+tempY[i];
#endif
	return true;
}
