/*!\file GaussTetra.c
 * \brief: implementation of the GaussTetra object
 */

#include <math.h>
#include "./GaussTetra.h"
#include "../../shared/io/Print/Print.h"
#include "../../shared/Exceptions/exceptions.h"
#include "../../shared/MemOps/MemOps.h"
#include "../../shared/Enum/Enum.h"
#include "../../shared/Numerics/GaussPoints.h"
#include "../../shared/Numerics/constants.h"

/*GaussTetra constructors and destructors:*/
GaussTetra::GaussTetra(){/*{{{*/

	numgauss=-1;

	weights=NULL;
	coords1=NULL;
	coords2=NULL;
	coords3=NULL;
	coords4=NULL;

	weight=UNDEF;
	coord1=UNDEF;
	coord2=UNDEF;
	coord3=UNDEF;
	coord4=UNDEF;
}
/*}}}*/
GaussTetra::GaussTetra(int order){/*{{{*/

	/*Get gauss points*/
	GaussLegendreTetra(&numgauss,&coords1,&coords2,&coords3,&coords4,&weights,order);

	/*Rescale weights if necessary*/
	IssmDouble sumweights = 0.;
	for(int i=0;i<numgauss;i++) sumweights += this->weights[i];
	if(sumweights==1.){
		for(int i=0;i<numgauss;i++) this->weights[i] = this->weights[i]/6.; /*rescale volume to 1/6*/
	}

	/*Check final weights in debugging mode*/
	#ifdef _ISSM_DEBUG_
	sumweights = 0.; for(int i=0;i<numgauss;i++) sumweights += this->weights[i];
	_assert_(sumweights>1./6.-1e-10);
	_assert_(sumweights<1./6.+1e-10);
	#endif

	/*Initialize static fields as undefinite*/
	weight=UNDEF;
	coord1=UNDEF;
	coord2=UNDEF;
	coord3=UNDEF;
}
/*}}}*/
GaussTetra::GaussTetra(int index1,int index2,int index3,int order){/*{{{*/

	/*Basal Tria*/
	if(index1==0 && index2==1 && index3==2){
		GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
		coords4=xNew<IssmDouble>(numgauss);
		for(int i=0;i<numgauss;i++) coords4[i]=0.;
	}
	else if(index1==0 && index2==3 && index3==1){
		GaussLegendreTria(&numgauss,&coords1,&coords2,&coords4,&weights,order);
		coords3=xNew<IssmDouble>(numgauss);
		for(int i=0;i<numgauss;i++) coords3[i]=0.;
	}
	else if(index1==1 && index2==3 && index3==2){
		GaussLegendreTria(&numgauss,&coords2,&coords3,&coords4,&weights,order);
		coords1=xNew<IssmDouble>(numgauss);
		for(int i=0;i<numgauss;i++) coords1[i]=0.;
	}
	else if(index1==0 && index2==2 && index3==3){
		GaussLegendreTria(&numgauss,&coords1,&coords3,&coords4,&weights,order);
		coords2=xNew<IssmDouble>(numgauss);
		for(int i=0;i<numgauss;i++) coords2[i]=0.;
	}
	else{
		_error_(index1 <<" "<<index2 <<" "<<index3 <<" Not supported yet");
	}
}
/*}}}*/
GaussTetra::~GaussTetra(){/*{{{*/
	xDelete<IssmDouble>(weights);
	xDelete<IssmDouble>(coords1);
	xDelete<IssmDouble>(coords2);
	xDelete<IssmDouble>(coords3);
	xDelete<IssmDouble>(coords4);
}
/*}}}*/

