Changeset 17797

Timestamp:

04/21/14 20:25:06 (11 years ago)

Author:

seroussi

Message:

NEW: added sigma_nn in outputs

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Element.cpp (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Penta.h (modified) (1 diff)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/Elements/Tetra.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (1 diff)
• classes/Elements/Tria.h (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -827 +827 @@
                                 input=this->inputs->GetInput(output_enum);
                                 break;
+                        case SigmaNNEnum: 
+                                this->ComputeSigmaNN();
+                                input=this->inputs->GetInput(output_enum);
+                                break;
                         default:
                                 _error_("input "<<EnumToStringx(output_enum)<<" not found in element");

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -200 +200 @@
                 virtual IssmDouble SurfaceArea(void)=0;
                 virtual void   InputDepthAverageAtBase(int enum_type,int average_enum_type)=0;
+                virtual void   ComputeSigmaNN(void)=0;
                 virtual void   ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
                 virtual void   ComputeStrainRate(Vector<IssmDouble>* eps)=0;

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -56 +56 @@
                 void   ComputeBasalStress(Vector<IssmDouble>* sigma_b);
                 void   ComputeStrainRate(Vector<IssmDouble>* eps);
+                void   ComputeSigmaNN(){_error_("not implemented yet");};
                 void   ComputeStressTensor();
                 void   ComputeDeviatoricStressTensor();

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -54 +54 @@
                 IssmDouble  CharacteristicLength(void);
                 void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
+                void        ComputeSigmaNN(){_error_("not implemented yet");};
                 void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
                 void        ComputeStressTensor(){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Elements/Tetra.h

@@ -54 +54 @@
                 IssmDouble  CharacteristicLength(void){_error_("not implemented yet");};
                 void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
+                void        ComputeSigmaNN(){_error_("not implemented yet");};
                 void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
                 void        ComputeStressTensor(){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -145 +145 @@
 void  Tria::ComputeBasalStress(Vector<IssmDouble>* 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 " <<EnumToStringx(domaintype)<<", extrude mesh or call ComputeDeviatoricStressTensor");
+                Input* pressure_input=inputs->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<IssmDouble>(xyz_list);
+                delete gauss;
+        }
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -52 +52 @@
                 IssmDouble  CharacteristicLength(void);
                 void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
+                void        ComputeSigmaNN();
                 void        ComputeStrainRate(Vector<IssmDouble>* eps);
                 void        ComputeStressTensor();

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -539 +539 @@
         HydrologyWaterVxEnum,
         HydrologyWaterVyEnum,
+        SigmaNNEnum,
         StressTensorEnum,
         StressTensorxxEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -528 +528 @@
                 case HydrologyWaterVxEnum : return "HydrologyWaterVx";
                 case HydrologyWaterVyEnum : return "HydrologyWaterVy";
+                case SigmaNNEnum : return "SigmaNN";
                 case StressTensorEnum : return "StressTensor";
                 case StressTensorxxEnum : return "StressTensorxx";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -540 +540 @@
               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 +629 @@
               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;