Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 18612) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 18613) @@ -213,4 +213,5 @@ ./shared/Elements/elements.h\ ./shared/Elements/Cuffey.cpp\ + ./shared/Elements/StressIntensityIntegralWeight.cpp\ ./shared/Elements/Paterson.cpp\ ./shared/Elements/Arrhenius.cpp\ Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 18612) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 18613) @@ -65,5 +65,7 @@ /*parameters: */ bool save_results; + int stabilization=1; femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + //femmodel->parameters->FindParam(&stabilization,LevelsetStabilizationEnum); /*activate formulation: */ @@ -71,5 +73,10 @@ if(VerboseSolution()) _printf0_("call computational core:\n"); - solutionsequence_linear(femmodel); + if(stabilization==4){ + solutionsequence_fct(femmodel); + } + else{ + solutionsequence_linear(femmodel); + } if(save_results){ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18612) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18613) @@ -1105,4 +1105,5 @@ name==SedimentHeadOldEnum || name==EplHeadOldEnum || + name==StressIntensityFactorEnum || name==HydrologydcEplThicknessOldEnum || name==HydrologydcEplInitialThicknessEnum || @@ -1203,4 +1204,8 @@ input=this->inputs->GetInput(output_enum); break; + case StressIntensityFactorEnum: + this->StressIntensityFactor(); + input=this->inputs->GetInput(output_enum); + break; default: _error_("input "<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();ivend();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(profweight*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); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18612) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18613) @@ -57,4 +57,5 @@ 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); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18612) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18613) @@ -57,4 +57,5 @@ 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");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18612) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18613) @@ -57,4 +57,5 @@ 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); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18612) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18613) @@ -56,4 +56,5 @@ 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); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18612) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18613) @@ -1636,4 +1636,16 @@ } /*}}}*/ +void FemModel::StressIntensityFactorx(){/*{{{*/ + + /*Initialize input as 0*/ + InputUpdateFromConstantx(this,0.,StressIntensityFactorEnum); + + /*Update input for basal element only*/ + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->StressIntensityFactor(); + } +} + /*}}}*/ #ifdef _HAVE_DAKOTA_ void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18612) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18613) @@ -81,4 +81,5 @@ 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); Index: /issm/trunk-jpl/src/c/shared/Elements/StressIntensityIntegralWeight.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/StressIntensityIntegralWeight.cpp (revision 18613) +++ /issm/trunk-jpl/src/c/shared/Elements/StressIntensityIntegralWeight.cpp (revision 18613) @@ -0,0 +1,31 @@ +/* \file StressIntensityIntegralWeight.cpp + * \Weight to integrate the stress along the ice flow direction to compute the stress intensity factor : Linear Fracture Mechanics (see "Combining damage and fracture mechanics to model calving",Krug 2014 in the appendix) + */ + +#include + +#include "../Numerics/types.h" + +IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness){ + + /*output: */ + IssmDouble beta; + + /*intermediaries: */ + IssmDouble M1,M2,M3,x,y,d; + const double pi = 3.141592653589793; + x = water_depth/thickness; + y = depth; + d = water_depth; + + M1 = 0.0719768-1.513476*x-61.1001*pow(x,2)+1554.95*pow(x,3)-14583.8*pow(x,4)+71590.7*pow(x,5)-205384*pow(x,6)+356469*pow(x,7)-368270*pow(x,8)+208233*pow(x,9)-49544*pow(x,10); + //printf("M1 : %g",M1); + M2 = 0.246984+6.47583*x+176.456*pow(x,2)-4058.76*pow(x,3)+37303.8*pow(x,4)-181755*pow(x,5)+520551*pow(x,6)-904370*pow(x,7)+936863*pow(x,8)-531940*pow(x,9)+127291*pow(x,10); + //printf("M2 : %g",M2); + M3 = 0.529659-22.3235*x+532.074*pow(x,2)-5479.53*pow(x,3)+28592.2*pow(x,4)-81388.6*pow(x,5)+128746*pow(x,6)-106246*pow(x,7)+35780.7*pow(x,8); + //printf("M3 : %g",M3); + + beta = 2/sqrt(2*pi*(d-y))*(1+M1*sqrt(1-y/d)+M2*(1-y/d)+M3*pow((1-y/d),1.5)); + + return beta; +} Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 18612) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 18613) @@ -22,4 +22,5 @@ IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.); +IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness); /*Print arrays*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18612) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18613) @@ -204,4 +204,5 @@ DamageEnum, NewDamageEnum, + StressIntensityFactorEnum, MaterialsRhoIceEnum, MaterialsRhoSeawaterEnum, @@ -400,4 +401,5 @@ MeshdeformationAnalysisEnum, LevelsetAnalysisEnum, + LevelsetStabilizationEnum, ExtrapolationAnalysisEnum, LsfReinitializationAnalysisEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18612) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18613) @@ -212,4 +212,5 @@ case DamageEnum : return "Damage"; case NewDamageEnum : return "NewDamage"; + case StressIntensityFactorEnum : return "StressIntensityFactor"; case MaterialsRhoIceEnum : return "MaterialsRhoIce"; case MaterialsRhoSeawaterEnum : return "MaterialsRhoSeawater"; @@ -403,4 +404,5 @@ case MeshdeformationAnalysisEnum : return "MeshdeformationAnalysis"; case LevelsetAnalysisEnum : return "LevelsetAnalysis"; + case LevelsetStabilizationEnum : return "LevelsetStabilization"; case ExtrapolationAnalysisEnum : return "ExtrapolationAnalysis"; case LsfReinitializationAnalysisEnum : return "LsfReinitializationAnalysis"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18612) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18613) @@ -215,4 +215,5 @@ 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,9 +260,9 @@ 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,9 +383,9 @@ 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,4 +413,5 @@ 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,10 +506,10 @@ 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,10 +629,10 @@ 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,10 +752,10 @@ 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; Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18612) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18613) @@ -204,4 +204,5 @@ 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,4 +396,5 @@ 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] Index: /issm/trunk-jpl/src/m/enum/LevelsetStabilizationEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/LevelsetStabilizationEnum.m (revision 18613) +++ /issm/trunk-jpl/src/m/enum/LevelsetStabilizationEnum.m (revision 18613) @@ -0,0 +1,11 @@ +function macro=LevelsetStabilizationEnum() +%LEVELSETSTABILIZATIONENUM - Enum of LevelsetStabilization +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=LevelsetStabilizationEnum() + +macro=StringToEnum('LevelsetStabilization'); Index: /issm/trunk-jpl/src/m/enum/StressIntensityFactorEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/StressIntensityFactorEnum.m (revision 18613) +++ /issm/trunk-jpl/src/m/enum/StressIntensityFactorEnum.m (revision 18613) @@ -0,0 +1,11 @@ +function macro=StressIntensityFactorEnum() +%STRESSINTENSITYFACTORENUM - Enum of StressIntensityFactor +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=StressIntensityFactorEnum() + +macro=StringToEnum('StressIntensityFactor');