/*Methods*/
int GaussTetra::begin(void){/*{{{*/

	/*Check that this has been initialized*/
	_assert_(numgauss>0);
	_assert_(weights);
	_assert_(coords1);
	_assert_(coords2);
	_assert_(coords3);
	_assert_(coords4);

	/*return first gauss index*/
	return 0;
}
/*}}}*/
void GaussTetra::Echo(void){/*{{{*/

	_printf_("GaussTetra:\n");
	_printf_("   numgauss: " << numgauss << "\n");

	if (weights){
	 _printf_("   weights = ["); 
	 for(int i=0;i<numgauss;i++) _printf_(" " << weights[i] << "\n");
	 _printf_("]\n");
	}
	else _printf_("weights = NULL\n");
	if (coords1){
	 _printf_("   coords1 = ["); 
	 for(int i=0;i<numgauss;i++) _printf_(" " << coords1[i] << "\n");
	 _printf_("]\n");
	}
	else _printf_("coords1 = NULL\n");
	if (coords2){
	 _printf_("   coords2 = ["); 
	 for(int i=0;i<numgauss;i++) _printf_(" " << coords2[i] << "\n");
	 _printf_("]\n");
	}
	else _printf_("coords2 = NULL\n");
	if (coords3){
	 _printf_("   coords3 = ["); 
	 for(int i=0;i<numgauss;i++) _printf_(" " << coords3[i] << "\n");
	 _printf_("]\n");
	}
	else _printf_("coords3 = NULL\n");
	if (coords4){
		_printf_("   coords4 = ["); 
		for(int i=0;i<numgauss;i++) _printf_(" " << coords4[i] << "\n");
		_printf_("]\n");
	}
	else _printf_("coords4 = NULL\n");

	_printf_("   weight = " << weight << "\n");
	_printf_("   coord1 = " << coord1 << "\n");
	_printf_("   coord2 = " << coord2 << "\n");
	_printf_("   coord3 = " << coord3 << "\n");
	_printf_("   coord4 = " << coord4 << "\n");

}
/*}}}*/
int GaussTetra::end(void){/*{{{*/

	/*Check that this has been initialized*/
	_assert_(numgauss>0);
	_assert_(weights);
	_assert_(coords1);
	_assert_(coords2);
	_assert_(coords3);
	_assert_(coords4);

	/*return last gauss index +1*/
	return numgauss;
}
/*}}}*/
int GaussTetra::Enum(void){/*{{{*/
	return GaussTetraEnum;
}
/*}}}*/
void GaussTetra::GaussPoint(int ig){/*{{{*/

	/*Check input in debugging mode*/
	 _assert_(ig>=0 && ig< numgauss);

	 /*update static arrays*/
	 weight=weights[ig];
	 coord1=coords1[ig];
	 coord2=coords2[ig];
	 coord3=coords3[ig];
	 coord4=coords4[ig];

}
/*}}}*/
void GaussTetra::GaussNode(int finiteelement,int iv){/*{{{*/

	/*in debugging mode: check that the default constructor has been called*/
	_assert_(numgauss==-1);

	/*update static arrays*/
	switch(finiteelement){
		case P1Enum: case P1DGEnum:
			switch(iv){
				case 0: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
				case 1: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
				case 2: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
				case 3: coord1=0.; coord2=0.; coord3=0.; coord4=1.; break;
				default: _error_("node index should be in [0 3]");
			}
			break;
		case P1bubbleEnum: case P1bubblecondensedEnum:
			switch(iv){
				case 0: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
				case 1: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
				case 2: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
				case 3: coord1=0.; coord2=0.; coord3=0.; coord4=1.; break;
				case 4: coord1=1./4.; coord2=1./4.; coord3=1./4.; coord4=1./4.; break;
				default: _error_("node index should be in [0 4]");
			}
			break;
		case P2Enum:
			switch(iv){
				case 0: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
				case 1: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
				case 2: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
				case 3: coord1=0.; coord2=0.; coord3=0.; coord4=1.; break;

				case 4: coord1=0.; coord2=.5; coord3=.5; coord4=0.; break;
				case 5: coord1=.5; coord2=0.; coord3=.5; coord4=0.; break;
				case 6: coord1=.5; coord2=.5; coord3=0.; coord4=0.; break;
				case 7: coord1=.5; coord2=0.; coord3=0.; coord4=.5; break;
				case 8: coord1=0.; coord2=.5; coord3=0.; coord4=.5; break;
				case 9: coord1=0.; coord2=0.; coord3=.5; coord4=.5; break;
				default: _error_("node index should be in [0 9]");
			}
			break;
		default: _error_("Finite element "<<EnumToStringx(finiteelement)<<" not supported");
	}

}
/*}}}*/
void GaussTetra::GaussVertex(int iv){/*{{{*/

	/*in debugging mode: check that the default constructor has been called*/
	_assert_(numgauss==-1);

	/*update static arrays*/
	switch(iv){
		case 0: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
		case 1: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
		case 2: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
		case 3: coord1=0.; coord2=0.; coord3=0.; coord4=1.; break;
		default: _error_("vertex index should be in [0 3]");

	}

}
/*}}}*/
void GaussTetra::SynchronizeGaussBase(Gauss* gauss){/*{{{*/

	_error_("not supported");
}
/*}}}*/
