source: issm/oecreview/Archive/23390-24306/ISSM-23887-23888.diff@ 24308

Last change on this file since 24308 was 24307, checked in by Mathieu Morlighem, 5 years ago

NEW: adding Archive/23390-24306
File size: 15.2 KB

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

     
    15451545        *pM=total_mass_flux;
    15461546
    15471547}/*}}}*/
     1548void FemModel::GroundinglineMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/
     1549
     1550        IssmDouble local_mass_flux = 0;
     1551        IssmDouble total_mass_flux;
     1552
     1553        for(int i=0;i<this->elements->Size();i++){
     1554                Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
     1555                local_mass_flux+=element->GroundinglineMassFlux(scaled);
     1556        }
     1557        ISSM_MPI_Reduce(&local_mass_flux,&total_mass_flux,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     1558        ISSM_MPI_Bcast(&total_mass_flux,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     1559
     1560        /*Assign output pointers: */
     1561        *pM=total_mass_flux;
     1562
     1563}/*}}}*/
    15481564void FemModel::IceMassx(IssmDouble* pM, bool scaled){/*{{{*/
    15491565
    15501566        IssmDouble local_ice_mass = 0;
    … …  
    22052221                                        case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(&double_result,true);  break;
    22062222                                        case GroundedAreaEnum:                   this->GroundedAreax(&double_result,false);             break;
    22072223                                        case GroundedAreaScaledEnum:             this->GroundedAreax(&double_result,true);              break;
     2224                                        case GroundinglineMassFluxEnum:          this->GroundinglineMassFluxx(&double_result,false);    break;
    22082225                                        case FloatingAreaEnum:                   this->FloatingAreax(&double_result,false);             break;
    22092226                                        case FloatingAreaScaledEnum:             this->FloatingAreax(&double_result,true);              break;
    22102227                                        case MinVelEnum:                         this->MinVelx(&double_result);                         break;
    … …  
    24022419                case IceVolumeScaledEnum:                this->IceVolumex(responses, true); break;
    24032420                case IceVolumeAboveFloatationEnum:       this->IceVolumeAboveFloatationx(responses, false); break;
    24042421                case IceVolumeAboveFloatationScaledEnum: this->IceVolumeAboveFloatationx(responses, true); break;
    2405                 case IcefrontMassFluxEnum:               this->IcefrontMassFluxx(responses, true); break;
     2422                case IcefrontMassFluxEnum:               this->IcefrontMassFluxx(responses, false); break;
    24062423                case GroundedAreaEnum:                   this->GroundedAreax(responses, false); break;
    24072424                case GroundedAreaScaledEnum:             this->GroundedAreax(responses, true); break;
    24082425                case FloatingAreaEnum:                   this->FloatingAreax(responses, false); break;

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

     
    9494                bool        HasEdgeOnSurface();
    9595                IssmDouble  IceVolume(bool scaled);
    9696                IssmDouble  IceVolumeAboveFloatation(bool scaled);
     97                IssmDouble  IcefrontMassFlux(bool scaled);
     98                IssmDouble  GroundinglineMassFlux(bool scaled);
    9799                void        Ismip6FloatingiceMeltingRate(void);
    98100                void        InputDepthAverageAtBase(int enum_type,int average_enum_type);
    99101                void        InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};

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

     
    10791079        IssmDouble* H=NULL;;
    10801080        int nrfrontbed,numiceverts;
    10811081
     1082        if(!this->IsOnBase()) return 0;
     1083        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
     1084
    10821085        /*Retrieve all inputs and parameters*/
    10831086        GetInputListOnVertices(&bed[0],BedEnum);
    10841087        GetInputListOnVertices(&surfaces[0],SurfaceEnum);
    … …  
    10851088        GetInputListOnVertices(&bases[0],BaseEnum);
    10861089        GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
    10871090
    1088         if(!this->IsOnBase()) return 0;
    1089         if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
    1090 
    10911091        nrfrontbed=0;
    10921092        for(int i=0;i<NUMVERTICES2D;i++){
    10931093                /*Find if bed<0*/

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

     
    240240                virtual IssmDouble IceVolume(bool scaled)=0;
    241241                virtual IssmDouble IceVolumeAboveFloatation(bool scaled)=0;
    242242                virtual IssmDouble IcefrontMassFlux(bool scaled){_error_("not implemented");};
     243                virtual IssmDouble GroundinglineMassFlux(bool scaled){_error_("not implemented");};
    243244                virtual void       InputDepthAverageAtBase(int enum_type,int average_enum_type)=0;
    244245                virtual void       InputExtrude(int input_enum,int start)=0;
    245246                virtual void       InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum)=0;

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

     
    14081408        IssmDouble* H=NULL;;
    14091409        int nrfrontbed,numiceverts;
    14101410
     1411        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
     1412
    14111413        /*Retrieve all inputs and parameters*/
    14121414        GetInputListOnVertices(&bed[0],BedEnum);
    14131415        GetInputListOnVertices(&surfaces[0],SurfaceEnum);
    … …  
    14141416        GetInputListOnVertices(&bases[0],BaseEnum);
    14151417        GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
    14161418
    1417         if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
    1418 
    14191419        nrfrontbed=0;
    14201420        for(int i=0;i<NUMVERTICES;i++){
    14211421                /*Find if bed<0*/
    … …  
    19461946        }
    19471947}
    19481948/*}}}*/
     1949IssmDouble Tria::IcefrontMassFlux(bool scaled){/*{{{*/
     1950
     1951        /*Make sure there is an ice front here*/
     1952        if(!IsIceInElement() || !IsIcefront()) return 0;
     1953
     1954        /*Scaled not implemented yet...*/
     1955        _assert_(!scaled);
     1956
     1957        int               domaintype,index1,index2;
     1958        const IssmPDouble epsilon = 1.e-15;
     1959        IssmDouble        s1,s2;
     1960        IssmDouble        gl[NUMVERTICES];
     1961        IssmDouble        xyz_front[2][3];
     1962
     1963
     1964        IssmDouble *xyz_list = NULL;
     1965        this->GetVerticesCoordinates(&xyz_list);
     1966
     1967        /*Recover parameters and values*/
     1968        parameters->FindParam(&domaintype,DomainTypeEnum);
     1969        GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
     1970
     1971        /*Be sure that values are not zero*/
     1972        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     1973        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     1974        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     1975
     1976        if(domaintype==Domain2DverticalEnum){
     1977                _error_("not implemented");
     1978        }
     1979        else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
     1980                int pt1 = 0;
     1981                int pt2 = 1;
     1982                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     1983
     1984                        /*Portion of the segments*/
     1985                        s1=gl[2]/(gl[2]-gl[1]);
     1986                        s2=gl[2]/(gl[2]-gl[0]);
     1987                        if(gl[2]<0.){
     1988                                pt1 = 1; pt2 = 0;
     1989                        }
     1990                        xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
     1991                        xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
     1992                        xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
     1993                        xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
     1994                        xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
     1995                        xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
     1996                }
     1997                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
     1998
     1999                        /*Portion of the segments*/
     2000                        s1=gl[0]/(gl[0]-gl[1]);
     2001                        s2=gl[0]/(gl[0]-gl[2]);
     2002                        if(gl[0]<0.){
     2003                                pt1 = 1; pt2 = 0;
     2004                        }
     2005
     2006                        xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
     2007                        xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
     2008                        xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
     2009                        xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
     2010                        xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
     2011                        xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
     2012                }
     2013                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
     2014
     2015                        /*Portion of the segments*/
     2016                        s1=gl[1]/(gl[1]-gl[0]);
     2017                        s2=gl[1]/(gl[1]-gl[2]);
     2018                        if(gl[1]<0.){
     2019                                pt1 = 1; pt2 = 0;
     2020                        }
     2021
     2022                        xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
     2023                        xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
     2024                        xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
     2025                        xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
     2026                        xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
     2027                        xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
     2028                }
     2029                else{
     2030                        _error_("case not possible");
     2031                }
     2032
     2033        }
     2034        else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
     2035
     2036        /*Some checks in debugging mode*/
     2037        _assert_(s1>=0 && s1<=1.);
     2038        _assert_(s2>=0 && s2<=1.);
     2039
     2040        /*Get normal vector*/
     2041        IssmDouble normal[3];
     2042        this->NormalSection(&normal[0],&xyz_front[0][0]);
     2043        normal[0] = -normal[0];
     2044        normal[1] = -normal[1];
     2045
     2046        /*Get inputs*/
     2047        IssmDouble flux = 0.;
     2048        IssmDouble vx,vy,thickness,Jdet;
     2049        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     2050        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     2051        Input* vx_input=NULL;
     2052        Input* vy_input=NULL;
     2053        if(domaintype==Domain2DhorizontalEnum){
     2054                vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     2055                vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     2056        }
     2057        else{
     2058                vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
     2059                vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
     2060        }
     2061
     2062        /*Start looping on Gaussian points*/
     2063        Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
     2064        for(int ig=gauss->begin();ig<gauss->end();ig++){
     2065
     2066                gauss->GaussPoint(ig);
     2067                thickness_input->GetInputValue(&thickness,gauss);
     2068                vx_input->GetInputValue(&vx,gauss);
     2069                vy_input->GetInputValue(&vy,gauss);
     2070                this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
     2071
     2072                flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
     2073        }
     2074
     2075        return flux;
     2076}
     2077/*}}}*/
     2078IssmDouble Tria::GroundinglineMassFlux(bool scaled){/*{{{*/
     2079
     2080        /*Make sure there is an ice front here*/
     2081        if(!IsIceInElement()) return 0;
     2082        if(!IsZeroLevelset(MaskGroundediceLevelsetEnum)) return 0;
     2083
     2084        /*Scaled not implemented yet...*/
     2085        _assert_(!scaled);
     2086
     2087        int               domaintype,index1,index2;
     2088        const IssmPDouble epsilon = 1.e-15;
     2089        IssmDouble        s1,s2;
     2090        IssmDouble        gl[NUMVERTICES];
     2091        IssmDouble        xyz_front[2][3];
     2092
     2093        IssmDouble *xyz_list = NULL;
     2094        this->GetVerticesCoordinates(&xyz_list);
     2095
     2096        /*Recover parameters and values*/
     2097        parameters->FindParam(&domaintype,DomainTypeEnum);
     2098        GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
     2099
     2100        /*Be sure that values are not zero*/
     2101        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     2102        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     2103        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     2104
     2105        if(domaintype==Domain2DverticalEnum){
     2106                _error_("not implemented");
     2107        }
     2108        else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
     2109                int pt1 = 0;
     2110                int pt2 = 1;
     2111                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     2112
     2113                        /*Portion of the segments*/
     2114                        s1=gl[2]/(gl[2]-gl[1]);
     2115                        s2=gl[2]/(gl[2]-gl[0]);
     2116                        if(gl[2]<0.){
     2117                                pt1 = 1; pt2 = 0;
     2118                        }
     2119                        xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
     2120                        xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
     2121                        xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
     2122                        xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
     2123                        xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
     2124                        xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
     2125                }
     2126                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
     2127
     2128                        /*Portion of the segments*/
     2129                        s1=gl[0]/(gl[0]-gl[1]);
     2130                        s2=gl[0]/(gl[0]-gl[2]);
     2131                        if(gl[0]<0.){
     2132                                pt1 = 1; pt2 = 0;
     2133                        }
     2134
     2135                        xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
     2136                        xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
     2137                        xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
     2138                        xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
     2139                        xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
     2140                        xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
     2141                }
     2142                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
     2143
     2144                        /*Portion of the segments*/
     2145                        s1=gl[1]/(gl[1]-gl[0]);
     2146                        s2=gl[1]/(gl[1]-gl[2]);
     2147                        if(gl[1]<0.){
     2148                                pt1 = 1; pt2 = 0;
     2149                        }
     2150
     2151                        xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
     2152                        xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
     2153                        xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
     2154                        xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
     2155                        xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
     2156                        xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
     2157                }
     2158                else{
     2159                        _error_("case not possible");
     2160                }
     2161
     2162        }
     2163        else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
     2164
     2165        /*Some checks in debugging mode*/
     2166        _assert_(s1>=0 && s1<=1.);
     2167        _assert_(s2>=0 && s2<=1.);
     2168
     2169        /*Get normal vector*/
     2170        IssmDouble normal[3];
     2171        this->NormalSection(&normal[0],&xyz_front[0][0]);
     2172
     2173        /*Get inputs*/
     2174        IssmDouble flux = 0.;
     2175        IssmDouble vx,vy,thickness,Jdet;
     2176        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     2177        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     2178        Input* vx_input=NULL;
     2179        Input* vy_input=NULL;
     2180        if(domaintype==Domain2DhorizontalEnum){
     2181                vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     2182                vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     2183        }
     2184        else{
     2185                vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
     2186                vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
     2187        }
     2188
     2189        /*Start looping on Gaussian points*/
     2190        Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
     2191        for(int ig=gauss->begin();ig<gauss->end();ig++){
     2192
     2193                gauss->GaussPoint(ig);
     2194                thickness_input->GetInputValue(&thickness,gauss);
     2195                vx_input->GetInputValue(&vx,gauss);
     2196                vy_input->GetInputValue(&vy,gauss);
     2197                this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
     2198
     2199                flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
     2200        }
     2201
     2202        return flux;
     2203}
     2204/*}}}*/
    19492205IssmDouble Tria::IceVolume(bool scaled){/*{{{*/
    19502206
    19512207        /*The volume of a truncated prism is area_base * 1/numedges sum(length of edges)*/

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

     
    101101                void GroundedAreax(IssmDouble* pV, bool scaled);
    102102                void IcefrontAreax();
    103103                void IcefrontMassFluxx(IssmDouble* presponse, bool scaled);
     104                void GroundinglineMassFluxx(IssmDouble* presponse, bool scaled);
    104105                void IceMassx(IssmDouble* pV, bool scaled);
    105106                void IceVolumex(IssmDouble* pV, bool scaled);
    106107                void IceVolumeAboveFloatationx(IssmDouble* pV, bool scaled);
Note: See TracBrowser for help on using the repository browser.