Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 21388) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 21389) @@ -175,4 +175,8 @@ iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum); switch(materialstype){ + case MatenhancediceEnum: + iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_E",MaterialsRheologyEEnum); + break; case MatdamageiceEnum: iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum); Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 21388) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 21389) @@ -119,4 +119,8 @@ iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum); switch(materialstype){ + case MatenhancediceEnum: + iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_E",MaterialsRheologyEEnum); + break; case MatdamageiceEnum: iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum); Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 21388) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 21389) @@ -31,4 +31,5 @@ this->element=NULL; this->isdamaged=false; + this->isenhanced=false; return; } @@ -50,8 +51,17 @@ int materialtype; iomodel->FindConstant(&materialtype,"md.materials.type"); - if(materialtype==MatdamageiceEnum) this->isdamaged = true; - else if(materialtype==MaticeEnum) this->isdamaged = false; + if(materialtype==MatdamageiceEnum){ + this->isdamaged = true; + this->isenhanced = false; + } + else if(materialtype==MaticeEnum){ + this->isdamaged = false; + this->isenhanced = false; + } + else if(materialtype==MatenhancediceEnum){ + this->isdamaged = false; + this->isenhanced = true; + } else _error_("Material type not recognized"); - return; @@ -78,4 +88,5 @@ matice->element =(Element*)this->helement->delivers(); matice->isdamaged = this->isdamaged; + matice->isenhanced = this->isenhanced; return matice; @@ -95,4 +106,5 @@ matice->element =element_in; matice->isdamaged = this->isdamaged; + matice->isenhanced = this->isenhanced; return matice; @@ -124,4 +136,5 @@ MARSHALLING(mid); MARSHALLING(isdamaged); + MARSHALLING(isenhanced); this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); this->element=(Element*)this->helement->delivers(); @@ -207,4 +220,22 @@ } /*}}}*/ +IssmDouble Matice::GetE(){/*{{{*/ + + _assert_(this->isenhanced); + /*Output*/ + IssmDouble E; + if(this->isenhanced)element->inputs->GetInputAverage(&E,MaterialsRheologyEEnum); + return E; +} +/*}}}*/ +IssmDouble Matice::GetEbar(){/*{{{*/ + + _assert_(this->isenhanced); + /*Output*/ + IssmDouble Ebar; + if(this->isenhanced)element->inputs->GetInputAverage(&Ebar,MaterialsRheologyEbarEnum); + return Ebar; +} +/*}}}*/ IssmDouble Matice::GetN(){/*{{{*/ @@ -219,4 +250,9 @@ return this->isdamaged; +} +/*}}}*/ +bool Matice::IsEnhanced(){/*{{{*/ + + return this->isenhanced; } /*}}}*/ @@ -225,8 +261,8 @@ (1-D) B viscosity= ------------------------- - 2 eps_eff ^[(n-1)/n] - - where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate - and n the flow law exponent. + 2 E^[1/n] eps_eff ^[(n-1)/n] + + where viscosity is the viscosity, B the flow law parameter , eps_eff is the effective strain rate, + n the flow law exponent, and E is the enhancement factor. If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we @@ -238,5 +274,5 @@ /*Intermediary: */ - IssmDouble B,D=0.,n; + IssmDouble B,D=0.,E=1.,n; /*Get B and n*/ @@ -247,8 +283,12 @@ _assert_(D>=0. && D<1.); } + if(this->isenhanced){ + E=GetE(); + _assert_(E>0.); + } if (n==1.){ - /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */ - viscosity=(1.-D)*B/2.; + /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2E: */ + viscosity=(1.-D)*B/(2.*E); } else{ @@ -262,5 +302,5 @@ else{ - viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n)); + viscosity=(1.-D)*B/(2.*pow(E,1./n)*pow(eps_eff,(n-1.)/n)); } } @@ -277,7 +317,8 @@ (1-D) B viscosity= ------------------------- - 2 eps_eff ^[(n-1)/n] - - where B the flow law parameter, eps_eff is the effective strain rate and n the flow law exponent. + 2 E^[1/n] eps_eff ^[(n-1)/n] + + where B the flow law parameter, eps_eff is the effective strain rate, n the flow law exponent, + and E is the enhancement factor. If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we @@ -289,5 +330,5 @@ /*Intermediary: */ - IssmDouble B,D=0.,n; + IssmDouble B,D=0.,E=1.,n; /*Get B and n*/ @@ -298,8 +339,12 @@ _assert_(D>=0. && D<1.); } + if(this->isenhanced){ + E=GetEbar(); + _assert_(E>0.); + } if (n==1.){ - /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */ - viscosity=(1.-D)*B/2.; + /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2E: */ + viscosity=(1.-D)*B/(2.*E); } else{ @@ -312,5 +357,5 @@ else{ - viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n)); + viscosity=(1.-D)*B/(2.*pow(E,1./n)*pow(eps_eff,(n-1.)/n)); } } @@ -475,5 +520,5 @@ /*Intermediary: */ - IssmDouble D=0.,n; + IssmDouble D=0.,E=1.,n; /*Get B and n*/ @@ -483,12 +528,16 @@ _assert_(D>=0. && D<1.); } + if(this->isenhanced){ + E=GetE(); + _assert_(E>0.); + } if(n==1.){ - /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */ - dmudB=(1.-D)/2.; + /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2E: */ + dmudB=(1.-D)/(2.*E); } else{ if(eps_eff==0.) dmudB = 0.; - else dmudB = (1.-D)/(2.*pow(eps_eff,(n-1.)/n)); + else dmudB = (1.-D)/(2.*pow(E,1./n)*pow(eps_eff,(n-1.)/n)); } @@ -503,5 +552,5 @@ /*Intermediary: */ - IssmDouble n,B; + IssmDouble n,B,E=1.; /*Get B and n*/ @@ -509,12 +558,16 @@ B=GetBbar(); _assert_(this->isdamaged); + if(this->isenhanced){ + E=GetE(); + _assert_(E>0.); + } if(n==1.){ - /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */ - dmudD=-B/2.; + /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2E: */ + dmudD=-B/(2.*E); } else{ if(eps_eff==0.) dmudD = 0.; - else dmudD = -B/(2.*pow(eps_eff,(n-1.)/n)); + else dmudD = -B/(2.*pow(E,1./n)*pow(eps_eff,(n-1.)/n)); } @@ -636,5 +689,5 @@ IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ IssmDouble epsilon2d[2];/* epsilon=[exx,exy]; */ - IssmDouble eps_eff; + IssmDouble eps_eff,E=1.0; if(dim==3){ Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 21388) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 21389) @@ -27,4 +27,5 @@ int mid; bool isdamaged; + bool isenhanced; Hook *helement; Element *element; @@ -70,6 +71,9 @@ IssmDouble GetD(); IssmDouble GetDbar(); + IssmDouble GetE(); + IssmDouble GetEbar(); IssmDouble GetN(); bool IsDamage(); + bool IsEnhanced(); void ResetHooks(); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 21388) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 21389) @@ -82,4 +82,5 @@ case MaticeEnum: case MatdamageiceEnum: + case MatenhancediceEnum: case MatestarEnum: iomodel->FindConstant(&this->rho_ice,"md.materials.rho_ice"); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 21388) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 21389) @@ -56,4 +56,20 @@ case 2: elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); + break; + case 3: + break; + default: + _error_("Mesh not supported yet"); + } + break; + case MatenhancediceEnum: + iomodel->FetchDataToInput(elements,"md.materials.rheology_B",MaterialsRheologyBEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_n",MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,"md.materials.rheology_E",MaterialsRheologyEEnum); + for (i=0;inumberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel)); + switch(iomodel->domaindim){ + case 2: + elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); + elements->InputDuplicate(MaterialsRheologyEEnum,MaterialsRheologyEbarEnum); break; case 3: Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21388) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21389) @@ -214,4 +214,6 @@ MaterialsRheologyLawEnum, MaterialsRheologyNEnum, + MaterialsRheologyEEnum, + MaterialsRheologyEbarEnum, MaterialsRheologyEcEnum, MaterialsRheologyEcbarEnum, @@ -828,4 +830,5 @@ MaticeEnum, MatdamageiceEnum, + MatenhancediceEnum, MatestarEnum, MatparEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21388) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21389) @@ -220,4 +220,6 @@ case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw"; case MaterialsRheologyNEnum : return "MaterialsRheologyN"; + case MaterialsRheologyEEnum : return "MaterialsRheologyE"; + case MaterialsRheologyEbarEnum : return "MaterialsRheologyEbar"; case MaterialsRheologyEcEnum : return "MaterialsRheologyEc"; case MaterialsRheologyEcbarEnum : return "MaterialsRheologyEcbar"; @@ -803,4 +805,5 @@ case MaticeEnum : return "Matice"; case MatdamageiceEnum : return "Matdamageice"; + case MatenhancediceEnum : return "Matenhancedice"; case MatestarEnum : return "Matestar"; case MatparEnum : return "Matpar"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21388) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21389) @@ -223,4 +223,6 @@ else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum; else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum; + else if (strcmp(name,"MaterialsRheologyE")==0) return MaterialsRheologyEEnum; + else if (strcmp(name,"MaterialsRheologyEbar")==0) return MaterialsRheologyEbarEnum; else if (strcmp(name,"MaterialsRheologyEc")==0) return MaterialsRheologyEcEnum; else if (strcmp(name,"MaterialsRheologyEcbar")==0) return MaterialsRheologyEcbarEnum; @@ -258,10 +260,10 @@ else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum; else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum; - else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum; - else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum; + if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum; + else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum; + else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum; else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum; else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; @@ -381,10 +383,10 @@ else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; - else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; - else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum; + if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; + else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; + else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum; else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; else if (strcmp(name,"SmbV")==0) return SmbVEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; else if (strcmp(name,"Vx")==0) return VxEnum; - else if (strcmp(name,"VxPicard")==0) return VxPicardEnum; - else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"Vy")==0) return VyEnum; + if (strcmp(name,"VxPicard")==0) return VxPicardEnum; + else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; + else if (strcmp(name,"Vy")==0) return VyEnum; else if (strcmp(name,"VyPicard")==0) return VyPicardEnum; else if (strcmp(name,"Vz")==0) return VzEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; - else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; - else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; else stage=6; } if(stage==6){ - if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; + if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; + else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; + else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum; @@ -750,10 +752,10 @@ else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum; - else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; - else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; + if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; + else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; + else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; else if (strcmp(name,"LevelsetReinitFrequency")==0) return LevelsetReinitFrequencyEnum; @@ -821,4 +823,5 @@ else if (strcmp(name,"Matice")==0) return MaticeEnum; else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; + else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; else if (strcmp(name,"Matestar")==0) return MatestarEnum; else if (strcmp(name,"Matpar")==0) return MatparEnum; @@ -872,11 +875,11 @@ else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum; else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; - else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum; - else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; - else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; + if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum; + else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; + else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; + else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum; else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 21388) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 21389) @@ -86,4 +86,5 @@ case 2: return MatestarEnum; case 3: return MaticeEnum; + case 4: return MatenhancediceEnum; default: _error_("Marshalled materials code \""<