Index: ../trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp	(revision 18945)
+++ ../trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp	(revision 18946)
@@ -82,47 +82,49 @@
 	/*Get faces and elements*/
 	CreateFaces(iomodel);
 
-	/*!All elements have been partitioned above, only create elements for this CPU: */
-	for(int i=0;i<iomodel->numberoffaces;i++){
+	if(iomodel->domaintype==Domain2DhorizontalEnum){
+		/*!All elements have been partitioned above, only create elements for this CPU: */
+		for(int i=0;i<iomodel->numberoffaces;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;
+				/*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");
+				}
 			}
-			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: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 18945)
+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 18946)
@@ -84,8 +84,8 @@
 		element->GetInputListOnVertices(s,SurfaceEnum);
 
 		/*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
 		for(v=0;v<numvertices;v++){
Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 18945)
+++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 18946)
@@ -63,32 +63,8 @@
 /*}}}*/
 
 /*Object virtual functions definitions:*/
-void Matice::Echo(void){/*{{{*/
+Object*   Matice::copy() {/*{{{*/
 
-	_printf_("Matice:\n");
-	_printf_("   mid: " << mid << "\n");
-	_printf_("   element:\n");
-	helement->Echo();
-}
-/*}}}*/
-void Matice::DeepEcho(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;
 
@@ -121,7 +97,31 @@
 	return matice;
 }
 /*}}}*/
+void      Matice::DeepEcho(void){/*{{{*/
 
+	_printf_("Matice:\n");
+	_printf_("   mid: " << mid << "\n");
+	_printf_("   element:\n");
+	helement->Echo();
+}		
+/*}}}*/
+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;
+
+}
+/*}}}*/
+
 /*Matice management*/
 void  Matice::Configure(Elements* elementsin){/*{{{*/
 
Index: ../trunk-jpl/src/c/classes/Materials/Material.h
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Material.h	(revision 18945)
+++ ../trunk-jpl/src/c/classes/Materials/Material.h	(revision 18946)
@@ -24,21 +24,21 @@
 		/*Numerics*/
 		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: ../trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 18945)
+++ ../trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 18946)
@@ -303,107 +303,31 @@
 /*}}}*/
 
 /*Matpar management: */
-void  Matpar::Configure(Elements* elementsin){/*{{{*/
+void       Matpar::Configure(Elements* elementsin){/*{{{*/
 
 	/*nothing done yet!*/
 
 }
 /*}}}*/
-IssmDouble Matpar::GetMaterialParameter(int enum_in){/*{{{*/
+void       Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
 
-	switch(enum_in){
-		case MaterialsRhoIceEnum:                    return this->rho_ice;
-		case MaterialsRhoSeawaterEnum:                  return this->rho_water;
-		case MaterialsRhoFreshwaterEnum:             return this->rho_freshwater;
-		case MaterialsMuWaterEnum:                   return this->mu_water;
-		case MaterialsHeatcapacityEnum:              return this->heatcapacity;
-		case MaterialsThermalconductivityEnum:       return this->thermalconductivity;
-		case MaterialsTemperateiceconductivityEnum:  return this->temperateiceconductivity;
-		case MaterialsLatentheatEnum:                return this->latentheat;
-		case MaterialsBetaEnum:                      return this->beta;
-		case MaterialsMeltingpointEnum:              return this->meltingpoint;
-		case ConstantsReferencetemperatureEnum:      return this->referencetemperature;
-		case MaterialsMixedLayerCapacityEnum:        return this->mixed_layer_capacity;
-		case MaterialsThermalExchangeVelocityEnum:   return this->thermal_exchange_velocity;
-		case HydrologydcSedimentPorosityEnum:        return this->sediment_porosity;
-		case HydrologydcSedimentThicknessEnum:       return this->sediment_thickness;
-		case HydrologydcSedimentCompressibilityEnum: return this->sediment_compressibility;
-		case HydrologydcEplPorosityEnum:             return this->epl_porosity;
-		case HydrologydcEplCompressibilityEnum:      return this->epl_compressibility;
-		case HydrologydcEplConductivityEnum:         return this->epl_conductivity;
-		case HydrologydcEplInitialThicknessEnum:     return this->epl_init_thickness;
-		case HydrologydcEplColapseThicknessEnum:     return this->epl_colapse_thickness;
-		case HydrologydcEplMaxThicknessEnum:         return this->epl_max_thickness;
-		case HydrologydcWaterCompressibilityEnum:    return this->water_compressibility;
-		case ConstantsGEnum:                         return this->g;
-		case MaterialsPoissonEnum:                   return this->poisson;
-		case MaterialsYoungModulusEnum:              return this->young_modulus;
-		case MaterialsRidgingExponentEnum:				return this->ridging_exponent;
-		case MaterialsCohesionEnum:		            return this->cohesion;
-		case MaterialsInternalFrictionCoefEnum:		return this->internal_friction_coef;
-		case MaterialsCompressionCoefEnum:				return this->compression_coef;
-		case MaterialsTractionCoefEnum:					return this->traction_coef;
-		case ConstantsOmegaEnum:							return this->omega;
-		case MaterialsTimeRelaxationStressEnum:		return this->time_relaxation_stress;
-		case MaterialsTimeRelaxationDamageEnum:		return this->time_relaxation_damage;
-		default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");
+	/*Ouput*/
+	IssmDouble temperature,waterfraction;
+
+	if(enthalpy<PureIceEnthalpy(pressure)){
+		temperature=referencetemperature+enthalpy/heatcapacity;
+		waterfraction=0.;
 	}
+	else{
+		temperature=TMeltingPoint(pressure);
+		waterfraction=(enthalpy-PureIceEnthalpy(pressure))/latentheat;
+	}
 
+	/*Assign output pointers:*/
+	*pwaterfraction=waterfraction;
+	*ptemperature=temperature;
 }
 /*}}}*/
-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::TMeltingPoint(IssmDouble pressure){/*{{{*/
-	return meltingpoint-beta*pressure;
-}
-/*}}}*/
-IssmDouble Matpar::PureIceEnthalpy(IssmDouble pressure){/*{{{*/
-	return heatcapacity*(TMeltingPoint(pressure)-referencetemperature);
-}
-/*}}}*/
 IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
 	if (enthalpy<PureIceEnthalpy(pressure))
 		return thermalconductivity/heatcapacity;
@@ -457,27 +381,69 @@
 	return kappa;
 }
 /*}}}*/
-void Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
+IssmDouble Matpar::GetMaterialParameter(int enum_in){/*{{{*/
 
-	/*Ouput*/
-	IssmDouble temperature,waterfraction;
-
-	if(enthalpy<PureIceEnthalpy(pressure)){
-		temperature=referencetemperature+enthalpy/heatcapacity;
-		waterfraction=0.;
+	switch(enum_in){
+		case MaterialsRhoIceEnum:                    return this->rho_ice;
+		case MaterialsRhoSeawaterEnum:               return this->rho_water;
+		case MaterialsRhoFreshwaterEnum:             return this->rho_freshwater;
+		case MaterialsMuWaterEnum:                   return this->mu_water;
+		case MaterialsHeatcapacityEnum:              return this->heatcapacity;
+		case MaterialsThermalconductivityEnum:       return this->thermalconductivity;
+		case MaterialsTemperateiceconductivityEnum:  return this->temperateiceconductivity;
+		case MaterialsLatentheatEnum:                return this->latentheat;
+		case MaterialsBetaEnum:                      return this->beta;
+		case MaterialsMeltingpointEnum:              return this->meltingpoint;
+		case ConstantsReferencetemperatureEnum:      return this->referencetemperature;
+		case MaterialsMixedLayerCapacityEnum:        return this->mixed_layer_capacity;
+		case MaterialsThermalExchangeVelocityEnum:   return this->thermal_exchange_velocity;
+		case HydrologydcSedimentPorosityEnum:        return this->sediment_porosity;
+		case HydrologydcSedimentThicknessEnum:       return this->sediment_thickness;
+		case HydrologydcSedimentCompressibilityEnum: return this->sediment_compressibility;
+		case HydrologydcEplPorosityEnum:             return this->epl_porosity;
+		case HydrologydcEplCompressibilityEnum:      return this->epl_compressibility;
+		case HydrologydcEplConductivityEnum:         return this->epl_conductivity;
+		case HydrologydcEplInitialThicknessEnum:     return this->epl_init_thickness;
+		case HydrologydcEplColapseThicknessEnum:     return this->epl_colapse_thickness;
+		case HydrologydcEplMaxThicknessEnum:         return this->epl_max_thickness;
+		case HydrologydcWaterCompressibilityEnum:    return this->water_compressibility;
+		case ConstantsGEnum:                         return this->g;
+		case MaterialsPoissonEnum:                   return this->poisson;
+		case MaterialsYoungModulusEnum:              return this->young_modulus;
+		case MaterialsRidgingExponentEnum:				return this->ridging_exponent;
+		case MaterialsCohesionEnum:		            return this->cohesion;
+		case MaterialsInternalFrictionCoefEnum:		return this->internal_friction_coef;
+		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;
+		case MaterialsTimeRelaxationDamageEnum:		return this->time_relaxation_damage;
+		default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");
 	}
-	else{
-		temperature=TMeltingPoint(pressure);
-		waterfraction=(enthalpy-PureIceEnthalpy(pressure))/latentheat;
-	}
 
-	/*Assign output pointers:*/
-	*pwaterfraction=waterfraction;
-	*ptemperature=temperature;
 }
 /*}}}*/
-void Matpar::ThermalToEnthalpy(IssmDouble * penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
+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;
 
@@ -492,25 +458,7 @@
 	*penthalpy=enthalpy;
 }
 /*}}}*/
-IssmDouble Matpar::GetLithosphereShearModulus(){		 /*{{{*/
-	return lithosphere_shear_modulus;			 
-}		 
-/*}}}*/ 
-IssmDouble Matpar::GetLithosphereDensity(){		 /*{{{*/
-	return lithosphere_density;			 
-}		 
-/*}}}*/ 
-IssmDouble Matpar::GetMantleDensity(){		 /*{{{*/
-	return mantle_density;			 
-}		 
-/*}}}*/ 
-IssmDouble Matpar::GetMantleShearModulus(){		 /*{{{*/
-	return mantle_shear_modulus;			 
-}		 
-/*}}}*/ 
-void  Matpar::ResetHooks(){/*{{{*/
-
-	//Nothing to be done
-	return;
+IssmDouble Matpar::TMeltingPoint(IssmDouble pressure){/*{{{*/
+	return meltingpoint-beta*pressure;
 }
 /*}}}*/
Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h
===================================================================
--- ../trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 18945)
+++ ../trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 18946)
@@ -109,28 +109,13 @@
 		void       ResetHooks();
 		/*}}}*/
 		/*Numerics: {{{*/
-		IssmDouble GetG();
-		IssmDouble GetRhoIce();
-		IssmDouble GetRhoWater();
-		IssmDouble GetRhoFreshwater();
-		IssmDouble GetBeta();
-		IssmDouble GetMeltingPoint();
-		IssmDouble TMeltingPoint(IssmDouble pressure);
-		IssmDouble PureIceEnthalpy(IssmDouble pressure);
+		void       EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
 		IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
 		IssmDouble GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
-		IssmDouble GetLithosphereShearModulus();
-		IssmDouble GetLithosphereDensity();
-		IssmDouble GetMantleShearModulus();
-		IssmDouble GetMantleDensity();
-		void       EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
+		IssmDouble GetMaterialParameter(int in_enum); 
+		IssmDouble PureIceEnthalpy(IssmDouble pressure);
 		void       ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
-		IssmDouble GetDesFac();
-		IssmDouble GetS0p(); 
-		IssmDouble GetS0t(); 
-		IssmDouble GetRlaps(); 
-		IssmDouble GetRlapslgm(); 
-		IssmDouble GetMaterialParameter(int in_enum); 
+		IssmDouble TMeltingPoint(IssmDouble pressure);
 		/*}}}*/
 
 };
Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 18945)
+++ ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 18946)
@@ -575,7 +575,7 @@
 	parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
 
 	/*Compute pressure melting point*/
