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

/*Include files: {{{1*/
#include "./../objects.h"
/*}}}*/

/*GaussPenta constructors and destructors:*/
/*FUNCTION GaussPenta::GaussPenta() {{{1*/
GaussPenta::GaussPenta(){

	numgauss=-1;

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

	weight=UNDEF;
	coord1=UNDEF;
	coord2=UNDEF;
	coord3=UNDEF;
	coord4=UNDEF;
}
/*}}}*/
/*FUNCTION GaussPenta::GaussPenta(int order_horiz,int order_vert) {{{1*/
GaussPenta::GaussPenta(int order_horiz,int order_vert){

	/*Intermediaries*/
	int     ighoriz,igvert;
	int     numgauss_horiz;
	int     numgauss_vert;
	double *coords1_horiz = NULL;
	double *coords2_horiz = NULL;
	double *coords3_horiz = NULL;
	double *weights_horiz  = NULL;
	double *coords_vert = NULL;
	double *weights_vert   = NULL;

	/*Get gauss points*/
	GaussLegendreTria(&numgauss_horiz,&coords1_horiz,&coords2_horiz,&coords3_horiz,&weights_horiz,order_horiz);
	GaussLegendreLinear(&coords_vert,&weights_vert,order_vert);
	numgauss_vert=order_vert;

	/*Allocate GaussPenta fields*/
	numgauss=numgauss_horiz*numgauss_vert;
	coords1=(double*)xmalloc(numgauss*sizeof(double));
	coords2=(double*)xmalloc(numgauss*sizeof(double));
	coords3=(double*)xmalloc(numgauss*sizeof(double));
	coords4=(double*)xmalloc(numgauss*sizeof(double));
	weights=(double*)xmalloc(numgauss*sizeof(double));

	/*Combine Horizontal and vertical points*/
	for (ighoriz=0; ighoriz<numgauss_horiz; ighoriz++){
		for (igvert=0; igvert<numgauss_vert; igvert++){
			coords1[numgauss_vert*ighoriz+igvert]=coords1_horiz[ighoriz];
			coords2[numgauss_vert*ighoriz+igvert]=coords2_horiz[ighoriz];
			coords3[numgauss_vert*ighoriz+igvert]=coords3_horiz[ighoriz];
			coords4[numgauss_vert*ighoriz+igvert]=coords_vert[igvert];
			weights[numgauss_vert*ighoriz+igvert]=weights_horiz[ighoriz]*weights_vert[igvert];
		}
	}

	/*Initialize static fields as undefinite*/
	weight=UNDEF;
	coord1=UNDEF;
	coord2=UNDEF;
	coord3=UNDEF;
	coord4=UNDEF;

	/*Clean up*/
	xfree((void**)&coords1_horiz);
	xfree((void**)&coords2_horiz);
	xfree((void**)&coords3_horiz);
	xfree((void**)&coords_vert);
	xfree((void**)&weights_horiz);
	xfree((void**)&weights_vert);
}
/*}}}*/
/*FUNCTION GaussPenta::GaussPenta(int index1, int index2, int order){{{1*/
GaussPenta::GaussPenta(int index1, int index2,int order){

	/*Intermediaties*/
	double *seg_coords  = NULL;
	double *seg_weights = NULL;
	int     i;

	/*Get Segment gauss points*/
	numgauss=order;
	GaussLegendreLinear(&seg_coords,&seg_weights,numgauss);

	/*Allocate GaussPenta fields*/
	coords1=(double*)xmalloc(numgauss*sizeof(double));
	coords2=(double*)xmalloc(numgauss*sizeof(double));
	coords3=(double*)xmalloc(numgauss*sizeof(double));
	coords4=(double*)xmalloc(numgauss*sizeof(double));
	weights=(double*)xmalloc(numgauss*sizeof(double));

	if(index1==0 && index2==3){
		for(i=0;i<numgauss;i++) coords1[i]=1.0;
		for(i=0;i<numgauss;i++) coords2[i]=0.0;
		for(i=0;i<numgauss;i++) coords3[i]=0.0;
		for(i=0;i<numgauss;i++) coords4[i]=seg_coords[i];
		for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
	}
	else if (index1==1 && index2==4){
		for(i=0;i<numgauss;i++) coords1[i]=0.0;
		for(i=0;i<numgauss;i++) coords2[i]=1.0;
		for(i=0;i<numgauss;i++) coords3[i]=0.0;
		for(i=0;i<numgauss;i++) coords4[i]=seg_coords[i];
		for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
	}
	else if (index1==2 && index2==5){
		for(i=0;i<numgauss;i++) coords1[i]=0.0;
		for(i=0;i<numgauss;i++) coords2[i]=0.0;
		for(i=0;i<numgauss;i++) coords3[i]=1.0;
		for(i=0;i<numgauss;i++) coords4[i]=seg_coords[i];
		for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
	}
	else{
		_error_("Penta not supported yet");
	}

	/*Initialize static fields as undefined*/
	weight=UNDEF;
	coord1=UNDEF;
	coord2=UNDEF;
	coord3=UNDEF;
	coord4=UNDEF;

	/*clean up*/
	xfree((void**)&seg_coords);
	xfree((void**)&seg_weights);

}
/*}}}*/
/*FUNCTION GaussPenta::GaussPenta(int index1, int index2, int index3, int order){{{1*/
GaussPenta::GaussPenta(int index1, int index2, int index3, int order){

	/*Basal Tria*/
	if(index1==0 && index2==1 && index3==2){

		/*Get GaussTria*/
		GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);

		/*compute z coordinate*/
		coords4=(double*)xmalloc(numgauss*sizeof(double));
		for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
	}
	else{
		_error_("Tria not supported yet");
	}

}
/*}}}*/
/*FUNCTION GaussPenta::GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert){{{1*/
GaussPenta::GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert){

	/*Intermediaties*/
	double *seg_horiz_coords  = NULL;
	double *seg_horiz_weights = NULL;
	double *seg_vert_coords   = NULL;
	double *seg_vert_weights  = NULL;
	int     i,j;

	/*get the gauss points using the product of two line rules*/
	GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
	GaussLegendreLinear(&seg_vert_coords, &seg_vert_weights, order_vert);

	/*Allocate GaussPenta fields*/
	numgauss=order_horiz*order_vert;
	coords1=(double*)xmalloc(numgauss*sizeof(double));
	coords2=(double*)xmalloc(numgauss*sizeof(double));
	coords3=(double*)xmalloc(numgauss*sizeof(double));
	coords4=(double*)xmalloc(numgauss*sizeof(double));
	weights=(double*)xmalloc(numgauss*sizeof(double));

	/*Quads: get the gauss points using the product of two line rules  */
	if(index1==0 && index2==1 && index3==4 && index4==3){
		for(i=0;i<order_horiz;i++){
			for(j=0;j<order_vert;j++){
				coords1[i*order_vert+j]=  0.5*(1-seg_horiz_coords[i]);
				coords2[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
				coords3[i*order_vert+j]=0.0;
				coords4[i*order_vert+j]=seg_vert_coords[j];
				weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
			}
		}
	}
	else if(index1==1 && index2==2 && index3==5 && index4==4){
		for(i=0;i<order_horiz;i++){
			for(j=0;j<order_vert;j++){
				coords1[i*order_vert+j]=0.0;
				coords2[i*order_vert+j]=  0.5*(1-seg_horiz_coords[i]);
				coords3[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
				coords4[i*order_vert+j]=seg_vert_coords[j];
				weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
			}
		}
	}
	else if(index1==2 && index2==0 && index3==3 && index4==5){
		for(i=0;i<order_horiz;i++){
			for(j=0;j<order_vert;j++){
				coords1[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
				coords2[i*order_vert+j]=0.0;
				coords3[i*order_vert+j]=  0.5*(1-seg_horiz_coords[i]);
				coords4[i*order_vert+j]=seg_vert_coords[j];
				weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
			}
		}
	}
	else{
		_error_("Tria not supported yet (user provided indices %i %i %i %i)",index1,index2,index3,index4);
	}

	/*clean-up*/
	xfree((void**)&seg_horiz_coords);
	xfree((void**)&seg_horiz_weights);
	xfree((void**)&seg_vert_coords);
	xfree((void**)&seg_vert_weights);
}
/*}}}*/
/*FUNCTION GaussPenta::~GaussPenta(){{{1*/
GaussPenta::~GaussPenta(){
	xfree((void**)&weights);
	xfree((void**)&coords1);
	xfree((void**)&coords2);
	xfree((void**)&coords3);
	xfree((void**)&coords4);
}
/*}}}*/

/*Methods*/
/*FUNCTION GaussPenta::Echo{{{1*/
void GaussPenta::Echo(void){

	printf("GaussPenta:\n");
	printf("   numgauss: %i\n",numgauss);

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

	printf("   weight = %g\n",weight);
	printf("   coord1 = %g\n",coord1);
	printf("   coord2 = %g\n",coord2);
	printf("   coord3 = %g\n",coord3);
	printf("   coord4 = %g\n",coord4);

}
/*}}}*/
/*FUNCTION GaussPenta::GaussCenter{{{1*/
void GaussPenta::GaussCenter(void){

	/*update static arrays*/
	coord1=ONETHIRD;
	coord2=ONETHIRD;
	coord3=ONETHIRD;
	coord4=0.0;

}
/*}}}*/
/*FUNCTION GaussPenta::GaussPoint{{{1*/
void GaussPenta::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];

}
/*}}}*/
/*FUNCTION GaussPenta::GaussVertex{{{1*/
void GaussPenta::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= -1;
			break;
		case 1:
			coord1=0; coord2=1; coord3=0; coord4= -1;
			break;
		case 2:
			coord1=0; coord2=0; coord3=1; coord4= -1;
			break;
		case 3:
			coord1=1; coord2=0; coord3=0; coord4= +1;
			break;
		case 4:
			coord1=0; coord2=1; coord3=0; coord4= +1;
			break;
		case 5:
			coord1=0; coord2=0; coord3=1; coord4= +1;
			break;
		default:
			_error_("vertex index should be in [0 5]");

	}

}
/*}}}*/
/*FUNCTION GaussPenta::GaussFaceTria{{{1*/
void GaussPenta::GaussFaceTria(int index1, int index2, int index3, int order){

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

	/*Basal Tria*/
	if(index1==0 && index2==1 && index3==2){
		GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
		coords4=(double*)xmalloc(numgauss*sizeof(double));
		for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
	}
	else{
		_error_("Tria not supported yet");
	}

}
/*}}}*/
/*FUNCTION GaussPenta::begin{{{1*/
int GaussPenta::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;
}
/*}}}*/
/*FUNCTION GaussPenta::end{{{1*/
int GaussPenta::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;
}
/*}}}*/
/*FUNCTION GaussPenta::SynchronizeGaussTria{{{1*/
void GaussPenta::SynchronizeGaussTria(GaussTria* gauss_tria){

	gauss_tria->coord1=this->coord1;
	gauss_tria->coord2=this->coord2;
	gauss_tria->coord3=this->coord3;
	gauss_tria->weight=UNDEF;
}
/*}}}*/