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

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

#include "./Matestar.h"
#include "./Materials.h"
#include "../Inputs/Input.h"
#include "../Inputs/Inputs.h"
#include "../Inputs/TriaInput.h"
#include "../Inputs/PentaInput.h"
#include "../Inputs/ControlInput.h"
#include "../Elements/Element.h"
#include "../Elements/Tria.h"
#include "../Elements/Penta.h"
#include "../Params/Parameters.h"
#include "../Vertex.h"
#include "../Hook.h"
#include "../Node.h"
#include "../IoModel.h"
#include "../../shared/shared.h"

/*Matestar constructors and destructor*/
Matestar::Matestar(){/*{{{*/
	this->helement=NULL;
	this->element=NULL;
	return;
}
/*}}}*/
Matestar::Matestar(int matestar_mid,int index, IoModel* iomodel){/*{{{*/

	/*Intermediaries:*/
	int    matestar_eid;

	/*Initialize id*/
	this->mid=matestar_mid;

	/*Hooks: */
	matestar_eid=index+1;
	this->helement=new Hook(&matestar_eid,1);
	this->element=NULL;

	return;
}
/*}}}*/
Matestar::~Matestar(){/*{{{*/
	delete helement;
	return;
}
/*}}}*/

/*Object virtual functions definitions:*/
Object*   Matestar::copy() {/*{{{*/

	/*Output*/
	Matestar* matestar=NULL;

	/*Initialize output*/
	matestar=new Matestar();

	/*copy fields: */
	matestar->mid=this->mid;
	matestar->helement=(Hook*)this->helement->copy();
	matestar->element =(Element*)this->helement->delivers();

	return matestar;
}
/*}}}*/
Material* Matestar::copy2(Element* element_in) {/*{{{*/

	/*Output*/
	Matestar* matestar=NULL;

	/*Initialize output*/
	matestar=new Matestar();

	/*copy fields: */
	matestar->mid=this->mid;
	matestar->helement=(Hook*)this->helement->copy();
	matestar->element =element_in;

	return matestar;
}
/*}}}*/
void      Matestar::DeepEcho(void){/*{{{*/

	_printf_("Matestar:\n");
	_printf_("   mid: " << mid << "\n");
	_printf_("   element:\n");
	helement->Echo();
}		
/*}}}*/
void      Matestar::Echo(void){/*{{{*/

	_printf_("Matestar:\n");
	_printf_("   mid: " << mid << "\n");
	_printf_("   element:\n");
	helement->Echo();
}
/*}}}*/
int       Matestar::Id(void){ return mid; }/*{{{*/
/*}}}*/
int       Matestar::ObjectEnum(void){/*{{{*/

	return MatestarEnum;

}
/*}}}*/
void      Matestar::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/

	if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook(); 
	
	MARSHALLING_ENUM(MatestarEnum);
	MARSHALLING(mid);
	this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
	this->element=(Element*)this->helement->delivers();

}
/*}}}*/

