/*!\file PentaP1Input.c
 * \brief: implementation of the PentaP1Input 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 "../objects.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
#include "../../shared/shared.h"
#include "../../Container/Container.h"
#include "../../include/include.h"

/*PentaP1Input constructors and destructor*/
/*FUNCTION PentaP1Input::PentaP1Input(){{{1*/
PentaP1Input::PentaP1Input(){
	return;
}
/*}}}*/
/*FUNCTION PentaP1Input::PentaP1Input(int in_enum_type,double* values){{{1*/
PentaP1Input::PentaP1Input(int in_enum_type,double* in_values)
		:PentaRef(1)
{

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

	enum_type=in_enum_type;
	values[0]=in_values[0];
	values[1]=in_values[1];
	values[2]=in_values[2];
	values[3]=in_values[3];
	values[4]=in_values[4];
	values[5]=in_values[5];
}
/*}}}*/
/*FUNCTION PentaP1Input::~PentaP1Input(){{{1*/
PentaP1Input::~PentaP1Input(){
	return;
}
/*}}}*/

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

	printf("PentaP1Input:\n");
	printf("   enum: %i (%s)\n",this->enum_type,EnumToStringx(this->enum_type));
	printf("   values: [%g %g %g %g %g %g]\n",this->values[0],this->values[1],this->values[2],this->values[3],this->values[4],this->values[5]);
}
/*}}}*/
/*FUNCTION PentaP1Input::Id{{{1*/
int    PentaP1Input::Id(void){ return -1; }
/*}}}*/
/*FUNCTION PentaP1Input::MyRank{{{1*/
int    PentaP1Input::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
/*}}}*/
#ifdef _SERIAL_
/*FUNCTION PentaP1Input::Marshall{{{1*/
void  PentaP1Input::Marshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   enum_value=0;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*get enum value of PentaP1Input: */
	enum_value=PentaP1InputEnum;
	
	/*marshall enum: */
	memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value);
	
	/*marshall PentaP1Input data: */
	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	memcpy(marshalled_dataset,&values,sizeof(values));marshalled_dataset+=sizeof(values);

	*pmarshalled_dataset=marshalled_dataset;
}
/*}}}*/
/*FUNCTION PentaP1Input::MarshallSize{{{1*/
int   PentaP1Input::MarshallSize(){
	
	return sizeof(values)+
		+sizeof(enum_type)+
		+sizeof(int); //sizeof(int) for enum value
}
/*}}}*/
/*FUNCTION PentaP1Input::Demarshall{{{1*/
void  PentaP1Input::Demarshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   i;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
	 *object data (thanks to DataSet::Demarshall):*/
	memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	memcpy(&values,marshalled_dataset,sizeof(values));marshalled_dataset+=sizeof(values);

	/*return: */
	*pmarshalled_dataset=marshalled_dataset;
	return;
}
/*}}}*/
#endif
/*FUNCTION PentaP1Input::ObjectEnum{{{1*/
int PentaP1Input::ObjectEnum(void){

	return PentaP1InputEnum;

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

}
/*}}}*/
/*FUNCTION PentaP1Input::InstanceEnum{{{1*/
int PentaP1Input::InstanceEnum(void){

	return this->enum_type;

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

	/*output*/
	TriaP1Input* outinput=NULL;
	double newvalues[3];

	/*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 TriaP1Input(this->enum_type,&newvalues[0]);

	/*Assign output*/
	return outinput;

}
/*}}}*/
/*FUNCTION PentaP1Input::SpawnResult{{{1*/
ElementResult* PentaP1Input::SpawnResult(int step, double time){

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

}
/*}}}*/

