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

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "../classes.h"
#include "../../shared/shared.h"

/*PentaInput constructors and destructor*/
/*FUNCTION PentaInput::PentaInput(){{{*/
PentaInput::PentaInput(){
	values = NULL;
}
/*}}}*/
/*FUNCTION PentaInput::PentaInput(int in_enum_type,IssmDouble* values,int element_type){{{*/
PentaInput::PentaInput(int in_enum_type,IssmDouble* in_values,int element_type_in)
		:PentaRef(1)
{

	/*Set PentaRef*/
	this->SetElementType(element_type_in,0);
	this->element_type=element_type_in;

	/*Set Enum*/
	enum_type=in_enum_type;

	/*Set values*/
	this->values=xNew<IssmDouble>(this->NumberofNodes());
	for(int i=0;i<this->NumberofNodes();i++) values[i]=in_values[i];
}
/*}}}*/
/*FUNCTION PentaInput::~PentaInput(){{{*/
PentaInput::~PentaInput(){
	xDelete<IssmDouble>(this->values);
}
/*}}}*/

/*Object virtual functions definitions:*/
/*FUNCTION PentaInput::Echo {{{*/
void PentaInput::Echo(void){
	this->DeepEcho();
}
/*}}}*/
/*FUNCTION PentaInput::DeepEcho{{{*/
void PentaInput::DeepEcho(void){

	_printf_("PentaInput:\n");
	_printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
	_printf_("   values: [");
	for(int i=0;i<this->NumberofNodes();i++) _printf_(" "<<this->values[i]);
	_printf_("]\n");
}
/*}}}*/
/*FUNCTION PentaInput::Id{{{*/
int    PentaInput::Id(void){ return -1; }
/*}}}*/
/*FUNCTION PentaInput::ObjectEnum{{{*/
int PentaInput::ObjectEnum(void){

	return PentaInputEnum;

}
/*}}}*/
/*FUNCTION PentaInput::copy{{{*/
Object* PentaInput::copy() {

	return new PentaInput(this->enum_type,this->values,this->element_type);

}
/*}}}*/

/*PentaInput management*/
/*FUNCTION PentaInput::InstanceEnum{{{*/
int PentaInput::InstanceEnum(void){

	return this->enum_type;

}
/*}}}*/
/*FUNCTION PentaInput::SpawnTriaInput{{{*/
Input* PentaInput::SpawnTriaInput(int location){

	/*output*/
	TriaInput* outinput=NULL;
	IssmDouble newvalues[3]; //Assume P1 interpolation only for now

	/*Create arrow of indices depending on location (0=base 1=surface)*/
	int indices[3];
	switch(location){
		case 0:
			indices[0] = 0;
			indices[1] = 1;
			indices[2] = 2;
			break;
		case 1:
			indices[0] = 3;
			indices[1] = 4;
			indices[2] = 5;
			break;
		default:
			_error_("case "<<location<<" not supported");
	}

	/*Loop over the new indices*/
	for(int i=0;i<3;i++){

		/*Check index value*/
		_assert_(indices[i]>=0 && indices[i]<6);

		/*Assign value to new input*/
		newvalues[i]=this->values[indices[i]];
	}

	/*Create new Tria input*/
	outinput=new TriaInput(this->enum_type,&newvalues[0],P1Enum);

	/*Assign output*/
	return outinput;
}
/*}}}*/
/*FUNCTION PentaInput::SpawnResult{{{*/
ElementResult* PentaInput::SpawnResult(int step, IssmDouble time){

	return new PentaP1ElementResult(this->enum_type,this->values,step,time);

}
/*}}}*/

