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

/*Headers:*/
/*{{{*/
#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 "../../../Container/Container.h"
#include "../../../EnumDefinitions/EnumDefinitions.h"
#include "../../../shared/shared.h"
#include "../../../include/include.h"
/*}}}*/	

/*Constructors/destructors*/
/*FUNCTION Friction::Friction() {{{*/
Friction::Friction(){
	this->element_type=NULL;
	this->inputs=NULL;
	this->matpar=NULL;
}
/*}}}*/
/*FUNCTION Friction::Friction(const char* element_type, Inputs* inputs,Matpar* matpar,int analysis_type){{{*/
Friction::Friction(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int in_analysis_type){

	this->analysis_type=in_analysis_type;
	this->inputs=inputs_in;
	this->element_type=xNew<char>(strlen(element_type_in)+1);
	xMemCpy<char>(this->element_type,element_type_in,(strlen(element_type_in)+1));

	this->matpar=matpar_in;
}
/*}}}*/
/*FUNCTION Friction::~Friction() {{{*/
Friction::~Friction(){
	xDelete<char>(element_type);
}
/*}}}*/

/*methods: */
/*FUNCTION Friction::Echo {{{*/
void Friction::Echo(void){
	_printLine_("Friction:");
	_printLine_("   analysis_type: " << EnumToStringx(analysis_type));
	_printLine_("   element_type: " << this->element_type);
	inputs->Echo();
	matpar->Echo();
}
/*}}}*/
/*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{*/
void Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){

	/*This routine calculates the basal friction coefficient 
	  alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/

	/*diverse: */
	IssmDouble  r,s;
	IssmDouble  drag_p, drag_q;
	IssmDouble  gravity,rho_ice,rho_water;
	IssmDouble  Neff;
	IssmDouble  thickness,bed;
	IssmDouble  vx,vy,vz,vmag;
	IssmDouble  drag_coefficient;
	IssmDouble  alpha2;

	/*Recover parameters: */
	inputs->GetInputValue(&drag_p,FrictionPEnum);
	inputs->GetInputValue(&drag_q,FrictionQEnum);
	this->GetInputValue(&thickness, gauss,ThicknessEnum);
	this->GetInputValue(&bed, gauss,BedEnum);
	this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);

	/*Get material parameters: */
	gravity=matpar->GetG();
	rho_ice=matpar->GetRhoIce();
	rho_water=matpar->GetRhoWater();

	//compute r and q coefficients: */
	r=drag_q/drag_p;
	s=1./drag_p;

	//From bed and thickness, compute effective pressure when drag is viscous:
	Neff=gravity*(rho_ice*thickness+rho_water*bed);

	/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
	  the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
	  for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
	  replace it by Neff=0 (ie, equival it to an ice shelf)*/
	if (Neff<0)Neff=0;

	if(strcmp(element_type,"2d")==0){
		this->GetInputValue(&vx, gauss,vxenum);
		this->GetInputValue(&vy, gauss,vyenum);
		vmag=sqrt(pow(vx,2)+pow(vy,2));
	}
	else if (strcmp(element_type,"3d")==0){
		this->GetInputValue(&vx, gauss,vxenum);
		this->GetInputValue(&vy, gauss,vyenum);
		this->GetInputValue(&vz, gauss,vzenum);
		vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
	}
	else _error_("element_type "<< element_type << " not supported yet");

	/*Checks that s-1>0 if v=0*/
        if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf");

	alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
	_assert_(!xIsNan<IssmDouble>(alpha2));

	/*Assign output pointers:*/
	*palpha2=alpha2;
}
/*}}}*/
/*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{*/
void Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){

	/*This routine calculates the basal friction coefficient 
	  alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/

	/*diverse: */
	IssmDouble  r,s;
	IssmDouble  drag_p, drag_q;
	IssmDouble  gravity,rho_ice,rho_water;
	IssmDouble  Neff;
	IssmDouble  thickness,bed;
	IssmDouble  vx,vy,vz,vmag;
	IssmDouble  drag_coefficient;
	IssmDouble  alpha2;

	/*Recover parameters: */
	inputs->GetInputValue(&drag_p,FrictionPEnum);
	inputs->GetInputValue(&drag_q,FrictionQEnum);
	this->GetInputValue(&thickness, gauss,ThicknessEnum);
	this->GetInputValue(&bed, gauss,BedEnum);
	this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);

	/*Get material parameters: */
	gravity=matpar->GetG();
	rho_ice=matpar->GetRhoIce();
	rho_water=matpar->GetRhoWater();

	//compute r and q coefficients: */
	r=drag_q/drag_p;
	s=1./drag_p;

	//From bed and thickness, compute effective pressure when drag is viscous:
	Neff=gravity*(rho_ice*thickness+rho_water*bed);

	/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
	  the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
	  for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
	  replace it by Neff=0 (ie, equival it to an ice shelf)*/
	if (Neff<0)Neff=0;

	if(strcmp(element_type,"2d")==0){
		this->GetInputValue(&vx, gauss,vxenum);
		this->GetInputValue(&vy, gauss,vyenum);
		vmag=sqrt(pow(vx,2)+pow(vy,2));
	}
	else if (strcmp(element_type,"3d")==0){
		this->GetInputValue(&vx, gauss,vxenum);
		this->GetInputValue(&vy, gauss,vyenum);
		this->GetInputValue(&vz, gauss,vzenum);
		vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
	}
	else _error_("element_type "<< element_type << " not supported yet");

	/*Checks that s-1>0 if v=0*/
	if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf");

	alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
	_assert_(!xIsNan<IssmDouble>(alpha2));

	/*Assign output pointers:*/
	*palpha2=alpha2;
}
/*}}}*/
/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum) {{{*/
void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum){

	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
	 * FrictionGetAlphaComplement is used in control methods on drag, and it computes: 
	 * alpha_complement= Neff ^r * vel ^s*/

	/*diverse: */
	IssmDouble  r,s;
	IssmDouble  vx,vy,vz,vmag;
	IssmDouble  drag_p,drag_q;
	IssmDouble  Neff;
	IssmDouble  drag_coefficient;
	IssmDouble  bed,thickness;
	IssmDouble  gravity,rho_ice,rho_water;
	IssmDouble  alpha_complement;

	/*Recover parameters: */
	inputs->GetInputValue(&drag_p,FrictionPEnum);
	inputs->GetInputValue(&drag_q,FrictionQEnum);
	this->GetInputValue(&thickness, gauss,ThicknessEnum);
	this->GetInputValue(&bed, gauss,BedEnum);
	this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);

	/*Get material parameters: */
	gravity=matpar->GetG();
	rho_ice=matpar->GetRhoIce();
	rho_water=matpar->GetRhoWater();

	//compute r and q coefficients: */
	r=drag_q/drag_p;
	s=1./drag_p;

	//From bed and thickness, compute effective pressure when drag is viscous:
	Neff=gravity*(rho_ice*thickness+rho_water*bed);

	/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
	  the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
	  for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
	  replace it by Neff=0 (ie, equival it to an ice shelf)*/
	if (Neff<0)Neff=0;

	//We need the velocity magnitude to evaluate the basal stress:
	if(strcmp(element_type,"2d")==0){
		this->GetInputValue(&vx, gauss,vxenum);
		this->GetInputValue(&vy, gauss,vyenum);
		vmag=sqrt(pow(vx,2)+pow(vy,2));
	}
	else if (strcmp(element_type,"3d")==0){
		this->GetInputValue(&vx, gauss,vxenum);
		this->GetInputValue(&vy, gauss,vyenum);
		this->GetInputValue(&vz, gauss,vzenum);
		vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
	}
	else _error_("element_type "<< element_type << " not supported yet");

	/*Checks that s-1>0 if v=0*/
	if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf");

	alpha_complement=pow(Neff,r)*pow(vmag,(s-1));            _assert_(!xIsNan<IssmDouble>(alpha_complement));

	/*Assign output pointers:*/
	*palpha_complement=alpha_complement;
}
/*}}}*/
/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum) {{{*/
void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){

	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
	 * FrictionGetAlphaComplement is used in control methods on drag, and it computes: 
	 * alpha_complement= Neff ^r * vel ^s*/

	/*diverse: */
	IssmDouble  r,s;
	IssmDouble  vx,vy,vz,vmag;
	IssmDouble  drag_p,drag_q;
	IssmDouble  Neff;
	IssmDouble  drag_coefficient;
	IssmDouble  bed,thickness;
	IssmDouble  gravity,rho_ice,rho_water;
	IssmDouble  alpha_complement;

	/*Recover parameters: */
	inputs->GetInputValue(&drag_p,FrictionPEnum);
	inputs->GetInputValue(&drag_q,FrictionQEnum);
	this->GetInputValue(&thickness, gauss,ThicknessEnum);
	this->GetInputValue(&bed, gauss,BedEnum);
	this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);

	/*Get material parameters: */
	gravity=matpar->GetG();
	rho_ice=matpar->GetRhoIce();
	rho_water=matpar->GetRhoWater();

	//compute r and q coefficients: */
	r=drag_q/drag_p;
	s=1./drag_p;

	//From bed and thickness, compute effective pressure when drag is viscous:
	Neff=gravity*(rho_ice*thickness+rho_water*bed);

	/*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
	  the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
	  for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
	  replace it by Neff=0 (ie, equival it to an ice shelf)*/
	if (Neff<0)Neff=0;

	//We need the velocity magnitude to evaluate the basal stress:
	if(strcmp(element_type,"2d")==0){
		this->GetInputValue(&vx, gauss,vxenum);
		this->GetInputValue(&vy, gauss,vyenum);
		vmag=sqrt(pow(vx,2)+pow(vy,2));
	}
	else if (strcmp(element_type,"3d")==0){
		this->GetInputValue(&vx, gauss,vxenum);
		this->GetInputValue(&vy, gauss,vyenum);
		this->GetInputValue(&vz, gauss,vzenum);
		vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
	}
	else _error_("element_type "<< element_type << " not supported yet");

	/*Checks that s-1>0 if v=0*/
	if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf");

	alpha_complement=pow(Neff,r)*pow(vmag,(s-1));            _assert_(!xIsNan<IssmDouble>(alpha_complement));

	/*Assign output pointers:*/
	*palpha_complement=alpha_complement;
}
/*}}}*/
/*FUNCTION Friction::GetInputValue{{{*/
void Friction::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type){

	Input* input=inputs->GetInput(enum_type);
	if(!input) _error_("input " << EnumToStringx(enum_type) << " not found");
	input->GetInputValue(pvalue,gauss);

}
/*}}}*/
/*FUNCTION Friction::GetInputValue{{{*/
void Friction::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type){

	Input* input=inputs->GetInput(enum_type);
	if(!input) _error_("input " << EnumToStringx(enum_type) << " not found");
	input->GetInputValue(pvalue,gauss);

}
/*}}}*/
