Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 19032) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 19033) @@ -68,4 +68,8 @@ iomodel->FetchDataToInput(elements,CalvingpiMeltingrateEnum); break; + case CalvingDevEnum: + iomodel->FetchDataToInput(elements,CalvingpiCoeffEnum); + iomodel->FetchDataToInput(elements,CalvingMeltingrateEnum); + break; default: _error_("Calving law "<GetInput(VxEnum); _assert_(vx_input); - vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input); - } - else{ - if(dim==1){ + switch(domaintype){ + case Domain2DverticalEnum: vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input); - } - if(dim==2){ + break; + case Domain2DhorizontalEnum: + vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input); + vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input); + break; + case Domain3DEnum: vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input); vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input); - } + break; + default: _error_("mesh "<GetInput(CalvingratexEnum); _assert_(calvingratex_input); - calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); - } - else{ - if(dim==1){ + switch(domaintype){ + case Domain2DverticalEnum: calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input); - } - if(dim==2){ + break; + case Domain2DhorizontalEnum: + calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input); + calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); + break; + case Domain3DEnum: calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input); calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input); - } + break; + default: _error_("mesh "<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){ + switch(domaintype){ + case Domain2DverticalEnum: calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input); - } - if(dim==2){ + break; + case Domain2DhorizontalEnum: + calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input); + calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); + break; + case Domain3DEnum: calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input); calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input); - } + break; + default: _error_("mesh "<GetInput(CalvingpiMeltingrateEnum); _assert_(meltingrate_input); + break; + case CalvingDevEnum: + switch(domaintype){ + case Domain2DverticalEnum: + calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input); + break; + case Domain2DhorizontalEnum: + calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input); + calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); + break; + case Domain3DEnum: + calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input); + calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input); + break; + default: _error_("mesh "<GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input); break; default: @@ -255,5 +279,5 @@ GetB(B,basalelement,xyz_list,gauss); GetBprime(Bprime,basalelement,xyz_list,gauss); - vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity + vx_input->GetInputValue(&v[0],gauss); vy_input->GetInputValue(&v[1],gauss); @@ -273,39 +297,22 @@ if(norm_dlsf>1.e-10) for(i=0;iGetInputValue(&c[0],gauss); // in 3D case, add mesh velocity + case CalvingPiEnum: + case CalvingDevEnum: + calvingratex_input->GetInputValue(&c[0],gauss); if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss); meltingrate_input->GetInputValue(&meltingrate,gauss); - norm_calving=0.; for(i=0;iGetInputValue(&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;iNodalFunctions(basis, gauss); - thickness_input->GetInputValue(&thickness,gauss); + thickness_input->GetInputValue(&thickness,gauss); _assert_(thickness>0); surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19032) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19033) @@ -171,4 +171,5 @@ virtual void CalvingRateLevermann(void)=0; virtual void CalvingRatePi(void)=0; + virtual void CalvingRateDev(void){_error_("not implemented yet");}; virtual IssmDouble CharacteristicLength(void)=0; virtual void ComputeBasalStress(Vector* sigma_b)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 19032) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 19033) @@ -326,4 +326,61 @@ delete gauss; +} +/*}}}*/ +void Tria::CalvingRateDev(){/*{{{*/ + + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + IssmDouble calvingratex[NUMVERTICES]; + IssmDouble calvingratey[NUMVERTICES]; + IssmDouble calvingrate[NUMVERTICES]; + IssmDouble lambda1,lambda2,ex,ey,vx,vy,vel; + + /* 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); + + /* Start looping on the number of vertices: */ + GaussTria* gauss=new GaussTria(); + for(int iv=0;ivGaussVertex(iv); + + /*Get velocity components and thickness*/ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel=sqrt(vx*vx+vy*vy)+1.e-14; + + /*Compute strain rate and viscosity: */ + this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + + /*Get Eigen values*/ + Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]); + + /*Process Eigen values*/ + lambda1>0? lambda1 = pow(lambda1,.3) : lambda1=0.; + lambda2>0? lambda2 = pow(lambda2,.3) : lambda2=0.; + lambda1 = lambda1*5.e-2; + lambda2 = lambda2*5.e-2; + _assert_(!xIsNan(lambda1)); + _assert_(!xIsNan(lambda2)); + + /*Assign values*/ + //calvingratex[iv]=ex*lambda1 - ey*lambda2; + //calvingratey[iv]=ey*lambda1 + ex*lambda2; + calvingratex[iv]=vx/vel*(lambda1 + lambda2); + calvingratey[iv]=vy/vel*(lambda1 + lambda2); + calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]); + } + + /*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; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 19032) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 19033) @@ -52,4 +52,6 @@ void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); void CalvingRateLevermann(); + void CalvingRatePi(); + void CalvingRateDev(); IssmDouble CharacteristicLength(void); void ComputeBasalStress(Vector* sigma_b); @@ -58,5 +60,4 @@ 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 19032) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19033) @@ -669,4 +669,13 @@ /*Now, update degrees of freedoms: */ NodesDofx(nodes,parameters,analysis_type); + +} +/*}}}*/ +void FemModel::ElementOperationx(void (Element::*function)(void)){ /*{{{*/ + + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + (element->*function)(); + } } @@ -1743,4 +1752,12 @@ } /*}}}*/ +void FemModel::CalvingRateDevx(){/*{{{*/ + + for(int i=0;iSize();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + element->CalvingRateDev(); + } +} +/*}}}*/ void FemModel::StrainRateparallelx(){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 19032) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 19033) @@ -62,4 +62,5 @@ /*Modules*/ + void ElementOperationx(void (Element::*function)(void)); void GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,IssmDouble* ug); void MassFluxx(IssmDouble* presponse); @@ -88,4 +89,5 @@ void CalvingRateLevermannx(); void CalvingRatePix(); + void CalvingRateDevx(); #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 19032) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 19033) @@ -103,17 +103,29 @@ if(islevelset){ - if(iscalving && calvinglaw==CalvingLevermannEnum){ - if(VerboseSolution()) _printf0_(" computing calving rate\n"); - femmodel->StrainRateparallelx(); - femmodel->StrainRateperpendicularx(); - femmodel->CalvingRateLevermannx(); + if(iscalving){ + switch(calvinglaw){ + case DefaultCalvingEnum: + break; + case CalvingLevermannEnum: + if(VerboseSolution()) _printf0_(" computing Levermann's calving rate\n"); + femmodel->StrainRateparallelx(); + femmodel->StrainRateperpendicularx(); + femmodel->CalvingRateLevermannx(); + break; + case CalvingPiEnum: + if(VerboseSolution()) _printf0_(" computing Pi's calving rate\n"); + femmodel->StrainRateparallelx(); + femmodel->DeviatoricStressx(); + femmodel->CalvingRatePix(); + break; + case CalvingDevEnum: + femmodel->CalvingRateDevx(); + femmodel->ElementOperationx(&Element::CalvingRateDev); + break; + default: + _error_("Caving law "<StrainRateparallelx(); - femmodel->DeviatoricStressx(); - femmodel->CalvingRatePix(); - } - if(VerboseSolution()) _printf0_(" computing movement of ice boundaries\n"); + if(VerboseSolution()) _printf0_(" computing levelset transport\n"); /* smoothen slope of lsf for computation of normal on ice domain*/ levelsetfunctionslope_core(femmodel);