/*!\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 "stdio.h"
#include "./Matice.h"
#include <string.h>
#include "../shared/shared.h"
#include "../EnumDefinitions/EnumDefinitions.h"

		
Matice::Matice(){
	return;
}
Matice::Matice(int matice_mid,double matice_B,double matice_n){
	mid=matice_mid;
	B=matice_B;
	n=matice_n;
	return;
}

Matice::~Matice(){
	return;
}
void Matice::Echo(void){

	printf("Matice:\n");
	printf("   mid: %i\n",mid);
	printf("   B: %g\n",B);
	printf("   n: %g\n",n);
	return;
}
		
void  Matice::Marshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   enum_type=0;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*get enum type of Matice: */
	enum_type=MaticeEnum();
	
	/*marshall enum: */
	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	
	/*marshall Matice data: */
	memcpy(marshalled_dataset,&mid,sizeof(mid));marshalled_dataset+=sizeof(mid);
	memcpy(marshalled_dataset,&B,sizeof(B));marshalled_dataset+=sizeof(B);
	memcpy(marshalled_dataset,&n,sizeof(n));marshalled_dataset+=sizeof(n);

	*pmarshalled_dataset=marshalled_dataset;
	return;
}
		
int   Matice::MarshallSize(){

	return sizeof(mid)+sizeof(B)+sizeof(n)+sizeof(int); //sizeof(int) for enum type
}

char* Matice::GetName(void){
	return "matice";
}
		

void  Matice::Demarshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
	 *object data (thanks to DataSet::Demarshall):*/

	memcpy(&mid,marshalled_dataset,sizeof(mid));marshalled_dataset+=sizeof(mid);
	memcpy(&B,marshalled_dataset,sizeof(B));marshalled_dataset+=sizeof(B);
	memcpy(&n,marshalled_dataset,sizeof(n));marshalled_dataset+=sizeof(n);

	/*return: */
	*pmarshalled_dataset=marshalled_dataset;
	return;
}

int Matice::Enum(void){

	return MaticeEnum();

}

int    Matice::GetId(void){ return mid; }

int    Matice::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
void  Matice::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}

#undef __FUNCT__ 
#define __FUNCT__ "Matice::UpdateFromInputs"
void  Matice::UpdateFromInputs(ParameterInputs* inputs){
	
	throw ErrorException(__FUNCT__," not supported yet!");

}
		
#undef __FUNCT__ 
#define __FUNCT__ "Matice::SetB"
void  Matice::SetB(double B_param){
	B=B_param;
}
		


#undef __FUNCT__ 
#define __FUNCT__ "MatIce::GetViscosity2d"
void  Matice::GetViscosity2d(double* pviscosity, double* epsilon){


	/*From a string tensor and a material object, return viscosity, using Glen's flow law.
	                                  2*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: */
	double viscosity;

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

	/*Intermediary value A and exponent e: */
	double A,e;

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

			/*Build viscosity: viscosity=2*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=4.5*pow(10,17);
			}
			else{
				e=(n-1)/2/n;
				viscosity=2*B/(2*pow(A,e));
			}
		}
	}
	#ifdef _DEBUG_
	_printf_("Viscosity %lf\n",viscosity);
	#endif

	/*Return: */
	*pviscosity=viscosity;
}

Object* Matice::copy() {
	return new Matice(*this); 
}


#undef __FUNCT__ 
#define __FUNCT__ "MatIce::GetViscosity2"
void  Matice::GetViscosity2(double* pviscosity2, double* epsilon){

	/*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 
	 *
	 *                                  2* (1-n)/2n
	 * mu2= -------------------------------------------------------------------
	 *     2[ (du/dx)^2+(dv/dy)^2+1/4*(du/dy+dv/dx)^2+du/dx*dv/dy ]^[(3n-1)/2n]
	 *
	 *     where mu2 is the second viscosity, (u,v) the velocity 
	 *     vector, and n the flow law exponent.
	 *
	 * If epsilon is NULL, it means this is the first time Gradjb is being run, and we 
	 * return mu20, initial viscosity.
	 */
	
	/*output: */
	double viscosity2;

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

	/*Intermediary value A and exponent e: */
	double A,e;

	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 viscosity2 for 0 shear areas: */
			viscosity2=4.5*pow(10,17);
		}
		else{
			e=(3*n-1)/2/n;
		
			viscosity2=2*(1-n)/2/n/(2*pow(A,e));
		}
	}
	else{
		viscosity2=4.5*pow(10,17);
	}
		
	#ifdef _DEBUG_
	_printf_("Viscosity2 %lf\n",mu2);
	#endif

	/*Return: */
	*pviscosity2=viscosity2;
}
