Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17796) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17797) @@ -827,4 +827,8 @@ input=this->inputs->GetInput(output_enum); break; + case SigmaNNEnum: + this->ComputeSigmaNN(); + input=this->inputs->GetInput(output_enum); + break; default: _error_("input "<* sigma_b)=0; virtual void ComputeStrainRate(Vector* eps)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17796) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17797) @@ -56,4 +56,5 @@ void ComputeBasalStress(Vector* sigma_b); void ComputeStrainRate(Vector* eps); + void ComputeSigmaNN(){_error_("not implemented yet");}; void ComputeStressTensor(); void ComputeDeviatoricStressTensor(); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17796) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17797) @@ -54,4 +54,5 @@ IssmDouble CharacteristicLength(void); void ComputeBasalStress(Vector* sigma_b){_error_("not implemented yet");}; + void ComputeSigmaNN(){_error_("not implemented yet");}; void ComputeStrainRate(Vector* eps){_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 17796) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17797) @@ -54,4 +54,5 @@ IssmDouble CharacteristicLength(void){_error_("not implemented yet");}; void ComputeBasalStress(Vector* sigma_b){_error_("not implemented yet");}; + void ComputeSigmaNN(){_error_("not implemented yet");}; void ComputeStrainRate(Vector* eps){_error_("not implemented yet");}; void ComputeStressTensor(){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17796) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17797) @@ -145,4 +145,59 @@ void Tria::ComputeBasalStress(Vector* eps){ _error_("Not Implemented yet"); +} +/*}}}*/ +/*FUNCTION Tria::ComputeSigmaNN {{{*/ +void Tria::ComputeSigmaNN(){ + + if(~IsOnBase()){ + IssmDouble sigma_nn=0; + this->inputs->AddInput(new TriaInput(SigmaNNEnum,&sigma_nn,P0Enum)); + return; + } + else{ + IssmDouble* xyz_list=NULL; + IssmDouble pressure,viscosity,sigma_nn; + IssmDouble sigma_xx,sigma_xy,sigma_yy; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble base_normal[2]; + int domaintype,dim=2; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will be needing: */ + this->FindParam(&domaintype,DomainTypeEnum); + if(domaintype==Domain2DhorizontalEnum) _error_("stress tensor calculation not supported for mesh of type " <GetInput(PressureEnum); _assert_(pressure_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of vertices: */ + Gauss* gauss=NewGaussBase(1); + gauss->GaussPoint(0); + + /*Compute strain rate viscosity and pressure: */ + this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + this->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,NULL); + pressure_input->GetInputValue(&pressure,gauss); + + /*Compute Stress*/ + sigma_xx=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure + sigma_yy=2*viscosity*epsilon[1]-pressure; + sigma_xy=2*viscosity*epsilon[2]; + + /*Get normal vector to the bed */ + NormalBase(&base_normal[0],xyz_list); + + /*Compute sigma_nn*/ + sigma_nn=sigma_xx*base_normal[0]*base_normal[0] + sigma_xy*base_normal[0]*base_normal[1] + sigma_yy*base_normal[1]*base_normal[1]; + + /*Add Stress tensor components into inputs*/ + this->inputs->AddInput(new TriaInput(SigmaNNEnum,&sigma_nn,P0Enum)); + + /*Clean up and return*/ + xDelete(xyz_list); + delete gauss; + } } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17796) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17797) @@ -52,4 +52,5 @@ IssmDouble CharacteristicLength(void); void ComputeBasalStress(Vector* sigma_b); + void ComputeSigmaNN(); void ComputeStrainRate(Vector* eps); void ComputeStressTensor(); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17796) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17797) @@ -539,4 +539,5 @@ HydrologyWaterVxEnum, HydrologyWaterVyEnum, + SigmaNNEnum, StressTensorEnum, StressTensorxxEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17796) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17797) @@ -528,4 +528,5 @@ case HydrologyWaterVxEnum : return "HydrologyWaterVx"; case HydrologyWaterVyEnum : return "HydrologyWaterVy"; + case SigmaNNEnum : return "SigmaNN"; case StressTensorEnum : return "StressTensor"; case StressTensorxxEnum : return "StressTensorxx"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17796) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17797) @@ -540,4 +540,5 @@ else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; + else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum; else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum; else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; - else if (strcmp(name,"SubelementMigration2")==0) return SubelementMigration2Enum; else stage=6; } if(stage==6){ - if (strcmp(name,"Contact")==0) return ContactEnum; + if (strcmp(name,"SubelementMigration2")==0) return SubelementMigration2Enum; + else if (strcmp(name,"Contact")==0) return ContactEnum; else if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum; else if (strcmp(name,"QmuMaskGroundediceLevelset")==0) return QmuMaskGroundediceLevelsetEnum;