/*!\file TransientInput.c
 * \brief: implementation of the TransientInput object
 */
/*Headers{{{1*/
#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"
/*}}}*/

/*TransientInput constructors and destructor*/
/*FUNCTION TransientInput::TransientInput(){{{1*/
TransientInput::TransientInput(){

	enum_type=UNDEF;
	inputs=NULL;
	this->numtimesteps=0;
	this->parameters=NULL;
	this->timesteps=NULL;

}
/*}}}*/
/*FUNCTION TransientInput::TransientInput(int in_enum_type){{{1*/
TransientInput::TransientInput(int in_enum_type)
{
	/*Set Enum*/
	enum_type=in_enum_type;

	/*Allocate values and timesteps, and copy: */
	this->numtimesteps=0;
	this->timesteps=NULL;
	inputs = new Inputs();
	this->parameters=NULL;

}
/*}}}*/
/*FUNCTION TransientInput::~TransientInput{{{1*/
TransientInput::~TransientInput(){
	xfree((void**)&this->timesteps);
	this->timesteps=NULL;
	this->numtimesteps=0;
	parameters=NULL;
	delete this->inputs;
	return;
}
/*}}}*/
/*FUNCTION void TransientInput::AddTimeInput(Input* input,double time){{{1*/
void TransientInput::AddTimeInput(Input* input,double time){

	/*insert values at time step: */
	if (this->numtimesteps>0 && time<=this->timesteps[this->numtimesteps-1]) _assert_("timestep values must increase sequentially");

	//copy timesteps, add the new time, delete previous timesteps, and add the new input: inputs->AddObject(input);
   double* old_timesteps=NULL;

	if (this->numtimesteps > 0){
		old_timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
		memcpy(old_timesteps,this->timesteps,this->numtimesteps*sizeof(double));
		xfree((void**)&this->timesteps); 
	}

	this->numtimesteps=this->numtimesteps+1;
   this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));

	if (this->numtimesteps > 1){
		memcpy(this->timesteps,old_timesteps,(this->numtimesteps-1)*sizeof(double));
		xfree((void**)&old_timesteps);
	}

	/*go ahead and plug: */
	this->timesteps[this->numtimesteps-1]=time;
	inputs->AddObject(input);

}
/*}}}*/

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

	int i;

	printf("TransientInput:\n");
	printf("   enum: %i (%s)\n",this->enum_type,EnumToStringx(this->enum_type));
	printf("   numtimesteps: %i\n",this->numtimesteps);
	printf("---inputs: \n"); 
	for(i=0;i<this->numtimesteps;i++){
		printf("   time: %g  \n",this->timesteps[i]);
		((Input*)this->inputs->GetObjectByOffset(i))->Echo();
	}
}
/*}}}*/
/*FUNCTION TransientInput::Id{{{1*/
int    TransientInput::Id(void){ return -1; }
/*}}}*/
/*FUNCTION TransientInput::MyRank{{{1*/
int    TransientInput::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
/*}}}*/
#ifdef _SERIAL_
/*FUNCTION TransientInput::Marshall{{{1*/
void  TransientInput::Marshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	char* marshalled_inputs=NULL;
	int   marshalled_inputs_size;
	int   enum_value=0;

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

	/*get enum value of TransientInput: */
	enum_value=TransientInputEnum;
	
	/*marshall enum: */
	memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value);

	/*marshall TransientInput data: */
	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	memcpy(marshalled_dataset,&numtimesteps,sizeof(numtimesteps));marshalled_dataset+=sizeof(numtimesteps);
	memcpy(marshalled_dataset,timesteps,numtimesteps*sizeof(double));marshalled_dataset+=numtimesteps*sizeof(double);

	/*marshal inputs*/
	marshalled_inputs_size=inputs->MarshallSize();
	marshalled_inputs=inputs->Marshall();
	memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
	marshalled_dataset+=marshalled_inputs_size;

	/*clean up and assign output pointer*/
	xfree((void**)&marshalled_inputs);
	*pmarshalled_dataset=marshalled_dataset;

}
/*}}}*
/*FUNCTION TransientInput::MarshallSize{{{1*/
int   TransientInput::MarshallSize(){

	return 
		+sizeof(enum_type)+
		+sizeof(numtimesteps)+
		+inputs->MarshallSize()
		+numtimesteps*sizeof(double)+
		+sizeof(int); //sizeof(int) for enum value
}
/*}}}*/
/*FUNCTION TransientInput::Demarshall{{{1*/
void  TransientInput::Demarshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;

	/*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(&numtimesteps,marshalled_dataset,sizeof(numtimesteps));marshalled_dataset+=sizeof(numtimesteps);

	/*allocate: */
	timesteps=(double*)xmalloc(numtimesteps*sizeof(double));
	memcpy(timesteps,marshalled_dataset,numtimesteps*sizeof(double));marshalled_dataset+=numtimesteps*sizeof(double);

	/*Demarshal values*/
	inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset);

	/*return: */
	*pmarshalled_dataset=marshalled_dataset;
	return;

}
/*}}}*/
#endif
/*FUNCTION TransientInput::ObjectEnum{{{1*/
int TransientInput::ObjectEnum(void){

	return TransientInputEnum;

}
/*}}}*/
/*FUNCTION TransientInput::copy{{{1*/
Object* TransientInput::copy() {

	TransientInput* output=NULL;

	output = new TransientInput();
	output->enum_type=this->enum_type;
	output->numtimesteps=this->numtimesteps;
	output->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
   memcpy(output->timesteps,this->timesteps,this->numtimesteps*sizeof(double));
	output->inputs=(Inputs*)this->inputs->Copy();
	output->parameters=this->parameters;

	return output;

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

	return this->enum_type;

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

	/*output*/
	TransientInput* outinput=NULL;

	/*Create new Transientinput (copy of current input)*/
	outinput=new TransientInput();
	outinput->enum_type=this->enum_type;
	outinput->numtimesteps=this->numtimesteps;
	outinput->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
	memcpy(outinput->timesteps,this->timesteps,this->numtimesteps*sizeof(double));
	outinput->inputs=(Inputs*)this->inputs->SpawnTriaInputs(indices);
	outinput->parameters=this->parameters;

	/*Assign output*/
	return outinput;

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

	ElementResult* elementresult=NULL;

	/*Ok, we want to spawn an ElementResult. We have the time, just get 
	 *the correct values: */
	Input* input=GetTimeInput(time);

	elementresult=input->SpawnResult(step,time);

   delete input;

	return elementresult;
}
/*}}}*/

