/*!\file Matdamageice.c
 * \brief: implementation of the Matdamageice 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 "../../../include/include.h"

/*Matdamageice constructors and destructor*/
/*FUNCTION Matdamageice::Matdamageice(){{{*/
Matdamageice::Matdamageice(){
	this->inputs=NULL;
	this->helement=NULL;
	return;
}
/*}}}*/
/*FUNCTION Matdamageice::Matdamageice(int id, int index, IoModel* iomodel, int num_vertices){{{*/
Matdamageice::Matdamageice(int matice_mid,int index, IoModel* iomodel){

	/*Intermediaries:*/
	int    matice_eid;

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

	/*Initialize inputs*/
	this->inputs=new Inputs();

	/*Initialize inputs from IoModel*/
	this->InputUpdateFromIoModel(index,iomodel);

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

	return;

}
/*}}}*/
/*FUNCTION Matdamageice::~Matdamageice(){{{*/
Matdamageice::~Matdamageice(){
	delete helement;
	delete inputs;
	return;
}
/*}}}*/

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

	_printLine_("Matdamageice:");
	_printLine_("   mid: " << mid);
	_printLine_("   inputs:");
	inputs->Echo();
	_printLine_("   element:");
	helement->Echo();
}
/*}}}*/
/*FUNCTION Matdamageice::DeepEcho {{{*/
void Matdamageice::DeepEcho(void){

	_printLine_("Matdamageice:");
	_printLine_("   mid: " << mid);
	_printLine_("   inputs:");
	inputs->DeepEcho();
	_printLine_("   element:");
	helement->Echo();
}		
/*}}}*/
/*FUNCTION Matdamageice::Id {{{*/
int    Matdamageice::Id(void){ return mid; }
/*}}}*/
/*FUNCTION Matdamageice::ObjectEnum{{{*/
int Matdamageice::ObjectEnum(void){

	return MatdamageiceEnum;

}
/*}}}*/
/*FUNCTION Matdamageice::copy {{{*/
Object* Matdamageice::copy() {

	/*Output*/
	Matdamageice* matice=NULL;

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

	/*copy fields: */
	matice->mid=this->mid;
	matice->helement=(Hook*)this->helement->copy();
	if(this->inputs) matice->inputs=(Inputs*)this->inputs->Copy();
	else  matice->inputs=new Inputs();

	return matice;
}
/*}}}*/

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

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

	IssmDouble B,n;

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

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

	/*Output*/
	IssmDouble B;

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

	/*Output*/
	IssmDouble Bbar;

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

	/*Output*/
	IssmDouble n;

	inputs->GetInputAverage(&n,MaterialsRheologyNEnum);
	return n;
}
/*}}}*/
/*FUNCTION Matdamageice::GetZ {{{*/
IssmDouble Matdamageice::GetZ(){

	/*Output*/
	IssmDouble Z;

	inputs->GetInputAverage(&Z,MaterialsRheologyZEnum);
	return Z;
}
/*}}}*/
/*FUNCTION Matdamageice::GetZbar {{{*/
IssmDouble Matdamageice::GetZbar(){

	/*Output*/
	IssmDouble Zbar;

	inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum);
	return Zbar;
}
/*}}}*/
/*FUNCTION Matdamageice::GetVectorFromInputs{{{*/
void  Matdamageice::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){

	/*Intermediaries*/
	Element *element= NULL;

	/*Recover element*/
	element=(Element*)helement->delivers();

	/*Check that input_enum is a material input*/
	if (!IsInput(input_enum)) return;

	switch(element->ObjectEnum()){

		case TriaEnum:{

			/*Prepare index list*/
			int doflist1[3];
			for(int i=0;i<3;i++) doflist1[i]=((Tria*)element)->nodes[i]->GetVertexPid();

			/*Get input (either in element or material)*/
			Input* input=inputs->GetInput(input_enum);
			if(!input) _error_("Input " << EnumToStringx(input_enum) << " not found in material");

			/*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
			input->GetVectorFromInputs(vector,&doflist1[0]);}
			break;

		default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet");
	}
}
/*}}}*/
/*FUNCTION Matdamageice::GetViscosity2d {{{*/
void  Matdamageice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){
	/*From a string tensor and a material object, return viscosity, using Glen's flow law.
												   Z * B
	  viscosity= -------------------------------------------------------------------
						  2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]

	  where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity 
	  vector, and n the flow law exponent.

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

	/*output: */
	IssmDouble viscosity;

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

	/*Intermediary: */
	IssmDouble A,e;
	IssmDouble Btmp,B,n,Z;

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

	if (n==1){
		/*Viscous behaviour! viscosity=B: */
		viscosity=B/2;
	}
	else{
		if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
			viscosity=0.5*pow((IssmDouble)10,(IssmDouble)14);
		}
		else{
			/*Retrive strain rate components: */
			exx=*(epsilon+0);
			eyy=*(epsilon+1);
			exy=*(epsilon+2);

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

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

	/*Return: */
	*pviscosity=viscosity;
}
/*}}}*/
/*FUNCTION Matdamageice::GetViscosity3d {{{*/
void  Matdamageice::GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* epsilon){

	/*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 
	 *
	 *                                               B
	 * viscosity3d= -------------------------------------------------------------------
	 *                      2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
	 *
	 *     where mu is the viscotiy, B the flow law parameter , (u,v) the velocity 
	 *     vector, and n the flow law exponent.
	 *
	 * If epsilon is NULL, it means this is the first time Emg is being run, and we 
	 * return g, initial viscosity.
	 */

	/*output: */
	IssmDouble viscosity3d;

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

	/*Intermediaries: */
	IssmDouble A,e;
	IssmDouble B,n,Z;

	/*Get B, Z and n*/
	n=GetN();
	Z=GetZ();
	B=Z*GetB();

	if (n==1){
		/*Viscous behaviour! viscosity3d=B: */
		viscosity3d=B/2;
	}
	else{
		if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
				(epsilon[3]==0) && (epsilon[4]==0)){
			viscosity3d=0.5*pow((IssmDouble)10,(IssmDouble)14);
		}
		else{

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

			/*Build viscosity: viscosity3d=2*B/(2*A^e) */
			A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy;
			if(A==0){
				/*Maxiviscosity3dm viscosity for 0 shear areas: */
				viscosity3d=2.25*pow((IssmDouble)10,(IssmDouble)17);
			}
			else{
				e=(n-1)/2/n;

				viscosity3d=B/(2*pow(A,e));
			}
		}
	}

	/*Checks in debugging mode*/
	if(viscosity3d<=0) _error_("Negative viscosity");
	_assert_(B>0);
	_assert_(n>0);

	/*Assign output pointers:*/
	*pviscosity3d=viscosity3d;
}
/*}}}*/
/*FUNCTION Matdamageice::GetViscosity3dStokes {{{*/
void  Matdamageice::GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon){
	/*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 
	 *
	 *                                          B
	 * viscosity3d= -------------------------------------------------------------------
	 *                   2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
	 *
	 *     where mu is the viscotiy, B the flow law parameter , (u,v) the velocity 
	 *     vector, and n the flow law exponent.
	 *
	 * If epsilon is NULL, it means this is the first time Emg is being run, and we 
	 * return g, initial viscosity.
	 */

	/*output: */
	IssmDouble viscosity3d;

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

	/*Intermediaries: */
	IssmDouble A,e;
	IssmDouble B,n,Z;
	IssmDouble eps0;

	/*Get B and n*/
	eps0=pow((IssmDouble)10,(IssmDouble)-27);
	n=GetN();
	Z=GetZ();
	B=Z*GetB();

	if (n==1){
		/*Viscous behaviour! viscosity3d=B: */
		viscosity3d=B/2;
	}
	else{
		if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
				(epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){
			viscosity3d=0.5*pow((IssmDouble)10,(IssmDouble)14);
		}
		else{

			/*Retrive strain rate components: */
			exx=*(epsilon+0);
			eyy=*(epsilon+1);
			ezz=*(epsilon+2); //not used
			exy=*(epsilon+3);
			exz=*(epsilon+4);
			eyz=*(epsilon+5);

			/*Build viscosity: viscosity3d=B/(2*A^e) */
			A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy+pow(eps0,2);
			if(A==0){
				/*Maxiviscosity3dm viscosity for 0 shear areas: */
				viscosity3d=2.25*pow((IssmDouble)10,(IssmDouble)17);
			}
			else{
				e=(n-1)/2/n;
				viscosity3d=B/(2*pow(A,e));
			}
		}
	}

	/*Checks in debugging mode*/
	if(viscosity3d<=0) _error_("Negative viscosity");
	_assert_(B>0);
	_assert_(n>0);

	/*Assign output pointers:*/
	*pviscosity3d=viscosity3d;
}
/*}}}*/
/*FUNCTION Matdamageice::GetViscosityComplement {{{*/
void  Matdamageice::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){
	/*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 
	 *
	 *  										                1
	 * 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 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((IssmDouble)10,(IssmDouble)17);
		}
		else{
			e=(n-1)/(2*n);

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

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

	/*Return: */
	*pviscosity_complement=viscosity_complement;
}
/*}}}*/
/*FUNCTION Matdamageice::GetViscosityZComplement {{{*/
void  Matdamageice::GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){
	/*Return viscosity derivative for control method d(mu)/dZ: 
	 *
	 *  										               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((IssmDouble)10,(IssmDouble)17);
		}
		else{
			e=(n-1)/(2*n);

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

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

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

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

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

	/*Get visocisty and n*/
	GetViscosity3d(&mu,epsilon);
	n=GetN();

	if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
				(epsilon[3]==0) && (epsilon[4]==0)){
		mu_prime=0.5*pow((IssmDouble)10,(IssmDouble)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;

		mu_prime=(1-n)/(2*n) * mu/eff2;
	}

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

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

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

	/*Get visocisty and n*/
	GetViscosity2d(&mu,epsilon);
	n=GetN();

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

		mu_prime=(1-n)/(2*n) * mu/eff2;
	}

	/*Assign output pointers:*/
	*pmu_prime=mu_prime;
}
/*}}}*/
/*FUNCTION Matdamageice::InputDuplicate{{{*/
void  Matdamageice::InputDuplicate(int original_enum,int new_enum){

	/*Call inputs method*/
	if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);

}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromVector(IssmDouble* vector, int name, int type) {{{*/
void  Matdamageice::InputUpdateFromVector(IssmDouble* vector, int name, int type){

	/*Intermediaries*/
	Element *element      = NULL;

	/*Recover element*/
	element=(Element*)helement->delivers();

	/*Check that name is an element input*/
	if (!IsInput(name)) return;

	switch(type){

		case VertexEnum:

			switch(element->ObjectEnum()){

				case TriaEnum: {
					IssmDouble values[3];
					for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->nodes[i]->GetVertexPid()];
					this->inputs->AddInput(new TriaP1Input(name,values));
					return;
				}
				default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet");
			}
		default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
	}
}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromVector(int* vector, int name, int type) {{{*/
void  Matdamageice::InputUpdateFromVector(int* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromVector(bool* vector, int name, int type) {{{*/
void  Matdamageice::InputUpdateFromVector(bool* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type) {{{*/
void  Matdamageice::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){

	/*Intermediaries*/
	Element *element      = NULL;
	Parameters* parameters= NULL;
	int         dim;

	/*Recover element*/
	element=(Element*)helement->delivers();

	/*Check that name is an element input*/
	if (!IsInput(name)) return;

	switch(type){

		case VertexEnum:

			switch(element->ObjectEnum()){

				case TriaEnum: {
					IssmDouble values[3];
					for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->nodes[i]->GetVertexSid()]; //use sid list, to index into serial oriented vector 
					this->inputs->AddInput(new TriaP1Input(name,values));
					/*Special case for rheology B in 2D: Pourave land for this solution{{{*/
					if(name==MaterialsRheologyBEnum){
						/*Are we in 2D?:*/
						if(element->ObjectEnum()==TriaEnum){
							parameters=((Tria*)(element))->parameters;
						}
						else{
							parameters=((Penta*)(element))->parameters;
						}
						parameters->FindParam(&dim,MeshDimensionEnum);
						if(dim==2){
							/*Dupliacte rheology input: */
							this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,values));
						}
					}
					/*}}}*/
					return;
				}
				default: _error_("element " << EnumToStringx(element->ObjectEnum()) << " not implemented yet");
			}
		default: _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
	}

}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromMatrixDakota(int* vector, int name, int type) {{{*/
void  Matdamageice::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols,int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromVectorDakota(int* vector, int name, int type) {{{*/
void  Matdamageice::InputUpdateFromVectorDakota(int* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromVectorDakota(bool* vector, int name, int type) {{{*/
void  Matdamageice::InputUpdateFromVectorDakota(bool* vector, int name, int type){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromConstant(IssmDouble constant, int name) {{{*/
void  Matdamageice::InputUpdateFromConstant(IssmDouble constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromConstant(int constant, int name) {{{*/
void  Matdamageice::InputUpdateFromConstant(int constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromConstant(bool constant, int name) {{{*/
void  Matdamageice::InputUpdateFromConstant(bool constant, int name){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromSolution{{{*/
void  Matdamageice::InputUpdateFromSolution(IssmDouble* solution){
	/*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Matdamageice::InputUpdateFromIoModel{{{*/
void Matdamageice::InputUpdateFromIoModel(int index, IoModel* iomodel){

	int i,j;

	int    dim;
	bool   control_analysis;
	int    num_control_type;

	/*Fetch parameters: */
	iomodel->Constant(&dim,MeshDimensionEnum);
	iomodel->Constant(&control_analysis,InversionIscontrolEnum);
	if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);

	/*if 2d*/
	if(dim==2){

		/*Intermediaries*/
		const int num_vertices = 3; //Tria has 3 vertices
		IssmDouble    nodeinputs[num_vertices];
		IssmDouble    cmmininputs[num_vertices];
		IssmDouble    cmmaxinputs[num_vertices];

		/*Get B*/
		if (iomodel->Data(MaterialsRheologyBEnum)) {
			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
			this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,nodeinputs));
		}

		/*Get n*/
		if (iomodel->Data(MaterialsRheologyNEnum)) {
			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
			this->inputs->AddInput(new TriaP1Input(MaterialsRheologyNEnum,nodeinputs));
		}

		/*Get Z*/
		if (iomodel->Data(MaterialsRheologyZEnum)) {
			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
			this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs));
		}

		/*Control Inputs*/
		#ifdef _HAVE_CONTROL_
		if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
			for(i=0;i<num_control_type;i++){
				switch(reCast<int>(iomodel->Data(InversionControlParametersEnum)[i])){
					case MaterialsRheologyBbarEnum:
						if (iomodel->Data(MaterialsRheologyBEnum)){
							_assert_(iomodel->Data(MaterialsRheologyBEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
							for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
							for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
							for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
							this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
						}
						break;
					case MaterialsRheologyZbarEnum:
						if (iomodel->Data(MaterialsRheologyZEnum)){
							_assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
							for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
							for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
							for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
							this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
						}
						break;

				}
			}
		}
		#endif
	}

	/*if 3d*/
	#ifdef _HAVE_3D_
	else if(dim==3){

		/*Intermediaries*/
		const int num_vertices = 6; //Penta has 6 vertices
		IssmDouble    nodeinputs[num_vertices];
		IssmDouble    cmmininputs[num_vertices];
		IssmDouble    cmmaxinputs[num_vertices];

		/*Get B*/
		if (iomodel->Data(MaterialsRheologyBEnum)) {
			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
			this->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,nodeinputs));
		}

		/*Get n*/
		if (iomodel->Data(MaterialsRheologyNEnum)) {
			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
			this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs));
		}

		/*Get Z*/
		if (iomodel->Data(MaterialsRheologyZEnum)) {
			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
			this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs));
		}

		/*Control Inputs*/
		#ifdef _HAVE_CONTROL_
		if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
			for(i=0;i<num_control_type;i++){
				switch(reCast<int>(iomodel->Data(InversionControlParametersEnum)[i])){
					case MaterialsRheologyBbarEnum:
						if (iomodel->Data(MaterialsRheologyBEnum)){
							_assert_(iomodel->Data(MaterialsRheologyBEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
							for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
							for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
							for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
							this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
						}
						break;
					case MaterialsRheologyZbarEnum:
						if (iomodel->Data(MaterialsRheologyZEnum)){
							_assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
							for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
							for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
							for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
							this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
						}
						break;
				}
			}
		}
		#endif
	}
	#endif
	else{
		_error_("Mesh type not supported yet!");
	}

	return;
}
/*}}}*/
/*FUNCTION Matdamageice::IsInput{{{*/
bool Matdamageice::IsInput(int name){
	if (
				name==MaterialsRheologyBEnum ||
				name==MaterialsRheologyBbarEnum ||
				name==MaterialsRheologyNEnum ||
				name==MaterialsRheologyZEnum ||
				name==MaterialsRheologyZbarEnum
		){
		return true;
	}
	else return false;
}
/*}}}*/
