Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17679) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17680) @@ -430,5 +430,5 @@ IssmDouble s_xx,s_xy,s_yy,s1,s2,stmp; IssmDouble J2s,Chi,Psi,PosPsi,NegPsi; - IssmDouble damage,sigma_xx,sigma_xy,sigma_yy; + IssmDouble damage,tau_xx,tau_xy,tau_yy; int equivstress; @@ -445,10 +445,10 @@ /*Compute stress tensor: */ - element->ComputeStressTensor(); + element->ComputeDeviatoricStressTensor(); /*retrieve what we need: */ - Input* sigma_xx_input = element->GetInput(StressTensorxxEnum); _assert_(sigma_xx_input); - Input* sigma_xy_input = element->GetInput(StressTensorxyEnum); _assert_(sigma_xy_input); - Input* sigma_yy_input = element->GetInput(StressTensoryyEnum); _assert_(sigma_yy_input); + Input* tau_xx_input = element->GetInput(DeviatoricStressxxEnum); _assert_(tau_xx_input); + Input* tau_xy_input = element->GetInput(DeviatoricStressxyEnum); _assert_(tau_xy_input); + Input* tau_yy_input = element->GetInput(DeviatoricStressyyEnum); _assert_(tau_yy_input); Input* damage_input = element->GetInput(DamageDEnum); _assert_(damage_input); @@ -462,12 +462,12 @@ damage_input->GetInputValue(&damage,gauss); - sigma_xx_input->GetInputValue(&sigma_xx,gauss); - sigma_xy_input->GetInputValue(&sigma_xy,gauss); - sigma_yy_input->GetInputValue(&sigma_yy,gauss); + tau_xx_input->GetInputValue(&tau_xx,gauss); + tau_xy_input->GetInputValue(&tau_xy,gauss); + tau_yy_input->GetInputValue(&tau_yy,gauss); /*Calculate effective stress components*/ - s_xx=sigma_xx/(1.-damage); - s_xy=sigma_xy/(1.-damage); - s_yy=sigma_yy/(1.-damage); + s_xx=tau_xx/(1.-damage); + s_xy=tau_xy/(1.-damage); + s_yy=tau_yy/(1.-damage); /*Calculate principal effective stresses*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17679) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17680) @@ -203,4 +203,5 @@ virtual void ComputeStrainRate(Vector* eps)=0; virtual void ComputeStressTensor(void)=0; + virtual void ComputeDeviatoricStressTensor(void)=0; virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17679) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17680) @@ -332,4 +332,56 @@ this->inputs->AddInput(new PentaInput(StressTensoryzEnum,&sigma_yz[0],P1Enum)); this->inputs->AddInput(new PentaInput(StressTensorzzEnum,&sigma_zz[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Penta::ComputeDeviatoricStressTensor {{{*/ +void Penta::ComputeDeviatoricStressTensor(){ + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble viscosity; + IssmDouble epsilon[6]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble tau_xx[NUMVERTICES]; + IssmDouble tau_yy[NUMVERTICES]; + IssmDouble tau_zz[NUMVERTICES]; + IssmDouble tau_xy[NUMVERTICES]; + IssmDouble tau_xz[NUMVERTICES]; + IssmDouble tau_yz[NUMVERTICES]; + GaussPenta* 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=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + + /* Start looping on the number of vertices: */ + gauss=new GaussPenta(); + for (int iv=0;ivGaussVertex(iv); + + /*Compute strain rate viscosity and pressure: */ + this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + + /*Compute Stress*/ + tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps + tau_yy[iv]=2*viscosity*epsilon[1]; + tau_zz[iv]=2*viscosity*epsilon[2]; + tau_xy[iv]=2*viscosity*epsilon[3]; + tau_xz[iv]=2*viscosity*epsilon[4]; + tau_yz[iv]=2*viscosity*epsilon[5]; + } + + /*Add Stress tensor components into inputs*/ + this->inputs->AddInput(new PentaInput(StressTensorxxEnum,&tau_xx[0],P1Enum)); + this->inputs->AddInput(new PentaInput(StressTensorxyEnum,&tau_xy[0],P1Enum)); + this->inputs->AddInput(new PentaInput(StressTensorxzEnum,&tau_xz[0],P1Enum)); + this->inputs->AddInput(new PentaInput(StressTensoryyEnum,&tau_yy[0],P1Enum)); + this->inputs->AddInput(new PentaInput(StressTensoryzEnum,&tau_yz[0],P1Enum)); + this->inputs->AddInput(new PentaInput(StressTensorzzEnum,&tau_zz[0],P1Enum)); /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17679) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17680) @@ -57,4 +57,5 @@ void ComputeStrainRate(Vector* eps); void ComputeStressTensor(); + void ComputeDeviatoricStressTensor(); void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17679) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17680) @@ -56,4 +56,5 @@ void ComputeStrainRate(Vector* eps){_error_("not implemented yet");}; void ComputeStressTensor(){_error_("not implemented yet");}; + void ComputeDeviatoricStressTensor(){_error_("not implemented yet");}; void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17679) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17680) @@ -56,4 +56,5 @@ void ComputeStrainRate(Vector* eps){_error_("not implemented yet");}; void ComputeStressTensor(){_error_("not implemented yet");}; + void ComputeDeviatoricStressTensor(){_error_("not implemented yet");}; void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17679) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17680) @@ -165,5 +165,5 @@ IssmDouble sigma_yz[NUMVERTICES]={0,0,0}; GaussTria* gauss=NULL; - int dim=2; + int meshtype,dim=2; /* Get node coordinates and dof list: */ @@ -171,4 +171,6 @@ /*Retrieve all inputs we will be needing: */ + this->FindParam(&meshtype,MeshTypeEnum); + if(meshtype==Mesh2DhorizontalEnum) _error_("stress tensor calculation not supported for mesh of type " <GetInput(PressureEnum); _assert_(pressure_input); Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); @@ -198,4 +200,55 @@ this->inputs->AddInput(new TriaInput(StressTensoryzEnum,&sigma_yz[0],P1Enum)); this->inputs->AddInput(new TriaInput(StressTensorzzEnum,&sigma_zz[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Tria::ComputeDeviatoricStressTensor {{{*/ +void Tria::ComputeDeviatoricStressTensor(){ + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble viscosity; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble tau_xx[NUMVERTICES]; + IssmDouble tau_yy[NUMVERTICES]; + IssmDouble tau_zz[NUMVERTICES]={0,0,0}; + IssmDouble tau_xy[NUMVERTICES]; + IssmDouble tau_xz[NUMVERTICES]={0,0,0}; + IssmDouble tau_yz[NUMVERTICES]={0,0,0}; + GaussTria* gauss=NULL; + int meshtype,dim=2; + + /* Get node coordinates and dof list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Retrieve all inputs we will be needing: */ + this->FindParam(&meshtype,MeshTypeEnum); + if(meshtype!=Mesh2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of vertices: */ + gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + /*Compute strain rate and viscosity: */ + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input); + + /*Compute Stress*/ + tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps + tau_yy[iv]=2*viscosity*epsilon[1]; + tau_xy[iv]=2*viscosity*epsilon[2]; + } + + /*Add Stress tensor components into inputs*/ + this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum)); + this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum)); + this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum)); + this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum)); + this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum)); + this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[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 17679) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17680) @@ -54,4 +54,5 @@ void ComputeStrainRate(Vector* eps); void ComputeStressTensor(); + void ComputeDeviatoricStressTensor(); void ComputeSurfaceNormalVelocity(); void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);