-	t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
+	t_pmp=matpar->GetMaterialParameter(MaterialsMeltingpointEnum)-matpar->GetMaterialParameter(MaterialsBetaEnum)*pressure;
 
 	/*Add penalty load*/
 	if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be  constrained to 0 when T<Tpmp, instead of using spcs, use penalties
@@ -650,7 +650,7 @@
 	parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
 
 	/*Compute pressure melting point*/
-	t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure;
+	t_pmp=matpar->GetMaterialParameter(MaterialsMeltingpointEnum)-matpar->GetMaterialParameter(MaterialsBetaEnum)*pressure;
 
 	/*Add penalty load
 	  This time, the penalty must have the same value as the one used for the thermal computation
@@ -685,7 +685,7 @@
 	parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum);
 
 	/*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: ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp	(revision 18945)
+++ ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp	(revision 18946)
@@ -504,9 +504,9 @@
 	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
 
 	/*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);
 	if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts");
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18945)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18946)
@@ -757,9 +757,9 @@
 			sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
 
 			/*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;
 
 			/*Compare basal stress to water pressure and determine whether it should ground*/
@@ -1351,8 +1351,8 @@
 
 	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);
 
 	/*First calculate the area of the base (cross section triangle)
@@ -1769,7 +1769,7 @@
 	thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
 	
 	/*Retrieve material parameters: */
