Index: issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- issm/trunk-jpl/src/c/Makefile.am	(revision 20686)
+++ issm/trunk-jpl/src/c/Makefile.am	(revision 20687)
@@ -61,5 +61,5 @@
 					./classes/Materials/Materials.cpp\
 					./classes/Materials/Matice.cpp\
-					./classes/Materials/Matearl.cpp\
+					./classes/Materials/Matestar.cpp\
 					./classes/Materials/Matpar.cpp\
 					./classes/Constraints/Constraints.cpp\
@@ -128,5 +128,5 @@
 					./shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp\
 					./shared/Elements/DrainageFunctionWaterfraction.cpp\
-					./shared/Elements/EarlComponents.cpp\
+					./shared/Elements/EstarComponents.cpp\
 					./shared/String/DescriptorIndex.cpp\
 					./toolkits/issm/IssmToolkitUtils.cpp\
Index: issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 20686)
+++ issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 20687)
@@ -283,6 +283,6 @@
 			dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
 		}
-		EarlStrainrateQuantities(&epso,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]);
-		lambdas[iv]=EarlLambdaS(epso,epsprime);
+		EstarStrainrateQuantities(&epso,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]);
+		lambdas[iv]=EstarLambdaS(epso,epsprime);
 		vorticityx[iv]=dvz[1]-dvy[2];
 		vorticityy[iv]=dvx[2]-dvz[0];
Index: issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 20686)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 20687)
@@ -2261,5 +2261,5 @@
 			if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
 			break;
-		case MatearlEnum:
+		case MatestarEnum:
 			this->InputDepthAverageAtBase(MaterialsRheologyKoEnum,MaterialsRheologyKobarEnum);
 			this->InputDepthAverageAtBase(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum);
@@ -2279,5 +2279,5 @@
 			this->inputs->DeleteInput(DamageDbarEnum);
 			break;
-		case MatearlEnum:
+		case MatestarEnum:
 			break;
 		default:
