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

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

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

	numgauss=-1;

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

	weight=UNDEF;
	coord1=UNDEF;
	coord2=UNDEF;
	coord3=UNDEF;
}
/*}}}*/
/*FUNCTION GaussTria::GaussTria(int order) {{{1*/
GaussTria::GaussTria(int order){

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

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

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

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

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

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

	/*Reverse index1 and 2 if necessary*/
	if (index1>index2){
		index3=index1; index1=index2; index2=index3;
		for(i=0;i<numgauss;i++) seg_coords[i]=-seg_coords[i];
	}

	/*Build Triangle Gauss point*/
	if (index1==0 && index2==1){
		for(i=0;i<numgauss;i++) coords1[i]=  0.5*(1-seg_coords[i]);
		for(i=0;i<numgauss;i++) coords2[i]=1-0.5*(1.-seg_coords[i]);
		for(i=0;i<numgauss;i++) coords3[i]=0;
		for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
	}
	else if (index1==0 && index2==2){
		for(i=0;i<numgauss;i++) coords1[i]=  0.5*(1-seg_coords[i]);
		for(i=0;i<numgauss;i++) coords2[i]=0;
		for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]);
		for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
	}
	else if (index1==1 && index2==2){
		for(i=0;i<numgauss;i++) coords1[i]=0;
		for(i=0;i<numgauss;i++) coords2[i]=  0.5*(1-seg_coords[i]);
		for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]);
		for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
	}
	else
	 ISSMERROR("The 2 indices provided are not supported yet (user provided %i and %i)",index1,index2);

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

	/*clean up*/
	xfree((void**)&seg_coords);
	xfree((void**)&seg_weights);
}
/*}}}*/
/*FUNCTION GaussTria::~GaussTria(){{{1*/
GaussTria::~GaussTria(){
	xfree((void**)&weights);
	xfree((void**)&coords1);
	xfree((void**)&coords2);
	xfree((void**)&coords3);
}
/*}}}*/

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

	printf("GaussTria:\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");

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

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

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

}
/*}}}*/
/*FUNCTION GaussTria::GaussEdgeCenter{{{1*/
void GaussTria::GaussEdgeCenter(int index1,int index2){

	int     index3;

	/*Reverse index1 and 2 if necessary*/
	if (index1>index2){
		index3=index1; index1=index2; index2=index3;
	}

	/*update static arrays*/
	if (index1==0 && index2==1){
		coord1=0.5;
		coord2=0.5;
		coord3=0.0;
	}
	else if (index1==0 && index2==2){
		coord1=0.5;
		coord2=0.0;
		coord3=0.5;
	}
	else if (index1==1 && index2==2){
		coord1=0.0;
		coord2=0.5;
		coord3=0.5;
	}
	else
	 ISSMERROR("The 2 indices provided are not supported yet (user provided %i and %i)",index1,index2);

}
/*}}}*/
/*FUNCTION GaussTria::GaussPoint{{{1*/
void GaussTria::GaussPoint(int ig){

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

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

}
/*}}}*/
/*FUNCTION GaussTria::GaussFromCoords{{{1*/
void GaussTria::GaussFromCoords(double x,double y,double* xyz_list){

	/*Intermediaries*/
	double    area = 0;
	double    x1,y1,x2,y2,x3,y3;

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

	x1=*(xyz_list+3*0+0); y1=*(xyz_list+3*0+1);
	x2=*(xyz_list+3*1+0); y2=*(xyz_list+3*1+1);
	x3=*(xyz_list+3*2+0); y3=*(xyz_list+3*2+1);

	area=(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;

	/*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
	coord1=((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area;

	/*Get second area coordinate = det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
	coord2=((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area;

	/*Get third  area coordinate 1-area1-area2: */
	coord3=1-coord1-coord2;

}
/*}}}*/
/*FUNCTION GaussTria::GaussVertex{{{1*/
void GaussTria::GaussVertex(int iv){

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

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

	}

}
/*}}}*/
/*FUNCTION GaussTria::begin{{{1*/
int GaussTria::begin(void){

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

	/*return first gauss index*/
	return 0;
}
/*}}}*/
/*FUNCTION GaussTria::end{{{1*/
int GaussTria::end(void){

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

	/*return last gauss index +1*/
	return numgauss;
}
/*}}}*/