Changeset 18613

Timestamp:

10/10/14 13:04:07 (10 years ago)

Author:

srebuffi

Message:

NEW: added stress intensity factor calculation

Location:

issm/trunk-jpl/src

Files:

• c/Makefile.am (modified) (1 diff)
• c/analyses/LevelsetAnalysis.cpp (modified) (2 diffs)
• c/classes/Elements/Element.cpp (modified) (2 diffs)
• c/classes/Elements/Element.h (modified) (1 diff)
• c/classes/Elements/Penta.cpp (modified) (1 diff)
• c/classes/Elements/Penta.h (modified) (1 diff)
• c/classes/Elements/Seg.h (modified) (1 diff)
• c/classes/Elements/Tetra.h (modified) (1 diff)
• c/classes/Elements/Tria.h (modified) (1 diff)
• c/classes/FemModel.cpp (modified) (1 diff)
• c/classes/FemModel.h (modified) (1 diff)
• c/shared/Elements/StressIntensityIntegralWeight.cpp (added)
• c/shared/Elements/elements.h (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (2 diffs)
• c/shared/Enum/EnumToStringx.cpp (modified) (2 diffs)
• c/shared/Enum/StringToEnumx.cpp (modified) (7 diffs)
• m/enum/EnumDefinitions.py (modified) (2 diffs)
• m/enum/LevelsetStabilizationEnum.m (added)
• m/enum/StressIntensityFactorEnum.m (added)

issm/trunk-jpl/src/c/Makefile.am

@@ -213 +213 @@
                                         ./shared/Elements/elements.h\
                                         ./shared/Elements/Cuffey.cpp\
+                                        ./shared/Elements/StressIntensityIntegralWeight.cpp\
                                         ./shared/Elements/Paterson.cpp\
                                         ./shared/Elements/Arrhenius.cpp\

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

@@ -65 +65 @@
         /*parameters: */
         bool save_results;
+        int  stabilization=1;
         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+        //femmodel->parameters->FindParam(&stabilization,LevelsetStabilizationEnum);
 
         /*activate formulation: */
@@ -71 +73 @@
 
         if(VerboseSolution()) _printf0_("call computational core:\n");
-        solutionsequence_linear(femmodel);
+        if(stabilization==4){
+                solutionsequence_fct(femmodel);
+        }
+        else{
+                solutionsequence_linear(femmodel);
+        }
 
         if(save_results){

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

@@ -1105 +1105 @@
                                 name==SedimentHeadOldEnum ||
                                 name==EplHeadOldEnum ||
+                                name==StressIntensityFactorEnum ||
                                 name==HydrologydcEplThicknessOldEnum ||
                                 name==HydrologydcEplInitialThicknessEnum ||
@@ -1203 +1204 @@
                                 input=this->inputs->GetInput(output_enum);
                                 break;
+                        case StressIntensityFactorEnum:
+                                this->StressIntensityFactor();
+                                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

@@ -219 +219 @@
                 virtual void   ComputeStressTensor(void)=0;
                 virtual void   ComputeDeviatoricStressTensor(void)=0;
+                virtual void    StressIntensityFactor(void)=0;
 
                 virtual void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0;

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

@@ -345 +345 @@
         /*Clean up and return*/
         delete gauss;
+}
+/*}}}*/
+void       Penta::StressIntensityFactor(){/*{{{*/
+
+        /* Check if we are on the base */
+        if(!IsOnBase()) return;
+
+        IssmDouble  ki[6]={0.};
+        IssmDouble  const_grav=9.81;
+        IssmDouble  rho_ice=900;
+        IssmDouble  rho_water=1000;
+        IssmDouble  Jdet[3];
+        IssmDouble  pressure,vx,vy,vel,deviaxx,deviaxy,deviayy,water_depth,prof,stress_xx,thickness;
+
+        Penta* penta=this;
+        for(;;){
+        
+                IssmDouble  xyz_list[NUMVERTICES][3];
+                /* Get node coordinates and dof list: */
+                ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
+
+                ///*Compute the Jacobian for the vertical integration*/
+                Jdet[0]=(xyz_list[3][2]-xyz_list[0][2])*0.5;
+                Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5;
+                Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5;
+        
+                /*Retrieve all inputs we will need*/
+                Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
+                Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
+                Input* vel_input=inputs->GetInput(VelEnum);                                _assert_(vel_input);
+                Input* pressure_input=inputs->GetInput(PressureEnum);                      _assert_(pressure_input);
+                Input* deviaxx_input=inputs->GetInput(DeviatoricStressxxEnum);             _assert_(deviaxx_input);
+                Input* deviaxy_input=inputs->GetInput(DeviatoricStressxyEnum);             _assert_(deviaxy_input);
+                Input* deviayy_input=inputs->GetInput(DeviatoricStressyyEnum);             _assert_(deviayy_input);
+                Input* surface_input=inputs->GetInput(SurfaceEnum);                                                             _assert_(surface_input);
+                Input* thickness_input=inputs->GetInput(ThicknessEnum);                                                 _assert_(thickness_input);
+                
+                /* Start looping on the number of 2D vertices: */
+                for(int ig=0;ig<3;ig++){
+                        GaussPenta* gauss=new GaussPenta(ig,3+ig,11);
+                        for (int iv=gauss->begin();iv<gauss->end();iv++){
+                                gauss->GaussPoint(iv);
+
+                                /* Get the value we need*/
+                                pressure_input->GetInputValue(&pressure,gauss);
+                                vx_input->GetInputValue(&vx,gauss);
+                                vy_input->GetInputValue(&vy,gauss);
+                                vel_input->GetInputValue(&vel,gauss);
+                                deviaxx_input->GetInputValue(&deviaxx,gauss);
+                                deviaxy_input->GetInputValue(&deviaxy,gauss);
+                                deviayy_input->GetInputValue(&deviayy,gauss);
+                                surface_input->GetInputValue(&water_depth,gauss);
+                                thickness_input->GetInputValue(&thickness,gauss);
+                                prof=water_depth-penta->GetZcoord(&xyz_list[0][0],gauss);
+
+                                /*stress_xx= Deviatoric stress along the ice flow direction plus cryostatic pressure */
+                                stress_xx=(vx*vx*(deviaxx)+vy*vy*(deviayy)+2*vy*vx*deviaxy)/(vel*vel+1.e-6);
+
+                                if(prof<water_depth&prof<thickness){
+                                        /* Compute the local stress intensity factor*/ 
+                                        ki[ig]+=Jdet[ig]*gauss->weight*stress_xx*StressIntensityIntegralWeight(prof,min(water_depth,thickness),thickness);
+                                }
+                        }
+                        delete gauss;
+                }
+                        
+                /*Stop if we have reached the surface/base*/
+                if(penta->IsOnSurface()) break;
+                
+                /*get upper Penta*/
+                penta=penta->GetUpperPenta();
+                _assert_(penta->Id()!=this->id);
+        }
+
+        /*Add input*/
+        this->inputs->AddInput(new PentaInput(StressIntensityFactorEnum,&ki[0],P1Enum));
+        this->InputExtrude(StressIntensityFactorEnum,-1);
 }
 /*}}}*/

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

@@ -57 +57 @@
                 void   ComputeStressTensor();
                 void   ComputeDeviatoricStressTensor();
+                void   StressIntensityFactor();
                 void   Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
                 void   ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);

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

