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: sm/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 -#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: sm/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 +#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;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matearl(i+1,i,iomodel)); + for(i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matestar(i+1,i,iomodel)); switch(iomodel->domaindim){ case 2: Index: sm/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 -#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 +#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;