Changeset 24119


Ignore:
Timestamp:
08/09/19 14:20:22 (6 years ago)
Author:
youngmc3
Message:

CHG: (Total)CalvingFluxLevelset, (Total)CalvingMeltingFluxLevelset in 3D

Location:
issm/trunk-jpl/src/c/classes
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r24097 r24119  
    329329        delete gauss;
    330330
     331}
     332/*}}}*/
     333void       Penta::CalvingFluxLevelset(){/*{{{*/
     334
     335        /*Make sure there is an ice front here*/
     336        if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)){
     337                IssmDouble flux_per_area=0;
     338                this->inputs->AddInput(new PentaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
     339        }
     340        else{
     341                int               domaintype,index1,index2;
     342                const IssmPDouble epsilon = 1.e-15;
     343                IssmDouble        s1,s2;
     344                IssmDouble        gl[NUMVERTICES];
     345                IssmDouble        xyz_front[2][3];
     346
     347                IssmDouble *xyz_list = NULL;
     348                this->GetVerticesCoordinates(&xyz_list);
     349
     350                /*Recover parameters and values*/
     351                GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
     352
     353                /*Be sure that values are not zero*/
     354                if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     355                if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     356                if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     357
     358                int pt1 = 0;
     359                int pt2 = 1;
     360                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     361
     362                        /*Portion of the segments*/
     363                        s1=gl[2]/(gl[2]-gl[1]);
     364                        s2=gl[2]/(gl[2]-gl[0]);
     365                        if(gl[2]<0.){
     366                                pt1 = 1; pt2 = 0;
     367                        }
     368                        xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
     369                        xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
     370                        xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
     371                        xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
     372                        xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
     373                        xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
     374                }
     375                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
     376
     377                        /*Portion of the segments*/
     378                        s1=gl[0]/(gl[0]-gl[1]);
     379                        s2=gl[0]/(gl[0]-gl[2]);
     380                        if(gl[0]<0.){
     381                                pt1 = 1; pt2 = 0;
     382                        }
     383
     384                        xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
     385                        xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
     386                        xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
     387                        xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
     388                        xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
     389                        xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
     390                }
     391                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
     392
     393                        /*Portion of the segments*/
     394                        s1=gl[1]/(gl[1]-gl[0]);
     395                        s2=gl[1]/(gl[1]-gl[2]);
     396                        if(gl[1]<0.){
     397                                pt1 = 1; pt2 = 0;
     398                        }
     399
     400                        xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
     401                        xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
     402                        xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
     403                        xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
     404                        xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
     405                        xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
     406                }
     407                else{
     408                        _error_("case not possible");
     409                }
     410
     411                /*Some checks in debugging mode*/
     412                _assert_(s1>=0 && s1<=1.);
     413                _assert_(s2>=0 && s2<=1.);
     414
     415                /*Get normal vector*/
     416                IssmDouble normal[3];
     417                this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
     418                normal[0] = -normal[0];
     419                normal[1] = -normal[1];
     420
     421                /*Get inputs*/
     422                IssmDouble flux = 0.;
     423                IssmDouble area = 0.;
     424                IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
     425                IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     426                Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     427                Input* calvingratex_input=NULL;
     428                Input* calvingratey_input=NULL;
     429                calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     430                calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     431
     432                /*Start looping on Gaussian points*/
     433                Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
     434                for(int ig=gauss->begin();ig<gauss->end();ig++){
     435
     436                        gauss->GaussPoint(ig);
     437                        thickness_input->GetInputValue(&thickness,gauss);
     438                        calvingratex_input->GetInputValue(&calvingratex,gauss);
     439                        calvingratey_input->GetInputValue(&calvingratey,gauss);
     440                        this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
     441
     442                        flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
     443                        area += Jdet*gauss->weight*thickness;
     444
     445                        flux_per_area=flux/area;
     446                }
     447
     448                this->inputs->AddInput(new PentaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
     449               
     450                /*Clean up and return*/
     451                delete gauss;
     452        }
     453}
     454/*}}}*/
     455void       Penta::CalvingMeltingFluxLevelset(){/*{{{*/
     456
     457        /*Make sure there is an ice front here*/
     458        if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)){
     459                IssmDouble flux_per_area=0;
     460                this->inputs->AddInput(new PentaInput(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum));
     461        }
     462        else{
     463                int               domaintype,index1,index2;
     464                const IssmPDouble epsilon = 1.e-15;
     465                IssmDouble        s1,s2;
     466                IssmDouble        gl[NUMVERTICES];
     467                IssmDouble        xyz_front[2][3];
     468
     469                IssmDouble *xyz_list = NULL;
     470                this->GetVerticesCoordinates(&xyz_list);
     471
     472                /*Recover parameters and values*/
     473                GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
     474
     475                /*Be sure that values are not zero*/
     476                if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     477                if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     478                if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     479
     480                int pt1 = 0;
     481                int pt2 = 1;
     482                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     483
     484                        /*Portion of the segments*/
     485                        s1=gl[2]/(gl[2]-gl[1]);
     486                        s2=gl[2]/(gl[2]-gl[0]);
     487                        if(gl[2]<0.){
     488                                pt1 = 1; pt2 = 0;
     489                        }
     490                        xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
     491                        xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
     492                        xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
     493                        xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
     494                        xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
     495                        xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
     496                }
     497                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
     498
     499                        /*Portion of the segments*/
     500                        s1=gl[0]/(gl[0]-gl[1]);
     501                        s2=gl[0]/(gl[0]-gl[2]);
     502                        if(gl[0]<0.){
     503                                pt1 = 1; pt2 = 0;
     504                        }
     505
     506                        xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
     507                        xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
     508                        xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
     509                        xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
     510                        xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
     511                        xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
     512                }
     513                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
     514
     515                        /*Portion of the segments*/
     516                        s1=gl[1]/(gl[1]-gl[0]);
     517                        s2=gl[1]/(gl[1]-gl[2]);
     518                        if(gl[1]<0.){
     519                                pt1 = 1; pt2 = 0;
     520                        }
     521
     522                        xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
     523                        xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
     524                        xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
     525                        xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
     526                        xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
     527                        xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
     528                }
     529                else{
     530                        _error_("case not possible");
     531                }
     532
     533                /*Some checks in debugging mode*/
     534                _assert_(s1>=0 && s1<=1.);
     535                _assert_(s2>=0 && s2<=1.);
     536
     537                /*Get normal vector*/
     538                IssmDouble normal[3];
     539                this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
     540                normal[0] = -normal[0];
     541                normal[1] = -normal[1];
     542
     543                /*Get inputs*/
     544                IssmDouble flux = 0.;
     545                IssmDouble area = 0.;
     546                IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
     547                IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     548                Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     549                Input* calvingratex_input=NULL;
     550                Input* calvingratey_input=NULL;
     551                Input* vx_input=NULL;
     552                Input* vy_input=NULL;
     553                Input* meltingrate_input=NULL;
     554                calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     555                calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     556                vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     557                vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     558                meltingrate_input=inputs->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input);
     559
     560                /*Start looping on Gaussian points*/
     561                Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
     562                for(int ig=gauss->begin();ig<gauss->end();ig++){
     563
     564                        gauss->GaussPoint(ig);
     565                        thickness_input->GetInputValue(&thickness,gauss);
     566                        calvingratex_input->GetInputValue(&calvingratex,gauss);
     567                        calvingratey_input->GetInputValue(&calvingratey,gauss);
     568                        vx_input->GetInputValue(&vx,gauss);
     569                        vy_input->GetInputValue(&vy,gauss);
     570                        vel=vx*vx+vy*vy;
     571                        meltingrate_input->GetInputValue(&meltingrate,gauss);   
     572                        meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
     573                        meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
     574                        this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
     575
     576                        flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
     577                        area += Jdet*gauss->weight*thickness;
     578
     579                        flux_per_area=flux/area;
     580                }
     581
     582                this->inputs->AddInput(new PentaInput(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum));
     583               
     584                /*Clean up and return*/
     585                delete gauss;
     586        }
    331587}
    332588/*}}}*/
    … …  
    15041760        /*Clean up and return*/
    15051761        return groundedarea;
     1762}
     1763/*}}}*/
     1764IssmDouble Penta::IcefrontMassFluxLevelset(bool scaled){/*{{{*/
     1765
     1766        /*Make sure there is an ice front here*/
     1767        if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum) || !IsOnBase()) return 0;
     1768
     1769        /*Scaled not implemented yet...*/
     1770        _assert_(!scaled);
     1771
     1772        int               domaintype,index1,index2;
     1773        const IssmPDouble epsilon = 1.e-15;
     1774        IssmDouble        s1,s2;
     1775        IssmDouble        gl[NUMVERTICES];
     1776        IssmDouble        xyz_front[2][3];
     1777
     1778        IssmDouble *xyz_list = NULL;
     1779        this->GetVerticesCoordinates(&xyz_list);
     1780
     1781        /*Recover parameters and values*/
     1782        GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
     1783
     1784        /*Be sure that values are not zero*/
     1785        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     1786        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     1787        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     1788
     1789        int pt1 = 0;
     1790        int pt2 = 1;
     1791        if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     1792
     1793                /*Portion of the segments*/
     1794                s1=gl[2]/(gl[2]-gl[1]);
     1795                s2=gl[2]/(gl[2]-gl[0]);
     1796                if(gl[2]<0.){
     1797                        pt1 = 1; pt2 = 0;
     1798                }
     1799                xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
     1800                xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
     1801                xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
     1802                xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
     1803                xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
     1804                xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
     1805        }
     1806        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
     1807
     1808                /*Portion of the segments*/
     1809                s1=gl[0]/(gl[0]-gl[1]);
     1810                s2=gl[0]/(gl[0]-gl[2]);
     1811                if(gl[0]<0.){
     1812                        pt1 = 1; pt2 = 0;
     1813                }
     1814
     1815                xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
     1816                xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
     1817                xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
     1818                xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
     1819                xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
     1820                xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
     1821        }
     1822        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
     1823
     1824                /*Portion of the segments*/
     1825                s1=gl[1]/(gl[1]-gl[0]);
     1826                s2=gl[1]/(gl[1]-gl[2]);
     1827                if(gl[1]<0.){
     1828                        pt1 = 1; pt2 = 0;
     1829                }
     1830
     1831                xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
     1832                xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
     1833                xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
     1834                xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
     1835                xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
     1836                xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
     1837        }
     1838        else{
     1839                _error_("case not possible");
     1840        }
     1841
     1842
     1843        /*Some checks in debugging mode*/
     1844        _assert_(s1>=0 && s1<=1.);
     1845        _assert_(s2>=0 && s2<=1.);
     1846
     1847        /*Get normal vector*/
     1848        IssmDouble normal[3];
     1849        this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
     1850        normal[0] = -normal[0];
     1851        normal[1] = -normal[1];
     1852       
     1853        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
     1854        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
     1855
     1856        /*Get inputs*/
     1857        IssmDouble flux = 0.;
     1858        IssmDouble vx,vy,thickness,Jdet;
     1859        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     1860        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     1861        Input* vx_input=NULL;
     1862        Input* vy_input=NULL;
     1863        vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
     1864        vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
     1865
     1866        /*Start looping on Gaussian points*/
     1867        Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
     1868        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1869
     1870                gauss->GaussPoint(ig);
     1871                thickness_input->GetInputValue(&thickness,gauss);
     1872                vx_input->GetInputValue(&vx,gauss);
     1873                vy_input->GetInputValue(&vy,gauss);
     1874                this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
     1875
     1876                flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
     1877        }
     1878
     1879        return flux;
    15061880}
    15071881/*}}}*/
    … …  
    21482522}
    21492523/*}}}*/
     2524Gauss*     Penta::NewGaussBase(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz){/*{{{*/
     2525
     2526        IssmDouble  area_coordinates[2][3];
     2527
     2528        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
     2529
     2530        return new GaussPenta(area_coordinates,order_horiz);
     2531}
     2532/*}}}*/
    21502533Gauss*     Penta::NewGaussLine(int vertex1,int vertex2,int order){/*{{{*/
    21512534        return new GaussPenta(vertex1,vertex2,order);
    … …  
    22972680
    22982681        for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
     2682}
     2683/*}}}*/
     2684void       Penta::NormalSectionBase(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
     2685
     2686        /*Build unit outward pointing vector*/
     2687        IssmDouble vector[2];
     2688        IssmDouble norm;
     2689
     2690        vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
     2691        vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
     2692
     2693        norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
     2694
     2695        normal[0]= + vector[1]/norm;
     2696        normal[1]= - vector[0]/norm;
    22992697}
    23002698/*}}}*/
    … …  
    25032901void       Penta::RignotMeltParameterization(){/*{{{*/
    25042902
    2505    IssmDouble A, B, alpha, beta;
     2903        if(!this->IsOnBase()) return;
     2904   
     2905        IssmDouble A, B, alpha, beta;
    25062906        IssmDouble bed,qsg,qsg_basin,TF,yts;
    25072907        int numbasins;
    … …  
    25432943
    25442944                        /* calculate melt rates */
    2545                         meltrates[iv]=(A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta)/86400; //[m/s]
     2945                        meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s]
    25462946                }       
    25472947
    … …  
    26853085        if(this->inputs->GetInput(CalvingratexEnum)) this->InputDepthAverageAtBase(CalvingratexEnum,CalvingratexAverageEnum);
    26863086        if(this->inputs->GetInput(CalvingrateyEnum)) this->InputDepthAverageAtBase(CalvingrateyEnum,CalvingrateyAverageEnum);
     3087       
    26873088        Tria* tria=(Tria*)SpawnTria(0,1,2);
    26883089        switch(this->material->ObjectEnum()){
    … …  
    30063407        return dt;
    30073408}/*}}}*/
     3409IssmDouble Penta::TotalCalvingFluxLevelset(bool scaled){/*{{{*/
     3410
     3411        /*Make sure there is an ice front here*/
     3412        if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum) || !IsOnBase()) return 0;
     3413
     3414        /*Scaled not implemented yet...*/
     3415        _assert_(!scaled);
     3416
     3417        int               domaintype,index1,index2;
     3418        const IssmPDouble epsilon = 1.e-15;
     3419        IssmDouble        s1,s2;
     3420        IssmDouble        gl[NUMVERTICES];
     3421        IssmDouble        xyz_front[2][3];
     3422
     3423        IssmDouble *xyz_list = NULL;
     3424        this->GetVerticesCoordinates(&xyz_list);
     3425
     3426        /*Recover parameters and values*/
     3427        GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
     3428
     3429        /*Be sure that values are not zero*/
     3430        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     3431        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     3432        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     3433
     3434        int pt1 = 0;
     3435        int pt2 = 1;
     3436        if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     3437
     3438                /*Portion of the segments*/
     3439                s1=gl[2]/(gl[2]-gl[1]);
     3440                s2=gl[2]/(gl[2]-gl[0]);
     3441                if(gl[2]<0.){
     3442                        pt1 = 1; pt2 = 0;
     3443                }
     3444                xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
     3445                xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
     3446                xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
     3447                xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
     3448                xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
     3449                xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
     3450        }
     3451        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
     3452
     3453                /*Portion of the segments*/
     3454                s1=gl[0]/(gl[0]-gl[1]);
     3455                s2=gl[0]/(gl[0]-gl[2]);
     3456                if(gl[0]<0.){
     3457                        pt1 = 1; pt2 = 0;
     3458                }
     3459
     3460                xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
     3461                xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
     3462                xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
     3463                xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
     3464                xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
     3465                xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
     3466        }
     3467        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
     3468
     3469                /*Portion of the segments*/
     3470                s1=gl[1]/(gl[1]-gl[0]);
     3471                s2=gl[1]/(gl[1]-gl[2]);
     3472                if(gl[1]<0.){
     3473                        pt1 = 1; pt2 = 0;
     3474                }
     3475
     3476                xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
     3477                xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
     3478                xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
     3479                xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
     3480                xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
     3481                xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
     3482        }
     3483        else{
     3484                _error_("case not possible");
     3485        }
     3486
     3487
     3488        /*Some checks in debugging mode*/
     3489        _assert_(s1>=0 && s1<=1.);
     3490        _assert_(s2>=0 && s2<=1.);
     3491
     3492        /*Get normal vector*/
     3493        IssmDouble normal[3];
     3494        this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
     3495        normal[0] = -normal[0];
     3496        normal[1] = -normal[1];
     3497
     3498        /*Get inputs*/
     3499        IssmDouble flux = 0.;
     3500        IssmDouble calvingratex,calvingratey,thickness,Jdet;
     3501        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     3502        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     3503        Input* calvingratex_input=NULL;
     3504        Input* calvingratey_input=NULL;
     3505        calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     3506        calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     3507
     3508        /*Start looping on Gaussian points*/
     3509        Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
     3510        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3511
     3512                gauss->GaussPoint(ig);
     3513                thickness_input->GetInputValue(&thickness,gauss);
     3514                calvingratex_input->GetInputValue(&calvingratex,gauss);
     3515                calvingratey_input->GetInputValue(&calvingratey,gauss);
     3516                this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
     3517
     3518                flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
     3519        }
     3520
     3521        return flux;
     3522               
     3523        /*Clean up and return*/
     3524        delete gauss;
     3525}
     3526/*}}}*/
     3527IssmDouble Penta::TotalCalvingMeltingFluxLevelset(bool scaled){/*{{{*/
     3528
     3529        /*Make sure there is an ice front here*/
     3530        if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum) || !IsOnBase()) return 0;
     3531
     3532        /*Scaled not implemented yet...*/
     3533        _assert_(!scaled);
     3534
     3535        int               domaintype,index1,index2;
     3536        const IssmPDouble epsilon = 1.e-15;
     3537        IssmDouble        s1,s2;
     3538        IssmDouble        gl[NUMVERTICES];
     3539        IssmDouble        xyz_front[2][3];
     3540
     3541        IssmDouble *xyz_list = NULL;
     3542        this->GetVerticesCoordinates(&xyz_list);
     3543
     3544        /*Recover parameters and values*/
     3545        GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
     3546
     3547        /*Be sure that values are not zero*/
     3548        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     3549        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     3550        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     3551
     3552        int pt1 = 0;
     3553        int pt2 = 1;
     3554        if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     3555
     3556                /*Portion of the segments*/
     3557                s1=gl[2]/(gl[2]-gl[1]);
     3558                s2=gl[2]/(gl[2]-gl[0]);
     3559                if(gl[2]<0.){
     3560                        pt1 = 1; pt2 = 0;
     3561                }
     3562                xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
     3563                xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
     3564                xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
     3565                xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
     3566                xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
     3567                xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
     3568        }
     3569        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
     3570
     3571                /*Portion of the segments*/
     3572                s1=gl[0]/(gl[0]-gl[1]);
     3573                s2=gl[0]/(gl[0]-gl[2]);
     3574                if(gl[0]<0.){
     3575                        pt1 = 1; pt2 = 0;
     3576                }
     3577
     3578                xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
     3579                xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
     3580                xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
     3581                xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
     3582                xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
     3583                xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
     3584        }
     3585        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
     3586
     3587                /*Portion of the segments*/
     3588                s1=gl[1]/(gl[1]-gl[0]);
     3589                s2=gl[1]/(gl[1]-gl[2]);
     3590                if(gl[1]<0.){
     3591                        pt1 = 1; pt2 = 0;
     3592                }
     3593
     3594                xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
     3595                xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
     3596                xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
     3597                xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
     3598                xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
     3599                xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
     3600        }
     3601        else{
     3602                _error_("case not possible");
     3603        }
     3604
     3605
     3606        /*Some checks in debugging mode*/
     3607        _assert_(s1>=0 && s1<=1.);
     3608        _assert_(s2>=0 && s2<=1.);
     3609
     3610        /*Get normal vector*/
     3611        IssmDouble normal[3];
     3612        this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
     3613        normal[0] = -normal[0];
     3614        normal[1] = -normal[1];
     3615       
     3616        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
     3617        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
     3618
     3619        /*Get inputs*/
     3620        IssmDouble flux = 0.;
     3621        IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet;
     3622        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     3623        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     3624        Input* calvingratex_input=NULL;
     3625        Input* calvingratey_input=NULL;
     3626        Input* vx_input=NULL;
     3627        Input* vy_input=NULL;
     3628        Input* meltingrate_input=NULL;
     3629        calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     3630        calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     3631        vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
     3632        vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
     3633        meltingrate_input=inputs->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input);
     3634
     3635        /*Start looping on Gaussian points*/
     3636        Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
     3637        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3638
     3639                gauss->GaussPoint(ig);
     3640                thickness_input->GetInputValue(&thickness,gauss);
     3641                calvingratex_input->GetInputValue(&calvingratex,gauss);
     3642                calvingratey_input->GetInputValue(&calvingratey,gauss);
     3643                vx_input->GetInputValue(&vx,gauss);
     3644                vy_input->GetInputValue(&vy,gauss);
     3645                vel=vx*vx+vy*vy;
     3646                meltingrate_input->GetInputValue(&meltingrate,gauss);   
     3647                meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
     3648                meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
     3649                this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
     3650
     3651                flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
     3652        }
     3653
     3654        return flux;
     3655       
     3656        /*Clean up and return*/
     3657        delete gauss;
     3658}
     3659/*}}}*/
    30083660IssmDouble Penta::TotalFloatingBmb(bool scaled){/*{{{*/
    30093661

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

    r24089 r24119  
    5050                void           CalvingRateVonmises();
    5151                void           CalvingRateLevermann();
     52                void           CalvingFluxLevelset();
     53                void           CalvingMeltingFluxLevelset();
    5254                IssmDouble     CharacteristicLength(void){_error_("not implemented yet");};
    5355                void           ComputeBasalStress(void);
    … …  
    9395                void           GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
    9496                IssmDouble     GroundedArea(bool scaled);
     97                IssmDouble              IcefrontMassFluxLevelset(bool scaled);
    9598                IssmDouble     IceVolume(bool scaled);
    9699                IssmDouble     IceVolumeAboveFloatation(bool scaled);
    … …  
    122125                Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
    123126                Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
     127                Gauss*         NewGaussBase(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz);
    124128                Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order);
    125129                Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order){_error_("not implemented yet");};
    … …  
    139143                void             NormalBase(IssmDouble* bed_normal, IssmDouble* xyz_list);
    140144                void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
     145                void           NormalSectionBase(IssmDouble* normal,IssmDouble* xyz_list);
    141146                void             NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list);
    142147                int            NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
    … …  
    163168                int            TensorInterpolation(){_error_("not implemented yet");};
    164169                IssmDouble     TimeAdapt();
     170                IssmDouble              TotalCalvingFluxLevelset(bool scaled);
     171                IssmDouble              TotalCalvingMeltingFluxLevelset(bool scaled);
    165172                IssmDouble     TotalFloatingBmb(bool scaled);
    166173                IssmDouble     TotalGroundedBmb(bool scaled);

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

    r24118 r24119  
    492492        }
    493493        else{
    494         int               domaintype,index1,index2;
    495         const IssmPDouble epsilon = 1.e-15;
    496         IssmDouble        s1,s2;
    497         IssmDouble        gl[NUMVERTICES];
    498         IssmDouble        xyz_front[2][3];
    499 
    500 
    501         IssmDouble *xyz_list = NULL;
    502         this->GetVerticesCoordinates(&xyz_list);
    503 
    504         /*Recover parameters and values*/
    505         parameters->FindParam(&domaintype,DomainTypeEnum);
    506         GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
    507 
    508         /*Be sure that values are not zero*/
    509         if(gl[0]==0.) gl[0]=gl[0]+epsilon;
    510         if(gl[1]==0.) gl[1]=gl[1]+epsilon;
    511         if(gl[2]==0.) gl[2]=gl[2]+epsilon;
    512 
    513         if(domaintype==Domain2DverticalEnum){
    514                 _error_("not implemented");
    515         }
    516         else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
    517                 int pt1 = 0;
    518                 int pt2 = 1;
    519                 if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
    520 
    521                         /*Portion of the segments*/
    522                         s1=gl[2]/(gl[2]-gl[1]);
    523                         s2=gl[2]/(gl[2]-gl[0]);
    524                         if(gl[2]<0.){
    525                                 pt1 = 1; pt2 = 0;
     494                int               domaintype,index1,index2;
     495                const IssmPDouble epsilon = 1.e-15;
     496                IssmDouble        s1,s2;
     497                IssmDouble        gl[NUMVERTICES];
     498                IssmDouble        xyz_front[2][3];
     499
     500                IssmDouble *xyz_list = NULL;
     501                this->GetVerticesCoordinates(&xyz_list);
     502
     503                /*Recover parameters and values*/
     504                parameters->FindParam(&domaintype,DomainTypeEnum);
     505                GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
     506
     507                /*Be sure that values are not zero*/
     508                if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     509                if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     510                if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     511
     512                if(domaintype==Domain2DverticalEnum){
     513                        _error_("not implemented");
     514                }
     515                else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
     516                        int pt1 = 0;
     517                        int pt2 = 1;
     518                        if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     519
     520                                /*Portion of the segments*/
     521                                s1=gl[2]/(gl[2]-gl[1]);
     522                                s2=gl[2]/(gl[2]-gl[0]);
     523                                if(gl[2]<0.){
     524                                        pt1 = 1; pt2 = 0;
     525                                }
     526                                xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
     527                                xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
     528                                xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
     529                                xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
     530                                xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
     531                                xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
    526532                        }
    527                         xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
    528                         xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
    529                         xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
    530                         xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
    531                         xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
    532                         xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
    533                 }
    534                 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
    535 
    536                         /*Portion of the segments*/
    537                         s1=gl[0]/(gl[0]-gl[1]);
    538                         s2=gl[0]/(gl[0]-gl[2]);
    539                         if(gl[0]<0.){
    540                                 pt1 = 1; pt2 = 0;
     533                        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
     534
     535                                /*Portion of the segments*/
     536                                s1=gl[0]/(gl[0]-gl[1]);
     537                                s2=gl[0]/(gl[0]-gl[2]);
     538                                if(gl[0]<0.){
     539                                        pt1 = 1; pt2 = 0;
     540                                }
     541
     542                                xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
     543                                xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
     544                                xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
     545                                xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
     546                                xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
     547                                xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
    541548                        }
    542 
    543                         xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
    544                         xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
    545                         xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
    546                         xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
    547                         xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
    548                         xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
    549                 }
    550                 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
    551 
    552                         /*Portion of the segments*/
    553                         s1=gl[1]/(gl[1]-gl[0]);
    554                         s2=gl[1]/(gl[1]-gl[2]);
    555                         if(gl[1]<0.){
    556                                 pt1 = 1; pt2 = 0;
     549                        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
     550
     551                                /*Portion of the segments*/
     552                                s1=gl[1]/(gl[1]-gl[0]);
     553                                s2=gl[1]/(gl[1]-gl[2]);
     554                                if(gl[1]<0.){
     555                                        pt1 = 1; pt2 = 0;
     556                                }
     557
     558                                xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
     559                                xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
     560                                xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
     561                                xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
     562                                xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
     563                                xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
    557564                        }
    558 
    559                         xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
    560                         xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
    561                         xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
    562                         xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
    563                         xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
    564                         xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
     565                        else{
     566                                _error_("case not possible");
     567                        }
     568
     569                }
     570                else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
     571
     572                /*Some checks in debugging mode*/
     573                _assert_(s1>=0 && s1<=1.);
     574                _assert_(s2>=0 && s2<=1.);
     575
     576                /*Get normal vector*/
     577                IssmDouble normal[3];
     578                this->NormalSection(&normal[0],&xyz_front[0][0]);
     579                normal[0] = -normal[0];
     580                normal[1] = -normal[1];
     581
     582                /*Get inputs*/
     583                IssmDouble flux = 0.;
     584                IssmDouble area = 0.;
     585                IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
     586                IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     587                Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     588                Input* calvingratex_input=NULL;
     589                Input* calvingratey_input=NULL;
     590                if(domaintype==Domain2DhorizontalEnum){
     591                        calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     592                        calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    565593                }
    566594                else{
    567                         _error_("case not possible");
    568                 }
    569 
    570         }
    571         else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
    572 
    573         /*Some checks in debugging mode*/
    574         _assert_(s1>=0 && s1<=1.);
    575         _assert_(s2>=0 && s2<=1.);
    576 
    577         /*Get normal vector*/
    578         IssmDouble normal[3];
    579         this->NormalSection(&normal[0],&xyz_front[0][0]);
    580         normal[0] = -normal[0];
    581         normal[1] = -normal[1];
    582 
    583         /*Get inputs*/
    584         IssmDouble flux = 0.;
    585         IssmDouble area = 0.;
    586         IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
    587         IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
    588         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    589         Input* calvingratex_input=NULL;
    590         Input* calvingratey_input=NULL;
    591         if(domaintype==Domain2DhorizontalEnum){
    592                 calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    593                 calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    594         }
    595         else{
    596                 calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
    597                 calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
    598         }
    599 
    600         /*Start looping on Gaussian points*/
    601         Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
    602         for(int ig=gauss->begin();ig<gauss->end();ig++){
    603 
    604                 gauss->GaussPoint(ig);
    605                 thickness_input->GetInputValue(&thickness,gauss);
    606                 calvingratex_input->GetInputValue(&calvingratex,gauss);
    607                 calvingratey_input->GetInputValue(&calvingratey,gauss);
    608                 this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
    609 
    610                 flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
    611                 area += Jdet*gauss->weight*thickness;
    612 
    613                 flux_per_area=flux/area;
    614         }
    615        
    616         this->inputs->AddInput(new TriaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
     595                        calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     596                        calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     597                }
     598
     599                /*Start looping on Gaussian points*/
     600                Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
     601                for(int ig=gauss->begin();ig<gauss->end();ig++){
     602
     603                        gauss->GaussPoint(ig);
     604                        thickness_input->GetInputValue(&thickness,gauss);
     605                        calvingratex_input->GetInputValue(&calvingratex,gauss);
     606                        calvingratey_input->GetInputValue(&calvingratey,gauss);
     607                        this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
     608
     609                        flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
     610                        area += Jdet*gauss->weight*thickness;
     611
     612                        flux_per_area=flux/area;
     613                }
     614
     615                this->inputs->AddInput(new TriaInput(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum));
     616               
     617                /*Clean up and return*/
     618                delete gauss;
    617619        }
    618620}
    … …  
    626628        }
    627629        else{
    628         int               domaintype,index1,index2;
    629         const IssmPDouble epsilon = 1.e-15;
    630         IssmDouble        s1,s2;
    631         IssmDouble        gl[NUMVERTICES];
    632         IssmDouble        xyz_front[2][3];
    633 
    634 
    635         IssmDouble *xyz_list = NULL;
    636         this->GetVerticesCoordinates(&xyz_list);
    637 
    638         /*Recover parameters and values*/
    639         parameters->FindParam(&domaintype,DomainTypeEnum);
    640         GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
    641 
    642         /*Be sure that values are not zero*/
    643         if(gl[0]==0.) gl[0]=gl[0]+epsilon;
    644         if(gl[1]==0.) gl[1]=gl[1]+epsilon;
    645         if(gl[2]==0.) gl[2]=gl[2]+epsilon;
    646 
    647         if(domaintype==Domain2DverticalEnum){
    648                 _error_("not implemented");
    649         }
    650         else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
    651                 int pt1 = 0;
    652                 int pt2 = 1;
    653                 if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
    654 
    655                         /*Portion of the segments*/
    656                         s1=gl[2]/(gl[2]-gl[1]);
    657                         s2=gl[2]/(gl[2]-gl[0]);
    658                         if(gl[2]<0.){
    659                                 pt1 = 1; pt2 = 0;
     630                int               domaintype,index1,index2;
     631                const IssmPDouble epsilon = 1.e-15;
     632                IssmDouble        s1,s2;
     633                IssmDouble        gl[NUMVERTICES];
     634                IssmDouble        xyz_front[2][3];
     635
     636
     637                IssmDouble *xyz_list = NULL;
     638                this->GetVerticesCoordinates(&xyz_list);
     639
     640                /*Recover parameters and values*/
     641                parameters->FindParam(&domaintype,DomainTypeEnum);
     642                GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
     643
     644                /*Be sure that values are not zero*/
     645                if(gl[0]==0.) gl[0]=gl[0]+epsilon;
     646                if(gl[1]==0.) gl[1]=gl[1]+epsilon;
     647                if(gl[2]==0.) gl[2]=gl[2]+epsilon;
     648
     649                if(domaintype==Domain2DverticalEnum){
     650                        _error_("not implemented");
     651                }
     652                else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
     653                        int pt1 = 0;
     654                        int pt2 = 1;
     655                        if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     656
     657                                /*Portion of the segments*/
     658                                s1=gl[2]/(gl[2]-gl[1]);
     659                                s2=gl[2]/(gl[2]-gl[0]);
     660                                if(gl[2]<0.){
     661                                        pt1 = 1; pt2 = 0;
     662                                }
     663                                xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
     664                                xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
     665                                xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
     666                                xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
     667                                xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
     668                                xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
    660669                        }
    661                         xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
    662                         xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
    663                         xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
    664                         xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
    665                         xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
    666                         xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
    667                 }
    668                 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
    669 
    670                         /*Portion of the segments*/
    671                         s1=gl[0]/(gl[0]-gl[1]);
    672                         s2=gl[0]/(gl[0]-gl[2]);
    673                         if(gl[0]<0.){
    674                                 pt1 = 1; pt2 = 0;
     670                        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
     671
     672                                /*Portion of the segments*/
     673                                s1=gl[0]/(gl[0]-gl[1]);
     674                                s2=gl[0]/(gl[0]-gl[2]);
     675                                if(gl[0]<0.){
     676                                        pt1 = 1; pt2 = 0;
     677                                }
     678
     679                                xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
     680                                xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
     681                                xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
     682                                xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
     683                                xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
     684                                xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
    675685                        }
    676 
    677                         xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
    678                         xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
    679                         xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
    680                         xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
    681                         xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
    682                         xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
    683                 }
    684                 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
    685 
    686                         /*Portion of the segments*/
    687                         s1=gl[1]/(gl[1]-gl[0]);
    688                         s2=gl[1]/(gl[1]-gl[2]);
    689                         if(gl[1]<0.){
    690                                 pt1 = 1; pt2 = 0;
     686                        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
     687
     688                                /*Portion of the segments*/
     689                                s1=gl[1]/(gl[1]-gl[0]);
     690                                s2=gl[1]/(gl[1]-gl[2]);
     691                                if(gl[1]<0.){
     692                                        pt1 = 1; pt2 = 0;
     693                                }
     694
     695                                xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
     696                                xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
     697                                xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
     698                                xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
     699                                xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
     700                                xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
    691701                        }
    692 
    693                         xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
    694                         xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
    695                         xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
    696                         xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
    697                         xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
    698                         xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
     702                        else{
     703                                _error_("case not possible");
     704                        }
     705
     706                }
     707                else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
     708
     709                /*Some checks in debugging mode*/
     710                _assert_(s1>=0 && s1<=1.);
     711                _assert_(s2>=0 && s2<=1.);
     712
     713                /*Get normal vector*/
     714                IssmDouble normal[3];
     715                this->NormalSection(&normal[0],&xyz_front[0][0]);
     716                normal[0] = -normal[0];
     717                normal[1] = -normal[1];
     718
     719                /*Get inputs*/
     720                IssmDouble flux = 0.;
     721                IssmDouble area = 0.;
     722                IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
     723                IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
     724                Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     725                Input* calvingratex_input=NULL;
     726                Input* calvingratey_input=NULL;
     727                Input* vx_input=NULL;
     728                Input* vy_input=NULL;
     729                Input* meltingrate_input=NULL;
     730                if(domaintype==Domain2DhorizontalEnum){
     731                        calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
     732                        calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
     733                        vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     734                        vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     735                        meltingrate_input=inputs->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input);
    699736                }
    700737                else{
    701                         _error_("case not possible");
    702                 }
    703 
    704         }
    705         else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
    706 
    707         /*Some checks in debugging mode*/
    708         _assert_(s1>=0 && s1<=1.);
    709         _assert_(s2>=0 && s2<=1.);
    710 
    711         /*Get normal vector*/
    712         IssmDouble normal[3];
    713         this->NormalSection(&normal[0],&xyz_front[0][0]);
    714         normal[0] = -normal[0];
    715         normal[1] = -normal[1];
    716 
    717         /*Get inputs*/
    718         IssmDouble flux = 0.;
    719         IssmDouble area = 0.;
    720         IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
    721         IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
    722         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    723         Input* calvingratex_input=NULL;
    724         Input* calvingratey_input=NULL;
    725         Input* vx_input=NULL;
    726         Input* vy_input=NULL;
    727         Input* meltingrate_input=NULL;
    728         if(domaintype==Domain2DhorizontalEnum){
    729                 calvingratex_input=inputs->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    730                 calvingratey_input=inputs->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    731                 vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    732                 vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    733                 meltingrate_input=inputs->GetInput(CalvingMeltingrateEnum); _assert_(meltingrate_input);
    734         }
    735         else{
    736                 calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
    737                 calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
    738         }
    739 
    740         /*Start looping on Gaussian points*/
    741         Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
    742         for(int ig=gauss->begin();ig<gauss->end();ig++){
    743 
    744                 gauss->GaussPoint(ig);
    745                 thickness_input->GetInputValue(&thickness,gauss);
    746                 calvingratex_input->GetInputValue(&calvingratex,gauss);
    747                 calvingratey_input->GetInputValue(&calvingratey,gauss);
    748                 vx_input->GetInputValue(&vx,gauss);
    749                 vy_input->GetInputValue(&vy,gauss);
    750                 vel=vx*vx+vy*vy;
    751                 meltingrate_input->GetInputValue(&meltingrate,gauss);   
    752                 meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
    753                 meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
    754                 this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
    755 
    756                 flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
    757                 area += Jdet*gauss->weight*thickness;
    758 
    759                 flux_per_area=flux/area;
    760         }
    761        
    762         this->inputs->AddInput(new TriaInput(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum));
     738                        calvingratex_input=inputs->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
     739                        calvingratey_input=inputs->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
     740                }
     741
     742                /*Start looping on Gaussian points*/
     743                Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
     744                for(int ig=gauss->begin();ig<gauss->end();ig++){
     745
     746                        gauss->GaussPoint(ig);
     747                        thickness_input->GetInputValue(&thickness,gauss);
     748                        calvingratex_input->GetInputValue(&calvingratex,gauss);
     749                        calvingratey_input->GetInputValue(&calvingratey,gauss);
     750                        vx_input->GetInputValue(&vx,gauss);
     751                        vy_input->GetInputValue(&vy,gauss);
     752                        vel=vx*vx+vy*vy;
     753                        meltingrate_input->GetInputValue(&meltingrate,gauss);   
     754                        meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
     755                        meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
     756                        this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);
     757
     758                        flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
     759                        area += Jdet*gauss->weight*thickness;
     760
     761                        flux_per_area=flux/area;
     762                }
     763
     764                this->inputs->AddInput(new TriaInput(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum));
     765               
     766                /*Clean up and return*/
     767                delete gauss;
    763768        }
    764769}
    … …  
    22992304                flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
    23002305        }
    2301 
     2306        delete gauss;
    23022307        return flux;
    23032308}
    … …  
    24282433                flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
    24292434        }
    2430 
     2435        delete gauss;
    24312436        return flux;
    24322437}
    … …  
    36733678                TF_input->GetInputValue(&TF,gauss);
    36743679
    3675                 if(basinid[iv]==0 || basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
     3680                if(basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
    36763681                else{
    36773682                        /* change the unit of qsg (m^3/d -> m/d) with ice front area */
    … …  
    36793684
    36803685                        /* calculate melt rates */
    3681                         meltrates[iv]=(A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta)/86400; //[m/s]
     3686                        meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s]
    36823687                }       
    36833688
    … …  
    42574262        /*Get inputs*/
    42584263        IssmDouble flux = 0.;
    4259         IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
     4264        IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet;
    42604265        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
    42614266        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

    r22075 r24119  
    433433}
    434434/*}}}*/
     435GaussPenta::GaussPenta(IssmDouble area_coordinates[2][3],int order_horiz){/*{{{*/
     436
     437        /*Intermediaties*/
     438        IssmPDouble *seg_horiz_coords  = NULL;
     439        IssmPDouble *seg_horiz_weights = NULL;
     440
     441        /*get the gauss points using the product of two line rules*/
     442        GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
     443
     444        /*Allocate GaussPenta fields*/
     445        numgauss=order_horiz;
     446        coords1=xNew<IssmDouble>(numgauss);
     447        coords2=xNew<IssmDouble>(numgauss);
     448        coords3=xNew<IssmDouble>(numgauss);
     449        coords4=xNew<IssmDouble>(numgauss);
     450        weights=xNew<IssmDouble>(numgauss);
     451
     452        /*Quads: get the gauss points using the product of two line rules  */
     453        for(int i=0;i<order_horiz;i++){
     454                coords1[i]=0.5*(area_coordinates[0][0]+area_coordinates[1][0]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][0]-area_coordinates[0][0]);
     455                coords2[i]=0.5*(area_coordinates[0][1]+area_coordinates[1][1]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][1]-area_coordinates[0][1]);
     456                coords3[i]=0.5*(area_coordinates[0][2]+area_coordinates[1][2]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][2]-area_coordinates[0][2]);
     457                coords4[i]=0.;
     458                weights[i]=seg_horiz_weights[i];
     459        }
     460
     461        /*clean-up*/
     462        xDelete<IssmPDouble>(seg_horiz_coords);
     463        xDelete<IssmPDouble>(seg_horiz_weights);
     464}
     465/*}}}*/
    435466GaussPenta::~GaussPenta(){/*{{{*/
    436467        xDelete<IssmDouble>(weights);

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h

    r20827 r24119  
    3737                GaussPenta(int index,IssmDouble r1, IssmDouble r2,bool maintlyfloating,int order);
    3838                GaussPenta(IssmDouble area_coordinates[4][3],int order_horiz,int order_vert);
     39                GaussPenta(IssmDouble area_coordinates[2][3],int order_horiz);
    3940                ~GaussPenta();
    4041
Note: See TracChangeset for help on using the changeset viewer.