@@ -57 +57 @@
                 void        ComputeStressTensor(){_error_("not implemented yet");};
                 void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
+                void        StressIntensityFactor(void){_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");};

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

@@ -57 +57 @@
                 void        ComputeStressTensor(){_error_("not implemented yet");};
                 void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
+                void        StressIntensityFactor(void){_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);

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

@@ -56 +56 @@
                 void        ComputeDeviatoricStressTensor();
                 void        ComputeSurfaceNormalVelocity();
+                void        StressIntensityFactor(void){_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);

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -1636 +1636 @@
 }
 /*}}}*/
+void FemModel::StressIntensityFactorx(){/*{{{*/
+
+        /*Initialize input as 0*/
+        InputUpdateFromConstantx(this,0.,StressIntensityFactorEnum);
+
+        /*Update input for basal element only*/
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+                element->StressIntensityFactor();
+        }
+}
+        /*}}}*/
 #ifdef  _HAVE_DAKOTA_
 void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -81 +81 @@
                 void ElementResponsex(IssmDouble* presponse,int response_enum);
                 void BalancethicknessMisfitx(IssmDouble* pV);
+                void StressIntensityFactorx();
                 #ifdef  _HAVE_DAKOTA_
                 void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);

issm/trunk-jpl/src/c/shared/Elements/elements.h

@@ -22 +22 @@
                                           IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
 IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.);
+IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness);
 
 /*Print arrays*/

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

@@ -204 +204 @@
         DamageEnum,
         NewDamageEnum,
+        StressIntensityFactorEnum,
         MaterialsRhoIceEnum,
         MaterialsRhoSeawaterEnum,
@@ -400 +401 @@
         MeshdeformationAnalysisEnum,
         LevelsetAnalysisEnum,