/*Object functions*/
/*FUNCTION PentaP1Input::GetInputValue(double* pvalue,GaussPenta* gauss){{{1*/
void PentaP1Input::GetInputValue(double* pvalue,GaussPenta* gauss){

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

}
/*}}}*/
/*FUNCTION PentaP1Input::GetInputDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){{{1*/
void PentaP1Input::GetInputDerivativeValue(double* p, double* xyz_list, GaussPenta* gauss){

	/*Call PentaRef function*/
	PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss);
}
/*}}}*/
/*FUNCTION PentaP1Input::GetVxStrainRate3d{{{1*/
void PentaP1Input::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, GaussPenta* gauss){
	int i,j;

	const int numnodes=6;
	const int DOFVELOCITY=3;
	double B[8][27];
	double B_reduced[6][DOFVELOCITY*numnodes];
	double velocity[numnodes][DOFVELOCITY];

	/*Get B matrix: */
	GetBStokes(&B[0][0], xyz_list, gauss);
	/*Create a reduced matrix of B to get rid of pressure */
	for (i=0;i<6;i++){
		for (j=0;j<3;j++){
			B_reduced[i][j]=B[i][j];
		}
		for (j=4;j<7;j++){
			B_reduced[i][j-1]=B[i][j];
		}
		for (j=8;j<11;j++){
			B_reduced[i][j-2]=B[i][j];
		}
		for (j=12;j<15;j++){
			B_reduced[i][j-3]=B[i][j];
		}
		for (j=16;j<19;j++){
			B_reduced[i][j-4]=B[i][j];
		}
		for (j=20;j<23;j++){
			B_reduced[i][j-5]=B[i][j];
		}
	}

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

}
/*}}}*/
/*FUNCTION PentaP1Input::GetVyStrainRate3d{{{1*/
void PentaP1Input::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, GaussPenta* gauss){
	int i,j;

	const int numnodes=6;
	const int DOFVELOCITY=3;
	double B[8][27];
	double B_reduced[6][DOFVELOCITY*numnodes];
	double velocity[numnodes][DOFVELOCITY];

	/*Get B matrix: */
	GetBStokes(&B[0][0], xyz_list, gauss);
	/*Create a reduced matrix of B to get rid of pressure */
	for (i=0;i<6;i++){
		for (j=0;j<3;j++){
			B_reduced[i][j]=B[i][j];
		}
		for (j=4;j<7;j++){
			B_reduced[i][j-1]=B[i][j];
		}
		for (j=8;j<11;j++){
			B_reduced[i][j-2]=B[i][j];
		}
		for (j=12;j<15;j++){
			B_reduced[i][j-3]=B[i][j];
		}
		for (j=16;j<19;j++){
			B_reduced[i][j-4]=B[i][j];
		}
		for (j=20;j<23;j++){
			B_reduced[i][j-5]=B[i][j];
		}
	}

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

}
/*}}}*/
/*FUNCTION PentaP1Input::GetVzStrainRate3d{{{1*/
void PentaP1Input::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, GaussPenta* gauss){
	int i,j;

	const int numnodes=6;
	const int DOFVELOCITY=3;
	double B[8][27];
	double B_reduced[6][DOFVELOCITY*numnodes];
	double velocity[numnodes][DOFVELOCITY];

	/*Get B matrix: */
	GetBStokes(&B[0][0], xyz_list, gauss);
	/*Create a reduced matrix of B to get rid of pressure */
	for (i=0;i<6;i++){
		for (j=0;j<3;j++){
			B_reduced[i][j]=B[i][j];
		}
		for (j=4;j<7;j++){
			B_reduced[i][j-1]=B[i][j];
		}
		for (j=8;j<11;j++){
			B_reduced[i][j-2]=B[i][j];
		}
		for (j=12;j<15;j++){
			B_reduced[i][j-3]=B[i][j];
		}
		for (j=16;j<19;j++){
			B_reduced[i][j-4]=B[i][j];
		}
		for (j=20;j<23;j++){
			B_reduced[i][j-5]=B[i][j];
		}
	}

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

	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvz,0);

}
/*}}}*/
/*FUNCTION PentaP1Input::GetVxStrainRate3dPattyn{{{1*/
void PentaP1Input::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, GaussPenta* gauss){

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

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

	/*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 PentaP1Input::GetVyStrainRate3dPattyn{{{1*/
void PentaP1Input::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, GaussPenta* gauss){

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

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

	/*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 PentaP1Input::ChangeEnum{{{1*/
void PentaP1Input::ChangeEnum(int newenumtype){
	this->enum_type=newenumtype;
}
/*}}}*/
/*FUNCTION PentaP1Input::GetInputAverage{{{1*/
void PentaP1Input::GetInputAverage(double* pvalue){
	*pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
}
/*}}}*/

/*Intermediary*/
/*FUNCTION PentaP1Input::SquareMin{{{1*/
void PentaP1Input::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){

	int i;
	const int numnodes=6;
	double valuescopy[numnodes];
	double 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 PentaP1Input::ConstrainMin{{{1*/
void PentaP1Input::ConstrainMin(double minimum){
	
	int i;
	const int numnodes=6;

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

	/*Output*/
	const int numnodes=6;
	double norm=0;

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

	const int numnodes=6;
	double    max=values[0];

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

	const int numnodes=6;
	double    max=fabs(values[0]);

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

	const int numnodes=6;
	double    min=values[0];

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

	const int numnodes=6;
	double    min=fabs(values[0]);

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

	for(i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
}
/*}}}*/
/*FUNCTION PentaP1Input::AXPY{{{1*/
void PentaP1Input::AXPY(Input* xinput,double scalar){

	int i;
	const int numnodes=6;

	/*xinput is of the same type, so cast it: */

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

		case PentaP1InputEnum:{
			PentaP1Input* cast_input=(PentaP1Input*)xinput;
			for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*(cast_input->values[i]);}
			return;
		case ControlInputEnum:{
			ControlInput* cont_input=(ControlInput*)xinput;
			if(cont_input->values->ObjectEnum()!=PentaP1InputEnum) _error_("not supported yet");
			PentaP1Input* cast_input=(PentaP1Input*)cont_input->values;
			for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*(cast_input->values[i]);}
			return;
		default:
			_error_("not implemented yet");
	}

}
/*}}}*/
/*FUNCTION PentaP1Input::Constrain{{{1*/
void PentaP1Input::Constrain(double cm_min, double cm_max){

	int i;
	const int numnodes=6;
		
	if(!isnan(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
	if(!isnan(cm_max)) for(i=0;i<numnodes;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;

}
/*}}}*/
/*FUNCTION PentaP1Input::Extrude{{{1*/
void PentaP1Input::Extrude(void){

	int i;

	/*First 3 values copied on 3 last values*/
	for(i=0;i<3;i++) this->values[3+i]=this->values[i];
}
/*}}}*/
/*FUNCTION PentaP1Input::VerticallyIntegrate{{{1*/
void PentaP1Input::VerticallyIntegrate(Input* thickness_input){

	/*Intermediaries*/
	int i;
	const int  numnodes = 6;
	int        num_thickness_values;
	double    *thickness_values = NULL;

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

	/*Get Thickness value pointer*/
	thickness_input->GetValuesPtr(&thickness_values,&num_thickness_values);

	/*vertically integrate depending on type:*/
	switch(thickness_input->ObjectEnum()){

		case PentaP1InputEnum:
			for(i=0;i<3;i++){
				this->values[i]=0.5*(this->values[i]+this->values[i+3]) * thickness_values[i];
				this->values[i+3]=this->values[i];
			}
			return;

		default:
			_error_("not implemented yet");
	}
}
/*}}}*/
/*FUNCTION PentaP1Input::PointwiseDivide{{{1*/
Input* PentaP1Input::PointwiseDivide(Input* inputB){

	/*Ouput*/
	PentaP1Input* outinput=NULL;

	/*Intermediaries*/
	int               i;
	PentaP1Input *xinputB     = NULL;
	int               B_numvalues;
	const int         numnodes    = 6;
	double            AdotBvalues[numnodes];

	/*Check that inputB is of the same type*/
	if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type %s",EnumToStringx(inputB->ObjectEnum()));
	xinputB=(PentaP1Input*)inputB;

	/*Create point wise sum*/
	for(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 PentaP1Input(this->enum_type,&AdotBvalues[0]);

	/*Return output pointer*/
	return outinput;

}
/*}}}*/
/*FUNCTION PentaP1Input::PointwiseMin{{{1*/
Input* PentaP1Input::PointwiseMin(Input* inputB){

	/*Ouput*/
	PentaP1Input* outinput=NULL;

	/*Intermediaries*/
	int               i;
	PentaP1Input *xinputB     = NULL;
	int               B_numvalues;
	const int         numnodes    = 6;
	double            minvalues[numnodes];

	/*Check that inputB is of the same type*/
	if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type %s",EnumToStringx(inputB->ObjectEnum()));
	xinputB=(PentaP1Input*)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 Penta vertex input (copy of current input)*/
	outinput=new PentaP1Input(this->enum_type,&minvalues[0]);

	/*Return output pointer*/
	return outinput;

}
/*}}}*/
/*FUNCTION PentaP1Input::PointwiseMax{{{1*/
Input* PentaP1Input::PointwiseMax(Input* inputB){

	/*Ouput*/
	PentaP1Input* outinput=NULL;

	/*Intermediaries*/
	int               i;
	PentaP1Input *xinputB     = NULL;
	int               B_numvalues;
	const int         numnodes    = 6;
	double            maxvalues[numnodes];

	/*Check that inputB is of the same type*/
	if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type %s",EnumToStringx(inputB->ObjectEnum()));
	xinputB=(PentaP1Input*)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 Penta vertex input (copy of current input)*/
	outinput=new PentaP1Input(this->enum_type,&maxvalues[0]);

	/*Return output pointer*/
	return outinput;

}
/*}}}*/
/*FUNCTION PentaP1Input::GetVectorFromInputs{{{1*/
void PentaP1Input::GetVectorFromInputs(Vec vector,int* doflist){

	const int numvertices=6;
	VecSetValues(vector,numvertices,doflist,(const double*)this->values,INSERT_VALUES);

} /*}}}*/
/*FUNCTION PentaP1Input::GetValuesPtr{{{1*/
void PentaP1Input::GetValuesPtr(double** pvalues,int* pnum_values){

	*pvalues=this->values;
	*pnum_values=6;

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