Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16274) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16275) @@ -96,4 +96,5 @@ virtual void ElementResponse(IssmDouble* presponse,int response_enum)=0; virtual IssmDouble IceVolume(void)=0; + virtual IssmDouble IceVolumeAboveFloatation(void)=0; virtual IssmDouble TotalSmb(void)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16274) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16275) @@ -3684,4 +3684,35 @@ } /*}}}*/ +/*FUNCTION Penta::IceVolumeAboveFloatation {{{*/ +IssmDouble Penta::IceVolumeAboveFloatation(void){ + + /*Volume above floatation: H + rho_water/rho_ice*bathymetry for nodes on the bed*/ + IssmDouble rho_ice,rho_water; + IssmDouble base,bed,surface,bathymetry; + IssmDouble xyz_list[NUMVERTICES][3]; + + if(NoIceInElement() || IsFloating() || !IsOnBed())return 0; + + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*First calculate the area of the base (cross section triangle) + * http://en.wikipedia.org/wiki/Pentangle + * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ + base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); + + /*Now get the average height above floatation*/ + Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); + Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); + Input* bathymetry_input = inputs->GetInput(BathymetryEnum); _assert_(bathymetry_input); + surface_input->GetInputAverage(&surface); + bed_input->GetInputAverage(&bed); + bathymetry_input->GetInputAverage(&bathymetry); + + /*Return: */ + return base*(surface - bed + rho_water/rho_ice * bathymetry); +} +/*}}}*/ /*FUNCTION Penta::MinVel{{{*/ void Penta::MinVel(IssmDouble* pminvel){ @@ -4849,5 +4880,5 @@ IssmDouble B_average,s_average; int *doflist = NULL; - //IssmDouble pressure[numdof]; + IssmDouble pressure[numdof]; /*Get dof list: */ @@ -4857,7 +4888,7 @@ for(i=0;imatpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]); - //if(values[i]TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50; + GetInputListOnVertices(&pressure[0],PressureEnum); + if(values[i]>matpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]); + if(values[i]TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50; /*Check solution*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16274) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16275) @@ -115,4 +115,5 @@ void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); IssmDouble IceVolume(void); + IssmDouble IceVolumeAboveFloatation(void); IssmDouble TotalSmb(void); void MinVel(IssmDouble* pminvel); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16274) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16275) @@ -2537,4 +2537,35 @@ /*Return: */ return base*(surface-bed); +} +/*}}}*/ +/*FUNCTION Tria::IceVolumeAboveFloatation {{{*/ +IssmDouble Tria::IceVolumeAboveFloatation(void){ + + /*The volume above floatation: H + rho_water/rho_ice * bathymetry */ + IssmDouble rho_ice,rho_water; + IssmDouble base,surface,bed,bathymetry; + IssmDouble xyz_list[NUMVERTICES][3]; + + if(NoIceInElement() || IsFloating())return 0; + + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*First calculate the area of the base (cross section triangle) + * http://en.wikipedia.org/wiki/Triangle + * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ + base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); + + /*Now get the average height and bathymetry*/ + Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); + Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); + Input* bathymetry_input = inputs->GetInput(BathymetryEnum); _assert_(bathymetry_input); + surface_input->GetInputAverage(&surface); + bed_input->GetInputAverage(&bed); + bathymetry_input->GetInputAverage(&bathymetry); + + /*Return: */ + return base*(surface-bed+rho_water/rho_ice*bathymetry); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16274) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16275) @@ -112,4 +112,5 @@ void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); IssmDouble IceVolume(void); + IssmDouble IceVolumeAboveFloatation(void); IssmDouble TotalSmb(void); void MinVel(IssmDouble* pminvel); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16274) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16275) @@ -430,35 +430,36 @@ #ifdef _HAVE_RESPONSES_ - case IceVolumeEnum: this->IceVolumex(responses); break; - case MinVelEnum: this->MinVelx(responses); break; - case MaxVelEnum: this->MaxVelx( responses); break; - case MinVxEnum: this->MinVxx(responses); break; - case MaxVxEnum: this->MaxVxx( responses); break; - case MaxAbsVxEnum: this->MaxAbsVxx( responses); break; - case MinVyEnum: this->MinVyx(responses); break; - case MaxVyEnum: this->MaxVyx( responses); break; - case MaxAbsVyEnum: this->MaxAbsVyx( responses); break; - case MinVzEnum: this->MinVzx(responses); break; - case MaxVzEnum: this->MaxVzx( responses); break; - case MaxAbsVzEnum: this->MaxAbsVzx( responses); break; - case MassFluxEnum: this->MassFluxx( responses); break; + case IceVolumeEnum: this->IceVolumex(responses); break; + case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(responses); break; + case MinVelEnum: this->MinVelx(responses); break; + case MaxVelEnum: this->MaxVelx( responses); break; + case MinVxEnum: this->MinVxx(responses); break; + case MaxVxEnum: this->MaxVxx( responses); break; + case MaxAbsVxEnum: this->MaxAbsVxx( responses); break; + case MinVyEnum: this->MinVyx(responses); break; + case MaxVyEnum: this->MaxVyx( responses); break; + case MaxAbsVyEnum: this->MaxAbsVyx( responses); break; + case MinVzEnum: this->MinVzx(responses); break; + case MaxVzEnum: this->MaxVzx( responses); break; + case MaxAbsVzEnum: this->MaxAbsVzx( responses); break; + case MassFluxEnum: this->MassFluxx( responses); break; #ifdef _HAVE_CONTROL_ - case SurfaceAbsVelMisfitEnum: SurfaceAbsVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; - case SurfaceRelVelMisfitEnum: SurfaceRelVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; - case SurfaceLogVelMisfitEnum: SurfaceLogVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; - case SurfaceLogVxVyMisfitEnum: SurfaceLogVxVyMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; - case SurfaceAverageVelMisfitEnum:SurfaceAverageVelMisfitx(responses,this,weight_index); break; - case ThicknessAbsMisfitEnum: ThicknessAbsMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; - case ThicknessAbsGradientEnum: this->ThicknessAbsGradientx(responses,weight_index); break; - case ThicknessAlongGradientEnum: ThicknessAlongGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; - case ThicknessAcrossGradientEnum:ThicknessAcrossGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; - case RheologyBbarAbsGradientEnum:RheologyBbarAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; + case SurfaceAbsVelMisfitEnum: SurfaceAbsVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; + case SurfaceRelVelMisfitEnum: SurfaceRelVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; + case SurfaceLogVelMisfitEnum: SurfaceLogVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; + case SurfaceLogVxVyMisfitEnum: SurfaceLogVxVyMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; + case SurfaceAverageVelMisfitEnum: SurfaceAverageVelMisfitx(responses,this,weight_index); break; + case ThicknessAbsMisfitEnum: ThicknessAbsMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; + case ThicknessAbsGradientEnum: this->ThicknessAbsGradientx(responses,weight_index); break; + case ThicknessAlongGradientEnum: ThicknessAlongGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; + case ThicknessAcrossGradientEnum: ThicknessAcrossGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; + case RheologyBbarAbsGradientEnum: RheologyBbarAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break; - case BalancethicknessMisfitEnum:BalancethicknessMisfitx(responses,weight_index); break; + case BalancethicknessMisfitEnum: BalancethicknessMisfitx(responses,weight_index); break; #endif - case TotalSmbEnum: this->TotalSmbx(responses); break; - case MaterialsRheologyBbarEnum: this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break; - case VelEnum: this->ElementResponsex(responses,VelEnum); break; - case FrictionCoefficientEnum:NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break; + case TotalSmbEnum: this->TotalSmbx(responses); break; + case MaterialsRheologyBbarEnum: this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break; + case VelEnum: this->ElementResponsex(responses,VelEnum); break; + case FrictionCoefficientEnum: NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break; default: _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); break; #else @@ -491,4 +492,8 @@ Responsex(&output_value,"IceVolume",0); results->AddObject(new GenericExternalResult(results->Size()+1,IceVolumeEnum,reCast(output_value),step,time)); + break; + case IceVolumeAboveFloatationEnum: + Responsex(&output_value,"IceVolumeAboveFloatation",0); + results->AddObject(new GenericExternalResult(results->Size()+1,IceVolumeAboveFloatationEnum,reCast(output_value),step,time)); break; case TotalSmbEnum: @@ -934,4 +939,20 @@ }/*}}}*/ +void FemModel::IceVolumeAboveFloatationx(IssmDouble* pV){/*{{{*/ + + IssmDouble local_ice_volume_af = 0; + IssmDouble total_ice_volume_af; + + for(int i=0;ielements->Size();i++){ + Element* element=dynamic_cast(this->elements->GetObjectByOffset(i)); + local_ice_volume_af+=element->IceVolumeAboveFloatation(); + } + ISSM_MPI_Reduce(&local_ice_volume_af,&total_ice_volume_af,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&total_ice_volume_af,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Assign output pointers: */ + *pV=total_ice_volume_af; + +}/*}}}*/ void FemModel::ElementResponsex(IssmDouble* presponse,int response_enum){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 16274) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 16275) @@ -75,4 +75,5 @@ void TotalSmbx(IssmDouble* pSmb); void IceVolumex(IssmDouble* pV); + void IceVolumeAboveFloatationx(IssmDouble* pV); void ElementResponsex(IssmDouble* presponse,int response_enum); void BalancethicknessMisfitx(IssmDouble* pV,int weight_index); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16274) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16275) @@ -553,4 +553,5 @@ MaxAbsVzEnum, IceVolumeEnum, + IceVolumeAboveFloatationEnum, TotalSmbEnum, /*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16274) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16275) @@ -541,4 +541,5 @@ case MaxAbsVzEnum : return "MaxAbsVz"; case IceVolumeEnum : return "IceVolume"; + case IceVolumeAboveFloatationEnum : return "IceVolumeAboveFloatation"; case TotalSmbEnum : return "TotalSmb"; case AbsoluteEnum : return "Absolute"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16274) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16275) @@ -553,4 +553,5 @@ else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum; else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; + else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;