/*!\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*/
/*FUNCTION Friction::Friction() {{{*/
Friction::Friction(){
	this->element=NULL;
	this->dim=0;
}
/*}}}*/
/*FUNCTION Friction::Friction(Element* element,int dim){{{*/
Friction::Friction(Element* element_in,int dim_in){

	this->element=element_in;
	this->dim=dim_in;
}
/*}}}*/
/*FUNCTION Friction::~Friction() {{{*/
Friction::~Friction(){
}
/*}}}*/

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

	/*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  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(MaterialsRhoWaterEnum);
	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:
			vx_input->GetInputValue(&vx,gauss);
			vmag=sqrt(vx*vx);
			break;
		case 2:
			vx_input->GetInputValue(&vx,gauss);
			vy_input->GetInputValue(&vy,gauss);
			vmag=sqrt(vx*vx+vy*vy);
			break;
		case 3:
			vx_input->GetInputValue(&vx,gauss);
			vy_input->GetInputValue(&vy,gauss);
			vz_input->GetInputValue(&vz,gauss);
			vmag=sqrt(vx*vx+vy*vy+vz*vz);
			break;
		default:
			_error_("not supported");
	}

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

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

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

	/* 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  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(MaterialsRhoWaterEnum);
	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:
			vx_input->GetInputValue(&vx,gauss);
			vmag=sqrt(vx*vx);
			break;
		case 2:
			vx_input->GetInputValue(&vx,gauss);
			vy_input->GetInputValue(&vy,gauss);
			vmag=sqrt(vx*vx+vy*vy);
			break;
		case 3:
			vx_input->GetInputValue(&vx,gauss);
			vy_input->GetInputValue(&vy,gauss);
			vz_input->GetInputValue(&vz,gauss);
			vmag=sqrt(vx*vx+vy*vy+vz*vz);
			break;
		default:
			_error_("not supported");
	}

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

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

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