Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 18952) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 18953) @@ -63,4 +63,8 @@ iomodel->FetchDataToInput(elements,CalvinglevermannCoeffEnum); iomodel->FetchDataToInput(elements,CalvinglevermannMeltingrateEnum); + break; + case CalvingPiEnum: + iomodel->FetchDataToInput(elements,CalvingpiCoeffEnum); + iomodel->FetchDataToInput(elements,CalvingpiMeltingrateEnum); break; default: @@ -208,4 +212,20 @@ meltingrate_input = basalelement->GetInput(CalvinglevermannMeltingrateEnum); _assert_(meltingrate_input); break; + case CalvingPiEnum: + if(domaintype==Domain2DhorizontalEnum){ + calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input); + calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); + } + else{ + if(dim==1){ + calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input); + } + if(dim==2){ + calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input); + calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input); + } + } + meltingrate_input = basalelement->GetInput(CalvingpiMeltingrateEnum); _assert_(meltingrate_input); + break; default: _error_("Calving law "<GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity @@ -274,4 +295,18 @@ break; + + case CalvingPiEnum: + calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity + if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss); + meltingrate_input->GetInputValue(&meltingrate,gauss); + + norm_calving=0.; + for(i=0;i(w); xDelete(c); - xDelete(m); xDelete(dlsf); delete gauss; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18952) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18953) @@ -1369,5 +1369,16 @@ this->StrainRateparallel(); this->StrainRateperpendicular(); - this->CalvingRateLevermann(); + int calvinglaw; + this->FindParam(&calvinglaw,CalvingLawEnum); + switch(calvinglaw){ + case CalvingLevermannEnum: + this->CalvingRateLevermann(); + break; + case CalvingPiEnum: + this->CalvingRatePi(); + break; + default: + _error_("Calving law "<inputs->GetInput(output_enum); break; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18952) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18953) @@ -170,4 +170,5 @@ virtual void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0; virtual void CalvingRateLevermann(void)=0; + virtual void CalvingRatePi(void)=0; virtual IssmDouble CharacteristicLength(void)=0; virtual void ComputeBasalStress(Vector* sigma_b)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18952) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18953) @@ -201,6 +201,81 @@ calvingrate[iv]=0; } - calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6); - calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6); + calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14); + calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14); + } + + /*Add input*/ + this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum)); + this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum)); + this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; + +} +/*}}}*/ +void Penta::CalvingRatePi(){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + GaussPenta* gauss=NULL; + IssmDouble vx,vy,vel; + IssmDouble strainparallel; + IssmDouble sxx; + IssmDouble sxy; + IssmDouble syy; + IssmDouble sigVM; + IssmDouble thickness; + IssmDouble base; + IssmDouble hAB; + IssmDouble propcoeff; + IssmDouble calvingratex[NUMVERTICES]; + IssmDouble calvingratey[NUMVERTICES]; + IssmDouble calvingrate[NUMVERTICES]; + + + /* Get node coordinates and dof list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Retrieve all inputs and parameters we will need*/ + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum); _assert_(strainparallel_input); + Input* sxx_input=inputs->GetInput(DeviatoricStressxxEnum); _assert_(sxx_input); + Input* sxy_input=inputs->GetInput(DeviatoricStressxyEnum); _assert_(sxy_input); + Input* syy_input=inputs->GetInput(DeviatoricStressyyEnum); _assert_(syy_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* base_input=inputs->GetInput(BaseEnum); _assert_(base_input); + Input* picoeff_input=inputs->GetInput(CalvingpiCoeffEnum); _assert_(picoeff_input); + + /* Start looping on the number of vertices: */ + gauss=new GaussPenta(); + for (int iv=0;ivGaussVertex(iv); + + /* Get the value we need*/ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel=vx*vx+vy*vy; + strainparallel_input->GetInputValue(&strainparallel,gauss); + sxx_input->GetInputValue(&sxx,gauss); + sxy_input->GetInputValue(&sxy,gauss); + syy_input->GetInputValue(&syy,gauss); + thickness_input->GetInputValue(&thickness,gauss); + base_input->GetInputValue(&base,gauss); + picoeff_input->GetInputValue(&propcoeff,gauss); + + /* Computing sigma Von Mises*/ + sigVM=sqrt(sxx*sxx+syy*syy+3*sxy*sxy-sxx*syy); + + /* Computing heigth above buoyancy*/ + hAB=thickness+1028/920*base; + + /*Calving rate for Pi criterion proportionnal to the product of the strain rate along the ice flow direction and the Von Mises stress and the square of the glacier width (hardcoded) divided by the height above buoyancy and the max of the ice velocity power 3 and the ice density (ignored here)*/ + calvingrate[iv]=propcoeff*strainparallel*sigVM*25.e6/hAB/1e9; + if(calvingrate[iv]<0){ + calvingrate[iv]=0; + } + calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-6); + calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-6); } @@ -2433,5 +2508,5 @@ /*strainparallel= Strain rate along the ice flow direction */ - strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6); + strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-14); } @@ -2480,5 +2555,5 @@ /*strainperpendicular= Strain rate perpendicular to the ice flow direction */ - strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6); + strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-14); } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18952) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18953) @@ -49,4 +49,5 @@ IssmDouble CharacteristicLength(void){_error_("not implemented yet");}; void CalvingRateLevermann(); + void CalvingRatePi(); void ComputeBasalStress(Vector* sigma_b); void ComputeDeviatoricStressTensor(); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18952) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18953) @@ -49,4 +49,5 @@ void ComputeSigmaNN(){_error_("not implemented yet");}; void ComputeStressTensor(){_error_("not implemented yet");}; + void CalvingRatePi(void){_error_("not implemented yet");}; void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18952) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18953) @@ -49,4 +49,5 @@ void ComputeStressTensor(){_error_("not implemented yet");}; void ComputeDeviatoricStressTensor(){_error_("not implemented yet");}; + void CalvingRatePi(void){_error_("not implemented yet");}; void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18952) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18953) @@ -234,6 +234,84 @@ calvingrate[iv]=0; } - calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6); - calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6); + calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14); + calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14); + } + + /*Add input*/ + this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum)); + this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum)); + this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; + +} +/*}}}*/ +void Tria::CalvingRatePi(){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + GaussTria* gauss=NULL; + IssmDouble vx,vy,vel; + IssmDouble strainparallel; + IssmDouble sxx; + IssmDouble sxy; + IssmDouble syy; + IssmDouble sigVM; + IssmDouble thickness; + IssmDouble base; + IssmDouble hAB; + IssmDouble propcoeff; + IssmDouble calvingratex[NUMVERTICES]; + IssmDouble calvingratey[NUMVERTICES]; + IssmDouble calvingrate[NUMVERTICES]; + + + /* Get node coordinates and dof list: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*Retrieve all inputs and parameters we will need*/ + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum); _assert_(strainparallel_input); + Input* sxx_input=inputs->GetInput(DeviatoricStressxxEnum); _assert_(sxx_input); + Input* sxy_input=inputs->GetInput(DeviatoricStressxyEnum); _assert_(sxy_input); + Input* syy_input=inputs->GetInput(DeviatoricStressyyEnum); _assert_(syy_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* base_input=inputs->GetInput(BaseEnum); _assert_(base_input); + Input* picoeff_input=inputs->GetInput(CalvingpiCoeffEnum); _assert_(picoeff_input); + + + + /* Start looping on the number of vertices: */ + gauss=new GaussTria(); + for (int iv=0;ivGaussVertex(iv); + + /* Get the value we need*/ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel=vx*vx+vy*vy; + strainparallel_input->GetInputValue(&strainparallel,gauss); + sxx_input->GetInputValue(&sxx,gauss); + sxy_input->GetInputValue(&sxy,gauss); + syy_input->GetInputValue(&syy,gauss); + thickness_input->GetInputValue(&thickness,gauss); + base_input->GetInputValue(&base,gauss); + picoeff_input->GetInputValue(&propcoeff,gauss); + + /* Computing sigma Von Mises*/ + sigVM=sqrt(sxx*sxx+syy*syy+3*sxy*sxy-sxx*syy); + + /* Computing heigth above buoyancy*/ + hAB=thickness+1028/920*base; + + /*Calving rate for Pi criterion proportionnal to the product of the strain rate along the ice flow direction and the Von Mises stress and the square of the glacier width (hardcoded) divided by the height above buoyancy and the max of the ice velocity power 3 and the ice density (ignored here)*/ + + calvingrate[iv]=propcoeff*strainparallel*sigVM*25.e6/hAB/1e9; + if(calvingrate[iv]<0){ + calvingrate[iv]=0; + } + calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14); + calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14); } @@ -2584,5 +2662,5 @@ /*strainparallel= Strain rate along the ice flow direction */ - strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6); + strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-14); } @@ -2630,5 +2708,5 @@ /*strainperpendicular= Strain rate perpendicular to the ice flow direction */ - strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6); + strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-14); } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18952) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18953) @@ -58,4 +58,5 @@ void ComputeStressTensor(); void ComputeSurfaceNormalVelocity(); + void CalvingRatePi(); void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18952) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18953) @@ -1735,4 +1735,12 @@ } /*}}}*/ +void FemModel::CalvingRatePix(){/*{{{*/ + + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->CalvingRatePi(); + } +} +/*}}}*/ void FemModel::StrainRateparallelx(){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18952) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18953) @@ -86,4 +86,5 @@ void StrainRateperpendicularx(); void CalvingRateLevermannx(); + void CalvingRatePix(); #ifdef _HAVE_DAKOTA_ void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses); Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 18952) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 18953) @@ -109,4 +109,9 @@ femmodel->CalvingRateLevermannx(); } + if(iscalving && calvinglaw==CalvingPiEnum){ + if(VerboseSolution()) _printf0_(" computing calving rate\n"); + femmodel->StrainRateparallelx(); + femmodel->CalvingRatePix(); + } if(VerboseSolution()) _printf0_(" computing movement of ice boundaries\n"); /* smoothen slope of lsf for computation of normal on ice domain*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18952) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18953) @@ -215,8 +215,11 @@ CalvingMeltingrateEnum, CalvingLevermannEnum, + CalvingPiEnum, DefaultCalvingEnum, CalvingRequestedOutputsEnum, CalvinglevermannCoeffEnum, CalvinglevermannMeltingrateEnum, + CalvingpiCoeffEnum, + CalvingpiMeltingrateEnum, CalvingratexEnum, CalvingrateyEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18952) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18953) @@ -223,8 +223,11 @@ case CalvingMeltingrateEnum : return "CalvingMeltingrate"; case CalvingLevermannEnum : return "CalvingLevermann"; + case CalvingPiEnum : return "CalvingPi"; case DefaultCalvingEnum : return "DefaultCalving"; case CalvingRequestedOutputsEnum : return "CalvingRequestedOutputs"; case CalvinglevermannCoeffEnum : return "CalvinglevermannCoeff"; case CalvinglevermannMeltingrateEnum : return "CalvinglevermannMeltingrate"; + case CalvingpiCoeffEnum : return "CalvingpiCoeff"; + case CalvingpiMeltingrateEnum : return "CalvingpiMeltingrate"; case CalvingratexEnum : return "Calvingratex"; case CalvingrateyEnum : return "Calvingratey"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18952) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18953) @@ -226,8 +226,11 @@ else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum; + else if (strcmp(name,"CalvingPi")==0) return CalvingPiEnum; else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; else if (strcmp(name,"CalvingRequestedOutputs")==0) return CalvingRequestedOutputsEnum; else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum; + else if (strcmp(name,"CalvingpiCoeff")==0) return CalvingpiCoeffEnum; + else if (strcmp(name,"CalvingpiMeltingrate")==0) return CalvingpiMeltingrateEnum; else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum; else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum; @@ -257,11 +260,11 @@ else if (strcmp(name,"MeshNumberofvertices")==0) return MeshNumberofverticesEnum; else if (strcmp(name,"MeshUpperelements")==0) return MeshUpperelementsEnum; - else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum; - else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum; - else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"MeshX")==0) return MeshXEnum; + if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum; + else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum; + else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum; + else if (strcmp(name,"MeshX")==0) return MeshXEnum; else if (strcmp(name,"MeshY")==0) return MeshYEnum; else if (strcmp(name,"MeshZ")==0) return MeshZEnum; @@ -380,11 +383,11 @@ else if (strcmp(name,"SurfaceforcingsMelt")==0) return SurfaceforcingsMeltEnum; else if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum; - else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; - else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum; - else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum; + if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; + else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum; + else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum; + else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum; else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum; else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; @@ -503,11 +506,11 @@ else if (strcmp(name,"MassconName")==0) return MassconNameEnum; else if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum; - else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum; - else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; - else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum; + if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum; + else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; + else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum; + else if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum; else if (strcmp(name,"MassconaxpbyNamex")==0) return MassconaxpbyNamexEnum; else if (strcmp(name,"MassconaxpbyNamey")==0) return MassconaxpbyNameyEnum; @@ -626,11 +629,11 @@ else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; - else if (strcmp(name,"StrainRate")==0) return StrainRateEnum; - else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; - else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; + if (strcmp(name,"StrainRate")==0) return StrainRateEnum; + else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; + else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; + else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum; @@ -749,11 +752,11 @@ else if (strcmp(name,"Index")==0) return IndexEnum; else if (strcmp(name,"Indexed")==0) return IndexedEnum; - else if (strcmp(name,"Intersect")==0) return IntersectEnum; - else if (strcmp(name,"Nodal")==0) return NodalEnum; - else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; + if (strcmp(name,"Intersect")==0) return IntersectEnum; + else if (strcmp(name,"Nodal")==0) return NodalEnum; + else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; + else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; else if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum; else if (strcmp(name,"RootPath")==0) return RootPathEnum; Index: /issm/trunk-jpl/src/m/classes/calvingpi.m =================================================================== --- /issm/trunk-jpl/src/m/classes/calvingpi.m (revision 18953) +++ /issm/trunk-jpl/src/m/classes/calvingpi.m (revision 18953) @@ -0,0 +1,53 @@ +%CALVINGPI class definition +% +% Usage: +% calvingpi=calvingpi(); + +classdef calvingpi + properties (SetAccess=public) + coeff = NaN; + meltingrate = NaN; + end + methods + function obj = calvingpi(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + case 1 + inputstruct=varargin{1}; + list1 = properties('calvingpi'); + list2 = fieldnames(inputstruct); + for i=1:length(list1) + fieldname = list1{i}; + if ismember(fieldname,list2), + obj.(fieldname) = inputstruct.(fieldname); + end + end + otherwise + error('constructor not supported'); + end + end % }}} + function obj = setdefaultparameters(obj) % {{{ + %Proportionality coefficient in Pi model + obj.coeff=2e13; + end % }}} + function md = checkconsistency(obj,md,solution,analyses) % {{{ + %Early return + if (solution~=TransientSolutionEnum() | md.transient.iscalving==0), return; end + md = checkfield(md,'fieldname','calving.coeff','>',0,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0); + end % }}} + function disp(obj) % {{{ + disp(sprintf(' Calving Pi parameters:')); + fielddisplay(obj,'coeff','proportionality coefficient in Pi model'); + fielddisplay(obj,'meltingrate','melting rate at given location [m/a]'); + + end % }}} + function marshall(obj,md,fid) % {{{ + yts=365.0*24.0*3600.0; + WriteData(fid,'enum',CalvingLawEnum(),'data',CalvingPiEnum(),'format','Integer'); + WriteData(fid,'enum',CalvingpiCoeffEnum(),'data',obj.coeff,'format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts); + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 18952) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 18953) @@ -836,4 +836,5 @@ if isa(md.calving,'calving'),md.calving.calvingrate=project3d(md,'vector',md.calving.calvingrate,'type','node');end; if isa(md.calving,'calvinglevermann'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end; + if isa(md.calving,'calvingpi'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end; % Hydrologydc variables Index: /issm/trunk-jpl/src/m/enum/CalvingPiEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/CalvingPiEnum.m (revision 18953) +++ /issm/trunk-jpl/src/m/enum/CalvingPiEnum.m (revision 18953) @@ -0,0 +1,11 @@ +function macro=CalvingPiEnum() +%CALVINGPIENUM - Enum of CalvingPi +% +% 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=CalvingPiEnum() + +macro=StringToEnum('CalvingPi'); Index: /issm/trunk-jpl/src/m/enum/CalvingpiCoeffEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/CalvingpiCoeffEnum.m (revision 18953) +++ /issm/trunk-jpl/src/m/enum/CalvingpiCoeffEnum.m (revision 18953) @@ -0,0 +1,11 @@ +function macro=CalvingpiCoeffEnum() +%CALVINGPICOEFFENUM - Enum of CalvingpiCoeff +% +% 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=CalvingpiCoeffEnum() + +macro=StringToEnum('CalvingpiCoeff'); Index: /issm/trunk-jpl/src/m/enum/CalvingpiMeltingrateEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/CalvingpiMeltingrateEnum.m (revision 18953) +++ /issm/trunk-jpl/src/m/enum/CalvingpiMeltingrateEnum.m (revision 18953) @@ -0,0 +1,11 @@ +function macro=CalvingpiMeltingrateEnum() +%CALVINGPIMELTINGRATEENUM - Enum of CalvingpiMeltingrate +% +% 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=CalvingpiMeltingrateEnum() + +macro=StringToEnum('CalvingpiMeltingrate'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18952) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18953) @@ -215,8 +215,11 @@ def CalvingMeltingrateEnum(): return StringToEnum("CalvingMeltingrate")[0] def CalvingLevermannEnum(): return StringToEnum("CalvingLevermann")[0] +def CalvingPiEnum(): return StringToEnum("CalvingPi")[0] def DefaultCalvingEnum(): return StringToEnum("DefaultCalving")[0] def CalvingRequestedOutputsEnum(): return StringToEnum("CalvingRequestedOutputs")[0] def CalvinglevermannCoeffEnum(): return StringToEnum("CalvinglevermannCoeff")[0] def CalvinglevermannMeltingrateEnum(): return StringToEnum("CalvinglevermannMeltingrate")[0] +def CalvingpiCoeffEnum(): return StringToEnum("CalvingpiCoeff")[0] +def CalvingpiMeltingrateEnum(): return StringToEnum("CalvingpiMeltingrate")[0] def CalvingratexEnum(): return StringToEnum("Calvingratex")[0] def CalvingrateyEnum(): return StringToEnum("Calvingratey")[0]