Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 18735) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 18736) @@ -48,6 +48,8 @@ int finiteelement; + bool islevelset; iomodel->Constant(&finiteelement,DamageElementinterpEnum); + iomodel->Constant(&islevelset,TransientIslevelsetEnum); /*Update elements: */ @@ -68,4 +70,8 @@ iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); iomodel->FetchDataToInput(elements,PressureEnum); + + if(islevelset){ + iomodel->FetchDataToInput(elements,IceMaskNodeActivationEnum); + } }/*}}}*/ @@ -441,4 +447,9 @@ void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ + bool islevelset; + femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum); + if(islevelset){ + SetActiveNodesLSMx(femmodel); + } return; }/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 18735) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 18736) @@ -105,6 +105,6 @@ IssmDouble Jdet, dt, D_scalar; IssmDouble h,hx,hy,hz; - IssmDouble vel, calvingrate; - IssmDouble norm_dlsf; + IssmDouble vel; +// IssmDouble norm_dlsf; IssmDouble* xyz_list = NULL; @@ -130,5 +130,5 @@ IssmDouble* w = xNew(dim); IssmDouble* c = xNew(dim); - IssmDouble* dlsf = xNew(dim); + //IssmDouble* dlsf = xNew(dim); /*Retrieve all inputs and parameters*/ @@ -137,21 +137,28 @@ Input* vx_input=NULL; Input* vy_input=NULL; + Input* calvingratex_input=NULL; + Input* calvingratey_input=NULL; if(domaintype==Domain2DhorizontalEnum){ vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input); vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input); + calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input); + calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); } else{ if(dim==1){ vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input); + calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input); } if(dim==2){ vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input); vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input); - } - } - - Input* lsf_slopex_input = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); - Input* lsf_slopey_input = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); - Input* calvingrate_input = basalelement->GetInput(MasstransportCalvingrateEnum); _assert_(calvingrate_input); + calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input); + calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input); + } + } + +// Input* lsf_slopex_input = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input); +// Input* lsf_slopey_input = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input); + //Input* calvingrate_input = basalelement->GetInput(MasstransportCalvingrateEnum); _assert_(calvingrate_input); /* Start looping on the number of gaussian points: */ @@ -178,16 +185,18 @@ vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity vy_input->GetInputValue(&v[1],gauss); - lsf_slopex_input->GetInputValue(&dlsf[0],gauss); - lsf_slopey_input->GetInputValue(&dlsf[1],gauss); - calvingrate_input->GetInputValue(&calvingrate,gauss); - - norm_dlsf=0.; - for(i=0;i1.e-10) - for(i=0;iGetInputValue(&c[0],gauss); // in 3D case, add mesh velocity + calvingratey_input->GetInputValue(&c[1],gauss); + //lsf_slopex_input->GetInputValue(&dlsf[0],gauss); + //lsf_slopey_input->GetInputValue(&dlsf[1],gauss); + //calvingrate_input->GetInputValue(&calvingrate,gauss); + + //norm_dlsf=0.; + //for(i=0;i1.e-10) + // for(i=0;iElementSizes(&hx,&hy,&hz); - h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); + h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) ); kappa=h*vel/2.; for(row=0;rowElementSizes(&hx,&hy,&hz); - h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); + h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) ); for(row=0;row(w); xDelete(c); - xDelete(dlsf); + //xDelete(dlsf); delete gauss; if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18735) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18736) @@ -1210,4 +1210,12 @@ input=this->inputs->GetInput(output_enum); break; + case CalvingratexEnum: + case CalvingrateyEnum: + case MasstransportCalvingrateEnum: + this->StrainRateparallel(); + this->StrainRateperpendicular(); + this->CalvingRateLevermann(); + input=this->inputs->GetInput(output_enum); + break; case StrainRateparallelEnum: this->StrainRateparallel(); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18735) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 18736) @@ -221,5 +221,6 @@ virtual void StressIntensityFactor(void)=0; virtual void StrainRateparallel(void)=0; - virtual void StrainRateperpendicular(void)=0; + virtual void StrainRateperpendicular(void)=0; + virtual void CalvingRateLevermann(void)=0; virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18735) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18736) @@ -345,4 +345,58 @@ /*Clean up and return*/ delete gauss; +} +/*}}}*/ +void Penta::CalvingRateLevermann(){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + GaussPenta* gauss=NULL; + IssmDouble vx,vy,vel; + IssmDouble strainparallel; + IssmDouble propcoeff; + IssmDouble strainperpendicular; + 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* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input); + this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum); + + /* 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); + strainperpendicular_input->GetInputValue(&strainperpendicular,gauss); + + /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */ + calvingrate[iv]=propcoeff*strainparallel*strainperpendicular; + if(calvingrate[iv]<0){ + calvingrate[iv]=0; + } + calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6); + calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6); + } + + /*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(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; + } /*}}}*/ @@ -426,5 +480,6 @@ void Penta::StrainRateparallel(){/*{{{*/ - IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble *xyz_list = NULL; + IssmDouble epsilon[6]; GaussPenta* gauss=NULL; IssmDouble vx,vy,vel; @@ -435,12 +490,10 @@ /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - + this->GetVerticesCoordinates(&xyz_list); + /*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* strainxx_input=inputs->GetInput(StrainRatexxEnum); _assert_(strainxx_input); - Input* strainxy_input=inputs->GetInput(StrainRatexyEnum); _assert_(strainxy_input); - Input* strainyy_input=inputs->GetInput(StrainRateyyEnum); _assert_(strainyy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); /* Start looping on the number of vertices: */ @@ -453,7 +506,10 @@ vy_input->GetInputValue(&vy,gauss); vel=vx*vx+vy*vy; - strainxx_input->GetInputValue(&strainxx,gauss); - strainxy_input->GetInputValue(&strainxy,gauss); - strainyy_input->GetInputValue(&strainyy,gauss); + + /*Compute strain rate and viscosity: */ + this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); + strainxx=epsilon[0]; + strainyy=epsilon[1]; + strainxy=epsilon[3]; /*strainparallel= Strain rate along the ice flow direction */ @@ -466,9 +522,11 @@ /*Clean up and return*/ delete gauss; + xDelete(xyz_list); } /*}}}*/ void Penta::StrainRateperpendicular(){/*{{{*/ - IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble *xyz_list = NULL; + IssmDouble epsilon[6]; GaussPenta* gauss=NULL; IssmDouble vx,vy,vel; @@ -479,13 +537,10 @@ /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + this->GetVerticesCoordinates(&xyz_list); /*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* strainxx_input=inputs->GetInput(StrainRatexxEnum); _assert_(strainxx_input); - Input* strainxy_input=inputs->GetInput(StrainRatexyEnum); _assert_(strainxy_input); - Input* strainyy_input=inputs->GetInput(StrainRateyyEnum); _assert_(strainyy_input); - + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); /* Start looping on the number of vertices: */ @@ -498,7 +553,10 @@ vy_input->GetInputValue(&vy,gauss); vel=vx*vx+vy*vy; - strainxx_input->GetInputValue(&strainxx,gauss); - strainxy_input->GetInputValue(&strainxy,gauss); - strainyy_input->GetInputValue(&strainyy,gauss); + + /*Compute strain rate and viscosity: */ + this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input); + strainxx=epsilon[0]; + strainyy=epsilon[1]; + strainxy=epsilon[3]; /*strainperpendicular= Strain rate perpendicular to the ice flow direction */ @@ -511,4 +569,5 @@ /*Clean up and return*/ delete gauss; + xDelete(xyz_list); } /*}}}*/ @@ -2176,4 +2235,6 @@ if(this->inputs->GetInput(VxEnum)) this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); if(this->inputs->GetInput(VyEnum)) this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); + if(this->inputs->GetInput(CalvingratexEnum)) this->InputDepthAverageAtBase(CalvingratexEnum,CalvingratexAverageEnum); + if(this->inputs->GetInput(CalvingrateyEnum)) this->InputDepthAverageAtBase(CalvingrateyEnum,CalvingrateyAverageEnum); Tria* tria=(Tria*)SpawnTria(0,1,2); this->inputs->DeleteInput(MaterialsRheologyBbarEnum); @@ -2181,4 +2242,6 @@ this->inputs->DeleteInput(VxAverageEnum); this->inputs->DeleteInput(VyAverageEnum); + this->inputs->DeleteInput(CalvingratexAverageEnum); + this->inputs->DeleteInput(CalvingrateyAverageEnum); return tria; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18735) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 18736) @@ -60,4 +60,5 @@ void StrainRateparallel(); void StrainRateperpendicular(); + void CalvingRateLevermann(); 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 18735) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 18736) @@ -58,6 +58,7 @@ void ComputeDeviatoricStressTensor(){_error_("not implemented yet");}; void StressIntensityFactor(void){_error_("not implemented yet");}; - void StrainRateparallel(void){_error_("not implemented yet");}; - void StrainRateperpendicular(void){_error_("not implemented yet");}; + void StrainRateparallel(void){_error_("not implemented yet");}; + void StrainRateperpendicular(void){_error_("not implemented yet");}; + void CalvingRateLevermann(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 18735) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18736) @@ -60,4 +60,5 @@ void StrainRateparallel(void){_error_("not implemented yet");}; void StrainRateperpendicular(void){_error_("not implemented yet");}; + void CalvingRateLevermann(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.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18735) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18736) @@ -276,93 +276,4 @@ } /*}}}*/ -void Tria::StrainRateperpendicular(){/*{{{*/ - - IssmDouble xyz_list[NUMVERTICES][3]; - GaussPenta* gauss=NULL; - IssmDouble vx,vy,vel; - IssmDouble strainxx; - IssmDouble strainxy; - IssmDouble strainyy; - IssmDouble strainperpendicular[NUMVERTICES]; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*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* strainxx_input=inputs->GetInput(StrainRatexxEnum); _assert_(strainxx_input); - Input* strainxy_input=inputs->GetInput(StrainRatexyEnum); _assert_(strainxy_input); - Input* strainyy_input=inputs->GetInput(StrainRateyyEnum); _assert_(strainyy_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; - strainxx_input->GetInputValue(&strainxx,gauss); - strainxy_input->GetInputValue(&strainxy,gauss); - strainyy_input->GetInputValue(&strainyy,gauss); - - /*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); - } - - /*Add input*/ - this->inputs->AddInput(new PentaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum)); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ -void Tria::StrainRateparallel(){/*{{{*/ - - IssmDouble xyz_list[NUMVERTICES][3]; - GaussPenta* gauss=NULL; - IssmDouble vx,vy,vel; - IssmDouble strainxx; - IssmDouble strainxy; - IssmDouble strainyy; - IssmDouble strainparallel[NUMVERTICES]; - - /* Get node coordinates and dof list: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*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* strainxx_input=inputs->GetInput(StrainRatexxEnum); _assert_(strainxx_input); - Input* strainxy_input=inputs->GetInput(StrainRatexyEnum); _assert_(strainxy_input); - Input* strainyy_input=inputs->GetInput(StrainRateyyEnum); _assert_(strainyy_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; - strainxx_input->GetInputValue(&strainxx,gauss); - strainxy_input->GetInputValue(&strainxy,gauss); - strainyy_input->GetInputValue(&strainyy,gauss); - - /*strainparallel= Strain rate along the ice flow direction */ - strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6); - } - - /*Add input*/ - this->inputs->AddInput(new PentaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum)); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ void Tria::ComputeDeviatoricStressTensor(){/*{{{*/ @@ -413,4 +324,150 @@ /*Clean up and return*/ delete gauss; +} +/*}}}*/ +void Tria::StrainRateparallel(){/*{{{*/ + + IssmDouble *xyz_list = NULL; + IssmDouble epsilon[3]; + GaussTria* gauss=NULL; + IssmDouble vx,vy,vel; + IssmDouble strainxx; + IssmDouble strainxy; + IssmDouble strainyy; + IssmDouble strainparallel[NUMVERTICES]; + + /* Get node coordinates and dof list: */ + this->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will need*/ + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_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; + + /*Compute strain rate viscosity and pressure: */ + this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + strainxx=epsilon[0]; + strainyy=epsilon[1]; + strainxy=epsilon[2]; + + /*strainparallel= Strain rate along the ice flow direction */ + strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6); + } + + /*Add input*/ + this->inputs->AddInput(new TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); +} +/*}}}*/ +void Tria::StrainRateperpendicular(){/*{{{*/ + + IssmDouble *xyz_list = NULL; + GaussTria* gauss=NULL; + IssmDouble epsilon[3]; + IssmDouble vx,vy,vel; + IssmDouble strainxx; + IssmDouble strainxy; + IssmDouble strainyy; + IssmDouble strainperpendicular[NUMVERTICES]; + + /* Get node coordinates and dof list: */ + this->GetVerticesCoordinates(&xyz_list); + + /*Retrieve all inputs we will need*/ + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_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; + + /*Compute strain rate viscosity and pressure: */ + this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + strainxx=epsilon[0]; + strainyy=epsilon[1]; + strainxy=epsilon[2]; + + /*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); + } + + /*Add input*/ + this->inputs->AddInput(new TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; + xDelete(xyz_list); +} +/*}}}*/ +void Tria::CalvingRateLevermann(){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + GaussTria* gauss=NULL; + IssmDouble vx,vy,vel; + IssmDouble strainparallel; + IssmDouble propcoeff; + IssmDouble strainperpendicular; + 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* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input); + this->parameters->FindParam(&propcoeff,MasstransportLevermannCalvingCoeffEnum); + + /* 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); + strainperpendicular_input->GetInputValue(&strainperpendicular,gauss); + + /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */ + calvingrate[iv]=propcoeff*strainparallel*strainperpendicular; + if(calvingrate[iv]<0){ + calvingrate[iv]=0; + } + calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6); + calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6); + } + + /*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(MasstransportCalvingrateEnum,&calvingrate[0],P1Enum)); + + /*Clean up and return*/ + delete gauss; + } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18735) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 18736) @@ -55,8 +55,9 @@ void ComputeStressTensor(); void ComputeDeviatoricStressTensor(); + void StrainRateparallel(); + void StrainRateperpendicular(); void ComputeSurfaceNormalVelocity(); void StressIntensityFactor(void){_error_("not implemented yet");}; - void StrainRateparallel(); - void StrainRateperpendicular(); + void CalvingRateLevermann(); 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 18735) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18736) @@ -1695,4 +1695,28 @@ } /*}}}*/ +void FemModel::CalvingRateLevermannx(){/*{{{*/ + + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->CalvingRateLevermann(); + } +} +/*}}}*/ +void FemModel::StrainRateparallelx(){/*{{{*/ + + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->StrainRateparallel(); + } +} +/*}}}*/ +void FemModel::StrainRateperpendicularx(){/*{{{*/ + + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->StrainRateperpendicular(); + } +} +/*}}}*/ #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 18735) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18736) @@ -82,4 +82,7 @@ void BalancethicknessMisfitx(IssmDouble* pV); void StressIntensityFactorx(); + void StrainRateparallelx(); + void StrainRateperpendicularx(); + void CalvingRateLevermannx(); #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 18735) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 18736) @@ -21,5 +21,5 @@ /*parameters: */ IssmDouble starttime,finaltime,dt,yts; - bool isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology; + bool isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalvingrate; bool save_results,dakota_analysis; bool time_adapt=false; @@ -53,4 +53,5 @@ femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum); femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); + femmodel->parameters->FindParam(&iscalvingrate,MasstransportIscalvingrateEnum); if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum); femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum); @@ -101,4 +102,10 @@ if(islevelset){ + if(iscalvingrate){ + if(VerboseSolution()) _printf0_(" computing calving rate\n"); + femmodel->StrainRateparallelx(); + femmodel->StrainRateperpendicularx(); + femmodel->CalvingRateLevermannx(); + } 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/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 18735) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 18736) @@ -83,4 +83,6 @@ parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum)); parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum)); + parameters->AddObject(iomodel->CopyConstantObject(MasstransportIscalvingrateEnum)); + parameters->AddObject(iomodel->CopyConstantObject(MasstransportLevermannCalvingCoeffEnum)); parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum)); parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum)); Index: /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 18735) +++ /issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 18736) @@ -59,5 +59,5 @@ analysis_type==BalancethicknessAnalysisEnum || analysis_type==HydrologyDCInefficientAnalysisEnum || - analysis_type==DamageEvolutionAnalysisEnum || + //analysis_type==DamageEvolutionAnalysisEnum || analysis_type==HydrologyDCEfficientAnalysisEnum || analysis_type==LevelsetAnalysisEnum || Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18735) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18736) @@ -208,4 +208,8 @@ NewDamageEnum, StressIntensityFactorEnum, + CalvingratexEnum, + CalvingrateyEnum, + CalvingratexAverageEnum, + CalvingrateyAverageEnum, StrainRateparallelEnum, StrainRateperpendicularEnum, @@ -245,4 +249,6 @@ MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script) MasstransportHydrostaticAdjustmentEnum, + MasstransportIscalvingrateEnum, + MasstransportLevermannCalvingCoeffEnum, MasstransportIsfreesurfaceEnum, MasstransportMinThicknessEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18735) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18736) @@ -216,4 +216,8 @@ case NewDamageEnum : return "NewDamage"; case StressIntensityFactorEnum : return "StressIntensityFactor"; + case CalvingratexEnum : return "Calvingratex"; + case CalvingrateyEnum : return "Calvingratey"; + case CalvingratexAverageEnum : return "CalvingratexAverage"; + case CalvingrateyAverageEnum : return "CalvingrateyAverage"; case StrainRateparallelEnum : return "StrainRateparallel"; case StrainRateperpendicularEnum : return "StrainRateperpendicular"; @@ -253,4 +257,6 @@ case MiscellaneousNameEnum : return "MiscellaneousName"; case MasstransportHydrostaticAdjustmentEnum : return "MasstransportHydrostaticAdjustment"; + case MasstransportIscalvingrateEnum : return "MasstransportIscalvingrate"; + case MasstransportLevermannCalvingCoeffEnum : return "MasstransportLevermannCalvingCoeff"; case MasstransportIsfreesurfaceEnum : return "MasstransportIsfreesurface"; case MasstransportMinThicknessEnum : return "MasstransportMinThickness"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18735) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18736) @@ -219,4 +219,8 @@ else if (strcmp(name,"NewDamage")==0) return NewDamageEnum; else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum; + else 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; else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; @@ -256,12 +260,14 @@ else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum; else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum; + else stage=3; + } + if(stage==3){ + if (strcmp(name,"MasstransportIscalvingrate")==0) return MasstransportIscalvingrateEnum; + else if (strcmp(name,"MasstransportLevermannCalvingCoeff")==0) return MasstransportLevermannCalvingCoeffEnum; else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum; else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum; else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum; - else stage=3; - } - if(stage==3){ - if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum; + else if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum; else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum; else if (strcmp(name,"MasstransportVertexPairing")==0) return MasstransportVertexPairingEnum; @@ -377,5 +383,8 @@ else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum; else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum; - else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum; else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum; else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum; @@ -383,8 +392,5 @@ else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum; else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum; + else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum; else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; @@ -500,5 +506,8 @@ else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; else if (strcmp(name,"Tetra")==0) return TetraEnum; - else if (strcmp(name,"TetraInput")==0) return TetraInputEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"TetraInput")==0) return TetraInputEnum; else if (strcmp(name,"Penta")==0) return PentaEnum; else if (strcmp(name,"PentaInput")==0) return PentaInputEnum; @@ -506,8 +515,5 @@ else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum; else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"Air")==0) return AirEnum; + else if (strcmp(name,"Air")==0) return AirEnum; else if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"Melange")==0) return MelangeEnum; @@ -623,5 +629,8 @@ else if (strcmp(name,"P2xP4")==0) return P2xP4Enum; else if (strcmp(name,"P1P1")==0) return P1P1Enum; - else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum; else if (strcmp(name,"MINI")==0) return MINIEnum; else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum; @@ -629,8 +638,5 @@ else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum; + else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum; else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum; else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum; @@ -746,5 +752,8 @@ else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; - else if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum; else if (strcmp(name,"Seaiceocean")==0) return SeaiceoceanEnum; else if (strcmp(name,"SeaiceThickness")==0) return SeaiceThicknessEnum; @@ -752,8 +761,5 @@ else if (strcmp(name,"SeaiceMinConcentration")==0) return SeaiceMinConcentrationEnum; else if (strcmp(name,"SeaiceMinThickness")==0) return SeaiceMinThicknessEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"SeaiceMaxThickness")==0) return SeaiceMaxThicknessEnum; + else if (strcmp(name,"SeaiceMaxThickness")==0) return SeaiceMaxThicknessEnum; else if (strcmp(name,"SeaiceSpcvx")==0) return SeaiceSpcvxEnum; else if (strcmp(name,"SeaiceSpcvy")==0) return SeaiceSpcvyEnum;