source: issm/oecreview/Archive/25834-26739/ISSM-26652-26653.diff

Last change on this file was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 9.4 KB

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

     
    12951295        else                  phi = f1*f2;
    12961296       
    12971297        /*Compute weights*/
    1298         gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
     1298        gauss = this->NewGauss(point,f1,f2,1-trapezeisnegative,2);
    12991299
    13001300        total_weight = 0.0;
    13011301        for(int i=0;i<NUMVERTICES2D;i++)weights[i] = 0;
    … …  
    13051305                total_weight += gauss->weight;
    13061306        }
    13071307
    1308         /*Normalize to phi*/
    1309         if(total_weight>0.){
    1310                 for(int i=0;i<NUMVERTICES2D;i++) weights[i] = weights[i]*phi/total_weight;
    1311         }
     1308        /*Normalizing to phi such that weights provide coefficients for integration over subelement (for averaging:phi*weights)*/
     1309   if(total_weight>0.) for(int i=0;i<NUMVERTICES;i++) weights[i] = weights[i]*phi/total_weight;
    13121310        else for(int i=0;i<NUMVERTICES2D;i++) weights[i] = 0.0;
    13131311
    13141312        /*Cleanup*/
    … …  
    22062204               
    22072205                IssmDouble basetot;
    22082206                height = 0.0;
    2209                 for(int i=0;i<NUMVERTICES2D;i++) height += weights[i]*heights[i];
     2207                for(int i=0;i<NUMVERTICES2D;i++) height += weights[i]/phi*heights[i];
    22102208                basetot = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    22112209                base    = basetot*phi; 
    22122210       
    … …  
    22162214                        Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
    22172215                        /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
    22182216                        scalefactor = 0.0;
    2219                         for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]*scalefactor_vertices[i];
     2217                        for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
    22202218                        base = base*scalefactor;
    22212219                        xDelete<IssmDouble>(scalefactor_vertices);
    22222220                }
    … …  
    43884386        /*The smb[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
    43894387        IssmDouble base,smb,rho_ice,scalefactor;
    43904388        IssmDouble Total_Smb=0;
     4389        IssmDouble lsf[NUMVERTICES];
    43914390        IssmDouble xyz_list[NUMVERTICES][3];
    43924391
    43934392        /*Get material parameters :*/
    … …  
    44034402        base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    44044403
    44054404        /*Now get the average SMB over the element*/
    4406         Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
     4405        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     4406   if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
     4407                /*Partially ice-covered element*/
     4408                bool mainlyice;
     4409      int point;
     4410      IssmDouble* smb_vertices = xNew<IssmDouble>(NUMVERTICES);
     4411                IssmDouble  weights[NUMVERTICES2D];
     4412                IssmDouble  lsf2d[NUMVERTICES2D];
     4413      IssmDouble f1,f2,phi;
     4414      Element::GetInputListOnVertices(&smb_vertices[0],SmbMassBalanceEnum);
     4415                for(int i=0;i<NUMVERTICES2D;i++) lsf2d[i] = lsf[i];
     4416                GetFractionGeometry2D(weights,&phi,&point,&f1,&f2,&mainlyice,lsf2d);
     4417                smb = 0.0;
     4418                for(int i=0;i<NUMVERTICES2D;i++) smb += weights[i]*smb_vertices[i];
    44074419
    4408         smb_input->GetInputAverage(&smb);
    4409         if(scaled==true){
    4410                 Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
    4411                 scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     4420                if(scaled==true){
     4421         IssmDouble* scalefactor_vertices   = xNew<IssmDouble>(NUMVERTICES);
     4422         Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
     4423         /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
     4424         scalefactor = 0.0;
     4425         for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
     4426         xDelete<IssmDouble>(scalefactor_vertices);
     4427                }
     4428                else scalefactor = 1.0;
     4429
     4430                /*Cleanup*/
     4431      xDelete<IssmDouble>(smb_vertices);
    44124432        }
     4433
    44134434        else{
    4414                 scalefactor=1.;
     4435                /*Fully ice-covered element*/
     4436                Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
     4437                smb_input->GetInputAverage(&smb);
     4438
     4439                if(scaled==true){
     4440                        Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     4441                        scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     4442                }
     4443                else scalefactor=1.0;
    44154444        }
     4445
    44164446        Total_Smb=rho_ice*base*smb*scalefactor;// smb on element in kg s-1
    44174447
    44184448        /*Return: */

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

     
    18161816        else phi=f1*f2;
    18171817
    18181818        /*Compute weights:*/
    1819         gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
     1819        gauss = this->NewGauss(point,f1,f2,1-trapezeisnegative,2); //VV correction (16Nov2021)
    18201820
    18211821        total_weight=0;
    18221822        for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
    … …  
    18251825                for(int i=0;i<NUMVERTICES;i++)weights[i]+=loadweights_g[i]*gauss->weight;
    18261826                total_weight+=gauss->weight;
    18271827        }
    1828         //normalize to phi.
     1828        /*Normalizing to phi such that weights provide coefficients for integration over subelement (for averaging:phi*weights)*/
    18291829        if(total_weight>0.) for(int i=0;i<NUMVERTICES;i++)weights[i]/=total_weight/phi;
    18301830        else for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
    18311831
    … …  
    35453545                Haverage/=IssmDouble(numthk);
    35463546        }*/
    35473547
    3548    IssmDouble* lsf = xNew<IssmDouble>(NUMVERTICES);
    3549    Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     3548   IssmDouble lsf[NUMVERTICES];
     3549        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
    35503550   /*Deal with partially ice-covered elements*/
    35513551        if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
    35523552                bool istrapneg;
    … …  
    35623562      GetFractionGeometry(weights,&phi,&point,&f1,&f2,&istrapneg,lsf);
    35633563      for(int i=0;i<NUMVERTICES;i++) Hice[i] = surfaces[i]-bases[i];
    35643564      Haverage = 0.0;
    3565       for(int i=0;i<NUMVERTICES;i++) Haverage += weights[i]*Hice[i];
     3565                /*Use weights[i]/phi to get average thickness over subelement*/
     3566      for(int i=0;i<NUMVERTICES;i++) Haverage += weights[i]/phi*Hice[i];
    35663567                /*Get back area of ice-covered base*/
    35673568                area_basetot = this->GetArea();
    35683569                area_base    = phi*area_basetot;
    … …  
    35723573                        IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
    35733574                        Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
    35743575                        scalefactor = 0.0;
    3575                         for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]*scalefactor_vertices[i];
     3576                        for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
    35763577                        area_base = area_base*scalefactor;
    35773578                        xDelete<IssmDouble>(scalefactor_vertices);
    35783579                }
    … …  
    36053606        xDelete<int>(indices);
    36063607        xDelete<IssmDouble>(H);
    36073608        xDelete<IssmDouble>(SF);
    3608         xDelete<IssmDouble>(lsf);
    36093609
    36103610        if(domaintype==Domain2DverticalEnum){
    36113611          return area_base;
    … …  
    54785478        /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
    54795479        IssmDouble base,smb,rho_ice,scalefactor;
    54805480        IssmDouble Total_Smb=0;
     5481        IssmDouble lsf[NUMVERTICES];
    54815482        IssmDouble xyz_list[NUMVERTICES][3];
    54825483
    54835484        /*Get material parameters :*/
    … …  
    54915492         * http://en.wikipedia.org/wiki/Triangle
    54925493         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    54935494        base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2
     5495       
     5496        /*Now get the average SMB over the element*/
     5497        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     5498        if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
     5499                /*Partially ice-covered element*/
     5500                bool mainlyice;
     5501      int point;
     5502      IssmDouble* weights       = xNew<IssmDouble>(NUMVERTICES);
     5503      IssmDouble* smb_vertices  = xNew<IssmDouble>(NUMVERTICES);
     5504      IssmDouble f1,f2,phi;
     5505               
     5506                Element::GetInputListOnVertices(&smb_vertices[0],SmbMassBalanceEnum);
     5507                GetFractionGeometry(weights,&phi,&point,&f1,&f2,&mainlyice,lsf);
     5508                smb = 0.0;
     5509                for(int i=0;i<NUMVERTICES;i++) smb += weights[i]*smb_vertices[i];
     5510       
     5511                if(scaled==true){
     5512         IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
     5513         Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
     5514         scalefactor = 0.0;
     5515         for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
     5516         xDelete<IssmDouble>(scalefactor_vertices);
     5517      }
     5518                else scalefactor = 1.0;
    54945519
    5495         /*Now get the average SMB over the element*/
    5496         Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
    5497         smb_input->GetInputAverage(&smb);       // average smb on element in m ice s-1
    5498         if(scaled==true){
    5499                 Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
    5500                 scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     5520                /*Cleanup*/
     5521      xDelete<IssmDouble>(weights);
     5522      xDelete<IssmDouble>(smb_vertices);
    55015523        }
    55025524        else{
    5503                 scalefactor=1.;
     5525                /*Fully ice-covered element*/
     5526                Input* smb_input = this->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
     5527                smb_input->GetInputAverage(&smb);   // average smb on element in m ice s-1
     5528
     5529                if(scaled==true){
     5530                        Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     5531                        scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     5532                }
     5533                else scalefactor=1.0;
    55045534        }
     5535       
    55055536   Total_Smb=rho_ice*base*smb*scalefactor;      // smb on element in kg s-1
    55065537
    55075538        /*Return: */
Note: See TracBrowser for help on using the repository browser.