/*Object functions*/
/*FUNCTION PentaInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){{{*/
void PentaInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){

	/*Call PentaRef function*/
	PentaRef::GetInputValue(pvalue,&values[0],gauss);

}
/*}}}*/
/*FUNCTION PentaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){{{*/
void PentaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){

	/*Call PentaRef function*/
	PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss);
}
/*}}}*/
/*FUNCTION PentaInput::GetVxStrainRate3d{{{*/
void PentaInput::GetVxStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){

	/*Intermediary*/
	int         numnodes=this->NumberofNodes();
	IssmDouble* B=xNew<IssmDouble>(numnodes*(NDOF3*numnodes));
	IssmDouble* velocity=xNew<IssmDouble>(numnodes*NDOF3);

	/*Get B matrix: */
	GetBFSstrainrate(B,xyz_list,gauss);

	/*Here, we are computing the strain rate of (vx,0,0)*/
	for(int i=0;i<numnodes;i++){
		velocity[NDOF3*i+0]=this->values[i];
		velocity[NDOF3*i+1]=0.;
		velocity[NDOF3*i+2]=0.;
	}
	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply(B,numnodes,NDOF3*numnodes,0,velocity,NDOF3*numnodes,1,0,epsilonvx,0);

	/*Clean-up*/
	xDelete<IssmDouble>(B);
	xDelete<IssmDouble>(velocity);

}
/*}}}*/
/*FUNCTION PentaInput::GetVyStrainRate3d{{{*/
void PentaInput::GetVyStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){

	/*Intermediary*/
	int         numnodes=this->NumberofNodes();
	IssmDouble* B=xNew<IssmDouble>(numnodes*(NDOF3*numnodes));
	IssmDouble* velocity=xNew<IssmDouble>(numnodes*NDOF3);

	/*Get B matrix: */
	GetBFSstrainrate(B,xyz_list,gauss);

	/*Here, we are computing the strain rate of (0,vy,0)*/
	for(int i=0;i<numnodes;i++){
		velocity[NDOF3*i+0]=0.;
		velocity[NDOF3*i+1]=this->values[i];
		velocity[NDOF3*i+2]=0.;
	}
	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply(B,numnodes,NDOF3*numnodes,0,velocity,NDOF3*numnodes,1,0,epsilonvx,0);

	/*Clean-up*/
	xDelete<IssmDouble>(B);
	xDelete<IssmDouble>(velocity);
}
/*}}}*/
/*FUNCTION PentaInput::GetVzStrainRate3d{{{*/
void PentaInput::GetVzStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){

	/*Intermediary*/
	int         numnodes=this->NumberofNodes();
	IssmDouble* B=xNew<IssmDouble>(numnodes*(NDOF3*numnodes));
	IssmDouble* velocity=xNew<IssmDouble>(numnodes*NDOF3);

	/*Get B matrix: */
	GetBFSstrainrate(B,xyz_list,gauss);

	/*Here, we are computing the strain rate of (0,0,vz)*/
	for(int i=0;i<numnodes;i++){
		velocity[NDOF3*i+0]=0.;
		velocity[NDOF3*i+1]=0.;
		velocity[NDOF3*i+2]=this->values[i];
	}
	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply(B,numnodes,NDOF3*numnodes,0,velocity,NDOF3*numnodes,1,0,epsilonvx,0);

	/*Clean-up*/
	xDelete<IssmDouble>(B);
	xDelete<IssmDouble>(velocity);
}
/*}}}*/
/*FUNCTION PentaInput::GetVxStrainRate3dHO{{{*/
void PentaInput::GetVxStrainRate3dHO(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){

	int i;
	const int numnodes=6;
	IssmDouble B[5][NDOF2*numnodes];
	IssmDouble velocity[numnodes][NDOF2];

	/*Get B matrix: */
	GetBHO(&B[0][0], xyz_list, gauss);

	_assert_(this->NumberofNodes()==6); //Check Tria too

	/*Here, we are computing the strain rate of (vx,0)*/
	for(i=0;i<numnodes;i++){
		velocity[i][0]=this->values[i];
		velocity[i][1]=0.0;
	}

	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply( &B[0][0],5,NDOF2*numnodes,0,
				&velocity[0][0],NDOF2*numnodes,1,0,
				epsilonvx,0);

}
/*}}}*/
/*FUNCTION PentaInput::GetVyStrainRate3dHO{{{*/
void PentaInput::GetVyStrainRate3dHO(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){

	int i;
	const int numnodes=6;
	IssmDouble B[5][NDOF2*numnodes];
	IssmDouble velocity[numnodes][NDOF2];

	/*Get B matrix: */
	GetBHO(&B[0][0], xyz_list, gauss);
	_assert_(this->NumberofNodes()==6); //Check Tria too

	/*Here, we are computing the strain rate of (0,vy)*/
	for(i=0;i<numnodes;i++){
		velocity[i][0]=0.0;
		velocity[i][1]=this->values[i];
	}

	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply( &B[0][0],5,NDOF2*numnodes,0,
				&velocity[0][0],NDOF2*numnodes,1,0,
				epsilonvy,0);

}
/*}}}*/
/*FUNCTION PentaInput::ChangeEnum{{{*/
void PentaInput::ChangeEnum(int newenumtype){
	this->enum_type=newenumtype;
}
/*}}}*/
/*FUNCTION PentaInput::GetInputAverage{{{*/
void PentaInput::GetInputAverage(IssmDouble* pvalue){

	int        numnodes  = this->NumberofNodes();
	IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
	IssmDouble value     = 0.;

	for(int i=0;i<numnodes;i++) value+=values[i];
	value = value/numnodesd;

	*pvalue=value;
}
/*}}}*/

