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

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

#include "./Matice.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"

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

	/*Intermediaries:*/
	int    matice_eid;

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

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

	 /*Other perporties*/
   int    materialtype;
   iomodel->Constant(&materialtype,MaterialsEnum);
   if(materialtype==MatdamageiceEnum) this->isdamaged = true;
   else if(materialtype==MaticeEnum) this->isdamaged = false;
   else _error_("Material type not recognized");

	return;

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

/*Object virtual functions definitions:*/
void Matice::Echo(void){/*{{{*/

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

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

	return MaticeEnum;

}
/*}}}*/
Object* Matice::copy() {/*{{{*/

	/*Output*/
	Matice* matice=NULL;

	/*Initialize output*/
	matice=new Matice();

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

	return matice;
}
/*}}}*/
Material* Matice::copy2(Element* element_in) {/*{{{*/

	/*Output*/
	Matice* matice=NULL;

	/*Initialize output*/
	matice=new Matice();

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

	return matice;
}
/*}}}*/

/*Matice management*/
void  Matice::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  Matice::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/

}
/*}}}*/
IssmDouble Matice::GetA(){/*{{{*/
	/*
	 * A = 1/B^n
	 */

	IssmDouble B,n;

	element->inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
	n=this->GetN();

	return pow(B,-n);
}
/*}}}*/
IssmDouble Matice::GetAbar(){/*{{{*/
	/*
	 * A = 1/B^n
	 */

	IssmDouble B,n;

	element->inputs->GetInputAverage(&B,MaterialsRheologyBbarEnum);
	n=this->GetN();

	return pow(B,-n);
}
/*}}}*/
IssmDouble Matice::GetB(){/*{{{*/

	/*Output*/
	IssmDouble B;

	element->inputs->GetInputAverage(&B,MaterialsRheologyBEnum);
	return B;
}
/*}}}*/
IssmDouble Matice::GetBbar(){/*{{{*/

	/*Output*/
	IssmDouble Bbar;

	element->inputs->GetInputAverage(&Bbar,MaterialsRheologyBbarEnum);
	return Bbar;
}
/*}}}*/
IssmDouble Matice::GetN(){/*{{{*/

	/*Output*/
	IssmDouble n;

	element->inputs->GetInputAverage(&n,MaterialsRheologyNEnum);
	return n;
}
/*}}}*/
IssmDouble Matice::GetD(){/*{{{*/

	_assert_(this->isdamaged);
	/*Output*/
	IssmDouble D;
	if(this->isdamaged)element->inputs->GetInputAverage(&D,DamageDEnum);
	return D;
}
/*}}}*/
IssmDouble Matice::GetDbar(){/*{{{*/

	_assert_(this->isdamaged);
	/*Output*/
	IssmDouble Dbar;
	if(this->isdamaged)element->inputs->GetInputAverage(&Dbar,DamageDbarEnum);
	return Dbar;
}
/*}}}*/
bool Matice::IsDamage(){/*{{{*/

	return this->isdamaged;
}
/*}}}*/
void  Matice::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
	/*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(this->isdamaged){
		D=GetD();
		_assert_(D>=0. && D<1.);
	}

	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  Matice::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/

	/*output: */
	IssmDouble dmudB;

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

	/*Get B and n*/
	n=GetN(); _assert_(n>0.);
	if(this->isdamaged){
		D=GetD();
		_assert_(D>=0. && D<1.);
	}

	if(n==1.){
		/*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */
		dmudB=(1.-D)/2.;
	}
	else{
		if(eps_eff==0.) dmudB = 0.;
		else            dmudB = (1.-D)/(2.*pow(eps_eff,(n-1.)/n));
	}

	/*Return: */
	*pdmudB=dmudB;
}
/*}}}*/
void  Matice::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){/*{{{*/

	/*output: */
	IssmDouble dmudD;

	/*Intermediary: */
	IssmDouble n,B;

	/*Get B and n*/
	n=GetN(); _assert_(n>0.);
	B=GetBbar();
	_assert_(this->isdamaged);

	if(n==1.){
		/*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */
		dmudD=-B/2.;
	}
	else{
		if(eps_eff==0.) dmudD = 0.;
		else            dmudD = -B/(2.*pow(eps_eff,(n-1.)/n));
	}

	/*Return: */
	*pdmudD=dmudD;
}
/*}}}*/
void  Matice::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
	/*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=GetBbar(); _assert_(B>0.);
	n=GetN();    _assert_(n>0.);
	if(this->isdamaged){
		D=GetDbar();
		_assert_(D>=0. && D<1.);
	}

	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=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  Matice::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
	/*Return viscosity accounting for steady state power law creep [Thomas and SSA, 1982]: 
	 *
	 *  										                (1-D)
	 * viscosity= -------------------------------------------------------------------
	 *  				  2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
	 *
	 * If epsilon is NULL, it means this is the first time Gradjb is being run, and we 
	 * return mu20, initial viscosity.
	 */

	/*output: */
	IssmDouble viscosity_complement;

	/*input strain rate: */
	IssmDouble exx,eyy,exy;

	/*Intermediary value A and exponent e: */
	IssmDouble A,e;
	IssmDouble D=0.,n;

	/*Get D and n*/
	if(this->isdamaged){
		D=GetDbar(); /* GetD()? */
		_assert_(D>=0. && D<1.);
	}
	n=GetN();

	if(epsilon){
		exx=*(epsilon+0);
		eyy=*(epsilon+1);
		exy=*(epsilon+2);

		/*Build viscosity: mu2=(1-D)/(2*A^e) */
		A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
		if(A==0){
			/*Maximum viscosity_complement for 0 shear areas: */
			viscosity_complement=2.25*pow(10.,17);
		}
		else{
			e=(n-1)/(2*n);

			viscosity_complement=(1-D)/(2*pow(A,e));
		}
	}
	else{
		viscosity_complement=4.5*pow(10.,17);
	}

	/*Checks in debugging mode*/
	_assert_(D>=0 && D<1);
	_assert_(n>0);
	_assert_(viscosity_complement>0);

	/*Return: */
	*pviscosity_complement=viscosity_complement;
}
/*}}}*/
void  Matice::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
	/*Return viscosity derivative for control method d(mu)/dD: 
	 *
	 *  										               B 
	 * dviscosity= - -------------------------------------------------------------------
	 *  				  2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
	 *
	 * If epsilon is NULL, it means this is the first time Gradjb is being run, and we 
	 * return mu20, initial viscosity.
	 */

	/*output: */
	IssmDouble viscosity_complement;

	/*input strain rate: */
	IssmDouble exx,eyy,exy;

	/*Intermediary value A and exponent e: */
	IssmDouble A,e;
	IssmDouble B,n;

	/*Get B and n*/
	B=GetBbar();
	n=GetN();

	if(epsilon){
		exx=*(epsilon+0);
		eyy=*(epsilon+1);
		exy=*(epsilon+2);

		/*Build viscosity: mu2=B/(2*A^e) */
		A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
		if(A==0){
			/*Maximum viscosity_complement for 0 shear areas: */
			viscosity_complement=- 2.25*pow(10.,17);
		}
		else{
			e=(n-1)/(2*n);

			viscosity_complement=- B/(2*pow(A,e));
		}
	}
	else{
		viscosity_complement=- 4.5*pow(10.,17);
	}

	/*Checks in debugging mode*/
	_assert_(B>0);
	_assert_(n>0);
	_assert_(viscosity_complement<0);

	/*Return: */
	*pviscosity_complement=viscosity_complement;
}
/*}}}*/
void  Matice::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/

	/*output: */
	IssmDouble mu_prime;
	IssmDouble mu,n,eff2;

	/*input strain rate: */
	IssmDouble exx,eyy,exy,exz,eyz;


	if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
				(epsilon[3]==0) && (epsilon[4]==0)){
		mu_prime=0.5*pow(10.,14);
	}
	else{

		/*Retrive strain rate components: */
		exx=epsilon[0];
		eyy=epsilon[1];
		exy=epsilon[2];
		exz=epsilon[3];
		eyz=epsilon[4];
		eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;

		GetViscosity(&mu,sqrt(eff2));
		n=GetN();
		mu_prime=(1.-n)/(2.*n) * mu/eff2;
	}

	/*Assign output pointers:*/
	*pmu_prime=mu_prime;
}
/*}}}*/
void  Matice::GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/

	/*output: */
	IssmDouble mu_prime;
	IssmDouble mu,n,eff2;

	/*input strain rate: */
	IssmDouble exx,eyy,exy;

	if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
		mu_prime=0.5*pow(10.,14);
	}
	else{
		/*Retrive strain rate components: */
		exx=epsilon[0];
		eyy=epsilon[1];
		exy=epsilon[2];
		eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy ;

		GetViscosityBar(&mu,sqrt(eff2));
		n=GetN();
		mu_prime=(1.-n)/(2.*n)*mu/eff2;
	}

	/*Assign output pointers:*/
	*pmu_prime=mu_prime;
}
/*}}}*/
void  Matice::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/

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

}
/*}}}*/
void  Matice::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols,int name, int type){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Matice::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Matice::InputUpdateFromConstant(int constant, int name){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Matice::InputUpdateFromConstant(bool constant, int name){/*{{{*/
	/*Nothing updated yet*/
}
/*}}}*/
void  Matice::ResetHooks(){/*{{{*/

	this->element=NULL;

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

}
/*}}}*/