-	rho_ice=matpar->GetRhoIce();
+	rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
 
 	/*Retrieve values of the levelset defining the masscon: */
 	values = xNew<IssmDouble>(NUMVERTICES);
@@ -1812,7 +1812,7 @@
 	GaussTria* gauss_2=NULL;
 
 	/*Get material parameters :*/
-	rho_ice=matpar->GetRhoIce();
+	rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
 
 	/*First off, check that this segment belongs to this element: */
 	if (segment_id!=this->id)_error_("error message: segment with id " << segment_id << " does not belong to element with id:" << this->id);
@@ -1875,7 +1875,7 @@
 	GaussTria* gauss_2=NULL;
 
 	/*Get material parameters :*/
-	rho_ice=matpar->GetRhoIce();
+	rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
 
 	/*First off, check that this segment belongs to this element: */
 	if (reCast<int>(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast<int>(*(segment+4)) << " does not belong to element with id:" << this->id);
@@ -2251,15 +2251,15 @@
   GetInputListOnVertices(&s[0],SurfaceEnum);
 
   /*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*/
    for (int iv = 0; iv<NUMVERTICES; iv++){
@@ -2284,8 +2284,8 @@
 	IssmDouble  rho_water,rho_ice,density;
 
 	/*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);
 	GetInputListOnVertices(&r[0],BedEnum);
@@ -2707,7 +2707,7 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 
 	/*Get material parameters :*/
-	rho_ice=matpar->GetRhoIce();
+	rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
 
    if(!IsIceInElement())return 0;
 
@@ -3100,11 +3100,11 @@
 	this->parameters->FindParam(&currenttime,TimeEnum);
 
 	/*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: */
 	thickness_input=inputs->GetInput(ThicknessEnum); 
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18945)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18946)
@@ -251,8 +251,8 @@
 	}
 
 	/*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: */
 	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
