Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18945) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18946) @@ -1226,6 +1226,6 @@ GetInputListOnVertices(&r[0],BedEnum); GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum); - rho_water = matpar->GetRhoWater(); - rho_ice = matpar->GetRhoIce(); + rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum); density = rho_ice/rho_water; @@ -2054,5 +2054,5 @@ surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); z=this->GetZcoord(xyz_list,gauss); - tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]); + tau_perp = matpar->GetMaterialParameter(MaterialsRhoIceEnum) * matpar->GetMaterialParameter(ConstantsGEnum) * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]); /* Get eps_b*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18945) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18946) @@ -252,6 +252,6 @@ /*recovre material parameters: */ - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + gravity=matpar->GetMaterialParameter(ConstantsGEnum); /* Get node coordinates and dof list: */ @@ -697,7 +697,7 @@ /*Compute water pressure*/ - IssmDouble rho_ice = matpar->GetRhoIce(); - IssmDouble rho_water = matpar->GetRhoWater(); - IssmDouble gravity = matpar->GetG(); + IssmDouble rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble gravity = matpar->GetMaterialParameter(ConstantsGEnum); water_pressure=gravity*rho_water*base; @@ -1177,6 +1177,6 @@ if(!IsIceInElement() || IsFloating() || !IsOnBase())return 0; - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -2049,13 +2049,13 @@ /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoFreshwater(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum); /*Get other pdd parameters*/ - desfac=matpar->GetDesFac(); - s0p=matpar->GetS0p(); - s0t=matpar->GetS0t(); - rlaps=matpar->GetRlaps(); - rlapslgm=matpar->GetRlapslgm(); + desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); + s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); + s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); + rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); + rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); /*measure the surface mass balance*/ @@ -2083,6 +2083,6 @@ /*material parameters: */ - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); density=rho_ice/rho_water; GetInputListOnVertices(&h[0],ThicknessEnum); @@ -2641,5 +2641,5 @@ /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); if(!IsIceInElement() || !IsOnSurface()) return 0.; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18945) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18946) @@ -758,7 +758,7 @@ /*Compute water pressure*/ - IssmDouble rho_ice = matpar->GetRhoIce(); - IssmDouble rho_water = matpar->GetRhoWater(); - IssmDouble gravity = matpar->GetG(); + IssmDouble rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble gravity = matpar->GetMaterialParameter(ConstantsGEnum); water_pressure=gravity*rho_water*base; @@ -1352,6 +1352,6 @@ if(!IsIceInElement() || IsFloating())return 0; - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -1770,5 +1770,5 @@ /*Retrieve material parameters: */ - rho_ice=matpar->GetRhoIce(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); /*Retrieve values of the levelset defining the masscon: */ @@ -1813,5 +1813,5 @@ /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); /*First off, check that this segment belongs to this element: */ @@ -1876,5 +1876,5 @@ /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); /*First off, check that this segment belongs to this element: */ @@ -2252,13 +2252,13 @@ /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoFreshwater(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); /*Get other pdd parameters*/ - desfac=matpar->GetDesFac(); - s0p=matpar->GetS0p(); - s0t=matpar->GetS0t(); - rlaps=matpar->GetRlaps(); - rlapslgm=matpar->GetRlapslgm(); + desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); + s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); + s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); + rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); + rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); /*measure the surface mass balance*/ @@ -2285,6 +2285,6 @@ /*material parameters: */ - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); density=rho_ice/rho_water; GetInputListOnVertices(&h[0],ThicknessEnum); @@ -2708,5 +2708,5 @@ /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); if(!IsIceInElement())return 0; @@ -3101,9 +3101,9 @@ /*recover material parameters: */ - lithosphere_shear_modulus=matpar->GetLithosphereShearModulus(); - lithosphere_density=matpar->GetLithosphereDensity(); - mantle_shear_modulus=matpar->GetMantleShearModulus(); - mantle_density=matpar->GetMantleDensity(); - rho_ice=matpar->GetRhoIce(); + lithosphere_shear_modulus=matpar->GetMaterialParameter(MaterialsLithosphereShearModulusEnum); + lithosphere_density=matpar->GetMaterialParameter(MaterialsLithosphereDensityEnum); + mantle_shear_modulus=matpar->GetMaterialParameter(MaterialsMantleShearModulusEnum); + mantle_density=matpar->GetMaterialParameter(MaterialsMantleDensityEnum); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); /*pull thickness averages: */ Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 18945) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 18946) @@ -576,5 +576,5 @@ /*Compute pressure melting point*/ - t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; + t_pmp=matpar->GetMaterialParameter(MaterialsMeltingpointEnum)-matpar->GetMaterialParameter(MaterialsBetaEnum)*pressure; /*Add penalty load*/ @@ -651,5 +651,5 @@ /*Compute pressure melting point*/ - t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; + t_pmp=matpar->GetMaterialParameter(MaterialsMeltingpointEnum)-matpar->GetMaterialParameter(MaterialsBetaEnum)*pressure; /*Add penalty load @@ -686,5 +686,5 @@ /*Compute pressure melting point*/ - t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; + t_pmp=matpar->GetMaterialParameter(MaterialsMeltingpointEnum)-matpar->GetMaterialParameter(MaterialsBetaEnum)*pressure; pe->values[0]=kmax*pow(10.,penalty_factor)*t_pmp; Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 18945) +++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 18946) @@ -505,7 +505,7 @@ /*Get some inputs: */ - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - gravity=matpar->GetG(); + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + gravity=matpar->GetMaterialParameter(ConstantsGEnum); tria1->GetInputValue(&h[0],nodes[0],ThicknessEnum); tria2->GetInputValue(&h[1],nodes[1],ThicknessEnum); Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 18945) +++ /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 18946) @@ -25,19 +25,19 @@ virtual Material* copy2(Element* element)=0; virtual void Configure(Elements* elements)=0; + virtual IssmDouble GetA()=0; + virtual IssmDouble GetAbar()=0; + virtual IssmDouble GetB()=0; + virtual IssmDouble GetBbar()=0; + virtual IssmDouble GetD()=0; + virtual IssmDouble GetDbar()=0; + virtual IssmDouble GetN()=0; virtual void GetViscosity(IssmDouble* pviscosity,IssmDouble epseff)=0; - virtual void GetViscosity_B(IssmDouble* pviscosity,IssmDouble epseff)=0; - virtual void GetViscosity_D(IssmDouble* pviscosity,IssmDouble epseff)=0; virtual void GetViscosityBar(IssmDouble* pviscosity,IssmDouble epseff)=0; virtual void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; virtual void GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; virtual void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; + virtual void GetViscosity_B(IssmDouble* pviscosity,IssmDouble epseff)=0; + virtual void GetViscosity_D(IssmDouble* pviscosity,IssmDouble epseff)=0; virtual void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; - virtual IssmDouble GetA()=0; - virtual IssmDouble GetAbar()=0; - virtual IssmDouble GetB()=0; - virtual IssmDouble GetBbar()=0; - virtual IssmDouble GetN()=0; - virtual IssmDouble GetD()=0; - virtual IssmDouble GetDbar()=0; virtual bool IsDamage()=0; virtual void ResetHooks()=0; Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18945) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18946) @@ -64,13 +64,39 @@ /*Object virtual functions definitions:*/ -void Matice::Echo(void){/*{{{*/ - - _printf_("Matice:\n"); - _printf_(" mid: " << mid << "\n"); - _printf_(" element:\n"); - helement->Echo(); -} -/*}}}*/ -void Matice::DeepEcho(void){/*{{{*/ +Object* Matice::copy() {/*{{{*/ + + /*Output*/ + Matice* matice=NULL; + + /*Initialize output*/ + matice=new Matice(); + + /*copy fields: */ + matice->mid=this->mid; + matice->helement=(Hook*)this->helement->copy(); + matice->element =(Element*)this->helement->delivers(); + matice->isdamaged = this->isdamaged; + + return matice; +} +/*}}}*/ +Material* Matice::copy2(Element* element_in) {/*{{{*/ + + /*Output*/ + Matice* matice=NULL; + + /*Initialize output*/ + matice=new Matice(); + + /*copy fields: */ + matice->mid=this->mid; + matice->helement=(Hook*)this->helement->copy(); + matice->element =element_in; + matice->isdamaged = this->isdamaged; + + return matice; +} +/*}}}*/ +void Matice::DeepEcho(void){/*{{{*/ _printf_("Matice:\n"); @@ -80,44 +106,18 @@ } /*}}}*/ -int Matice::Id(void){ return mid; }/*{{{*/ -/*}}}*/ -int Matice::ObjectEnum(void){/*{{{*/ +void Matice::Echo(void){/*{{{*/ + + _printf_("Matice:\n"); + _printf_(" mid: " << mid << "\n"); + _printf_(" element:\n"); + helement->Echo(); +} +/*}}}*/ +int Matice::Id(void){ return mid; }/*{{{*/ +/*}}}*/ +int Matice::ObjectEnum(void){/*{{{*/ return MaticeEnum; -} -/*}}}*/ -Object* Matice::copy() {/*{{{*/ - - /*Output*/ - Matice* matice=NULL; - - /*Initialize output*/ - matice=new Matice(); - - /*copy fields: */ - matice->mid=this->mid; - matice->helement=(Hook*)this->helement->copy(); - matice->element =(Element*)this->helement->delivers(); - matice->isdamaged = this->isdamaged; - - return matice; -} -/*}}}*/ -Material* Matice::copy2(Element* element_in) {/*{{{*/ - - /*Output*/ - Matice* matice=NULL; - - /*Initialize output*/ - matice=new Matice(); - - /*copy fields: */ - matice->mid=this->mid; - matice->helement=(Hook*)this->helement->copy(); - matice->element =element_in; - matice->isdamaged = this->isdamaged; - - return matice; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 18945) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 18946) @@ -304,8 +304,80 @@ /*Matpar management: */ -void Matpar::Configure(Elements* elementsin){/*{{{*/ +void Matpar::Configure(Elements* elementsin){/*{{{*/ /*nothing done yet!*/ +} +/*}}}*/ +void Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/ + + /*Ouput*/ + IssmDouble temperature,waterfraction; + + if(enthalpy(numvertices); + IssmDouble* dHpmp = xNew(numvertices); + + for(iv=0; iv=0)); + + if(allequalsign){ + kappa=GetEnthalpyDiffusionParameter(enthalpy[0], pressure[0]); + } + else { + /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice, + cf Patankar 1980, pp44 */ + kappa_c=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.); + kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.); + Hc=0.; Ht=0.; + for(iv=0; iv0.); + lambda = Hc/(Hc+Ht); + kappa = 1./(lambda/kappa_c + (1.-lambda)/kappa_t); + } + + /*Clean up and return*/ + xDelete(PIE); + xDelete(dHpmp); + return kappa; } /*}}}*/ @@ -314,5 +386,5 @@ switch(enum_in){ case MaterialsRhoIceEnum: return this->rho_ice; - case MaterialsRhoSeawaterEnum: return this->rho_water; + case MaterialsRhoSeawaterEnum: return this->rho_water; case MaterialsRhoFreshwaterEnum: return this->rho_freshwater; case MaterialsMuWaterEnum: return this->mu_water; @@ -344,4 +416,13 @@ case MaterialsCompressionCoefEnum: return this->compression_coef; case MaterialsTractionCoefEnum: return this->traction_coef; + case SurfaceforcingsDesfacEnum: return this->desfac; + case SurfaceforcingsS0pEnum: return this->s0p; + case SurfaceforcingsS0tEnum: return this->s0t; + case SurfaceforcingsRlapsEnum: return this->rlaps; + case SurfaceforcingsRlapslgmEnum: return this->rlapslgm; + case MaterialsLithosphereShearModulusEnum: return this->lithosphere_shear_modulus; + case MaterialsLithosphereDensityEnum: return this->lithosphere_density; + case MaterialsMantleDensityEnum: return this->mantle_density; + case MaterialsMantleShearModulusEnum: return this->mantle_shear_modulus; case ConstantsOmegaEnum: return this->omega; case MaterialsTimeRelaxationStressEnum: return this->time_relaxation_stress; @@ -352,47 +433,28 @@ } /*}}}*/ -IssmDouble Matpar::GetBeta(){/*{{{*/ - return beta; -} -/*}}}*/ -IssmDouble Matpar::GetG(){/*{{{*/ - return g; -} -/*}}}*/ -IssmDouble Matpar::GetMeltingPoint(){/*{{{*/ - return meltingpoint; -} -/*}}}*/ -IssmDouble Matpar::GetRhoIce(){/*{{{*/ - - return rho_ice; -} -/*}}}*/ -IssmDouble Matpar::GetRhoWater(){/*{{{*/ - return rho_water; -} -/*}}}*/ -IssmDouble Matpar::GetRhoFreshwater(){/*{{{*/ - return rho_freshwater; -} -/*}}}*/ -IssmDouble Matpar::GetDesFac(){/*{{{*/ - return desfac; -} -/*}}}*/ -IssmDouble Matpar::GetS0p(){/*{{{*/ - return s0p; -} -/*}}}*/ -IssmDouble Matpar::GetS0t(){/*{{{*/ - return s0t; -} -/*}}}*/ -IssmDouble Matpar::GetRlaps(){/*{{{*/ - return rlaps; -} -/*}}}*/ -IssmDouble Matpar::GetRlapslgm(){/*{{{*/ - return rlapslgm; +IssmDouble Matpar::PureIceEnthalpy(IssmDouble pressure){/*{{{*/ + return heatcapacity*(TMeltingPoint(pressure)-referencetemperature); +} +/*}}}*/ +void Matpar::ResetHooks(){/*{{{*/ + + //Nothing to be done + return; +} +/*}}}*/ +void Matpar::ThermalToEnthalpy(IssmDouble * penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/ + + /*Ouput*/ + IssmDouble enthalpy; + + if(temperature(numvertices); - IssmDouble* dHpmp = xNew(numvertices); - - for(iv=0; iv=0)); - - if(allequalsign){ - kappa=GetEnthalpyDiffusionParameter(enthalpy[0], pressure[0]); - } - else { - /* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice, - cf Patankar 1980, pp44 */ - kappa_c=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.); - kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.); - Hc=0.; Ht=0.; - for(iv=0; iv0.); - lambda = Hc/(Hc+Ht); - kappa = 1./(lambda/kappa_c + (1.-lambda)/kappa_t); - } - - /*Clean up and return*/ - xDelete(PIE); - xDelete(dHpmp); - return kappa; -} -/*}}}*/ -void Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/ - - /*Ouput*/ - IssmDouble temperature,waterfraction; - - if(enthalpynumberoffaces;i++){ + if(iomodel->domaintype==Domain2DhorizontalEnum){ + /*!All elements have been partitioned above, only create elements for this CPU: */ + for(int i=0;inumberoffaces;i++){ - /*Get left and right elements*/ - e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] - e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2] + /*Get left and right elements*/ + e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] + e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2] - /* 1) If the element e1 is in the current partition - * 2) and if the face of the element is shared by another element (internal face) - * 3) and if this element is not in the same partition: - * we must clone the nodes on this partition so that the loads (Numericalflux) - * will have access to their properties (dofs,...)*/ - if(my_elements[e1] && e2!=-2 && !my_elements[e2]){ + /* 1) If the element e1 is in the current partition + * 2) and if the face of the element is shared by another element (internal face) + * 3) and if this element is not in the same partition: + * we must clone the nodes on this partition so that the loads (Numericalflux) + * will have access to their properties (dofs,...)*/ + if(my_elements[e1] && e2!=-2 && !my_elements[e2]){ - /*1: Get vertices ids*/ - i1=iomodel->faces[4*i+0]; - i2=iomodel->faces[4*i+1]; + /*1: Get vertices ids*/ + i1=iomodel->faces[4*i+0]; + i2=iomodel->faces[4*i+1]; - /*2: Get the column where these ids are located in the index*/ - pos=UNDEF; - for(int j=0;j<3;j++){ - if(iomodel->elements[3*e2+j]==i1) pos=j; - } + /*2: Get the column where these ids are located in the index*/ + pos=UNDEF; + for(int j=0;j<3;j++){ + if(iomodel->elements[3*e2+j]==i1) pos=j; + } - /*3: We have the id of the elements and the position of the vertices in the index - * we can now create the corresponding nodes:*/ - if(pos==0){ - my_nodes[e2*3+0]=true; - my_nodes[e2*3+2]=true; - } - else if(pos==1){ - my_nodes[e2*3+1]=true; - my_nodes[e2*3+0]=true; - } - else if(pos==2){ - my_nodes[e2*3+2]=true; - my_nodes[e2*3+1]=true; - } - else{ - _error_("Problem in faces creation"); + /*3: We have the id of the elements and the position of the vertices in the index + * we can now create the corresponding nodes:*/ + if(pos==0){ + my_nodes[e2*3+0]=true; + my_nodes[e2*3+2]=true; + } + else if(pos==1){ + my_nodes[e2*3+1]=true; + my_nodes[e2*3+0]=true; + } + else if(pos==2){ + my_nodes[e2*3+2]=true; + my_nodes[e2*3+1]=true; + } + else{ + _error_("Problem in faces creation"); + } } } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 18945) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 18946) @@ -85,6 +85,6 @@ /*Get material parameters :*/ - rho_ice=element->matpar->GetRhoIce(); - rho_water=element->matpar->GetRhoFreshwater(); + rho_ice=element->matpar->GetMaterialParameter(MaterialsRhoIceEnum); + rho_water=element->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); // loop over all vertices