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

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

#include <stdio.h>
#include <string.h>
#include "../../classes.h"
#include "../../../EnumDefinitions/EnumDefinitions.h"
#include "../../../shared/shared.h"
#include "../../../Container/Container.h"
#include "../../../include/include.h"

/*TriaP1Input constructors and destructor*/
/*FUNCTION TriaP1Input::TriaP1Input(){{{*/
TriaP1Input::TriaP1Input(){
	return;
}
/*}}}*/
/*FUNCTION TriaP1Input::TriaP1Input(int in_enum_type,IssmDouble* values){{{*/
TriaP1Input::TriaP1Input(int in_enum_type,IssmDouble* in_values)
	:TriaRef(1)
{

	/*Set TriaRef*/
	this->SetElementType(P1Enum,0);
	this->element_type=P1Enum;

	/*Set Enum*/
	enum_type=in_enum_type;

	/*Set values*/
	values[0]=in_values[0];
	values[1]=in_values[1];
	values[2]=in_values[2];
}
/*}}}*/
/*FUNCTION TriaP1Input::~TriaP1Input(){{{*/
TriaP1Input::~TriaP1Input(){
	return;
}
/*}}}*/

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

	_printLine_("TriaP1Input:");
	_printLine_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")");
	_printLine_("   values: [" << this->values[0] << " " << this->values[1] << " " << this->values[2] << "]");
}
/*}}}*/
/*FUNCTION TriaP1Input::Id{{{*/
int    TriaP1Input::Id(void){ return -1; }
/*}}}*/
/*FUNCTION TriaP1Input::MyRank{{{*/
int    TriaP1Input::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
/*}}}*/
/*FUNCTION TriaP1Input::ObjectEnum{{{*/
int TriaP1Input::ObjectEnum(void){

	return TriaP1InputEnum;

}
/*}}}*/
/*FUNCTION TriaP1Input::copy{{{*/
Object* TriaP1Input::copy() {
	
	return new TriaP1Input(this->enum_type,this->values);

}
/*}}}*/
	
/*TriaP1Input management*/
/*FUNCTION TriaP1Input::InstanceEnum{{{*/
int TriaP1Input::InstanceEnum(void){

	return this->enum_type;

}
/*}}}*/
/*FUNCTION TriaP1Input::SpawnTriaInput{{{*/
Input* TriaP1Input::SpawnTriaInput(int* indices){

	/*output*/
	TriaP1Input* outinput=NULL;

	/*Create new Tria input (copy of current input)*/
	outinput=new TriaP1Input(this->enum_type,&this->values[0]);

	/*Assign output*/
	return outinput;

}
/*}}}*/
/*FUNCTION TriaP1Input::SpawnResult{{{*/
ElementResult* TriaP1Input::SpawnResult(int step, IssmDouble time){

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

}
/*}}}*/

/*Object functions*/
/*FUNCTION TriaP1Input::GetInputValue(IssmDouble* pvalue,GaussTria* gauss){{{*/
void TriaP1Input::GetInputValue(IssmDouble* pvalue,GaussTria* gauss){

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

}
/*}}}*/
/*FUNCTION TriaP1Input::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussTria* gauss){{{*/
void TriaP1Input::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussTria* gauss){

	/*Call TriaRef function*/
	TriaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss);
}
/*}}}*/
/*FUNCTION TriaP1Input::GetVxStrainRate2d{{{*/
void TriaP1Input::GetVxStrainRate2d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussTria* gauss){

	/*Intermediary*/
	int       i;
	const int numnodes=3;
	IssmDouble B[3][NDOF2*numnodes];
	IssmDouble velocity[3][NDOF2];

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

	/*Here, we are computing the strain rate of (vx,0)*/
	for(i=0;i<3;i++){
		velocity[i][0]=this->values[i];
		velocity[i][1]=0.0;
	}
	/*Get epsilon(vx) = B*velocity*/
	MatrixMultiply( &B[0][0],3,NDOF2*numnodes,0,
				&velocity[0][0],NDOF2*numnodes,1,0,
				epsilonvx,0);
}
/*}}}*/
/*FUNCTION TriaP1Input::GetVyStrainRate2d{{{*/
void TriaP1Input::GetVyStrainRate2d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussTria* gauss){

	/*Intermediary*/
	int       i;
	const int numnodes=3;
	IssmDouble B[3][NDOF2*numnodes];
	IssmDouble velocity[3][NDOF2];

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

	/*Here, we are computing the strain rate of (0,vy)*/
	for(i=0;i<3;i++){
		velocity[i][0]=0.0;
		velocity[i][1]=this->values[i];
	}
	/*Get epsilon(vy) = B*velocity*/
	MatrixMultiply( &B[0][0],3,NDOF2*numnodes,0,
				&velocity[0][0],NDOF2*numnodes,1,0,
				epsilonvy,0);
}
/*}}}*/
/*FUNCTION TriaP1Input::ChangeEnum{{{*/
void TriaP1Input::ChangeEnum(int newenumtype){
	this->enum_type=newenumtype;
}
/*}}}*/
/*FUNCTION TriaP1Input::GetInputAverage{{{*/
void TriaP1Input::GetInputAverage(IssmDouble* pvalue){
	*pvalue=1./3.*(values[0]+values[1]+values[2]);
}
/*}}}*/

