Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23899) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23900) @@ -241,4 +241,5 @@ virtual IssmDouble IceVolumeAboveFloatation(bool scaled)=0; virtual IssmDouble IcefrontMassFlux(bool scaled){_error_("not implemented");}; + virtual IssmDouble IcefrontMassFluxLevelset(bool scaled){_error_("not implemented");}; virtual IssmDouble GroundinglineMassFlux(bool scaled){_error_("not implemented");}; virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23899) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23900) @@ -1507,5 +1507,4 @@ /* Intermediaries */ - int i, dir,nrfrontnodes; IssmDouble levelset[NUMVERTICES]; int indicesfront[NUMVERTICES]; @@ -1515,13 +1514,12 @@ /* Get nodes where there is no ice */ - nrfrontnodes=0; - for(i=0;i=0.){ - indicesfront[nrfrontnodes]=i; - nrfrontnodes++; - } - } - - _assert_(nrfrontnodes==2); + indicesfront[num_frontnodes]=i; + num_frontnodes++; + } + } + _assert_(num_frontnodes==2); /* arrange order of frontnodes such that they are oriented counterclockwise */ @@ -1532,9 +1530,9 @@ } - IssmDouble* xyz_front = xNew(3*nrfrontnodes); + IssmDouble* xyz_front = xNew(3*2); /* Return nodes */ - for(i=0;iFindParam(&domaintype,DomainTypeEnum); + + /*Get ice front coordinates*/ + IssmDouble *xyz_list = NULL; + IssmDouble* xyz_front = NULL; + this->GetVerticesCoordinates(&xyz_list); + this->GetIcefrontCoordinates(&xyz_front,xyz_list,MaskIceLevelsetEnum); + + /*Get normal vector*/ + IssmDouble normal[3]; + this->NormalSection(&normal[0],xyz_front); + 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,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,gauss); + + flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]); + } + + return flux; +} +/*}}}*/ +IssmDouble Tria::IcefrontMassFluxLevelset(bool scaled){/*{{{*/ /*Make sure there is an ice front here*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 23899) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 23900) @@ -96,4 +96,5 @@ IssmDouble IceVolumeAboveFloatation(bool scaled); IssmDouble IcefrontMassFlux(bool scaled); + IssmDouble IcefrontMassFluxLevelset(bool scaled); IssmDouble GroundinglineMassFlux(bool scaled); void Ismip6FloatingiceMeltingRate(void); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23899) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 23900) @@ -1546,4 +1546,20 @@ }/*}}}*/ +void FemModel::IcefrontMassFluxLevelsetx(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->IcefrontMassFluxLevelset(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::GroundinglineMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/ @@ -2215,4 +2231,5 @@ case IceMassEnum: this->IceMassx(&double_result,false); break; case IcefrontMassFluxEnum: this->IcefrontMassFluxx(&double_result,false); break; + case IcefrontMassFluxLevelsetEnum: this->IcefrontMassFluxLevelsetx(&double_result,false); break; case IceMassScaledEnum: this->IceMassx(&double_result,true); break; case IceVolumeEnum: this->IceVolumex(&double_result,false); break; @@ -2421,4 +2438,5 @@ case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(responses, true); break; case IcefrontMassFluxEnum: this->IcefrontMassFluxx(responses, false); break; + case IcefrontMassFluxLevelsetEnum: this->IcefrontMassFluxLevelsetx(responses, false); break; case GroundedAreaEnum: this->GroundedAreax(responses, false); break; case GroundedAreaScaledEnum: this->GroundedAreax(responses, true); break; Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 23899) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 23900) @@ -102,4 +102,5 @@ void IcefrontAreax(); void IcefrontMassFluxx(IssmDouble* presponse, bool scaled); + void IcefrontMassFluxLevelsetx(IssmDouble* presponse, bool scaled); void GroundinglineMassFluxx(IssmDouble* presponse, bool scaled); void IceMassx(IssmDouble* pV, bool scaled); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23899) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23900) @@ -1056,4 +1056,5 @@ GroundinglineMassFluxEnum, IcefrontMassFluxEnum, + IcefrontMassFluxLevelsetEnum, MasstransportAnalysisEnum, MasstransportSolutionEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23899) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23900) @@ -1060,4 +1060,5 @@ case GroundinglineMassFluxEnum : return "GroundinglineMassFlux"; case IcefrontMassFluxEnum : return "IcefrontMassFlux"; + case IcefrontMassFluxLevelsetEnum : return "IcefrontMassFluxLevelset"; case MasstransportAnalysisEnum : return "MasstransportAnalysis"; case MasstransportSolutionEnum : return "MasstransportSolution"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23899) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23900) @@ -1084,4 +1084,5 @@ else if (strcmp(name,"GroundinglineMassFlux")==0) return GroundinglineMassFluxEnum; else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum; + else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum; else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"Nodal")==0) return NodalEnum; else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; - else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum; + if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; + else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum; else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum; else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"Vertices")==0) return VerticesEnum; else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; - else if (strcmp(name,"Water")==0) return WaterEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum; + if (strcmp(name,"Water")==0) return WaterEnum; + else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum; else if (strcmp(name,"XY")==0) return XYEnum; else if (strcmp(name,"XYZ")==0) return XYZEnum;