@@ -696,9 +696,9 @@
 			  + 2.*sigma_xy*bed_normal[0]*bed_normal[1] + 2.*sigma_xz*bed_normal[0]*bed_normal[2] + 2.*sigma_yz*bed_normal[1]*bed_normal[2];
 
 			/*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;
 
 			/*Compare basal stress to water pressure and determine whether it should ground*/
@@ -1176,8 +1176,8 @@
 
 	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);
 
 	/*First calculate the area of the base (cross section triangle)
@@ -2048,15 +2048,15 @@
   GetInputListOnVertices(&s[0],SurfaceEnum);
 
   /*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*/
    for (int iv = 0; iv < NUMVERTICES; iv++){
@@ -2082,8 +2082,8 @@
 	IssmDouble  rho_water,rho_ice,density;
 
 	/*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);
 	GetInputListOnVertices(&r[0],BedEnum);
@@ -2640,7 +2640,7 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 
 	/*Get material parameters :*/
-	rho_ice=matpar->GetRhoIce();
+	rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
 
 	if(!IsIceInElement() || !IsOnSurface()) return 0.;
 
Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18945)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18946)
@@ -1225,8 +1225,8 @@
 	GetInputListOnVertices(&b[0],BaseEnum);
 	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;
 
 	/*go through vertices, and update inputs, considering them to be TriaVertex type: */
@@ -2053,7 +2053,7 @@
 	surface_input->GetInputValue(&s,gauss);
 	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*/
 	this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