/*Intermediary*/
/*FUNCTION PentaInput::SquareMin{{{*/
void PentaInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){

	int        numnodes=this->NumberofNodes();
	IssmDouble squaremin;

	/*Now, figure out minimum of valuescopy: */
	squaremin=pow(this->values[0],2);
	for(int i=1;i<numnodes;i++){
		if(pow(this->values[i],2)<squaremin)squaremin=pow(this->values[i],2);
	}
	/*Assign output pointers:*/
	*psquaremin=squaremin;
}
/*}}}*/
/*FUNCTION PentaInput::ConstrainMin{{{*/
void PentaInput::ConstrainMin(IssmDouble minimum){

	int numnodes = this->NumberofNodes();
	for(int i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
}
/*}}}*/
/*FUNCTION PentaInput::InfinityNorm{{{*/
IssmDouble PentaInput::InfinityNorm(void){

	/*Output*/
	IssmDouble norm=0.;
	int numnodes=this->NumberofNodes();

	for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
	return norm;
}
/*}}}*/
/*FUNCTION PentaInput::Max{{{*/
IssmDouble PentaInput::Max(void){

	int  numnodes=this->NumberofNodes();
	IssmDouble max=values[0];

	for(int i=1;i<numnodes;i++){
		if(values[i]>max) max=values[i];
	}
	return max;
}
/*}}}*/
/*FUNCTION PentaInput::MaxAbs{{{*/
IssmDouble PentaInput::MaxAbs(void){

	int  numnodes=this->NumberofNodes();
	IssmDouble max=fabs(values[0]);

	for(int i=1;i<numnodes;i++){
		if(fabs(values[i])>max) max=fabs(values[i]);
	}
	return max;
}
/*}}}*/
/*FUNCTION PentaInput::Min{{{*/
IssmDouble PentaInput::Min(void){

	const int  numnodes=this->NumberofNodes();
	IssmDouble min=values[0];

	for(int i=1;i<numnodes;i++){
		if(values[i]<min) min=values[i];
	}
	return min;
}
/*}}}*/
/*FUNCTION PentaInput::MinAbs{{{*/
IssmDouble PentaInput::MinAbs(void){

	const int  numnodes=this->NumberofNodes();
	IssmDouble min=fabs(values[0]);

	for(int i=1;i<numnodes;i++){
		if(fabs(values[i])<min) min=fabs(values[i]);
	}
	return min;
}
/*}}}*/
/*FUNCTION PentaInput::Scale{{{*/
void PentaInput::Scale(IssmDouble scale_factor){

	const int numnodes=this->NumberofNodes();
	for(int i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
}
/*}}}*/
/*FUNCTION PentaInput::AXPY{{{*/
void PentaInput::AXPY(Input* xinput,IssmDouble scalar){

	const int numnodes=this->NumberofNodes();
	PentaInput* xpentainput=NULL;

	/*If xinput is a ControlInput, take its values directly*/
	if(xinput->ObjectEnum()==ControlInputEnum){
		xinput=((ControlInput*)xinput)->values;
	}

	/*xinput is of the same type, so cast it: */
	if(xinput->ObjectEnum()!=PentaInputEnum)
	  _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
	xpentainput=(PentaInput*)xinput;
	if(xpentainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xpentainput->element_type));

	/*Carry out the AXPY operation depending on type:*/
	for(int i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*xpentainput->values[i];

}
/*}}}*/
/*FUNCTION PentaInput::Constrain{{{*/
void PentaInput::Constrain(IssmDouble cm_min, IssmDouble cm_max){

	int i;
	const int numnodes=this->NumberofNodes();

	if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
	if(!xIsNan<IssmDouble>(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;

}
/*}}}*/
/*FUNCTION PentaInput::Extrude{{{*/
void PentaInput::Extrude(void){

	switch(this->element_type){
		case P1Enum:
			for(int i=0;i<3;i++) this->values[3+i]=this->values[i];
			break;
		default:
			_error_("not supported yet for type "<<EnumToStringx(this->element_type));
	}
}
/*}}}*/
/*FUNCTION PentaInput::VerticallyIntegrate{{{*/
void PentaInput::VerticallyIntegrate(Input* thickness_input){

	IssmDouble thickness;

	/*Check that input provided is a thickness*/
	if (thickness_input->InstanceEnum()!=ThicknessEnum) _error_("Input provided is not a Thickness (enum_type is " << EnumToStringx(thickness_input->InstanceEnum()) << ")");

	/*vertically integrate depending on type:*/
	switch(this->element_type){
		case P1Enum:{
			GaussPenta *gauss=new GaussPenta();
			for(int iv=0;iv<3;iv++){
				gauss->GaussVertex(iv);
				thickness_input->GetInputValue(&thickness,gauss);
				this->values[iv]=0.5*(this->values[iv]+this->values[iv+3]) * thickness;
				this->values[iv+3]=this->values[iv];
			}
			delete gauss;
			return; }
		default:
			_error_("not supported yet for type "<<EnumToStringx(this->element_type));
	}
}
/*}}}*/
/*FUNCTION PentaInput::PointwiseDivide{{{*/
Input* PentaInput::PointwiseDivide(Input* inputB){

	/*Ouput*/
	PentaInput* outinput=NULL;

	/*Intermediaries*/
	PentaInput *xinputB  = NULL;
	const int   numnodes = this->NumberofNodes();

	/*Check that inputB is of the same type*/
	if(inputB->ObjectEnum()!=PentaInputEnum)     _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
	xinputB=(PentaInput*)inputB;
	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));

	/*Allocate intermediary*/
	IssmDouble* AdotBvalues=xNew<IssmDouble>(numnodes);

	/*Create point wise sum*/
	for(int i=0;i<numnodes;i++){
		_assert_(xinputB->values[i]!=0);
		AdotBvalues[i]=this->values[i]/xinputB->values[i];
	}

	/*Create new Penta vertex input (copy of current input)*/
	outinput=new PentaInput(this->enum_type,AdotBvalues,this->element_type);

	/*Return output pointer*/
	xDelete<IssmDouble>(AdotBvalues);
	return outinput;

}
/*}}}*/
/*FUNCTION PentaInput::PointwiseMin{{{*/
Input* PentaInput::PointwiseMin(Input* inputB){

	/*Ouput*/
	PentaInput* outinput=NULL;

	/*Intermediaries*/
	int         i;
	PentaInput  *xinputB   = NULL;
	const int   numnodes  = this->NumberofNodes();
	IssmDouble *minvalues = xNew<IssmDouble>(numnodes);

	/*Check that inputB is of the same type*/
	if(inputB->ObjectEnum()!=PentaInputEnum)       _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
	xinputB=(PentaInput*)inputB;
	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));

	/*Create point wise min*/
	for(i=0;i<numnodes;i++){
		if(this->values[i] > xinputB->values[i]) minvalues[i]=xinputB->values[i];
		else minvalues[i]=this->values[i];
	}

	/*Create new Penta vertex input (copy of current input)*/
	outinput=new PentaInput(this->enum_type,&minvalues[0],this->element_type);

	/*Return output pointer*/
	xDelete<IssmDouble>(minvalues);
	return outinput;
}
/*}}}*/
/*FUNCTION PentaInput::PointwiseMax{{{*/
Input* PentaInput::PointwiseMax(Input* inputB){

	/*Ouput*/
	PentaInput* outinput=NULL;

	/*Intermediaries*/
	int         i;
	PentaInput  *xinputB   = NULL;
	const int   numnodes  = this->NumberofNodes();
	IssmDouble *maxvalues = xNew<IssmDouble>(numnodes);

	/*Check that inputB is of the same type*/
	if(inputB->ObjectEnum()!=PentaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
	xinputB=(PentaInput*)inputB;
	if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));

	/*Create point wise max*/
	for(i=0;i<numnodes;i++){
		if(this->values[i] < xinputB->values[i]) maxvalues[i]=xinputB->values[i];
		else maxvalues[i]=this->values[i];
	}

	/*Create new Penta vertex input (copy of current input)*/
	outinput=new PentaInput(this->enum_type,&maxvalues[0],this->element_type);

	/*Return output pointer*/
	xDelete<IssmDouble>(maxvalues);
	return outinput;
}
/*}}}*/
/*FUNCTION PentaInput::GetVectorFromInputs{{{*/
void PentaInput::GetVectorFromInputs(Vector<IssmDouble>* vector,int* doflist){

	const int numnodes=this->NumberofNodes();
	vector->SetValues(numnodes,doflist,this->values,INS_VAL);

} /*}}}*/
/*FUNCTION PentaInput::Configure{{{*/
void PentaInput::Configure(Parameters* parameters){
	/*do nothing: */
}
/*}}}*/