Index: issm/trunk-jpl/src/c/classes/Materials/Matearl.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Materials/Matearl.cpp	(revision 20686)
+++ 	(revision )
@@ -1,441 +1,0 @@
-/*!\file Matearl.c
- * \brief: implementation of the Matearl object
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include "./Matearl.h"
-#include "./Materials.h"
-#include "../Inputs/Input.h"
-#include "../Inputs/Inputs.h"
-#include "../Inputs/TriaInput.h"
-#include "../Inputs/PentaInput.h"
-#include "../Inputs/ControlInput.h"
-#include "../Elements/Element.h"
-#include "../Elements/Tria.h"
-#include "../Elements/Penta.h"
-#include "../Params/Parameters.h"
-#include "../Vertex.h"
-#include "../Hook.h"
-#include "../Node.h"
-#include "../IoModel.h"
-#include "../../shared/shared.h"
-
-/*Matearl constructors and destructor*/
-Matearl::Matearl(){/*{{{*/
-	this->helement=NULL;
-	this->element=NULL;
-	return;
-}
-/*}}}*/
-Matearl::Matearl(int matearl_mid,int index, IoModel* iomodel){/*{{{*/
-
-	/*Intermediaries:*/
-	int    matearl_eid;
-
-	/*Initialize id*/
-	this->mid=matearl_mid;
-
-	/*Hooks: */
-	matearl_eid=index+1;
-	this->helement=new Hook(&matearl_eid,1);
-	this->element=NULL;
-
-	return;
-}
-/*}}}*/
-Matearl::~Matearl(){/*{{{*/
-	delete helement;
-	return;
-}
-/*}}}*/
-
-/*Object virtual functions definitions:*/
-Object*   Matearl::copy() {/*{{{*/
-
-	/*Output*/
-	Matearl* matearl=NULL;
-
-	/*Initialize output*/
-	matearl=new Matearl();
-
-	/*copy fields: */
-	matearl->mid=this->mid;
-	matearl->helement=(Hook*)this->helement->copy();
-	matearl->element =(Element*)this->helement->delivers();
-
-	return matearl;
-}
-/*}}}*/
-Material* Matearl::copy2(Element* element_in) {/*{{{*/
-
-	/*Output*/
-	Matearl* matearl=NULL;
-
-	/*Initialize output*/
-	matearl=new Matearl();
-
-	/*copy fields: */
-	matearl->mid=this->mid;
-	matearl->helement=(Hook*)this->helement->copy();
-	matearl->element =element_in;
-
-	return matearl;
-}
-/*}}}*/
-void      Matearl::DeepEcho(void){/*{{{*/
-
-	_printf_("Matearl:\n");
-	_printf_("   mid: " << mid << "\n");
-	_printf_("   element:\n");
-	helement->Echo();
-}		
-/*}}}*/
-void      Matearl::Echo(void){/*{{{*/
-
-	_printf_("Matearl:\n");
-	_printf_("   mid: " << mid << "\n");
-	_printf_("   element:\n");
-	helement->Echo();
-}
-/*}}}*/
-int       Matearl::Id(void){ return mid; }/*{{{*/
-/*}}}*/
-int       Matearl::ObjectEnum(void){/*{{{*/
-
-	return MatearlEnum;
-
-}
-/*}}}*/
-void      Matearl::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
-
-	if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook(); 
-	
-	MARSHALLING_ENUM(MatearlEnum);
-	MARSHALLING(mid);
-	this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
-	this->element=(Element*)this->helement->delivers();
-
-}
-/*}}}*/
-
-/*Matearl management*/
-void  Matearl::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((DataSet*)elementsin);
-	this->element  = (Element*)helement->delivers();
-}
-/*}}}*/
-void  Matearl::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
-
-}
-/*}}}*/
-IssmDouble Matearl::GetA(){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-IssmDouble Matearl::GetAbar(){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-IssmDouble Matearl::GetB(){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-IssmDouble Matearl::GetBbar(){/*{{{*/
-
-	_error_("not implemented yet");
-}
-/*}}}*/
-IssmDouble Matearl::GetN(){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-IssmDouble Matearl::GetD(){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-IssmDouble Matearl::GetDbar(){/*{{{*/
-
-	_error_("not implemented yet");
-}
-/*}}}*/
-bool Matearl::IsDamage(){/*{{{*/
-
-	_error_("not implemented yet");
-}
-/*}}}*/
-void  Matearl::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
-	_error_("not implemented yet");
-	/*From a string tensor and a material object, return viscosity, using Glen's flow law.
-								(1-D) B
-	  viscosity= -------------------------
-						  2 eps_eff ^[(n-1)/n]
-
-	  where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
-	  and n the flow law exponent.
-
-	  If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we 
-	  return 10^14, initial viscosity.
-	  */
-
-	/*output: */
-	IssmDouble viscosity;
-
-	/*Intermediary: */
-	IssmDouble B,D=0.,n;
-
-	/*Get B and n*/
-	B=GetB(); _assert_(B>0.);
-	n=GetN(); _assert_(n>0.);
-
-	if (n==1.){
-		/*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
-		viscosity=(1.-D)*B/2.;
-	}
-	else{
-
-		/*if no strain rate, return maximum viscosity*/
-		if(eps_eff==0.){
-			viscosity = 1.e+14/2.;
-			//viscosity = B;
-			//viscosity=2.5*pow(10.,17);
-		}
-
-		else{
-			viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
-		}
-	}
-
-	/*Checks in debugging mode*/
-	if(viscosity<=0) _error_("Negative viscosity");
-
-	/*Return: */
-	*pviscosity=viscosity;
-}
-/*}}}*/
-void  Matearl::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-void  Matearl::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){/*{{{*/
-	 _error_("not implemented yet");
-}
-/*}}}*/
-void  Matearl::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-void  Matearl::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-void  Matearl::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-void  Matearl::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-void  Matearl::GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-void  Matearl::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
-
-}
-/*}}}*/
-void  Matearl::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
-
-}
-/*}}}*/
-void  Matearl::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols,int name, int type){/*{{{*/
-	/*Nothing updated yet*/
-}
-/*}}}*/
-void  Matearl::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/
-	/*Nothing updated yet*/
-}
-/*}}}*/
-void  Matearl::InputUpdateFromConstant(int constant, int name){/*{{{*/
-	/*Nothing updated yet*/
-}
-/*}}}*/
-void  Matearl::InputUpdateFromConstant(bool constant, int name){/*{{{*/
-	/*Nothing updated yet*/
-}
-/*}}}*/
-void  Matearl::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble vx,vy,vz;
-	IssmDouble dvx[3],dvy[3],dvz[3];
-	IssmDouble ko,Ec,Es;
-
-	/*Get velocity derivatives in all directions*/
-	_assert_(dim>1);
-	_assert_(vx_input);
-	vx_input->GetInputValue(&vx,gauss);
-	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-	_assert_(vy_input);
-	vy_input->GetInputValue(&vy,gauss);
-	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-	if(dim==3){
-		_assert_(vz_input);
-		vz_input->GetInputValue(&vz,gauss);
-		vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
-	}
-	else{
-		vz = 0.;
-		dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
-	}
-
-	/*Get enhancement factors and ko*/
-	Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
-	Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
-	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
-	ec_input->GetInputValue(&Ec,gauss);
-	es_input->GetInputValue(&Es,gauss);
-	ko_input->GetInputValue(&ko,gauss);
-
-	/*Compute viscosity*/
-	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
-}
-/*}}}*/
-void  Matearl::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-	this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
-}/*}}}*/
-void  Matearl::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble vx,vy,vz;
-	IssmDouble dvx[3],dvy[3],dvz[3];
-	IssmDouble ko,Ec,Es;
-
-	/*Get velocity derivatives in all directions*/
-	_assert_(dim==2 || dim==3);
-	_assert_(vx_input);
-	vx_input->GetInputValue(&vx,gauss);
-	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-	if(dim==3){
-		_assert_(vy_input);
-		vy_input->GetInputValue(&vy,gauss);
-		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-	}
-	else{
-		dvx[2] = 0.;
-		vy = 0.;
-		dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
-	}
-	vz = 0.;
-	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
-
-	/*Get enhancement factors and ko*/
-	Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
-	Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
-	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
-	ec_input->GetInputValue(&Ec,gauss);
-	es_input->GetInputValue(&Es,gauss);
-	ko_input->GetInputValue(&ko,gauss);
-
-	/*Compute viscosity*/
-	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
-}/*}}}*/
-void  Matearl::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
-void  Matearl::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
-void  Matearl::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
-	/*Intermediaries*/
-	IssmDouble vx,vy,vz;
-	IssmDouble dvx[3],dvy[3],dvz[3];
-	IssmDouble ko,Ec,Es;
-
-	/*Get velocity derivatives in all directions*/
-	_assert_(dim==1 || dim==2);
-	_assert_(vx_input);
-	vx_input->GetInputValue(&vx,gauss);
-	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-	if(dim==2){
-		_assert_(vy_input);
-		vy_input->GetInputValue(&vy,gauss);
-		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-	}
-	else{
-		dvx[1] = 0.;
-		dvx[2] = 0.;
-		vy = 0.;
-		dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
-	}
-	vz = 0.;
-	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
-
-	/*Get enhancement factors and ko*/
-	Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcbarEnum); _assert_(ec_input);
-	Input* es_input = element->inputs->GetInput(MaterialsRheologyEsbarEnum); _assert_(es_input);
-	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input);
-	ec_input->GetInputValue(&Ec,gauss);
-	es_input->GetInputValue(&Es,gauss);
-	ko_input->GetInputValue(&ko,gauss);
-
-	/*Compute viscosity*/
-	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
-}/*}}}*/
-void  Matearl::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
-void  Matearl::ResetHooks(){/*{{{*/
-
-	this->element=NULL;
-
-	/*Get Element type*/
-	this->helement->reset();
-
-}
-/*}}}*/
-IssmDouble Matearl::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble viscosity;
-	IssmDouble E,lambdas;
-	IssmDouble epso,epsprime_norm;
-	IssmDouble vmag,dvmag[3];
-
-	/*Calculate velocity magnitude and its derivative*/
-	vmag = sqrt(vx*vx+vy*vy+vz*vz);
-	if(vmag<1e-12){
-		dvmag[0]=0;
-		dvmag[1]=0;
-		dvmag[2]=0;
-	}
-	else{
-		dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
-		dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
-		dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
-	}
-
-	EarlStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
-	lambdas=EarlLambdaS(epso,epsprime_norm);
-
-	/*Get total enhancement factor E(lambdas)*/
-	E = Ec + (Es-Ec)*lambdas*lambdas;
-
-	/*Compute viscosity*/
-	_assert_(E>0.); 
-	_assert_(ko>0.); 
-	_assert_(epso>0.); 
-	viscosity = 1./(2*pow(ko*E*epso*epso,1./3.));
-
-	/*Assign output pointer*/
-	return viscosity;
-}
-/*}}}*/
Index: issm/trunk-jpl/src/c/classes/Materials/Matearl.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Materials/Matearl.h	(revision 20686)
+++ 	(revision )
@@ -1,87 +1,0 @@
-/*!\file Matearl.h
- * \brief: header file for matice object
- */
-
-#ifndef MATEARL_H_
-#define MATEARL_H_
-
-/*Headers:*/
-/*{{{*/
-#include "./Material.h"
-#include "../Hook.h"
-class IoModel;
-class Elements;
-class Element;
-class Loads;
-class Nodes;
-class Vertices;
-class Materials;
-class Parameters;
-class Gauss;
-class Input;
-/*}}}*/
-
-class Matearl: public Material{
-
-	private: 
-		int      mid;
-		Hook    *helement;
-		Element *element;
-
-	public:
-		/*Matearl constructors, destructors: {{{*/
-		Matearl();
-		Matearl(int mid,int i, IoModel* iomodel);
-		~Matearl();
-		/*}}}*/
-		/*Object virtual functions definitions:{{{ */
-		void  Echo();
-		void  DeepEcho();
-		int   Id(); 
-		int   ObjectEnum();
-		Object* copy();
-		void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
-		/*}}}*/
-		/*Update virtual funictions definitions: {{{*/
-		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
-		void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols, int name, int type);
-		void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
-		void  InputUpdateFromConstant(IssmDouble constant, int name);
-		void  InputUpdateFromConstant(int constant, int name);
-		void  InputUpdateFromConstant(bool constant, int name);
-		void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented");};
-		/*}}}*/
-		/*Material virtual functions resolution: {{{*/
-		void   Configure(Elements* elements);
-		Material*  copy2(Element* element);
-		void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
-		void       GetViscosity(IssmDouble* pviscosity, IssmDouble eps_eff);
-		void       GetViscosity_B(IssmDouble* pviscosity, IssmDouble eps_eff);
-		void       GetViscosity_D(IssmDouble* pviscosity, IssmDouble eps_eff);
-		void       GetViscosityBar(IssmDouble* pviscosity, IssmDouble eps_eff);
-		void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
-		void       GetViscosityDComplement(IssmDouble*, IssmDouble*);
-		void       GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
-		void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
-		IssmDouble GetA();
-		IssmDouble GetAbar();
-		IssmDouble GetB();
-		IssmDouble GetBbar();
-		IssmDouble GetD();
-		IssmDouble GetDbar();
-		IssmDouble GetN();
-		bool       IsDamage();
-		void       ResetHooks();
-
-		void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
-		void       ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
-		void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
-		void       ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
-		void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
-		void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
-		void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
-		/*}}}*/
-		IssmDouble GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz);
-};
-
-#endif  /* _MATEARL_H_ */
Index: issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp	(revision 20687)
+++ issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp	(revision 20687)
@@ -0,0 +1,441 @@
+/*!\file Matestar.c
+ * \brief: implementation of the Matestar object
+ */
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./Matestar.h"
+#include "./Materials.h"
+#include "../Inputs/Input.h"
+#include "../Inputs/Inputs.h"
+#include "../Inputs/TriaInput.h"
+#include "../Inputs/PentaInput.h"
+#include "../Inputs/ControlInput.h"
+#include "../Elements/Element.h"
+#include "../Elements/Tria.h"
+#include "../Elements/Penta.h"
+#include "../Params/Parameters.h"
+#include "../Vertex.h"
+#include "../Hook.h"
+#include "../Node.h"
+#include "../IoModel.h"
+#include "../../shared/shared.h"
+
+/*Matestar constructors and destructor*/
+Matestar::Matestar(){/*{{{*/
+	this->helement=NULL;
+	this->element=NULL;
+	return;
+}
+/*}}}*/
+Matestar::Matestar(int matestar_mid,int index, IoModel* iomodel){/*{{{*/
+
+	/*Intermediaries:*/
+	int    matestar_eid;
+
+	/*Initialize id*/
+	this->mid=matestar_mid;
+
+	/*Hooks: */
+	matestar_eid=index+1;
+	this->helement=new Hook(&matestar_eid,1);
+	this->element=NULL;
+
+	return;
+}
+/*}}}*/
+Matestar::~Matestar(){/*{{{*/
+	delete helement;
+	return;
+}
+/*}}}*/
+
+/*Object virtual functions definitions:*/
+Object*   Matestar::copy() {/*{{{*/
+
+	/*Output*/
+	Matestar* matestar=NULL;
+
+	/*Initialize output*/
+	matestar=new Matestar();
+
+	/*copy fields: */
+	matestar->mid=this->mid;
+	matestar->helement=(Hook*)this->helement->copy();
+	matestar->element =(Element*)this->helement->delivers();
+
+	return matestar;
+}
+/*}}}*/
+Material* Matestar::copy2(Element* element_in) {/*{{{*/
+
+	/*Output*/
+	Matestar* matestar=NULL;
+
+	/*Initialize output*/
+	matestar=new Matestar();
+
+	/*copy fields: */
+	matestar->mid=this->mid;
+	matestar->helement=(Hook*)this->helement->copy();
+	matestar->element =element_in;
+
+	return matestar;
+}
+/*}}}*/
+void      Matestar::DeepEcho(void){/*{{{*/
+
+	_printf_("Matestar:\n");
+	_printf_("   mid: " << mid << "\n");
+	_printf_("   element:\n");
+	helement->Echo();
+}		
+/*}}}*/
+void      Matestar::Echo(void){/*{{{*/
+
+	_printf_("Matestar:\n");
+	_printf_("   mid: " << mid << "\n");
+	_printf_("   element:\n");
+	helement->Echo();
+}
+/*}}}*/
+int       Matestar::Id(void){ return mid; }/*{{{*/
+/*}}}*/
+int       Matestar::ObjectEnum(void){/*{{{*/
+
+	return MatestarEnum;
+
+}
+/*}}}*/
+void      Matestar::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
+
+	if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook(); 
+	
+	MARSHALLING_ENUM(MatestarEnum);
+	MARSHALLING(mid);
+	this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
+	this->element=(Element*)this->helement->delivers();
+
+}
+/*}}}*/
+
+/*Matestar management*/
+void  Matestar::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((DataSet*)elementsin);
+	this->element  = (Element*)helement->delivers();
+}
+/*}}}*/
+void  Matestar::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
+
+}
+/*}}}*/
+IssmDouble Matestar::GetA(){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+IssmDouble Matestar::GetAbar(){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+IssmDouble Matestar::GetB(){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+IssmDouble Matestar::GetBbar(){/*{{{*/
+
+	_error_("not implemented yet");
+}
+/*}}}*/
+IssmDouble Matestar::GetN(){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+IssmDouble Matestar::GetD(){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+IssmDouble Matestar::GetDbar(){/*{{{*/
+
+	_error_("not implemented yet");
+}
+/*}}}*/
+bool Matestar::IsDamage(){/*{{{*/
+
+	_error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
+	_error_("not implemented yet");
+	/*From a string tensor and a material object, return viscosity, using Glen's flow law.
+								(1-D) B
+	  viscosity= -------------------------
+						  2 eps_eff ^[(n-1)/n]
+
+	  where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
+	  and n the flow law exponent.
+
+	  If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we 
+	  return 10^14, initial viscosity.
+	  */
+
+	/*output: */
+	IssmDouble viscosity;
+
+	/*Intermediary: */
+	IssmDouble B,D=0.,n;
+
+	/*Get B and n*/
+	B=GetB(); _assert_(B>0.);
+	n=GetN(); _assert_(n>0.);
+
+	if (n==1.){
+		/*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
+		viscosity=(1.-D)*B/2.;
+	}
+	else{
+
+		/*if no strain rate, return maximum viscosity*/
+		if(eps_eff==0.){
+			viscosity = 1.e+14/2.;
+			//viscosity = B;
+			//viscosity=2.5*pow(10.,17);
+		}
+
+		else{
+			viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
+		}
+	}
+
+	/*Checks in debugging mode*/
+	if(viscosity<=0) _error_("Negative viscosity");
+
+	/*Return: */
+	*pviscosity=viscosity;
+}
+/*}}}*/
+void  Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){/*{{{*/
+	 _error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
+
+}
+/*}}}*/
+void  Matestar::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
+
+}
+/*}}}*/
+void  Matestar::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols,int name, int type){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
+void  Matestar::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
+void  Matestar::InputUpdateFromConstant(int constant, int name){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
+void  Matestar::InputUpdateFromConstant(bool constant, int name){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
+void  Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble vx,vy,vz;
+	IssmDouble dvx[3],dvy[3],dvz[3];
+	IssmDouble ko,Ec,Es;
+
+	/*Get velocity derivatives in all directions*/
+	_assert_(dim>1);
+	_assert_(vx_input);
+	vx_input->GetInputValue(&vx,gauss);
+	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+	_assert_(vy_input);
+	vy_input->GetInputValue(&vy,gauss);
+	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+	if(dim==3){
+		_assert_(vz_input);
+		vz_input->GetInputValue(&vz,gauss);
+		vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+	}
+	else{
+		vz = 0.;
+		dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
+	}
+
+	/*Get enhancement factors and ko*/
+	Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
+	Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
+	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
+	ec_input->GetInputValue(&Ec,gauss);
+	es_input->GetInputValue(&Es,gauss);
+	ko_input->GetInputValue(&ko,gauss);
+
+	/*Compute viscosity*/
+	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
+}
+/*}}}*/
+void  Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
+	this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
+}/*}}}*/
+void  Matestar::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble vx,vy,vz;
+	IssmDouble dvx[3],dvy[3],dvz[3];
+	IssmDouble ko,Ec,Es;
+
+	/*Get velocity derivatives in all directions*/
+	_assert_(dim==2 || dim==3);
+	_assert_(vx_input);
+	vx_input->GetInputValue(&vx,gauss);
+	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+	if(dim==3){
+		_assert_(vy_input);
+		vy_input->GetInputValue(&vy,gauss);
+		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+	}
+	else{
+		dvx[2] = 0.;
+		vy = 0.;
+		dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
+	}
+	vz = 0.;
+	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
+
+	/*Get enhancement factors and ko*/
+	Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcEnum); _assert_(ec_input);
+	Input* es_input = element->inputs->GetInput(MaterialsRheologyEsEnum); _assert_(es_input);
+	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKoEnum); _assert_(ko_input);
+	ec_input->GetInputValue(&Ec,gauss);
+	es_input->GetInputValue(&Es,gauss);
+	ko_input->GetInputValue(&ko,gauss);
+
+	/*Compute viscosity*/
+	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
+}/*}}}*/
+void  Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+void  Matestar::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+void  Matestar::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+	/*Intermediaries*/
+	IssmDouble vx,vy,vz;
+	IssmDouble dvx[3],dvy[3],dvz[3];
+	IssmDouble ko,Ec,Es;
+
+	/*Get velocity derivatives in all directions*/
+	_assert_(dim==1 || dim==2);
+	_assert_(vx_input);
+	vx_input->GetInputValue(&vx,gauss);
+	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+	if(dim==2){
+		_assert_(vy_input);
+		vy_input->GetInputValue(&vy,gauss);
+		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+	}
+	else{
+		dvx[1] = 0.;
+		dvx[2] = 0.;
+		vy = 0.;
+		dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
+	}
+	vz = 0.;
+	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
+
+	/*Get enhancement factors and ko*/
+	Input* ec_input = element->inputs->GetInput(MaterialsRheologyEcbarEnum); _assert_(ec_input);
+	Input* es_input = element->inputs->GetInput(MaterialsRheologyEsbarEnum); _assert_(es_input);
+	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input);
+	ec_input->GetInputValue(&Ec,gauss);
+	es_input->GetInputValue(&Es,gauss);
+	ko_input->GetInputValue(&ko,gauss);
+
+	/*Compute viscosity*/
+	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
+}/*}}}*/
+void  Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+void  Matestar::ResetHooks(){/*{{{*/
+
+	this->element=NULL;
+
+	/*Get Element type*/
+	this->helement->reset();
+
+}
+/*}}}*/
+IssmDouble Matestar::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble viscosity;
+	IssmDouble E,lambdas;
+	IssmDouble epso,epsprime_norm;
+	IssmDouble vmag,dvmag[3];
+
+	/*Calculate velocity magnitude and its derivative*/
+	vmag = sqrt(vx*vx+vy*vy+vz*vz);
+	if(vmag<1e-12){
+		dvmag[0]=0;
+		dvmag[1]=0;
+		dvmag[2]=0;
+	}
+	else{
+		dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
+		dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
+		dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
+	}
+
+	EstarStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
+	lambdas=EstarLambdaS(epso,epsprime_norm);
+
+	/*Get total enhancement factor E(lambdas)*/
+	E = Ec + (Es-Ec)*lambdas*lambdas;
+
+	/*Compute viscosity*/
+	_assert_(E>0.); 
+	_assert_(ko>0.); 
+	_assert_(epso>0.); 
+	viscosity = 1./(2*pow(ko*E*epso*epso,1./3.));
+
+	/*Assign output pointer*/
+	return viscosity;
+}
+/*}}}*/
Index: issm/trunk-jpl/src/c/classes/Materials/Matestar.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Materials/Matestar.h	(revision 20687)
+++ issm/trunk-jpl/src/c/classes/Materials/Matestar.h	(revision 20687)
@@ -0,0 +1,87 @@
+/*!\file Matestar.h
+ * \brief: header file for matice object
+ */
+
+#ifndef MATEARL_H_
+#define MATEARL_H_
+
+/*Headers:*/
+/*{{{*/
+#include "./Material.h"
+#include "../Hook.h"
+class IoModel;
+class Elements;
+class Element;
+class Loads;
+class Nodes;
+class Vertices;
+class Materials;
+class Parameters;
+class Gauss;
+class Input;
+/*}}}*/
+
+class Matestar: public Material{
+
+	private: 
+		int      mid;
+		Hook    *helement;
+		Element *element;
+
+	public:
+		/*Matestar constructors, destructors: {{{*/
+		Matestar();
+		Matestar(int mid,int i, IoModel* iomodel);
+		~Matestar();
+		/*}}}*/
+		/*Object virtual functions definitions:{{{ */
+		void  Echo();
+		void  DeepEcho();
+		int   Id(); 
+		int   ObjectEnum();
+		Object* copy();
+		void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
+		/*}}}*/
+		/*Update virtual funictions definitions: {{{*/
+		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
+		void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols, int name, int type);
+		void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
+		void  InputUpdateFromConstant(IssmDouble constant, int name);
+		void  InputUpdateFromConstant(int constant, int name);
+		void  InputUpdateFromConstant(bool constant, int name);
+		void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented");};
+		/*}}}*/
+		/*Material virtual functions resolution: {{{*/
+		void   Configure(Elements* elements);
+		Material*  copy2(Element* element);
+		void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
+		void       GetViscosity(IssmDouble* pviscosity, IssmDouble eps_eff);
+		void       GetViscosity_B(IssmDouble* pviscosity, IssmDouble eps_eff);
+		void       GetViscosity_D(IssmDouble* pviscosity, IssmDouble eps_eff);
+		void       GetViscosityBar(IssmDouble* pviscosity, IssmDouble eps_eff);
+		void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
+		void       GetViscosityDComplement(IssmDouble*, IssmDouble*);
+		void       GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
+		void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
+		IssmDouble GetA();
+		IssmDouble GetAbar();
+		IssmDouble GetB();
+		IssmDouble GetBbar();
+		IssmDouble GetD();
+		IssmDouble GetDbar();
+		IssmDouble GetN();
+		bool       IsDamage();
+		void       ResetHooks();
+
+		void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
+		void       ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
+		void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
+		void       ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
+		void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
+		void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
+		void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
+		/*}}}*/
+		IssmDouble GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz);
+};
+
+#endif  /* _MATEARL_H_ */
Index: issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 20686)
+++ issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 20687)
@@ -82,5 +82,5 @@
 		case MaticeEnum:
 		case MatdamageiceEnum:
