Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21619) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21620) @@ -1561,4 +1561,8 @@ name==EsaNmotionEnum || name==EsaEmotionEnum || + name==EsaStrainratexxEnum || + name==EsaStrainratexyEnum || + name==EsaStrainrateyyEnum || + name==EsaRotationrateEnum || name==EsaDeltathicknessEnum || name==GiaWEnum || @@ -2121,4 +2125,8 @@ case DeviatoricStresszzEnum: case DeviatoricStresseffectiveEnum: this->ComputeDeviatoricStressTensor(); break; + case EsaStrainratexxEnum: + case EsaStrainratexyEnum: + case EsaStrainrateyyEnum: + case EsaRotationrateEnum: this->ComputeEsaStrainAndVorticity(); break; case SigmaNNEnum: this->ComputeSigmaNN(); break; case LambdaSEnum: this->ComputeLambdaS(); break; @@ -2673,4 +2681,24 @@ } /*}}}*/ +void Element::StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ + + /*Intermediaries*/ + IssmDouble dvx[3]; + IssmDouble dvy[3]; + + /*Check that both inputs have been found*/ + if(!vx_input || !vy_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + epsilon[0] = dvx[0]; // normal strain rate x-direction + epsilon[1] = dvy[1]; // normal strain rate y-direction + epsilon[2] = 0.5*(dvx[1] + dvy[0]); // shear strain rate + epsilon[3] = 0.5*(dvx[1] - dvy[0]); // rotation rate + +}/*}}}*/ void Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ /*Compute the 3d Strain Rate (6 components): Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 21619) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 21620) @@ -141,4 +141,5 @@ int Sid(); void SmbGemb(); + void StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); @@ -182,4 +183,5 @@ virtual void ComputeSigmaNN(void)=0; virtual void ComputeStressTensor(void)=0; + virtual void ComputeEsaStrainAndVorticity(void)=0; virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 21619) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 21620) @@ -53,4 +53,5 @@ void ComputeBasalStress(Vector* sigma_b); void ComputeDeviatoricStressTensor(); + void ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");}; void ComputeSigmaNN(){_error_("not implemented yet");}; void ComputeStressTensor(); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 21619) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 21620) @@ -48,4 +48,5 @@ void ComputeBasalStress(Vector* sigma_b){_error_("not implemented yet");}; void ComputeDeviatoricStressTensor(){_error_("not implemented yet");}; + void ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");}; void ComputeSigmaNN(){_error_("not implemented yet");}; void ComputeStressTensor(){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21619) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21620) @@ -50,4 +50,5 @@ void ComputeStressTensor(){_error_("not implemented yet");}; void ComputeDeviatoricStressTensor(){_error_("not implemented yet");}; + void ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");}; void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21619) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21620) @@ -410,4 +410,46 @@ this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum)); this->inputs->AddInput(new TriaInput(DeviatoricStresseffectiveEnum,&tau_e[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ +void Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[4]; /* epsilon=[exx,eyy,exy+ (shear),exy- (rotation)];*/ + IssmDouble strain_xx[NUMVERTICES]; + IssmDouble strain_yy[NUMVERTICES]; + IssmDouble strain_xy[NUMVERTICES]; + IssmDouble vorticity_xy[NUMVERTICES]; + GaussTria* gauss=NULL; + + /* Get node coordinates and dof list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Retrieve all inputs we will be needing: */ + Input* vx_input=this->GetInput(EsaEmotionEnum); _assert_(vx_input); + Input* vy_input=this->GetInput(EsaNmotionEnum); _assert_(vy_input); + + /* Start looping on the number of vertices: */ + gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + /*Compute strain rate and vorticity rate: */ + this->StrainRateESA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + + /*Compute Stress*/ + strain_xx[iv]=epsilon[0]; + strain_yy[iv]=epsilon[1]; + strain_xy[iv]=epsilon[2]; + vorticity_xy[iv]=epsilon[3]; + } + + /*Add Stress tensor components into inputs*/ + this->inputs->AddInput(new TriaInput(EsaStrainratexxEnum,&strain_xx[0],P1Enum)); + this->inputs->AddInput(new TriaInput(EsaStrainrateyyEnum,&strain_yy[0],P1Enum)); + this->inputs->AddInput(new TriaInput(EsaStrainratexyEnum,&strain_xy[0],P1Enum)); + this->inputs->AddInput(new TriaInput(EsaRotationrateEnum,&vorticity_xy[0],P1Enum)); /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 21619) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 21620) @@ -57,4 +57,5 @@ void ComputeBasalStress(Vector* sigma_b); void ComputeDeviatoricStressTensor(); + void ComputeEsaStrainAndVorticity(); void ComputeSigmaNN(); void ComputeStressTensor(); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21619) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21620) @@ -805,4 +805,8 @@ EsaNmotionEnum, EsaEmotionEnum, + EsaStrainratexxEnum, + EsaStrainratexyEnum, + EsaStrainrateyyEnum, + EsaRotationrateEnum, EsaDeltathicknessEnum, EsaUElasticEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21619) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21620) @@ -783,4 +783,8 @@ case EsaNmotionEnum : return "EsaNmotion"; case EsaEmotionEnum : return "EsaEmotion"; + case EsaStrainratexxEnum : return "EsaStrainratexx"; + case EsaStrainratexyEnum : return "EsaStrainratexy"; + case EsaStrainrateyyEnum : return "EsaStrainrateyy"; + case EsaRotationrateEnum : return "EsaRotationrate"; case EsaDeltathicknessEnum : return "EsaDeltathickness"; case EsaUElasticEnum : return "EsaUElastic"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21619) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21620) @@ -801,4 +801,8 @@ else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum; else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum; + else if (strcmp(name,"EsaStrainratexx")==0) return EsaStrainratexxEnum; + else if (strcmp(name,"EsaStrainratexy")==0) return EsaStrainratexyEnum; + else if (strcmp(name,"EsaStrainrateyy")==0) return EsaStrainrateyyEnum; + else if (strcmp(name,"EsaRotationrate")==0) return EsaRotationrateEnum; else if (strcmp(name,"EsaDeltathickness")==0) return EsaDeltathicknessEnum; else if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum; @@ -871,12 +875,12 @@ else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum; else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum; - else if (strcmp(name,"Option")==0) return OptionEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"Option")==0) return OptionEnum; else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; else if (strcmp(name,"OptionCell")==0) return OptionCellEnum; else if (strcmp(name,"OptionStruct")==0) return OptionStructEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; + else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum; else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; @@ -994,12 +998,12 @@ else if (strcmp(name,"Open")==0) return OpenEnum; else if (strcmp(name,"Air")==0) return AirEnum; - else if (strcmp(name,"Ice")==0) return IceEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"Melange")==0) return MelangeEnum; else if (strcmp(name,"Water")==0) return WaterEnum; else if (strcmp(name,"DataSet")==0) return DataSetEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Constraints")==0) return ConstraintsEnum; + else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; else if (strcmp(name,"Loads")==0) return LoadsEnum; else if (strcmp(name,"Materials")==0) return MaterialsEnum; Index: /issm/trunk-jpl/src/m/classes/esa.m =================================================================== --- /issm/trunk-jpl/src/m/classes/esa.m (revision 21619) +++ /issm/trunk-jpl/src/m/classes/esa.m (revision 21620) @@ -68,5 +68,5 @@ fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'); fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); - fielddisplay(self,'requested_outputs','additional outputs requested (default: EsaUmotion)'); + fielddisplay(self,'requested_outputs','additional outputs requested (e.g., EsaUmotion, EsaStrainratexx, EsaRotationrate)'); end % }}}