Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23887) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23888) @@ -1545,6 +1545,22 @@ *pM=total_mass_flux; }/*}}}*/ +void FemModel::GroundinglineMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/ + + IssmDouble local_mass_flux = 0; + IssmDouble total_mass_flux; + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + local_mass_flux+=element->GroundinglineMassFlux(scaled); + } + ISSM_MPI_Reduce(&local_mass_flux,&total_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_mass_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pM=total_mass_flux; + +}/*}}}*/ void FemModel::IceMassx(IssmDouble* pM, bool scaled){/*{{{*/ IssmDouble local_ice_mass = 0; @@ -2205,6 +2221,7 @@ case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(&double_result,true); break; case GroundedAreaEnum: this->GroundedAreax(&double_result,false); break; case GroundedAreaScaledEnum: this->GroundedAreax(&double_result,true); break; + case GroundinglineMassFluxEnum: this->GroundinglineMassFluxx(&double_result,false); break; case FloatingAreaEnum: this->FloatingAreax(&double_result,false); break; case FloatingAreaScaledEnum: this->FloatingAreax(&double_result,true); break; case MinVelEnum: this->MinVelx(&double_result); break; @@ -2402,7 +2419,7 @@ case IceVolumeScaledEnum: this->IceVolumex(responses, true); break; case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(responses, false); break; case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(responses, true); break; - case IcefrontMassFluxEnum: this->IcefrontMassFluxx(responses, true); break; + case IcefrontMassFluxEnum: this->IcefrontMassFluxx(responses, false); break; case GroundedAreaEnum: this->GroundedAreax(responses, false); break; case GroundedAreaScaledEnum: this->GroundedAreax(responses, true); break; case FloatingAreaEnum: this->FloatingAreax(responses, false); break; Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 23887) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 23888) @@ -94,6 +94,8 @@ bool HasEdgeOnSurface(); IssmDouble IceVolume(bool scaled); IssmDouble IceVolumeAboveFloatation(bool scaled); + IssmDouble IcefrontMassFlux(bool scaled); + IssmDouble GroundinglineMassFlux(bool scaled); void Ismip6FloatingiceMeltingRate(void); void InputDepthAverageAtBase(int enum_type,int average_enum_type); void InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/}; Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23887) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23888) @@ -1079,6 +1079,9 @@ IssmDouble* H=NULL;; int nrfrontbed,numiceverts; + if(!this->IsOnBase()) return 0; + if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0; + /*Retrieve all inputs and parameters*/ GetInputListOnVertices(&bed[0],BedEnum); GetInputListOnVertices(&surfaces[0],SurfaceEnum); @@ -1085,9 +1088,6 @@ GetInputListOnVertices(&bases[0],BaseEnum); GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum); - if(!this->IsOnBase()) return 0; - if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0; - nrfrontbed=0; for(int i=0;iGetVerticesCoordinates(&xyz_list); + + /*Recover parameters and values*/ + parameters->FindParam(&domaintype,DomainTypeEnum); + GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum); + + /*Be sure that values are not zero*/ + if(gl[0]==0.) gl[0]=gl[0]+epsilon; + if(gl[1]==0.) gl[1]=gl[1]+epsilon; + if(gl[2]==0.) gl[2]=gl[2]+epsilon; + + if(domaintype==Domain2DverticalEnum){ + _error_("not implemented"); + } + else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){ + int pt1 = 0; + int pt2 = 1; + if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + + /*Portion of the segments*/ + s1=gl[2]/(gl[2]-gl[1]); + s2=gl[2]/(gl[2]-gl[0]); + if(gl[2]<0.){ + pt1 = 1; pt2 = 0; + } + xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]); + xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]); + xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]); + xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]); + xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]); + xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]); + } + else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 + + /*Portion of the segments*/ + s1=gl[0]/(gl[0]-gl[1]); + s2=gl[0]/(gl[0]-gl[2]); + if(gl[0]<0.){ + pt1 = 1; pt2 = 0; + } + + xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]); + xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]); + xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]); + xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]); + xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]); + xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]); + } + else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 + + /*Portion of the segments*/ + s1=gl[1]/(gl[1]-gl[0]); + s2=gl[1]/(gl[1]-gl[2]); + if(gl[1]<0.){ + pt1 = 1; pt2 = 0; + } + + xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]); + xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]); + xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]); + xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]); + xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]); + xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]); + } + else{ + _error_("case not possible"); + } + + } + else _error_("mesh type "<=0 && s1<=1.); + _assert_(s2>=0 && s2<=1.); + + /*Get normal vector*/ + IssmDouble normal[3]; + this->NormalSection(&normal[0],&xyz_front[0][0]); + normal[0] = -normal[0]; + normal[1] = -normal[1]; + + /*Get inputs*/ + IssmDouble flux = 0.; + IssmDouble vx,vy,thickness,Jdet; + IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input=NULL; + Input* vy_input=NULL; + if(domaintype==Domain2DhorizontalEnum){ + vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + } + else{ + vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input); + vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input); + } + + /*Start looping on Gaussian points*/ + Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + thickness_input->GetInputValue(&thickness,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss); + + flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]); + } + + return flux; +} +/*}}}*/ +IssmDouble Tria::GroundinglineMassFlux(bool scaled){/*{{{*/ + + /*Make sure there is an ice front here*/ + if(!IsIceInElement()) return 0; + if(!IsZeroLevelset(MaskGroundediceLevelsetEnum)) return 0; + + /*Scaled not implemented yet...*/ + _assert_(!scaled); + + int domaintype,index1,index2; + const IssmPDouble epsilon = 1.e-15; + IssmDouble s1,s2; + IssmDouble gl[NUMVERTICES]; + IssmDouble xyz_front[2][3]; + + IssmDouble *xyz_list = NULL; + this->GetVerticesCoordinates(&xyz_list); + + /*Recover parameters and values*/ + parameters->FindParam(&domaintype,DomainTypeEnum); + GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum); + + /*Be sure that values are not zero*/ + if(gl[0]==0.) gl[0]=gl[0]+epsilon; + if(gl[1]==0.) gl[1]=gl[1]+epsilon; + if(gl[2]==0.) gl[2]=gl[2]+epsilon; + + if(domaintype==Domain2DverticalEnum){ + _error_("not implemented"); + } + else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){ + int pt1 = 0; + int pt2 = 1; + if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + + /*Portion of the segments*/ + s1=gl[2]/(gl[2]-gl[1]); + s2=gl[2]/(gl[2]-gl[0]); + if(gl[2]<0.){ + pt1 = 1; pt2 = 0; + } + xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]); + xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]); + xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]); + xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]); + xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]); + xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]); + } + else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 + + /*Portion of the segments*/ + s1=gl[0]/(gl[0]-gl[1]); + s2=gl[0]/(gl[0]-gl[2]); + if(gl[0]<0.){ + pt1 = 1; pt2 = 0; + } + + xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]); + xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]); + xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]); + xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]); + xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]); + xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]); + } + else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 + + /*Portion of the segments*/ + s1=gl[1]/(gl[1]-gl[0]); + s2=gl[1]/(gl[1]-gl[2]); + if(gl[1]<0.){ + pt1 = 1; pt2 = 0; + } + + xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]); + xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]); + xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]); + xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]); + xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]); + xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]); + } + else{ + _error_("case not possible"); + } + + } + else _error_("mesh type "<=0 && s1<=1.); + _assert_(s2>=0 && s2<=1.); + + /*Get normal vector*/ + IssmDouble normal[3]; + this->NormalSection(&normal[0],&xyz_front[0][0]); + + /*Get inputs*/ + IssmDouble flux = 0.; + IssmDouble vx,vy,thickness,Jdet; + IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input=NULL; + Input* vy_input=NULL; + if(domaintype==Domain2DhorizontalEnum){ + vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + } + else{ + vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input); + vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input); + } + + /*Start looping on Gaussian points*/ + Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + thickness_input->GetInputValue(&thickness,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss); + + flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]); + } + + return flux; +} +/*}}}*/ IssmDouble Tria::IceVolume(bool scaled){/*{{{*/ /*The volume of a truncated prism is area_base * 1/numedges sum(length of edges)*/ Index: ../trunk-jpl/src/c/classes/FemModel.h =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.h (revision 23887) +++ ../trunk-jpl/src/c/classes/FemModel.h (revision 23888) @@ -101,6 +101,7 @@ void GroundedAreax(IssmDouble* pV, bool scaled); void IcefrontAreax(); void IcefrontMassFluxx(IssmDouble* presponse, bool scaled); + void GroundinglineMassFluxx(IssmDouble* presponse, bool scaled); void IceMassx(IssmDouble* pV, bool scaled); void IceVolumex(IssmDouble* pV, bool scaled); void IceVolumeAboveFloatationx(IssmDouble* pV, bool scaled);