-		case MatearlEnum:
+		case MatestarEnum:
 			iomodel->FindConstant(&this->rho_ice,MaterialsRhoIceEnum);
 			iomodel->FindConstant(&this->rho_water,MaterialsRhoSeawaterEnum);
Index: issm/trunk-jpl/src/c/classes/classes.h
===================================================================
--- issm/trunk-jpl/src/c/classes/classes.h	(revision 20686)
+++ issm/trunk-jpl/src/c/classes/classes.h	(revision 20687)
@@ -83,5 +83,5 @@
 #include "./Materials/Material.h"
 #include "./Materials/Matice.h"
-#include "./Materials/Matearl.h"
+#include "./Materials/Matestar.h"
 #include "./Materials/Matpar.h"
 
Index: issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 20686)
+++ issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp	(revision 20687)
@@ -79,9 +79,9 @@
 			}
 			break;
-		case MatearlEnum:
+		case MatestarEnum:
 			iomodel->FetchDataToInput(elements,MaterialsRheologyKoEnum);
 			iomodel->FetchDataToInput(elements,MaterialsRheologyEcEnum);
 			iomodel->FetchDataToInput(elements,MaterialsRheologyEsEnum);
-			for(i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matearl(i+1,i,iomodel));
+			for(i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matestar(i+1,i,iomodel));
 			switch(iomodel->domaindim){
 				case 2:
Index: issm/trunk-jpl/src/c/shared/Elements/EarlComponents.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Elements/EarlComponents.cpp	(revision 20686)
+++ 	(revision )
@@ -1,123 +1,0 @@
-#include <math.h>
-#include "../Numerics/types.h"
-#include "../Exceptions/exceptions.h"
-#include "./elements.h"
-void EarlStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble omega[3],omega_norm;
-	IssmDouble nrsp[3],nrsp_norm;
-	IssmDouble eps[3][3],epso;
-	IssmDouble epsprime[3],epsprime_norm;
-
-	/*Get omega, correction for rigid body rotation*/
-	EarlOmega(&omega[0],vx,vy,vz,vmag,dvx,dvy,dvz,dvmag);
-
-	/*Non-rotating shear plane*/
-	nrsp[0] =  vy*omega[2] - vz*omega[1];
-	nrsp[1] =  vz*omega[0] - vx*omega[2];
-	nrsp[2] =  vx*omega[1] - vy*omega[0];
-
-	/*Normalize*/
-	nrsp_norm = sqrt(nrsp[0]*nrsp[0] + nrsp[1]*nrsp[1] + nrsp[2]*nrsp[2]);
-	if(nrsp_norm==0){
-		nrsp[0] = 0.;
-		nrsp[1] = 0.;
-		nrsp[2] = 1.;
-	}
-	else{
-		nrsp[0] =nrsp[0]/nrsp_norm;
-		nrsp[1] =nrsp[1]/nrsp_norm;
-		nrsp[2] =nrsp[2]/nrsp_norm;
-	}
-
-	/*Build strain rate tensor*/
-	eps[0][0] = dvx[0];             eps[0][1] = .5*(dvx[1]+dvy[0]); eps[0][2] = .5*(dvx[2]+dvz[0]);
-	eps[1][0] = .5*(dvx[1]+dvy[0]); eps[1][1] = dvy[1];             eps[1][2] = .5*(dvy[2]+dvz[1]);
-	eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2];
-
-	/*octahedral strain rate*/
-	epso = 0.;
-	for(int i=0;i<3;i++) for(int j=0;j<3;j++) epso += eps[i][j]*eps[i][j];
-	if(epso==0.) epso=1.e-14;
-	epso=sqrt(epso)/sqrt(3.);
-
-	/*Compute the shear strain rate on the non ratating shear plane*/
-	epsprime[0]=0.;
-	epsprime[1]=0.;
-	epsprime[2]=0.;
-	/*term #1: eps'.n */
-	for(int i=0;i<3;i++){
-		for(int j=0;j<3;j++){
-			epsprime[i] += eps[i][j]*nrsp[j];
-		}
-	}
-	/*term #2: ((eps'.n).n)n */
-	for(int i=0;i<3;i++){
-		for(int j=0;j<3;j++){
-			for(int k=0;k<3;k++){
-				epsprime[j] += -nrsp[i]*eps[i][k]*nrsp[k]*nrsp[j];
-			}
-		}
-	}
-	/*term #3: ((eps'.n).omega)omega */
-	for(int i=0;i<3;i++){
-		for(int j=0;j<3;j++){
-			for(int k=0;k<3;k++){
-				epsprime[j] += -nrsp[i]*eps[i][k]*omega[k]*omega[j];
-			}
-		}
-	}
-
-	/*Get norm of epsprime*/
-	epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]);
-	
-	/*Assign output pointers*/
-	*pepso=epso;
-	*pepsprime_norm=epsprime_norm;
-}/*}}}*/
-void EarlOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble omega_norm;
-
-	/*Create vorticity vector*/
-	_assert_(dvx && dvy && dvz && dvmag);
-	if(vmag<1e-12)vmag=1e-12;
-
-	/*Create vorticity vector, corrected for rigid body rotation
-	 * \overline{\omega} =\omega - \Omega
-	 *                   =-\nabla\times{\bf v} - 2*\nabla U\times{\bf v}/U;
-	 * check the magnitude of the second term -- if it is small, then the two
-	 * vorticities (omega and first term in omega) are approx. equal
-	 *
-	 * */
-	omega[0] = -(dvz[1] - dvy[2]) + 2*(dvmag[1]*vz - dvmag[2]*vy)/vmag;
-	omega[1] = -(dvx[2] - dvz[0]) + 2*(dvmag[2]*vx - dvmag[0]*vz)/vmag;
-	omega[2] = -(dvy[0] - dvx[1]) + 2*(dvmag[0]*vy - dvmag[1]*vx)/vmag;
-
-	/*Take out vorticity component aligned with the velocity*/
-	IssmDouble wdotv = vx/vmag*omega[0] + vy/vmag*omega[1] + vz/vmag*omega[2];
-	omega[0] = omega[0] - wdotv*vx/vmag;
-	omega[1] = omega[1] - wdotv*vy/vmag;
-	omega[2] = omega[2] - wdotv*vz/vmag;
-
-	/*Normalize*/
-	omega_norm = sqrt(omega[0]*omega[0] + omega[1]*omega[1] + omega[2]*omega[2]);
-	if(omega_norm==0){
-		omega[0] = 0.;
-		omega[1] = 0.;
-		omega[2] = 1.;
-	}
-	else{
-		omega[0] =omega[0]/omega_norm;
-		omega[1] =omega[1]/omega_norm;
-		omega[2] =omega[2]/omega_norm;
-	}
-
-}/*}}}*/
-IssmDouble EarlLambdaS(IssmDouble epso, IssmDouble epsprime_norm){/*{{{*/
-   _assert_(epso>0.);
-	_assert_(epsprime_norm>=0.); 
-	return sqrt(2*epsprime_norm*epsprime_norm/(3*epso*epso));
-}/*}}}*/
Index: issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp	(revision 20687)
+++ issm/trunk-jpl/src/c/shared/Elements/EstarComponents.cpp	(revision 20687)
@@ -0,0 +1,123 @@
+#include <math.h>
+#include "../Numerics/types.h"
+#include "../Exceptions/exceptions.h"
+#include "./elements.h"
+void EstarStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble omega[3],omega_norm;
+	IssmDouble nrsp[3],nrsp_norm;
+	IssmDouble eps[3][3],epso;
+	IssmDouble epsprime[3],epsprime_norm;
+
+	/*Get omega, correction for rigid body rotation*/
+	EstarOmega(&omega[0],vx,vy,vz,vmag,dvx,dvy,dvz,dvmag);
+
+	/*Non-rotating shear plane*/
+	nrsp[0] =  vy*omega[2] - vz*omega[1];
+	nrsp[1] =  vz*omega[0] - vx*omega[2];
+	nrsp[2] =  vx*omega[1] - vy*omega[0];
+
+	/*Normalize*/
+	nrsp_norm = sqrt(nrsp[0]*nrsp[0] + nrsp[1]*nrsp[1] + nrsp[2]*nrsp[2]);
+	if(nrsp_norm==0){
+		nrsp[0] = 0.;
+		nrsp[1] = 0.;
+		nrsp[2] = 1.;
+	}
+	else{
+		nrsp[0] =nrsp[0]/nrsp_norm;
+		nrsp[1] =nrsp[1]/nrsp_norm;
+		nrsp[2] =nrsp[2]/nrsp_norm;
+	}
+
+	/*Build strain rate tensor*/
+	eps[0][0] = dvx[0];             eps[0][1] = .5*(dvx[1]+dvy[0]); eps[0][2] = .5*(dvx[2]+dvz[0]);
+	eps[1][0] = .5*(dvx[1]+dvy[0]); eps[1][1] = dvy[1];             eps[1][2] = .5*(dvy[2]+dvz[1]);
+	eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2];
+
+	/*octahedral strain rate*/
+	epso = 0.;
+	for(int i=0;i<3;i++) for(int j=0;j<3;j++) epso += eps[i][j]*eps[i][j];
+	if(epso==0.) epso=1.e-14;
+	epso=sqrt(epso)/sqrt(3.);
+
+	/*Compute the shear strain rate on the non ratating shear plane*/
+	epsprime[0]=0.;
+	epsprime[1]=0.;
+	epsprime[2]=0.;
+	/*term #1: eps'.n */
+	for(int i=0;i<3;i++){
+		for(int j=0;j<3;j++){
+			epsprime[i] += eps[i][j]*nrsp[j];
+		}
+	}
+	/*term #2: ((eps'.n).n)n */
+	for(int i=0;i<3;i++){
+		for(int j=0;j<3;j++){
+			for(int k=0;k<3;k++){
+				epsprime[j] += -nrsp[i]*eps[i][k]*nrsp[k]*nrsp[j];
+			}
+		}
+	}
+	/*term #3: ((eps'.n).omega)omega */
+	for(int i=0;i<3;i++){
+		for(int j=0;j<3;j++){
+			for(int k=0;k<3;k++){
+				epsprime[j] += -nrsp[i]*eps[i][k]*omega[k]*omega[j];
+			}
+		}
+	}
+
+	/*Get norm of epsprime*/
+	epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]);
+	
+	/*Assign output pointers*/
+	*pepso=epso;
+	*pepsprime_norm=epsprime_norm;
+}/*}}}*/
+void EstarOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble omega_norm;
+
+	/*Create vorticity vector*/
+	_assert_(dvx && dvy && dvz && dvmag);
+	if(vmag<1e-12)vmag=1e-12;
+
+	/*Create vorticity vector, corrected for rigid body rotation
+	 * \overline{\omega} =\omega - \Omega
+	 *                   =-\nabla\times{\bf v} - 2*\nabla U\times{\bf v}/U;
+	 * check the magnitude of the second term -- if it is small, then the two
+	 * vorticities (omega and first term in omega) are approx. equal
+	 *
+	 * */
+	omega[0] = -(dvz[1] - dvy[2]) + 2*(dvmag[1]*vz - dvmag[2]*vy)/vmag;
+	omega[1] = -(dvx[2] - dvz[0]) + 2*(dvmag[2]*vx - dvmag[0]*vz)/vmag;
+	omega[2] = -(dvy[0] - dvx[1]) + 2*(dvmag[0]*vy - dvmag[1]*vx)/vmag;
+
+	/*Take out vorticity component aligned with the velocity*/
+	IssmDouble wdotv = vx/vmag*omega[0] + vy/vmag*omega[1] + vz/vmag*omega[2];
+	omega[0] = omega[0] - wdotv*vx/vmag;
+	omega[1] = omega[1] - wdotv*vy/vmag;
+	omega[2] = omega[2] - wdotv*vz/vmag;
+
+	/*Normalize*/
+	omega_norm = sqrt(omega[0]*omega[0] + omega[1]*omega[1] + omega[2]*omega[2]);
+	if(omega_norm==0){
+		omega[0] = 0.;
+		omega[1] = 0.;
+		omega[2] = 1.;
+	}
+	else{
+		omega[0] =omega[0]/omega_norm;
+		omega[1] =omega[1]/omega_norm;
+		omega[2] =omega[2]/omega_norm;
+	}
+
+}/*}}}*/
+IssmDouble EstarLambdaS(IssmDouble epso, IssmDouble epsprime_norm){/*{{{*/
+   _assert_(epso>0.);
+	_assert_(epsprime_norm>=0.); 
+	return sqrt(2*epsprime_norm*epsprime_norm/(3*epso*epso));
+}/*}}}*/
Index: issm/trunk-jpl/src/c/shared/Elements/elements.h
===================================================================
--- issm/trunk-jpl/src/c/shared/Elements/elements.h	(revision 20686)
+++ issm/trunk-jpl/src/c/shared/Elements/elements.h	(revision 20687)
@@ -15,7 +15,7 @@
 IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat);
 // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n);
