Changeset 18946


Ignore:
Timestamp:
12/04/14 21:24:24 (10 years ago)
Author:
seroussi
Message:

CHG: removed useless functions, call GetMaterialParameter instead

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

Legend:

Unmodified
Added
Removed

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

    r18920 r18946  
    12261226        GetInputListOnVertices(&r[0],BedEnum);
    12271227        GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
    1228         rho_water   = matpar->GetRhoWater();
    1229         rho_ice     = matpar->GetRhoIce();
     1228        rho_water   = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     1229        rho_ice     = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    12301230        density     = rho_ice/rho_water;
    12311231
    … …  
    20542054        surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
    20552055        z=this->GetZcoord(xyz_list,gauss);
    2056         tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
     2056        tau_perp = matpar->GetMaterialParameter(MaterialsRhoIceEnum) * matpar->GetMaterialParameter(ConstantsGEnum) * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
    20572057
    20582058        /* Get eps_b*/

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

    r18921 r18946  
    252252
    253253        /*recovre material parameters: */
    254         rho_ice=matpar->GetRhoIce();
    255         gravity=matpar->GetG();
     254        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     255        gravity=matpar->GetMaterialParameter(ConstantsGEnum);
    256256
    257257        /* Get node coordinates and dof list: */
    … …  
    697697
    698698                        /*Compute water pressure*/
    699                         IssmDouble rho_ice   = matpar->GetRhoIce();
    700                         IssmDouble rho_water = matpar->GetRhoWater();
    701                         IssmDouble gravity   = matpar->GetG();
     699                        IssmDouble rho_ice   = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     700                        IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     701                        IssmDouble gravity   = matpar->GetMaterialParameter(ConstantsGEnum);
    702702                        water_pressure=gravity*rho_water*base;
    703703
    … …  
    11771177        if(!IsIceInElement() || IsFloating() || !IsOnBase())return 0;
    11781178
    1179         rho_ice=matpar->GetRhoIce();
    1180         rho_water=matpar->GetRhoWater();
     1179        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     1180        rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    11811181        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    11821182
    … …  
    20492049
    20502050  /*Get material parameters :*/
    2051   rho_ice=matpar->GetRhoIce();
    2052   rho_water=matpar->GetRhoFreshwater();
     2051  rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     2052  rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    20532053
    20542054  /*Get other pdd parameters*/
    2055   desfac=matpar->GetDesFac();
    2056   s0p=matpar->GetS0p();
    2057   s0t=matpar->GetS0t();
    2058   rlaps=matpar->GetRlaps();
    2059   rlapslgm=matpar->GetRlapslgm();
     2055  desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
     2056  s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
     2057  s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
     2058  rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
     2059  rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
    20602060
    20612061   /*measure the surface mass balance*/
    … …  
    20832083
    20842084        /*material parameters: */
    2085         rho_water=matpar->GetRhoWater();
    2086         rho_ice=matpar->GetRhoIce();
     2085        rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     2086        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    20872087        density=rho_ice/rho_water;
    20882088        GetInputListOnVertices(&h[0],ThicknessEnum);
    … …  
    26412641
    26422642        /*Get material parameters :*/
    2643         rho_ice=matpar->GetRhoIce();
     2643        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    26442644
    26452645        if(!IsIceInElement() || !IsOnSurface()) return 0.;

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

    r18911 r18946  
    758758
    759759                        /*Compute water pressure*/
    760                         IssmDouble rho_ice   = matpar->GetRhoIce();
    761                         IssmDouble rho_water = matpar->GetRhoWater();
    762                         IssmDouble gravity   = matpar->GetG();
     760                        IssmDouble rho_ice   = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     761                        IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     762                        IssmDouble gravity   = matpar->GetMaterialParameter(ConstantsGEnum);
    763763                        water_pressure=gravity*rho_water*base;
    764764
    … …  
    13521352        if(!IsIceInElement() || IsFloating())return 0;
    13531353
    1354         rho_ice=matpar->GetRhoIce();
    1355         rho_water=matpar->GetRhoWater();
     1354        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     1355        rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    13561356        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    13571357
    … …  
    17701770       
    17711771        /*Retrieve material parameters: */
    1772         rho_ice=matpar->GetRhoIce();
     1772        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    17731773
    17741774        /*Retrieve values of the levelset defining the masscon: */
    … …  
    18131813
    18141814        /*Get material parameters :*/
    1815         rho_ice=matpar->GetRhoIce();
     1815        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    18161816
    18171817        /*First off, check that this segment belongs to this element: */
    … …  
    18761876
    18771877        /*Get material parameters :*/
    1878         rho_ice=matpar->GetRhoIce();
     1878        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    18791879
    18801880        /*First off, check that this segment belongs to this element: */
    … …  
    22522252
    22532253  /*Get material parameters :*/
    2254   rho_ice=matpar->GetRhoIce();
    2255   rho_water=matpar->GetRhoFreshwater();
     2254  rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     2255  rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    22562256
    22572257  /*Get other pdd parameters*/
    2258   desfac=matpar->GetDesFac();
    2259   s0p=matpar->GetS0p();
    2260   s0t=matpar->GetS0t();
    2261   rlaps=matpar->GetRlaps();
    2262   rlapslgm=matpar->GetRlapslgm();
     2258  desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
     2259  s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
     2260  s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
     2261  rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
     2262  rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
    22632263     
    22642264   /*measure the surface mass balance*/
    … …  
    22852285
    22862286        /*material parameters: */
    2287         rho_water=matpar->GetRhoWater();
    2288         rho_ice=matpar->GetRhoIce();
     2287        rho_water=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     2288        rho_ice=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    22892289        density=rho_ice/rho_water;
    22902290        GetInputListOnVertices(&h[0],ThicknessEnum);
    … …  
    27082708
    27092709        /*Get material parameters :*/
    2710         rho_ice=matpar->GetRhoIce();
     2710        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    27112711
    27122712   if(!IsIceInElement())return 0;
    … …  
    31013101
    31023102        /*recover material parameters: */
    3103         lithosphere_shear_modulus=matpar->GetLithosphereShearModulus();
    3104         lithosphere_density=matpar->GetLithosphereDensity();
    3105         mantle_shear_modulus=matpar->GetMantleShearModulus();
    3106         mantle_density=matpar->GetMantleDensity();
    3107         rho_ice=matpar->GetRhoIce();
     3103        lithosphere_shear_modulus=matpar->GetMaterialParameter(MaterialsLithosphereShearModulusEnum);
     3104        lithosphere_density=matpar->GetMaterialParameter(MaterialsLithosphereDensityEnum);
     3105        mantle_shear_modulus=matpar->GetMaterialParameter(MaterialsMantleShearModulusEnum);
     3106        mantle_density=matpar->GetMaterialParameter(MaterialsMantleDensityEnum);
     3107        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    31083108
    31093109        /*pull thickness averages: */

  • issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp

    r18926 r18946  
    576576
    577577        /*Compute pressure melting point*/
    578         t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
     578        t_pmp=matpar->GetMaterialParameter(MaterialsMeltingpointEnum)-matpar->GetMaterialParameter(MaterialsBetaEnum)*pressure;
    579579
    580580        /*Add penalty load*/
    … …  
    651651
    652652        /*Compute pressure melting point*/
    653         t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
     653        t_pmp=matpar->GetMaterialParameter(MaterialsMeltingpointEnum)-matpar->GetMaterialParameter(MaterialsBetaEnum)*pressure;
    654654
    655655        /*Add penalty load
    … …  
    686686
    687687        /*Compute pressure melting point*/
    688         t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
     688        t_pmp=matpar->GetMaterialParameter(MaterialsMeltingpointEnum)-matpar->GetMaterialParameter(MaterialsBetaEnum)*pressure;
    689689
    690690        pe->values[0]=kmax*pow(10.,penalty_factor)*t_pmp;

  • issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp

    r18926 r18946  
    505505
    506506        /*Get some inputs: */
    507         rho_ice=matpar->GetRhoIce();
    508         rho_water=matpar->GetRhoWater();
    509         gravity=matpar->GetG();
     507        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     508        rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     509        gravity=matpar->GetMaterialParameter(ConstantsGEnum);
    510510        tria1->GetInputValue(&h[0],nodes[0],ThicknessEnum);
    511511        tria2->GetInputValue(&h[1],nodes[1],ThicknessEnum);

  • issm/trunk-jpl/src/c/classes/Materials/Material.h

    r18237 r18946  
    2525                virtual Material*  copy2(Element* element)=0;
    2626                virtual void       Configure(Elements* elements)=0;
     27                virtual IssmDouble GetA()=0;
     28                virtual IssmDouble GetAbar()=0;
     29                virtual IssmDouble GetB()=0;
     30                virtual IssmDouble GetBbar()=0;
     31                virtual IssmDouble GetD()=0;
     32                virtual IssmDouble GetDbar()=0;
     33                virtual IssmDouble GetN()=0;
    2734                virtual void       GetViscosity(IssmDouble* pviscosity,IssmDouble epseff)=0;
    28                 virtual void       GetViscosity_B(IssmDouble* pviscosity,IssmDouble epseff)=0;
    29                 virtual void       GetViscosity_D(IssmDouble* pviscosity,IssmDouble epseff)=0;
    3035                virtual void       GetViscosityBar(IssmDouble* pviscosity,IssmDouble epseff)=0;
    3136                virtual void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;
    3237                virtual void       GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;
    3338                virtual void       GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0;
     39                virtual void       GetViscosity_B(IssmDouble* pviscosity,IssmDouble epseff)=0;
     40                virtual void       GetViscosity_D(IssmDouble* pviscosity,IssmDouble epseff)=0;
    3441                virtual void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0;
    35                 virtual IssmDouble GetA()=0;
    36                 virtual IssmDouble GetAbar()=0;
    37                 virtual IssmDouble GetB()=0;
    38                 virtual IssmDouble GetBbar()=0;
    39                 virtual IssmDouble GetN()=0;
    40                 virtual IssmDouble GetD()=0;
    41                 virtual IssmDouble GetDbar()=0;
    4242                virtual bool       IsDamage()=0;
    4343                virtual void       ResetHooks()=0;

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r18237 r18946  
    6464
    6565/*Object virtual functions definitions:*/
    66 void Matice::Echo(void){/*{{{*/
    67 
    68         _printf_("Matice:\n");
    69         _printf_("   mid: " << mid << "\n");
    70         _printf_("   element:\n");
    71         helement->Echo();
    72 }
    73 /*}}}*/
    74 void Matice::DeepEcho(void){/*{{{*/
     66Object*   Matice::copy() {/*{{{*/
     67
     68        /*Output*/
     69        Matice* matice=NULL;
     70
     71        /*Initialize output*/
     72        matice=new Matice();
     73
     74        /*copy fields: */
     75        matice->mid=this->mid;
     76        matice->helement=(Hook*)this->helement->copy();
     77        matice->element =(Element*)this->helement->delivers();
     78        matice->isdamaged = this->isdamaged;
     79
     80        return matice;
     81}
     82/*}}}*/
     83Material* Matice::copy2(Element* element_in) {/*{{{*/
     84
     85        /*Output*/
     86        Matice* matice=NULL;
     87
     88        /*Initialize output*/
     89        matice=new Matice();
     90
     91        /*copy fields: */
     92        matice->mid=this->mid;
     93        matice->helement=(Hook*)this->helement->copy();
     94        matice->element =element_in;
     95        matice->isdamaged = this->isdamaged;
     96
     97        return matice;
     98}
     99/*}}}*/
     100void      Matice::DeepEcho(void){/*{{{*/
    75101
    76102        _printf_("Matice:\n");
    … …  
    80106}               
    81107/*}}}*/
    82 int    Matice::Id(void){ return mid; }/*{{{*/
    83 /*}}}*/
    84 int Matice::ObjectEnum(void){/*{{{*/
     108void      Matice::Echo(void){/*{{{*/
     109
     110        _printf_("Matice:\n");
     111        _printf_("   mid: " << mid << "\n");
     112        _printf_("   element:\n");
     113        helement->Echo();
     114}
     115/*}}}*/
     116int       Matice::Id(void){ return mid; }/*{{{*/
     117/*}}}*/
     118int       Matice::ObjectEnum(void){/*{{{*/
    85119
    86120        return MaticeEnum;
    87121
    88 }
    89 /*}}}*/
    90 Object* Matice::copy() {/*{{{*/
    91 
    92         /*Output*/
    93         Matice* matice=NULL;
    94 
    95         /*Initialize output*/
    96         matice=new Matice();
    97 
    98         /*copy fields: */
    99         matice->mid=this->mid;
    100         matice->helement=(Hook*)this->helement->copy();
    101         matice->element =(Element*)this->helement->delivers();
    102         matice->isdamaged = this->isdamaged;
    103 
    104         return matice;
    105 }
    106 /*}}}*/
    107 Material* Matice::copy2(Element* element_in) {/*{{{*/
    108 
    109         /*Output*/
    110         Matice* matice=NULL;
    111 
    112         /*Initialize output*/
    113         matice=new Matice();
    114 
    115         /*copy fields: */
    116         matice->mid=this->mid;
    117         matice->helement=(Hook*)this->helement->copy();
    118         matice->element =element_in;
    119         matice->isdamaged = this->isdamaged;
    120 
    121         return matice;
    122122}
    123123/*}}}*/

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r18717 r18946  
    304304
    305305/*Matpar management: */
    306 void  Matpar::Configure(Elements* elementsin){/*{{{*/
     306void       Matpar::Configure(Elements* elementsin){/*{{{*/
    307307
    308308        /*nothing done yet!*/
    309309
     310}
     311/*}}}*/
     312void       Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
     313
     314        /*Ouput*/
     315        IssmDouble temperature,waterfraction;
     316
     317        if(enthalpy<PureIceEnthalpy(pressure)){
     318                temperature=referencetemperature+enthalpy/heatcapacity;
     319                waterfraction=0.;
     320        }
     321        else{
     322                temperature=TMeltingPoint(pressure);
     323                waterfraction=(enthalpy-PureIceEnthalpy(pressure))/latentheat;
     324        }
     325
     326        /*Assign output pointers:*/
     327        *pwaterfraction=waterfraction;
     328        *ptemperature=temperature;
     329}
     330/*}}}*/
     331IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
     332        if (enthalpy<PureIceEnthalpy(pressure))
     333                return thermalconductivity/heatcapacity;
     334        else
     335                return temperateiceconductivity/heatcapacity;
     336}
     337/*}}}*/
     338IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/
     339
     340        int         iv;
     341        IssmDouble  lambda;                 // fraction of cold ice
     342        IssmDouble  kappa,kappa_c,kappa_t;  //enthalpy conductivities
     343        IssmDouble  Hc,Ht;
     344        IssmDouble* PIE   = xNew<IssmDouble>(numvertices);
     345        IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
     346
     347        for(iv=0; iv<numvertices; iv++){
     348                PIE[iv]=PureIceEnthalpy(pressure[iv]);
     349                dHpmp[iv]=enthalpy[iv]-PIE[iv];
     350        }
     351
     352        bool allequalsign=true;
     353        if(dHpmp[0]<0)
     354                for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0));
     355        else
     356                for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0));
     357
     358        if(allequalsign){
     359                kappa=GetEnthalpyDiffusionParameter(enthalpy[0], pressure[0]);
     360        }
     361        else {
     362                /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
     363                 cf Patankar 1980, pp44 */
     364                kappa_c=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.);
     365                kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);
     366                Hc=0.; Ht=0.;
     367                for(iv=0; iv<numvertices;iv++){
     368                        if(enthalpy[iv]<PIE[iv])
     369                         Hc+=(PIE[iv]-enthalpy[iv]);
     370                        else
     371                         Ht+=(enthalpy[iv]-PIE[iv]);
     372                }
     373                _assert_((Hc+Ht)>0.);
     374                lambda = Hc/(Hc+Ht);
     375                kappa  = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
     376        }
     377
     378        /*Clean up and return*/
     379        xDelete<IssmDouble>(PIE);
     380        xDelete<IssmDouble>(dHpmp);
     381        return kappa;
    310382}
    311383/*}}}*/
    … …  
    314386        switch(enum_in){
    315387                case MaterialsRhoIceEnum:                    return this->rho_ice;
    316                 case MaterialsRhoSeawaterEnum:                  return this->rho_water;
     388                case MaterialsRhoSeawaterEnum:               return this->rho_water;
    317389                case MaterialsRhoFreshwaterEnum:             return this->rho_freshwater;
    318390                case MaterialsMuWaterEnum:                   return this->mu_water;
    … …  
    344416                case MaterialsCompressionCoefEnum:                              return this->compression_coef;
    345417                case MaterialsTractionCoefEnum:                                 return this->traction_coef;
     418                case SurfaceforcingsDesfacEnum:              return this->desfac;
     419                case SurfaceforcingsS0pEnum:                 return this->s0p;
     420                case SurfaceforcingsS0tEnum:                 return this->s0t;
     421                case SurfaceforcingsRlapsEnum:               return this->rlaps;
     422                case SurfaceforcingsRlapslgmEnum:            return this->rlapslgm;
     423                case MaterialsLithosphereShearModulusEnum:   return this->lithosphere_shear_modulus;
     424                case MaterialsLithosphereDensityEnum:        return this->lithosphere_density;
     425                case MaterialsMantleDensityEnum:             return this->mantle_density;
     426                case MaterialsMantleShearModulusEnum:        return this->mantle_shear_modulus;
    346427                case ConstantsOmegaEnum:                                                        return this->omega;
    347428                case MaterialsTimeRelaxationStressEnum:         return this->time_relaxation_stress;
    … …  
    352433}
    353434/*}}}*/
    354 IssmDouble Matpar::GetBeta(){/*{{{*/
    355         return beta;
    356 }
    357 /*}}}*/
    358 IssmDouble Matpar::GetG(){/*{{{*/
    359         return g;
    360 }
    361 /*}}}*/
    362 IssmDouble Matpar::GetMeltingPoint(){/*{{{*/
    363         return meltingpoint;
    364 }
    365 /*}}}*/
    366 IssmDouble Matpar::GetRhoIce(){/*{{{*/
    367 
    368         return rho_ice;
    369 }
    370 /*}}}*/
    371 IssmDouble Matpar::GetRhoWater(){/*{{{*/
    372         return rho_water;
    373 }
    374 /*}}}*/
    375 IssmDouble Matpar::GetRhoFreshwater(){/*{{{*/
    376         return rho_freshwater;
    377 }
    378 /*}}}*/
    379 IssmDouble Matpar::GetDesFac(){/*{{{*/
    380         return desfac;
    381 }
    382 /*}}}*/
    383 IssmDouble Matpar::GetS0p(){/*{{{*/
    384         return s0p;
    385 }
    386 /*}}}*/
    387 IssmDouble Matpar::GetS0t(){/*{{{*/
    388         return s0t;
    389 }
    390 /*}}}*/
    391 IssmDouble Matpar::GetRlaps(){/*{{{*/
    392         return rlaps;
    393 }
    394 /*}}}*/
    395 IssmDouble Matpar::GetRlapslgm(){/*{{{*/
    396         return rlapslgm;
     435IssmDouble Matpar::PureIceEnthalpy(IssmDouble pressure){/*{{{*/
     436        return heatcapacity*(TMeltingPoint(pressure)-referencetemperature);
     437}
     438/*}}}*/
     439void       Matpar::ResetHooks(){/*{{{*/
     440
     441        //Nothing to be done
     442        return;
     443}
     444/*}}}*/
     445void       Matpar::ThermalToEnthalpy(IssmDouble * penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
     446
     447        /*Ouput*/
     448        IssmDouble enthalpy;
     449
     450        if(temperature<TMeltingPoint(pressure)){
     451                enthalpy=heatcapacity*(temperature-referencetemperature);
     452        }
     453        else{
     454                enthalpy=PureIceEnthalpy(pressure)+latentheat*waterfraction;
     455        }
     456
     457        /*Assign output pointers:*/
     458        *penthalpy=enthalpy;
    397459}
    398460/*}}}*/
    … …  
    401463}
    402464/*}}}*/
    403 IssmDouble Matpar::PureIceEnthalpy(IssmDouble pressure){/*{{{*/
    404         return heatcapacity*(TMeltingPoint(pressure)-referencetemperature);
    405 }
    406 /*}}}*/
    407 IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
    408         if (enthalpy<PureIceEnthalpy(pressure))
    409                 return thermalconductivity/heatcapacity;
    410         else
    411                 return temperateiceconductivity/heatcapacity;
    412 }
    413 /*}}}*/
    414 IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/
    415 
    416         int         iv;
    417         IssmDouble  lambda;                 // fraction of cold ice
    418         IssmDouble  kappa,kappa_c,kappa_t;  //enthalpy conductivities
    419         IssmDouble  Hc,Ht;
    420         IssmDouble* PIE   = xNew<IssmDouble>(numvertices);
    421         IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
    422 
    423         for(iv=0; iv<numvertices; iv++){
    424                 PIE[iv]=PureIceEnthalpy(pressure[iv]);
    425                 dHpmp[iv]=enthalpy[iv]-PIE[iv];
    426         }
    427 
    428         bool allequalsign=true;
    429         if(dHpmp[0]<0)
    430                 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0));
    431         else
    432                 for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0));
    433 
    434         if(allequalsign){
    435                 kappa=GetEnthalpyDiffusionParameter(enthalpy[0], pressure[0]);
    436         }
    437         else {
    438                 /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
    439                  cf Patankar 1980, pp44 */
    440                 kappa_c=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.);
    441                 kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);
    442                 Hc=0.; Ht=0.;
    443                 for(iv=0; iv<numvertices;iv++){
    444                         if(enthalpy[iv]<PIE[iv])
    445                          Hc+=(PIE[iv]-enthalpy[iv]);
    446                         else
    447                          Ht+=(enthalpy[iv]-PIE[iv]);
    448                 }
    449                 _assert_((Hc+Ht)>0.);
    450                 lambda = Hc/(Hc+Ht);
    451                 kappa  = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
    452         }
    453 
    454         /*Clean up and return*/
    455         xDelete<IssmDouble>(PIE);
    456         xDelete<IssmDouble>(dHpmp);
    457         return kappa;
    458 }
    459 /*}}}*/
    460 void Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
    461 
    462         /*Ouput*/
    463         IssmDouble temperature,waterfraction;
    464 
    465         if(enthalpy<PureIceEnthalpy(pressure)){
    466                 temperature=referencetemperature+enthalpy/heatcapacity;
    467                 waterfraction=0.;
    468         }
    469         else{
    470                 temperature=TMeltingPoint(pressure);
    471                 waterfraction=(enthalpy-PureIceEnthalpy(pressure))/latentheat;
    472         }
    473 
    474         /*Assign output pointers:*/
    475         *pwaterfraction=waterfraction;
    476         *ptemperature=temperature;
    477 }
    478 /*}}}*/
    479 void Matpar::ThermalToEnthalpy(IssmDouble * penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
    480 
    481         /*Ouput*/
    482         IssmDouble enthalpy;
    483 
    484         if(temperature<TMeltingPoint(pressure)){
    485                 enthalpy=heatcapacity*(temperature-referencetemperature);
    486         }
    487         else{
    488                 enthalpy=PureIceEnthalpy(pressure)+latentheat*waterfraction;
    489         }
    490 
    491         /*Assign output pointers:*/
    492         *penthalpy=enthalpy;
    493 }
    494 /*}}}*/
    495 IssmDouble Matpar::GetLithosphereShearModulus(){                 /*{{{*/
    496         return lithosphere_shear_modulus;                       
    497 }               
    498 /*}}}*/
    499 IssmDouble Matpar::GetLithosphereDensity(){              /*{{{*/
    500         return lithosphere_density;                     
    501 }               
    502 /*}}}*/
    503 IssmDouble Matpar::GetMantleDensity(){           /*{{{*/
    504         return mantle_density;                   
    505 }               
    506 /*}}}*/
    507 IssmDouble Matpar::GetMantleShearModulus(){              /*{{{*/
    508         return mantle_shear_modulus;                     
    509 }               
    510 /*}}}*/
    511 void  Matpar::ResetHooks(){/*{{{*/
    512 
    513         //Nothing to be done
    514         return;
    515 }
    516 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r18717 r18946  
    110110                /*}}}*/
    111111                /*Numerics: {{{*/
    112                 IssmDouble GetG();
    113                 IssmDouble GetRhoIce();
    114                 IssmDouble GetRhoWater();
    115                 IssmDouble GetRhoFreshwater();
    116                 IssmDouble GetBeta();
    117                 IssmDouble GetMeltingPoint();
    118                 IssmDouble TMeltingPoint(IssmDouble pressure);
    119                 IssmDouble PureIceEnthalpy(IssmDouble pressure);
     112                void       EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
    120113                IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
    121114                IssmDouble GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
    122                 IssmDouble GetLithosphereShearModulus();
    123                 IssmDouble GetLithosphereDensity();
    124                 IssmDouble GetMantleShearModulus();
    125                 IssmDouble GetMantleDensity();
    126                 void       EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
     115                IssmDouble GetMaterialParameter(int in_enum);
     116                IssmDouble PureIceEnthalpy(IssmDouble pressure);
    127117                void       ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
    128                 IssmDouble GetDesFac();
    129                 IssmDouble GetS0p();
    130                 IssmDouble GetS0t();
    131                 IssmDouble GetRlaps();
    132                 IssmDouble GetRlapslgm();
    133                 IssmDouble GetMaterialParameter(int in_enum);
     118                IssmDouble TMeltingPoint(IssmDouble pressure);
    134119                /*}}}*/
    135120

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp

    r18532 r18946  
    8383        CreateFaces(iomodel);
    8484
    85         /*!All elements have been partitioned above, only create elements for this CPU: */
    86         for(int i=0;i<iomodel->numberoffaces;i++){
     85        if(iomodel->domaintype==Domain2DhorizontalEnum){
     86                /*!All elements have been partitioned above, only create elements for this CPU: */
     87                for(int i=0;i<iomodel->numberoffaces;i++){
    8788
    88                 /*Get left and right elements*/
    89                 e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
    90                 e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
     89                        /*Get left and right elements*/
     90                        e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
     91                        e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
    9192
    92                 /* 1) If the element e1 is in the current partition
    93                 * 2) and if the face of the element is shared by another element (internal face)
    94                 * 3) and if this element is not in the same partition:
    95                 * we must clone the nodes on this partition so that the loads (Numericalflux)
    96                 * will have access to their properties (dofs,...)*/
    97                 if(my_elements[e1] && e2!=-2 && !my_elements[e2]){
     93                        /* 1) If the element e1 is in the current partition
     94                        * 2) and if the face of the element is shared by another element (internal face)
     95                        * 3) and if this element is not in the same partition:
     96                        * we must clone the nodes on this partition so that the loads (Numericalflux)
     97                        * will have access to their properties (dofs,...)*/
     98                        if(my_elements[e1] && e2!=-2 && !my_elements[e2]){
    9899
    99                         /*1: Get vertices ids*/
    100                         i1=iomodel->faces[4*i+0];
    101                         i2=iomodel->faces[4*i+1];
     100                                /*1: Get vertices ids*/
     101                                i1=iomodel->faces[4*i+0];
     102                                i2=iomodel->faces[4*i+1];
    102103
    103                         /*2: Get the column where these ids are located in the index*/
    104                         pos=UNDEF;
    105                         for(int j=0;j<3;j++){
    106                                 if(iomodel->elements[3*e2+j]==i1) pos=j;
    107                         }
     104                                /*2: Get the column where these ids are located in the index*/
     105                                pos=UNDEF;
     106                                for(int j=0;j<3;j++){
     107                                        if(iomodel->elements[3*e2+j]==i1) pos=j;
     108                                }
    108109
    109                         /*3: We have the id of the elements and the position of the vertices in the index
    110                          * we can now create the corresponding nodes:*/
    111                         if(pos==0){
    112                                 my_nodes[e2*3+0]=true;
    113                                 my_nodes[e2*3+2]=true;
    114                         }
    115                         else if(pos==1){
    116                                 my_nodes[e2*3+1]=true;
    117                                 my_nodes[e2*3+0]=true;
    118                         }
    119                         else if(pos==2){
    120                                 my_nodes[e2*3+2]=true;
    121                                 my_nodes[e2*3+1]=true;
    122                         }
    123                         else{
    124                                 _error_("Problem in faces creation");
     110                                /*3: We have the id of the elements and the position of the vertices in the index
     111                                 * we can now create the corresponding nodes:*/
     112                                if(pos==0){
     113                                        my_nodes[e2*3+0]=true;
     114                                        my_nodes[e2*3+2]=true;
     115                                }
     116                                else if(pos==1){
     117                                        my_nodes[e2*3+1]=true;
     118                                        my_nodes[e2*3+0]=true;
     119                                }
     120                                else if(pos==2){
     121                                        my_nodes[e2*3+2]=true;
     122                                        my_nodes[e2*3+1]=true;
     123                                }
     124                                else{
     125                                        _error_("Problem in faces creation");
     126                                }
    125127                        }
    126128                }

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

    r18584 r18946  
    8585
    8686                /*Get material parameters :*/
    87                 rho_ice=element->matpar->GetRhoIce();
    88                 rho_water=element->matpar->GetRhoFreshwater();
     87                rho_ice=element->matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     88                rho_water=element->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    8989
    9090                // loop over all vertices
Note: See TracChangeset for help on using the changeset viewer.