Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26875) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26876) @@ -2873,77 +2873,40 @@ void Penta::MovingFrontalVelocity(void){/*{{{*/ if(!this->IsOnBase()) return; - int dim, domaintype, calvinglaw, i; + int domaintype, calvinglaw, i; IssmDouble v[3],w[3],c[3],m[3],dlsf[3]; IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice; IssmDouble migrationmax, calvinghaf, heaviside, haf_eps; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble movingfrontvx[NUMVERTICES]; - IssmDouble movingfrontvy[NUMVERTICES]; - IssmDouble vel; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble movingfrontvx[NUMVERTICES]; + IssmDouble movingfrontvy[NUMVERTICES]; + IssmDouble vel; + int dim=2; /* Get node coordinates and dof list: */ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vx_input = NULL; - Input* vy_input = NULL; Input* calvingratex_input = NULL; Input* calvingratey_input = NULL; - Input* lsf_slopex_input = NULL; - Input* lsf_slopey_input = NULL; Input* calvingrate_input = NULL; Input* meltingrate_input = NULL; - Input* gr_input = NULL; + + /*Load levelset function gradients*/ + Input *lsf_slopex_input = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); + Input *lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); + Input *vx_input = this->GetInput(VxEnum); _assert_(vx_input); + Input *vy_input = this->GetInput(VyEnum); _assert_(vy_input); + Input *gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); /*Get problem dimension and whether there is moving front or not*/ - this->FindParam(&domaintype,DomainTypeEnum); this->FindParam(&calvinglaw,CalvingLawEnum); - - switch(domaintype){ - case Domain2DverticalEnum: dim = 1; break; - case Domain2DhorizontalEnum: dim = 2; break; - case Domain3DEnum: dim = 2; break; - default: _error_("mesh "<GetInput(VxEnum); _assert_(vx_input); - break; - case Domain2DhorizontalEnum: - vx_input=this->GetInput(VxEnum); _assert_(vx_input); - vy_input=this->GetInput(VyEnum); _assert_(vy_input); - gr_input=this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); - break; - case Domain3DEnum: - vx_input=this->GetInput(VxAverageEnum); _assert_(vx_input); - vy_input=this->GetInput(VyAverageEnum); _assert_(vy_input); - gr_input=this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); - break; - default: _error_("mesh "<GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); - if(dim==2) lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); calvingrate_input = this->GetInput(CalvingCalvingrateEnum); _assert_(calvingrate_input); meltingrate_input = this->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); break; case CalvingLevermannEnum: - switch(domaintype){ - case Domain2DverticalEnum: - calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input); - break; - case Domain2DhorizontalEnum: - calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input); - calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); - break; - case Domain3DEnum: - calvingratex_input=this->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input); - calvingratey_input=this->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input); - break; - default: _error_("mesh "<GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input); + calvingratey_input=this->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input); meltingrate_input = this->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); break; @@ -2951,12 +2914,8 @@ case CalvingHabEnum: case CalvingCrevasseDepthEnum: - lsf_slopex_input = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); - if(dim==2) lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); meltingrate_input = this->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); break; case CalvingDev2Enum: this->FindParam(&calvinghaf,CalvingHeightAboveFloatationEnum); - lsf_slopex_input = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); - if(dim==2) lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); calvingrate_input = this->GetInput(CalvingCalvingrateEnum); _assert_(calvingrate_input); meltingrate_input = this->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); @@ -2975,4 +2934,7 @@ vy_input->GetInputValue(&v[1],&gauss); gr_input->GetInputValue(&groundedice,&gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + norm_dlsf=sqrt(dlsf[0]*dlsf[0] + dlsf[1]*dlsf[1]); /*Get calving speed*/ @@ -2980,14 +2942,8 @@ case DefaultCalvingEnum: case CalvingVonmisesEnum: - lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); calvingrate_input->GetInputValue(&calvingrate,&gauss); meltingrate_input->GetInputValue(&meltingrate,&gauss); if(groundedice<0) meltingrate = 0.; - - norm_dlsf=0.; - for(i=0;i1.e-10) @@ -3003,5 +2959,5 @@ case CalvingLevermannEnum: calvingratex_input->GetInputValue(&c[0],&gauss); - if(dim==2) calvingratey_input->GetInputValue(&c[1],&gauss); + calvingratey_input->GetInputValue(&c[1],&gauss); meltingrate_input->GetInputValue(&meltingrate,&gauss); norm_calving=0.; @@ -3014,12 +2970,6 @@ case CalvingHabEnum: case CalvingCrevasseDepthEnum: - lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); meltingrate_input->GetInputValue(&meltingrate,&gauss); if(groundedice<0) meltingrate = 0.; - - norm_dlsf=0.; - for(i=0;i1.e-10) @@ -3037,6 +2987,4 @@ case CalvingDev2Enum: { - lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); calvingrate_input->GetInputValue(&calvingrate,&gauss); meltingrate_input->GetInputValue(&meltingrate,&gauss); @@ -3062,8 +3010,4 @@ meltingrate=heaviside*meltingrate+0.; } - - norm_dlsf=0.; - for(i=0;i1.e-10) Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26875) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26876) @@ -4320,9 +4320,9 @@ Input* calvingratex_input = NULL; Input* calvingratey_input = NULL; - Input* lsf_slopex_input = NULL; Input* lsf_slopey_input = NULL; Input* calvingrate_input = NULL; Input* meltingrate_input = NULL; Input* gr_input = NULL; + /*Get problem dimension and whether there is moving front or not*/ @@ -4336,4 +4336,9 @@ default: _error_("mesh "<GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); + if(dim==2){lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);} + /*Load velocities*/ switch(domaintype){ @@ -4357,6 +4362,4 @@ case DefaultCalvingEnum: case CalvingVonmisesEnum: - lsf_slopex_input = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); - if(dim==2) lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); calvingrate_input = this->GetInput(CalvingCalvingrateEnum); _assert_(calvingrate_input); meltingrate_input = this->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); @@ -4382,12 +4385,7 @@ case CalvingHabEnum: case CalvingCrevasseDepthEnum: - lsf_slopex_input = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); - if(dim==2) lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); - meltingrate_input = this->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); break; case CalvingDev2Enum: this->FindParam(&calvinghaf,CalvingHeightAboveFloatationEnum); - lsf_slopex_input = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); - if(dim==2) lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); calvingrate_input = this->GetInput(CalvingCalvingrateEnum); _assert_(calvingrate_input); meltingrate_input = this->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); @@ -4398,6 +4396,4 @@ break; case CalvingParameterizationEnum: - lsf_slopex_input = this->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); - if(dim==2) lsf_slopey_input = this->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); calvingrate_input = this->GetInput(CalvingCalvingrateEnum); _assert_(calvingrate_input); meltingrate_input = this->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); @@ -4416,4 +4412,10 @@ vy_input->GetInputValue(&v[1],&gauss); gr_input->GetInputValue(&groundedice,&gauss); + lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); + if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); + + norm_dlsf=0.; + for(i=0;iGetInputValue(&dlsf[0],&gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); calvingrate_input->GetInputValue(&calvingrate,&gauss); meltingrate_input->GetInputValue(&meltingrate,&gauss); if(groundedice<0) meltingrate = 0.; - - norm_dlsf=0.; - for(i=0;i1.e-10) @@ -4454,11 +4450,5 @@ case CalvingHabEnum: case CalvingCrevasseDepthEnum: - lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); meltingrate_input->GetInputValue(&meltingrate,&gauss); - - norm_dlsf=0.; - for(i=0;i1.e-10) @@ -4476,6 +4466,4 @@ case CalvingDev2Enum: { - lsf_slopex_input->GetInputValue(&dlsf[0],&gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); calvingrate_input->GetInputValue(&calvingrate,&gauss); meltingrate_input->GetInputValue(&meltingrate,&gauss); @@ -4502,8 +4490,4 @@ } - norm_dlsf=0.; - for(i=0;i1.e-10) for(i=0;iGetInputValue(&dlsf[0],&gauss); - if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],&gauss); calvingrate_input->GetInputValue(&calvingrate,&gauss); meltingrate_input->GetInputValue(&meltingrate,&gauss); meltingrate = 0.; - - norm_dlsf=0.; - for(i=0;i1.e-10)