/*Object functions*/
/*FUNCTION TransientInput::GetInputValue(double* pvalue,GaussTria* gauss){{{1*/
void TransientInput::GetInputValue(double* pvalue,GaussTria* gauss){
	
	double time;

	/*First, recover current time from parameters: */
	this->parameters->FindParam(&time,TimeEnum);

	/*Retrieve interpolated values for this time step: */
	Input* input=GetTimeInput(time);

	/*Call input function*/
	input->GetInputValue(pvalue,gauss);

	delete input;

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

	double time;

	/*First, recover current time from parameters: */
	parameters->FindParam(&time,TimeEnum);

	/*Retrieve interpolated values for this time step: */
	Input* input=GetTimeInput(time);
		   
	/*Call input function*/
	input->GetInputDerivativeValue(p,xyz_list,gauss);

	delete input;

}
/*}}}*/
/*FUNCTION TransientInput::ChangeEnum{{{1*/
void TransientInput::ChangeEnum(int newenumtype){
	this->enum_type=newenumtype;
}
/*}}}*/
/*FUNCTION TransientInput::GetInputAverage{{{1*/
void TransientInput::GetInputAverage(double* pvalue){
	
	double time;

	/*First, recover current time from parameters: */
	parameters->FindParam(&time,TimeEnum);

	/*Retrieve interpolated values for this time step: */
	Input* input=GetTimeInput(time);

	/*Call input function*/
	input->GetInputAverage(pvalue);
			   
	delete input;

}
/*}}}*/

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

	double time;

	/*First, recover current time from parameters: */
	parameters->FindParam(&time,TimeEnum);

   /*Retrieve interpolated values for this time step: */
	Input* input=GetTimeInput(time);
		   
	/*Call input function*/
	input->SquareMin(psquaremin,process_units,parameters);
			   
	delete input;

}
/*}}}*/
/*FUNCTION TransientInput::InfinityNorm{{{1*/
double TransientInput::InfinityNorm(void){

	double time;
	double infnorm;

	/*First, recover current time from parameters: */
	parameters->FindParam(&time,TimeEnum);

   /*Retrieve interpolated values for this time step: */
	Input* input=GetTimeInput(time);

	/*Call input function*/
	infnorm=input->InfinityNorm();
			   
	/*Clean-up and return*/
	delete input;
	return infnorm;
}
/*}}}*/
/*FUNCTION TransientInput::Max{{{1*/
double TransientInput::Max(void){

	double time;
	double max;

	/*First, recover current time from parameters: */
	parameters->FindParam(&time,TimeEnum);

   /*Retrieve interpolated values for this time step: */
	Input* input=GetTimeInput(time);
		   
	/*Call input function*/
	max=input->Max();
			   
	delete input;

	return max;
}
/*}}}*/
/*FUNCTION TransientInput::MaxAbs{{{1*/
double TransientInput::MaxAbs(void){

	double time;
	double maxabs;

	/*First, recover current time from parameters: */
	parameters->FindParam(&time,TimeEnum);

	/*Retrieve interpolated values for this time step: */
	Input* input=GetTimeInput(time);

	/*Call input function*/
	maxabs=input->MaxAbs();

	/*Clean-up and return*/
	delete input;
	return maxabs;

}
/*}}}*/
/*FUNCTION TransientInput::Min{{{1*/
double TransientInput::Min(void){

	double time;
	double min;

	/*First, recover current time from parameters: */
	parameters->FindParam(&time,TimeEnum);

   /*Retrieve interpolated values for this time step: */
	Input* input=GetTimeInput(time);

	/*Call input function*/
	min=input->Min();

	/*Clean-up and return*/
	delete input;
	return min;

}
/*}}}*/
/*FUNCTION TransientInput::MinAbs{{{1*/
double TransientInput::MinAbs(void){

	double time;
	double minabs;

	/*First, recover current time from parameters: */
	parameters->FindParam(&time,TimeEnum);

	/*Retrieve interpolated values for this time step: */
	Input* input=GetTimeInput(time);

	/*Call input function*/
	minabs=input->MinAbs();
			   
	/*Clean-up and return*/
	delete input;
	return minabs;
}
/*}}}*/
/*FUNCTION TransientInput::GetVectorFromInputs{{{1*/
void TransientInput::GetVectorFromInputs(Vector* vector,int* doflist){

	double time;

	/*First, recover current time from parameters: */
	parameters->FindParam(&time,TimeEnum);

	/*Retrieve interpolated values for this time step: */
	Input* input=GetTimeInput(time);
		   
	/*Call input function*/
	input->GetVectorFromInputs(vector,doflist);
			   
	delete input;

} /*}}}*/
/*FUNCTION TransientInput::GetTimeInput{{{1*/
Input* TransientInput::GetTimeInput(double intime){

	int     i,j;
	double  deltat;
	double  alpha1,alpha2;
	bool    found=false;
	Input*  input=NULL;
	Input*  input1=NULL;
	Input*  input2=NULL;

	/*Ok, we have the time, go through the timesteps, and figure out which interval we 
	 *fall within. Then interpolate the values on this interval: */
	if(intime<this->timesteps[0]){
		/*get values for the first time: */
		input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
		found=true;
	}
	else if(intime>this->timesteps[this->numtimesteps-1]){
		/*get values for the last time: */
		input=(Input*)((Input*)this->inputs->GetObjectByOffset(numtimesteps-1))->copy();
		found=true;
	}
	else{
		/*Find which interval we fall within: */
		for(i=0;i<this->numtimesteps;i++){
			if(intime==this->timesteps[i]){
				/*We are right on one step time: */
				input=(Input*)((Input*)this->inputs->GetObjectByOffset(i))->copy();
				found=true;
				break; //we are done with the time interpolation.
			}
			else{
				if(this->timesteps[i]<intime && intime<this->timesteps[i+1]){
					/*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */
					deltat=this->timesteps[i+1]-this->timesteps[i];
					alpha2=(intime-this->timesteps[i])/deltat;
					alpha1=(1-alpha2);

					input1=(Input*)this->inputs->GetObjectByOffset(i); 
					input2=(Input*)this->inputs->GetObjectByOffset(i+1);

					input=(Input*)input1->copy();
					input->Scale(alpha1);
					input->AXPY(input2,alpha2);

					found=true;
					break;
				}
				else continue; //keep looking on the next interval
			}
		}
	}
	if(!found)_error_("did not find time interval on which to interpolate forcing values!");

	/*Assign output pointer*/
	return input;
}
/*}}}*/
/*FUNCTION TransientInput::Configure{{{1*/
void TransientInput::Configure(Parameters* parameters){
	this->parameters=parameters;
}
/*}}}*/