/*!\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 "../classes.h"
#include "shared/shared.h"
/*}}}*/	

/*Constructors/destructors*/
Friction::Friction(){/*{{{*/
	this->element=NULL;
	this->dim=0;
	this->law=0;

}
/*}}}*/
Friction::Friction(Element* element_in,int dim_in){/*{{{*/

	this->element=element_in;
	this->dim=dim_in;
	element_in->FindParam(&this->law,FrictionLawEnum);
}
/*}}}*/
Friction::~Friction(){/*{{{*/
}
/*}}}*/

/*methods: */
void Friction::Echo(void){/*{{{*/
	_printf_("Friction:\n");
	_printf_("   dim: " << this->dim<< "\n");
}
/*}}}*/
void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	switch(this->law){
		case 1:
			GetAlpha2Viscous(palpha2,gauss);
			break;
		case 2:
			GetAlpha2Weertman(palpha2,gauss);
			break;
		default:
			_error_("not supported");
	}

}/*}}}*/
void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/*This routine calculates the basal friction coefficient 
	  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), 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  Neff;
	IssmDouble  thickness,bed;
	IssmDouble  vx,vy,vz,vmag;
	IssmDouble  drag_coefficient;
	IssmDouble  alpha2;

	/*Recover parameters: */
	element->GetInputValue(&drag_p,FrictionPEnum);
	element->GetInputValue(&drag_q,FrictionQEnum);
	element->GetInputValue(&thickness, gauss,ThicknessEnum);
	element->GetInputValue(&bed, gauss,BaseEnum);
	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
	IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
	IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);

	//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(Neff<0)Neff=0;

	switch(dim){
		case 1:
			element->GetInputValue(&vx,gauss,VxEnum);
			vmag=sqrt(vx*vx);
			break;
		case 2:
			element->GetInputValue(&vx,gauss,VxEnum);
			element->GetInputValue(&vy,gauss,VyEnum);
			vmag=sqrt(vx*vx+vy*vy);
			break;
		case 3:
			element->GetInputValue(&vx,gauss,VxEnum);
			element->GetInputValue(&vy,gauss,VyEnum);
			element->GetInputValue(&vz,gauss,VzEnum);
			vmag=sqrt(vx*vx+vy*vy+vz*vz);
			break;
		default:
			_error_("not supported");
	}

	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0. && (s-1.)<0.) alpha2=0.;
	else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
	_assert_(!xIsNan<IssmDouble>(alpha2));

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */

	/*diverse: */
	IssmDouble  C,m;
	IssmDouble  vx,vy,vz,vmag;
	IssmDouble  alpha2;

	/*Recover parameters: */
	element->GetInputValue(&C,gauss,FrictionCEnum);
	element->GetInputValue(&m,FrictionMEnum);

	switch(dim){
		case 1:
			element->GetInputValue(&vx,gauss,VxEnum);
			vmag=sqrt(vx*vx);
			break;
		case 2:
			element->GetInputValue(&vx,gauss,VxEnum);
			element->GetInputValue(&vy,gauss,VyEnum);
			vmag=sqrt(vx*vx+vy*vy);
			break;
		case 3:
			element->GetInputValue(&vx,gauss,VxEnum);
			element->GetInputValue(&vy,gauss,VyEnum);
			element->GetInputValue(&vz,gauss,VzEnum);
			vmag=sqrt(vx*vx+vy*vy+vz*vz);
			break;
		default:
			_error_("not supported");
	}

	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
	else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));
	_assert_(!xIsNan<IssmDouble>(alpha2));

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/

	/* 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*/

	if(this->law!=1)_error_("not supported");

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

	/*Recover parameters: */
	element->GetInputValue(&drag_p,FrictionPEnum);
	element->GetInputValue(&drag_q,FrictionQEnum);
	element->GetInputValue(&thickness, gauss,ThicknessEnum);
	element->GetInputValue(&bed, gauss,BaseEnum);
	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
	IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
	IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);

	//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(Neff<0)Neff=0;

	//We need the velocity magnitude to evaluate the basal stress:
	switch(dim){
		case 1:
			element->GetInputValue(&vx,gauss,VxEnum);
			vmag=sqrt(vx*vx);
			break;
		case 2:
			element->GetInputValue(&vx,gauss,VxEnum);
			element->GetInputValue(&vy,gauss,VyEnum);
			vmag=sqrt(vx*vx+vy*vy);
			break;
		case 3:
			element->GetInputValue(&vx,gauss,VxEnum);
			element->GetInputValue(&vy,gauss,VyEnum);
			element->GetInputValue(&vz,gauss,VzEnum);
			vmag=sqrt(vx*vx+vy*vy+vz*vz);
			break;
		default:
			_error_("not supported");
	}

	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
	else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan<IssmDouble>(alpha_complement));

	/*Assign output pointers:*/
	*palpha_complement=alpha_complement;
}
/*}}}*/
