Index: ../trunk-jpl/src/c/classes/Materials/Matbafl.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Materials/Matbafl.cpp (revision 20650) +++ ../trunk-jpl/src/c/classes/Materials/Matbafl.cpp (revision 20651) @@ -384,6 +384,47 @@ 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 Matbafl::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ +void Matbafl::ResetHooks(){/*{{{*/ + + this->element=NULL; + + /*Get Element type*/ + this->helement->reset(); + +} +/*}}}*/ +IssmDouble Matbafl::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/ + + /*Intermediaries*/ + IssmDouble viscosity; + IssmDouble vorticity[3],vorticity_norm; + IssmDouble nrsp[3],nrsp_norm; + IssmDouble eps[3][3],epso; + IssmDouble epsprime[3],epsprime_norm; + IssmDouble E,lambdas; + + /*Create vorticity vector*/ + _assert_(dvx && dvy && dvz); + vorticity[0] = dvz[1] - dvy[2]; + vorticity[1] = dvx[2] - dvz[0]; + vorticity[2] = dvy[0] - dvx[1]; + + /*Normalize*/ + vorticity_norm = sqrt(vorticity[0]*vorticity[0] + vorticity[1]*vorticity[1] + vorticity[2]*vorticity[2]); + if(vorticity_norm==0){ + vorticity[0] = 0.; + vorticity[1] = 0.; + vorticity[2] = 1.; + } else{ vorticity[0] =vorticity[0]/vorticity_norm; vorticity[1] =vorticity[1]/vorticity_norm;