-IssmDouble EarlLambdaS(IssmDouble epso, IssmDouble epsprime_norm);
-void EarlOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag);
-void EarlStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag);
+IssmDouble EstarLambdaS(IssmDouble epso, IssmDouble epsprime_norm);
+void EstarOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag);
+void EstarStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag);
 IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec,
 				 IssmDouble* pdds, IssmDouble* pds, IssmDouble* melt, IssmDouble* accu, IssmDouble signorm, 
Index: issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 20686)
+++ issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 20687)
@@ -611,5 +611,5 @@
 	MaticeEnum,
 	MatdamageiceEnum,
-	MatearlEnum,
+	MatestarEnum,
 	MatparEnum,
 	NodeEnum,
Index: issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 20686)
+++ issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 20687)
@@ -601,5 +601,5 @@
 		case MaticeEnum : return "Matice";
 		case MatdamageiceEnum : return "Matdamageice";
-		case MatearlEnum : return "Matearl";
+		case MatestarEnum : return "Matestar";
 		case MatparEnum : return "Matpar";
 		case NodeEnum : return "Node";
Index: issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 20686)
+++ issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 20687)
@@ -613,5 +613,5 @@
 	      else if (strcmp(name,"Matice")==0) return MaticeEnum;
 	      else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
-	      else if (strcmp(name,"Matearl")==0) return MatearlEnum;
+	      else if (strcmp(name,"Matestar")==0) return MatestarEnum;
 	      else if (strcmp(name,"Matpar")==0) return MatparEnum;
 	      else if (strcmp(name,"Node")==0) return NodeEnum;