+        LevelsetStabilizationEnum,
         ExtrapolationAnalysisEnum,
         LsfReinitializationAnalysisEnum,

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

@@ -212 +212 @@
                 case DamageEnum : return "Damage";
                 case NewDamageEnum : return "NewDamage";
+                case StressIntensityFactorEnum : return "StressIntensityFactor";
                 case MaterialsRhoIceEnum : return "MaterialsRhoIce";
                 case MaterialsRhoSeawaterEnum : return "MaterialsRhoSeawater";
@@ -403 +404 @@
                 case MeshdeformationAnalysisEnum : return "MeshdeformationAnalysis";
                 case LevelsetAnalysisEnum : return "LevelsetAnalysis";
+                case LevelsetStabilizationEnum : return "LevelsetStabilization";
                 case ExtrapolationAnalysisEnum : return "ExtrapolationAnalysis";
                 case LsfReinitializationAnalysisEnum : return "LsfReinitializationAnalysis";

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

@@ -215 +215 @@
               else if (strcmp(name,"Damage")==0) return DamageEnum;
               else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
+              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
               else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
               else if (strcmp(name,"MaterialsRhoSeawater")==0) return MaterialsRhoSeawaterEnum;
@@ -259 +260 @@
               else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
               else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
-              else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
+              if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
+              else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
               else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
               else if (strcmp(name,"QmuMassFluxSegmentsPresent")==0) return QmuMassFluxSegmentsPresentEnum;
@@ -382 +383 @@
               else if (strcmp(name,"FlaimAnalysis")==0) return FlaimAnalysisEnum;
               else if (strcmp(name,"FlaimSolution")==0) return FlaimSolutionEnum;
-              else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
+              if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
+              else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
               else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
               else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
@@ -412 +413 @@
               else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;
               else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
+              else if (strcmp(name,"LevelsetStabilization")==0) return LevelsetStabilizationEnum;
               else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
               else if (strcmp(name,"LsfReinitializationAnalysis")==0) return LsfReinitializationAnalysisEnum;
@@ -504 +506 @@
               else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
               else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
-              else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
-              else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
+              if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
+              else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
+              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
               else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
               else if (strcmp(name,"Converged")==0) return ConvergedEnum;
@@ -627 +629 @@
               else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
               else if (strcmp(name,"Step")==0) return StepEnum;
-              else if (strcmp(name,"Time")==0) return TimeEnum;
-              else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
+              if (strcmp(name,"Time")==0) return TimeEnum;
+              else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
+              else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
               else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
               else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
@@ -750 +752 @@
               else if (strcmp(name,"BasalforcingsOceanSsh")==0) return BasalforcingsOceanSshEnum;
               else if (strcmp(name,"BasalforcingsOceanVx")==0) return BasalforcingsOceanVxEnum;
-              else if (strcmp(name,"BasalforcingsOceanVy")==0) return BasalforcingsOceanVyEnum;
-              else if (strcmp(name,"SurfaceforcingsRhoAir")==0) return SurfaceforcingsRhoAirEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum;
+              if (strcmp(name,"BasalforcingsOceanVy")==0) return BasalforcingsOceanVyEnum;
+              else if (strcmp(name,"SurfaceforcingsRhoAir")==0) return SurfaceforcingsRhoAirEnum;
+              else if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum;
               else if (strcmp(name,"SurfaceforcingsAirLinDragCoef")==0) return SurfaceforcingsAirLinDragCoefEnum;
               else if (strcmp(name,"SurfaceforcingsAirQuadDragCoef")==0) return SurfaceforcingsAirQuadDragCoefEnum;

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -204 +204 @@
 def DamageEnum(): return StringToEnum("Damage")[0]
 def NewDamageEnum(): return StringToEnum("NewDamage")[0]
+def StressIntensityFactorEnum(): return StringToEnum("StressIntensityFactor")[0]
 def MaterialsRhoIceEnum(): return StringToEnum("MaterialsRhoIce")[0]
 def MaterialsRhoSeawaterEnum(): return StringToEnum("MaterialsRhoSeawater")[0]
@@ -395 +396 @@
 def MeshdeformationAnalysisEnum(): return StringToEnum("MeshdeformationAnalysis")[0]
 def LevelsetAnalysisEnum(): return StringToEnum("LevelsetAnalysis")[0]
+def LevelsetStabilizationEnum(): return StringToEnum("LevelsetStabilization")[0]
 def ExtrapolationAnalysisEnum(): return StringToEnum("ExtrapolationAnalysis")[0]
 def LsfReinitializationAnalysisEnum(): return StringToEnum("LsfReinitializationAnalysis")[0]