Changeset 16275


Ignore:
Timestamp:
09/30/13 20:08:04 (11 years ago)
Author:
seroussi
Message:

NEW: added VolumeAboveFloatation calculation

Location:
issm/trunk-jpl/src/c
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r16272 r16275  
    9696                virtual void   ElementResponse(IssmDouble* presponse,int response_enum)=0;
    9797                virtual IssmDouble IceVolume(void)=0;
     98                virtual IssmDouble IceVolumeAboveFloatation(void)=0;
    9899                virtual IssmDouble TotalSmb(void)=0;
    99100                #endif

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r16272 r16275  
    36843684}
    36853685/*}}}*/
     3686/*FUNCTION Penta::IceVolumeAboveFloatation {{{*/
     3687IssmDouble Penta::IceVolumeAboveFloatation(void){
     3688
     3689        /*Volume above floatation: H + rho_water/rho_ice*bathymetry for nodes on the bed*/
     3690        IssmDouble rho_ice,rho_water;
     3691        IssmDouble base,bed,surface,bathymetry;
     3692        IssmDouble xyz_list[NUMVERTICES][3];
     3693
     3694        if(NoIceInElement() || IsFloating() || !IsOnBed())return 0;
     3695
     3696        rho_ice=matpar->GetRhoIce();
     3697        rho_water=matpar->GetRhoWater();
     3698        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     3699
     3700        /*First calculate the area of the base (cross section triangle)
     3701         * http://en.wikipedia.org/wiki/Pentangle
     3702         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     3703        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]));
     3704
     3705        /*Now get the average height above floatation*/
     3706        Input* surface_input    = inputs->GetInput(SurfaceEnum);    _assert_(surface_input);
     3707        Input* bed_input        = inputs->GetInput(BedEnum);        _assert_(bed_input);
     3708        Input* bathymetry_input = inputs->GetInput(BathymetryEnum); _assert_(bathymetry_input);
     3709        surface_input->GetInputAverage(&surface);
     3710        bed_input->GetInputAverage(&bed);
     3711        bathymetry_input->GetInputAverage(&bathymetry);
     3712
     3713        /*Return: */
     3714        return base*(surface - bed + rho_water/rho_ice * bathymetry);
     3715}
     3716/*}}}*/
    36863717/*FUNCTION Penta::MinVel{{{*/
    36873718void  Penta::MinVel(IssmDouble* pminvel){
    … …  
    48494880        IssmDouble  B_average,s_average;
    48504881        int        *doflist = NULL;
    4851         //IssmDouble  pressure[numdof];
     4882        IssmDouble  pressure[numdof];
    48524883
    48534884        /*Get dof list: */
    … …  
    48574888        for(i=0;i<numdof;i++){
    48584889                values[i]=solution[doflist[i]];
    4859                 //GetInputListOnVertices(&pressure[0],PressureEnum);
    4860                 //if(values[i]>matpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]);
    4861                 //if(values[i]<matpar->TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50;
     4890                GetInputListOnVertices(&pressure[0],PressureEnum);
     4891                if(values[i]>matpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]);
     4892                if(values[i]<matpar->TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50;
    48624893
    48634894                /*Check solution*/

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r16272 r16275  
    115115                void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
    116116                IssmDouble IceVolume(void);
     117                IssmDouble IceVolumeAboveFloatation(void);
    117118                IssmDouble TotalSmb(void);
    118119                void   MinVel(IssmDouble* pminvel);

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r16256 r16275  
    25372537        /*Return: */
    25382538        return base*(surface-bed);
     2539}
     2540/*}}}*/
     2541/*FUNCTION Tria::IceVolumeAboveFloatation {{{*/
     2542IssmDouble Tria::IceVolumeAboveFloatation(void){
     2543
     2544        /*The volume above floatation: H + rho_water/rho_ice * bathymetry */
     2545        IssmDouble rho_ice,rho_water;
     2546        IssmDouble base,surface,bed,bathymetry;
     2547        IssmDouble xyz_list[NUMVERTICES][3];
     2548
     2549        if(NoIceInElement() || IsFloating())return 0;
     2550
     2551        rho_ice=matpar->GetRhoIce();
     2552        rho_water=matpar->GetRhoWater();
     2553        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     2554
     2555        /*First calculate the area of the base (cross section triangle)
     2556         * http://en.wikipedia.org/wiki/Triangle
     2557         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     2558        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]));
     2559
     2560        /*Now get the average height and bathymetry*/
     2561        Input* surface_input    = inputs->GetInput(SurfaceEnum);    _assert_(surface_input);
     2562        Input* bed_input        = inputs->GetInput(BedEnum);        _assert_(bed_input);
     2563        Input* bathymetry_input = inputs->GetInput(BathymetryEnum); _assert_(bathymetry_input);
     2564        surface_input->GetInputAverage(&surface);
     2565        bed_input->GetInputAverage(&bed);
     2566        bathymetry_input->GetInputAverage(&bathymetry);
     2567
     2568        /*Return: */
     2569        return base*(surface-bed+rho_water/rho_ice*bathymetry);
    25392570}
    25402571/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r16272 r16275  
    112112                void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
    113113                IssmDouble IceVolume(void);
     114                IssmDouble IceVolumeAboveFloatation(void);
    114115                IssmDouble TotalSmb(void);
    115116                void       MinVel(IssmDouble* pminvel);

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r16272 r16275  
    430430
    431431                #ifdef _HAVE_RESPONSES_
    432                 case IceVolumeEnum:              this->IceVolumex(responses); break;
    433                 case MinVelEnum:                 this->MinVelx(responses); break;
    434                 case MaxVelEnum:                 this->MaxVelx(                  responses); break;
    435                 case MinVxEnum:                  this->MinVxx(responses); break;
    436                 case MaxVxEnum:                  this->MaxVxx(                   responses); break;
    437                 case MaxAbsVxEnum:               this->MaxAbsVxx(                responses); break;
    438                 case MinVyEnum:                  this->MinVyx(responses); break;
    439                 case MaxVyEnum:                  this->MaxVyx(                   responses); break;
    440                 case MaxAbsVyEnum:               this->MaxAbsVyx(                responses); break;
    441                 case MinVzEnum:                  this->MinVzx(responses); break;
    442                 case MaxVzEnum:                  this->MaxVzx(                   responses); break;
    443                 case MaxAbsVzEnum:               this->MaxAbsVzx(                responses); break;
    444                 case MassFluxEnum:               this->MassFluxx(         responses); break;
     432                case IceVolumeEnum:                this->IceVolumex(responses); break;
     433                case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(responses); break;
     434                case MinVelEnum:                   this->MinVelx(responses); break;
     435                case MaxVelEnum:                   this->MaxVelx(                  responses); break;
     436                case MinVxEnum:                    this->MinVxx(responses); break;
     437                case MaxVxEnum:                    this->MaxVxx(                   responses); break;
     438                case MaxAbsVxEnum:                 this->MaxAbsVxx(                responses); break;
     439                case MinVyEnum:                    this->MinVyx(responses); break;
     440                case MaxVyEnum:                    this->MaxVyx(                   responses); break;
     441                case MaxAbsVyEnum:                 this->MaxAbsVyx(                responses); break;
     442                case MinVzEnum:                    this->MinVzx(responses); break;
     443                case MaxVzEnum:                    this->MaxVzx(                   responses); break;
     444                case MaxAbsVzEnum:                 this->MaxAbsVzx(                responses); break;
     445                case MassFluxEnum:                 this->MassFluxx(         responses); break;
    445446                #ifdef _HAVE_CONTROL_
    446                 case SurfaceAbsVelMisfitEnum:    SurfaceAbsVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
    447                 case SurfaceRelVelMisfitEnum:    SurfaceRelVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
    448                 case SurfaceLogVelMisfitEnum:    SurfaceLogVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
    449                 case SurfaceLogVxVyMisfitEnum:   SurfaceLogVxVyMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
    450                 case SurfaceAverageVelMisfitEnum:SurfaceAverageVelMisfitx(responses,this,weight_index); break;
    451                 case ThicknessAbsMisfitEnum:     ThicknessAbsMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
    452                 case ThicknessAbsGradientEnum:   this->ThicknessAbsGradientx(responses,weight_index); break;
    453                 case ThicknessAlongGradientEnum: ThicknessAlongGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
    454                 case ThicknessAcrossGradientEnum:ThicknessAcrossGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
    455                 case RheologyBbarAbsGradientEnum:RheologyBbarAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
     447                case SurfaceAbsVelMisfitEnum:      SurfaceAbsVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
     448                case SurfaceRelVelMisfitEnum:      SurfaceRelVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
     449                case SurfaceLogVelMisfitEnum:      SurfaceLogVelMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
     450                case SurfaceLogVxVyMisfitEnum:     SurfaceLogVxVyMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
     451                case SurfaceAverageVelMisfitEnum:  SurfaceAverageVelMisfitx(responses,this,weight_index); break;
     452                case ThicknessAbsMisfitEnum:       ThicknessAbsMisfitx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
     453                case ThicknessAbsGradientEnum:     this->ThicknessAbsGradientx(responses,weight_index); break;
     454                case ThicknessAlongGradientEnum:   ThicknessAlongGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
     455                case ThicknessAcrossGradientEnum:  ThicknessAcrossGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
     456                case RheologyBbarAbsGradientEnum:  RheologyBbarAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
    456457                case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters,weight_index); break;
    457                 case BalancethicknessMisfitEnum:BalancethicknessMisfitx(responses,weight_index); break;
     458                case BalancethicknessMisfitEnum:   BalancethicknessMisfitx(responses,weight_index); break;
    458459                #endif
    459                 case TotalSmbEnum:                                      this->TotalSmbx(responses); break;
    460                 case MaterialsRheologyBbarEnum: this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break;
    461                 case VelEnum:                   this->ElementResponsex(responses,VelEnum); break;
    462                 case FrictionCoefficientEnum:NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break;
     460                case TotalSmbEnum:                                        this->TotalSmbx(responses); break;
     461                case MaterialsRheologyBbarEnum:    this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break;
     462                case VelEnum:                      this->ElementResponsex(responses,VelEnum); break;
     463                case FrictionCoefficientEnum:      NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break;
    463464                default: _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); break;
    464465                #else
    … …  
    491492                                        Responsex(&output_value,"IceVolume",0);
    492493                                        results->AddObject(new GenericExternalResult<double>(results->Size()+1,IceVolumeEnum,reCast<IssmPDouble>(output_value),step,time));
     494                                        break;
     495                                case IceVolumeAboveFloatationEnum:
     496                                        Responsex(&output_value,"IceVolumeAboveFloatation",0);
     497                                        results->AddObject(new GenericExternalResult<double>(results->Size()+1,IceVolumeAboveFloatationEnum,reCast<IssmPDouble>(output_value),step,time));
    493498                                        break;
    494499                                case TotalSmbEnum:
    … …  
    934939
    935940}/*}}}*/
     941void FemModel::IceVolumeAboveFloatationx(IssmDouble* pV){/*{{{*/
     942
     943        IssmDouble local_ice_volume_af = 0;
     944        IssmDouble total_ice_volume_af;
     945
     946        for(int i=0;i<this->elements->Size();i++){
     947                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
     948                local_ice_volume_af+=element->IceVolumeAboveFloatation();
     949        }
     950        ISSM_MPI_Reduce(&local_ice_volume_af,&total_ice_volume_af,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     951        ISSM_MPI_Bcast(&total_ice_volume_af,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     952
     953        /*Assign output pointers: */
     954        *pV=total_ice_volume_af;
     955
     956}/*}}}*/
    936957void FemModel::ElementResponsex(IssmDouble* presponse,int response_enum){/*{{{*/
    937958

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r16272 r16275  
    7575                void TotalSmbx(IssmDouble* pSmb);
    7676                void IceVolumex(IssmDouble* pV);
     77                void IceVolumeAboveFloatationx(IssmDouble* pV);
    7778                void ElementResponsex(IssmDouble* presponse,int response_enum);
    7879                void BalancethicknessMisfitx(IssmDouble* pV,int weight_index);

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r16252 r16275  
    553553        MaxAbsVzEnum,
    554554        IceVolumeEnum,
     555        IceVolumeAboveFloatationEnum,
    555556        TotalSmbEnum,
    556557        /*}}}*/

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r16252 r16275  
    541541                case MaxAbsVzEnum : return "MaxAbsVz";
    542542                case IceVolumeEnum : return "IceVolume";
     543                case IceVolumeAboveFloatationEnum : return "IceVolumeAboveFloatation";
    543544                case TotalSmbEnum : return "TotalSmb";
    544545                case AbsoluteEnum : return "Absolute";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r16252 r16275  
    553553              else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum;
    554554              else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
     555              else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
    555556              else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
    556557              else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
Note: See TracChangeset for help on using the changeset viewer.