Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18058) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18059) @@ -976,5 +976,5 @@ /*return if floating*/ - if(element->IsFloating())return; + if(element->IsFloating()) return; /*Intermediaries*/ @@ -1068,5 +1068,87 @@ }/*}}}*/ void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector* gradient,int control_index){/*{{{*/ - _error_("not implemented yet"); + + /*Intermediaries*/ + int domaintype,dim; + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + dim = 2; + break; + case Domain2DverticalEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + dim = 1; + break; + case Domain3DEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + dim = 2; + break; + default: _error_("mesh "<GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* ge = xNew(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->GradientIndexing(&vertexpidlist[0],control_index); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* rheologyb_input = element->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + + thickness_input->GetInputValue(&thickness,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss); + adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss); + + basalelement->dViscositydBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctionsP1(basis,gauss); + + /*Build gradient vector (actually -dJ/dB): */ + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ void AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector* gradient,int control_index){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18058) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18059) @@ -213,4 +213,31 @@ return divergence; }/*}}}*/ +void Element::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dmudB; + IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */ + IssmDouble epsilon1d; /* epsilon=[exx]; */ + IssmDouble eps_eff; + + if(dim==2){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ + this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); + } + else{ + /* eps_eff^2 = 1/2 exx^2*/ + this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); + eps_eff = sqrt(epsilon1d*epsilon1d/2.); + } + + /*Get viscosity*/ + material->GetViscosity_B(&dmudB,eps_eff); + + /*Assign output pointer*/ + *pdmudB=dmudB; + +} +/*}}}*/ void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/ matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18058) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18059) @@ -60,4 +60,5 @@ void DeepEcho(); void DeleteMaterials(void); + void dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); IssmDouble Divergence(void); void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure); Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 18058) +++ /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 18059) @@ -26,4 +26,5 @@ virtual void Configure(Elements* elements)=0; virtual void GetViscosity(IssmDouble* pviscosity,IssmDouble epseff)=0; + virtual void GetViscosity_B(IssmDouble* pviscosity,IssmDouble epseff)=0; virtual void GetViscosityBar(IssmDouble* pviscosity,IssmDouble epseff)=0; virtual void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18058) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18059) @@ -284,4 +284,54 @@ } /*}}}*/ +/*FUNCTION Matice::GetViscosity_B {{{*/ +void Matice::GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){ + /*From a string tensor and a material object, return viscosity, using Glen's flow law. + (1-D) + 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 D=0.,n; + + /*Get B and n*/ + n=GetN(); _assert_(n>0.); + if(this->isdamaged){ + D=GetD(); + _assert_(D>=0. && D<1.); + } + + if (n==1.){ + /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */ + viscosity=(1.-D)/2.; + } + else{ + + /*if no strain rate, return maximum viscosity*/ + if(eps_eff==0.){ + viscosity = 1.e+14/2.; + } + + else{ + viscosity=(1.-D)/(2.*pow(eps_eff,(n-1.)/n)); + } + } + + /*Checks in debugging mode*/ + if(viscosity<=0) _error_("Negative viscosity"); + + /*Return: */ + *pviscosity=viscosity; +} +/*}}}*/ /*FUNCTION Matice::GetViscosityBar {{{*/ void Matice::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){ Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 18058) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 18059) @@ -55,4 +55,5 @@ 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 GetViscosityBar(IssmDouble* pviscosity, IssmDouble eps_eff); void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 18058) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 18059) @@ -76,4 +76,5 @@ void Configure(Elements* elements); void GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");}; + void GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");}; void GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");}; void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};