/*Matestar management*/
void  Matestar::Configure(Elements* elementsin){/*{{{*/

	/*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 
	 * datasets, using internal ids and offsets hidden in hooks: */
	helement->configure((DataSet*)elementsin);
	this->element  = (Element*)helement->delivers();
}
/*}}}*/
void  Matestar::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/

}
/*}}}*/
IssmDouble Matestar::GetA(){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetAbar(){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetB(){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetBbar(){/*{{{*/

	_error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetN(){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetD(){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
IssmDouble Matestar::GetDbar(){/*{{{*/

	_error_("not implemented yet");
}
/*}}}*/
bool Matestar::IsDamage(){/*{{{*/

	_error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
	_error_("not implemented yet");
	/*From a string tensor and a material object, return viscosity, using Glen's flow law.
								(1-D) B
	  viscosity= -------------------------
						  2 eps_eff ^[(n-1)/n]

	  where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
	  and n the flow law exponent.

	  If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we 
	  return 10^14, initial viscosity.
	  */

	/*output: */
	IssmDouble viscosity;

	/*Intermediary: */
	IssmDouble B,D=0.,n;

	/*Get B and n*/
	B=GetB(); _assert_(B>0.);
	n=GetN(); _assert_(n>0.);

	if (n==1.){
		/*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
		viscosity=(1.-D)*B/2.;
	}
	else{

		/*if no strain rate, return maximum viscosity*/
		if(eps_eff==0.){
			viscosity = 1.e+14/2.;
			//viscosity = B;
			//viscosity=2.5*pow(10.,17);
		}

		else{
			viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
		}
	}

	/*Checks in debugging mode*/
	if(viscosity<=0) _error_("Negative viscosity");

	/*Return: */
	*pviscosity=viscosity;
}
/*}}}*/
void  Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){/*{{{*/
	 _error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
void  Matestar::GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
	_error_("not implemented yet");
}
/*}}}*/
void  Matestar::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/

}
/*}}}*/
void  Matestar::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/

}
/*}}}*/
void  Matestar::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols,int name, int type){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Matestar::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Matestar::InputUpdateFromConstant(int constant, int name){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Matestar::InputUpdateFromConstant(bool constant, int name){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble vx,vy,vz;
	IssmDouble dvx[3],dvy[3],dvz[3];
	IssmDouble ko,Ec,Es;

	/*Get velocity derivatives in all directions*/
	_assert_(dim>1);
	_assert_(vx_input);
	vx_input->GetInputValue(&vx,gauss);
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	_assert_(vy_input);
	vy_input->GetInputValue(&vy,gauss);
	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
	if(dim==3){
		_assert_(vz_input);
		vz_input->GetInputValue(&vz,gauss);
		vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
	}
	else{
		vz = 0.;
		dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
	}

	/*Get enhancement factors and ko*/
	Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
	Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
	ec_input->GetInputValue(&Ec,gauss);
	es_input->GetInputValue(&Es,gauss);
	ko_input->GetInputValue(&ko,gauss);

	/*Compute viscosity*/
	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
}
/*}}}*/
void  Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
	this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
}/*}}}*/
void  Matestar::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble vx,vy,vz;
	IssmDouble dvx[3],dvy[3],dvz[3];
	IssmDouble ko,Ec,Es;

	/*Get velocity derivatives in all directions*/
	_assert_(dim==2 || dim==3);
	_assert_(vx_input);
	vx_input->GetInputValue(&vx,gauss);
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	if(dim==3){
		_assert_(vy_input);
		vy_input->GetInputValue(&vy,gauss);
		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
	}
	else{
		dvx[2] = 0.;
		vy = 0.;
		dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
	}
	vz = 0.;
	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];

	/*Get enhancement factors and ko*/
	Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
	Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
	ec_input->GetInputValue(&Ec,gauss);
	es_input->GetInputValue(&Es,gauss);
	ko_input->GetInputValue(&ko,gauss);

	/*Compute viscosity*/
	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
}/*}}}*/
void  Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
	_error_("not implemented yet");
}/*}}}*/
void  Matestar::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
	_error_("not implemented yet");
}/*}}}*/
void  Matestar::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
	/*Intermediaries*/
	IssmDouble vx,vy,vz;
	IssmDouble dvx[3],dvy[3],dvz[3];
	IssmDouble ko,Ec,Es;

	/*Get velocity derivatives in all directions*/
	_assert_(dim==1 || dim==2);
	_assert_(vx_input);
	vx_input->GetInputValue(&vx,gauss);
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	if(dim==2){
		_assert_(vy_input);
		vy_input->GetInputValue(&vy,gauss);
		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
	}
	else{
		dvx[1] = 0.;
		dvx[2] = 0.;
		vy = 0.;
		dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
	}
	vz = 0.;
	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];

	/*Get enhancement factors and ko*/
	Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcbarEnum); _assert_(ec_input);
	Input* es_input = element->inputs->GetInput(MaterialsRheologyEsbarEnum); _assert_(es_input);
	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input);
	ec_input->GetInputValue(&Ec,gauss);
	es_input->GetInputValue(&Es,gauss);
	ko_input->GetInputValue(&ko,gauss);

	/*Compute viscosity*/
	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
}/*}}}*/
void  Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
	_error_("not implemented yet");
}/*}}}*/
void  Matestar::ResetHooks(){/*{{{*/

	this->element=NULL;

	/*Get Element type*/
	this->helement->reset();

}
/*}}}*/
IssmDouble Matestar::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/

	/*Intermediaries*/
	IssmDouble viscosity;
	IssmDouble E,lambdas;
	IssmDouble epso,epsprime_norm;
	IssmDouble vmag,dvmag[3];

	/*Calculate velocity magnitude and its derivative*/
	vmag = sqrt(vx*vx+vy*vy+vz*vz);
	if(vmag<1e-12){
		dvmag[0]=0;
		dvmag[1]=0;
		dvmag[2]=0;
	}
	else{
		dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
		dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
		dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
	}

	EstarStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
	lambdas=EstarLambdaS(epso,epsprime_norm);

	/*Get total enhancement factor E(lambdas)*/
	E = Ec + (Es-Ec)*lambdas*lambdas;

	/*Compute viscosity*/
	_assert_(E>0.); 
	_assert_(ko>0.); 
	_assert_(epso>0.); 
	viscosity = 1./(2*pow(ko*E*epso*epso,1./3.));

	/*Assign output pointer*/
	return viscosity;
}
/*}}}*/