/*Intermediary*/
/*FUNCTION TriaP1Input::SquareMin{{{*/
void TriaP1Input::SquareMin(IssmDouble* psquaremin, bool process_units,Parameters* parameters){

	int i;
	const int numnodes=3;
	IssmDouble valuescopy[numnodes];
	IssmDouble squaremin;

	/*First,  copy values, to process units if requested: */
	for(i=0;i<numnodes;i++)valuescopy[i]=this->values[i];

	/*Process units if requested: */
	if(process_units)UnitConversion(&valuescopy[0],numnodes,IuToExtEnum,enum_type);

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

	for(i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
}
/*}}}*/
/*FUNCTION TriaP1Input::InfinityNorm{{{*/
IssmDouble TriaP1Input::InfinityNorm(void){

	/*Output*/
	IssmDouble norm=0;
	const int numnodes=3;

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

	const int numnodes=3;
	IssmDouble    max=values[0];

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

	const int numnodes=3;
	IssmDouble    max=fabs(values[0]);

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

	const int numnodes=3;
	IssmDouble    min=values[0];

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

	const int numnodes=3;
	IssmDouble    min=fabs(values[0]);

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

	for(i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
}
/*}}}*/
/*FUNCTION TriaP1Input::ArtificialNoise{{{*/
void TriaP1Input::ArtificialNoise(IssmDouble min,IssmDouble max){

	int i;
	const int numnodes=3;
	IssmDouble noise;

	/*Compute random number between bounds:
	 * rand() outputs an integer in [0 RAND_MAX]
	 * (IssmDouble)rand()/RAND_MAX is in [0 1]
	 */
	 noise=min+(max-min)*(IssmDouble)rand()/RAND_MAX;

	for(i=0;i<numnodes;i++)values[i]=values[i]+noise;
}
/*}}}*/
/*FUNCTION TriaP1Input::AXPY{{{*/
void TriaP1Input::AXPY(Input* xinput,IssmDouble scalar){

	int i;
	const int numnodes=3;
	TriaP1Input*  xtriavertexinput=NULL;

	/*xinput is of the same type, so cast it: */
	xtriavertexinput=(TriaP1Input*)xinput;

	/*Carry out the AXPY operation depending on type:*/
	switch(xinput->ObjectEnum()){

		case TriaP1InputEnum :
			for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*xtriavertexinput->values[i];
			return;

		default :
			_error2_("not implemented yet");
	}

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

	int i;
	const int numnodes=3;
		
	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 TriaP1Input::GetVectorFromInputs{{{*/
void TriaP1Input::GetVectorFromInputs(Vector* vector,int* doflist){

	const int numvertices=3;
	vector->SetValues(numvertices,doflist,this->values,INS_VAL);

} /*}}}*/
/*FUNCTION TriaP1Input::GetValuesPtr{{{*/
void TriaP1Input::GetValuesPtr(IssmDouble** pvalues,int* pnum_values){

	*pvalues=this->values;
	if(pnum_values)*pnum_values=3;

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

	/*Ouput*/
	TriaP1Input* outinput=NULL;

	/*Intermediaries*/
	int               i;
	TriaP1Input *xinputB     = NULL;
	int               B_numvalues;
	const int         numnodes    = 3;
	IssmDouble            minvalues[numnodes];

	/*Check that inputB is of the same type*/
	if (inputB->ObjectEnum()!=TriaP1InputEnum) _error2_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
	xinputB=(TriaP1Input*)inputB;

	/*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 Tria vertex input (copy of current input)*/
	outinput=new TriaP1Input(this->enum_type,&minvalues[0]);

	/*Return output pointer*/
	return outinput;

}
/*}}}*/
/*FUNCTION TriaP1Input::PointwiseMax{{{*/
Input* TriaP1Input::PointwiseMax(Input* inputB){

	/*Ouput*/
	TriaP1Input* outinput=NULL;

	/*Intermediaries*/
	int               i;
	TriaP1Input *xinputB     = NULL;
	int               B_numvalues;
	const int         numnodes    = 3;
	IssmDouble            maxvalues[numnodes];

	/*Check that inputB is of the same type*/
	if (inputB->ObjectEnum()!=TriaP1InputEnum) _error2_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
	xinputB=(TriaP1Input*)inputB;

	/*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 Tria vertex input (copy of current input)*/
	outinput=new TriaP1Input(this->enum_type,&maxvalues[0]);

	/*Return output pointer*/
	return outinput;

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