Index: ../trunk-jpl/src/m/dev/issmversion.m
===================================================================
--- ../trunk-jpl/src/m/dev/issmversion.m	(revision 19047)
+++ ../trunk-jpl/src/m/dev/issmversion.m	(revision 19048)
@@ -16,7 +16,7 @@
 disp(['Build date: ' IssmConfig('PACKAGE_BUILD_DATE')]);
 disp(['Compiled on ' IssmConfig('HOST_VENDOR') ' ' IssmConfig('HOST_OS') ' ' IssmConfig('HOST_ARCH') ' by ' IssmConfig('USER_NAME')]);
 disp([' ']);
-disp(['Copyright (c) 2009-2014 California Institute of Technology']);
+disp(['Copyright (c) 2009-2015 California Institute of Technology']);
 disp([' ']);
 disp(['    to get started type: issmdoc']);
 disp([' ']);
Index: ../trunk-jpl/src/m/classes/initialization.m
===================================================================
--- ../trunk-jpl/src/m/classes/initialization.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/initialization.m	(revision 19048)
@@ -35,6 +35,21 @@
 			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="watercolumn" type="',class(self.watercolumn),'" default="',self.watercolumn,'">','     <section name="initialization" />','     <help> thickness of subglacial water [m] </help>','</parameter>');
 			fprintf(fid,'%s\n%s\n','</frame>');
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.vx=project3d(md,'vector',self.vx,'type','node');
+			self.vy=project3d(md,'vector',self.vy,'type','node');
+			self.vz=project3d(md,'vector',self.vz,'type','node');
+			self.vel=project3d(md,'vector',self.vel,'type','node');
+			self.temperature=project3d(md,'vector',self.temperature,'type','node');
+			self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node');
+			self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node','layer',1);
+			self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1);
+			self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1);
+			self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1);
+
+			%Lithostatic pressure by default
+			self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
+		end % }}}
 		function self = initialization(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/m1qn3inversion.py
===================================================================
--- ../trunk-jpl/src/m/classes/m1qn3inversion.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/m1qn3inversion.py	(revision 19048)
@@ -1,4 +1,5 @@
 import numpy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
 from StringToEnum import StringToEnum
@@ -93,6 +94,19 @@
 		string="%s\n%s"%(string,'   503: ThicknessAbsGradient')
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node')
+		self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node')
+		self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node')
+		self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node')
+		if not numpy.any(numpy.isnan(self.cost_functions_coefficients)):
+			self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node')
+		if not numpy.any(numpy.isnan(self.min_parameters)):
+			self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node')
+		if not numpy.any(numpy.isnan(self.max_parameters)):
+			self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		
 		#default is incomplete adjoint for now
Index: ../trunk-jpl/src/m/classes/thermal.py
===================================================================
--- ../trunk-jpl/src/m/classes/thermal.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/thermal.py	(revision 19048)
@@ -1,4 +1,5 @@
 import numpy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
 from checkfield import checkfield
@@ -42,6 +43,14 @@
 		string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.spctemperature=project3d(md,'vector',self.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',numpy.nan)
+		if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
+			self.spctemperature=numpy.nan*numpy.ones((md.mesh.numberofvertices,1))
+			pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
+			self.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
+		return self
+	#}}}
 	def defaultoutputs(self,md): # {{{
 
 		if self.isenthalpy:
Index: ../trunk-jpl/src/m/classes/gia.py
===================================================================
--- ../trunk-jpl/src/m/classes/gia.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/gia.py	(revision 19048)
@@ -1,4 +1,5 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from checkfield import checkfield
 from WriteData import WriteData
@@ -29,6 +30,11 @@
 		string="%s\n%s"%(string,fielddisplay(self,'cross_section_shape',"1: square-edged, 2: elliptical-edged surface"))
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.mantle_viscosity=project3d(md,'vector',self.mantle_viscosity,'type','node')
+		self.lithosphere_thickness=project3d(md,'vector',self.lithosphere_thickness,'type','node')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 
 		self.cross_section_shape=1; 
Index: ../trunk-jpl/src/m/classes/matdamageice.py
===================================================================
--- ../trunk-jpl/src/m/classes/matdamageice.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/matdamageice.py	(revision 19048)
@@ -1,4 +1,5 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import MaterialsEnum, MatdamageiceEnum, MaterialsRheologyLawEnum, MaterialsRhoSeawaterEnum
 from StringToEnum import StringToEnum
 from checkfield import checkfield
@@ -62,6 +63,11 @@
 
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node')
+		self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		#ice density (kg/m^3)
 		self.rho_ice=917.
Index: ../trunk-jpl/src/m/classes/frictionweertman.py
===================================================================
--- ../trunk-jpl/src/m/classes/frictionweertman.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/frictionweertman.py	(revision 19048)
@@ -1,4 +1,5 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from checkfield import checkfield
 from WriteData import WriteData
Index: ../trunk-jpl/src/m/classes/frictionwaterlayer.m
===================================================================
--- ../trunk-jpl/src/m/classes/frictionwaterlayer.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/frictionwaterlayer.m	(revision 19048)
@@ -35,6 +35,12 @@
 			md = checkfield(md,'fieldname','thermal.spctemperature','timeseries',1,'>=',0.);
 
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1);
+			self.p=project3d(md,'vector',self.p,'type','element');
+			self.q=project3d(md,'vector',self.q,'type','element');
+			self.water_layer=project3d(md,'vector',self.water_layer,'type','node','layer',1);
+		end % }}}
 		function disp(self) % {{{
 			disp(sprintf('Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b * 1/f(T)\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(bed+water_layer), r=q/p and s=1/p)'));
 			fielddisplay(self,'coefficient','frictiontemp coefficient [SI]');
Index: ../trunk-jpl/src/m/classes/calving.py
===================================================================
--- ../trunk-jpl/src/m/classes/calving.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/calving.py	(revision 19048)
@@ -1,4 +1,5 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from StringToEnum import StringToEnum
 from checkfield import checkfield
@@ -30,6 +31,11 @@
 
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.calvingrate=project3d(md,'vector',self.calvingrate,'type','node')
+		self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 
 		#stabilization = 2 by default
Index: ../trunk-jpl/src/m/classes/m1qn3inversion.m
===================================================================
--- ../trunk-jpl/src/m/classes/m1qn3inversion.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/m1qn3inversion.m	(revision 19048)
@@ -26,6 +26,15 @@
 
 	end
 	methods
+		function self = extrude(self,md) % {{{
+			self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
+			self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
+			self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
+			self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
+			if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
+			if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
+			if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
+		end % }}}
 		function self = m1qn3inversion(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/thermal.m
===================================================================
--- ../trunk-jpl/src/m/classes/thermal.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/thermal.m	(revision 19048)
@@ -17,30 +17,38 @@
 		requested_outputs = {};
 	end
 	methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '<!-- thermal -->\n');            
-                    
-            % thermal solution parameters
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Thermal solution parameters">','<section name="thermal" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="spctemperature" type="',class(self.spctemperature),'" default="',convert2str(self.spctemperature),'">','     <section name="thermal" />','     <help> temperature constraints (NaN means no constraint) [K] </help>','  </parameter>');
-                
-            % penalty_threshold drop-down (0, 1, or 2)
-            fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="penalty_threshold" type="alternative"  optional="false">','     <section name="thermal" />','     <help> 0: no, 1: artificial_diffusivity, 2: SUPG </help>');
-            fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n','       <option value="1" type="string" default="false"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="2" type="string" default="false"> </option>','</parameter>');
+		function createxml(self,fid) % {{{
+			fprintf(fid, '<!-- thermal -->\n');            
 
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="stabilization" type="',class(self.stabilization),'" default="',convert2str(self.stabilization),'">','     <section name="thermal" />','     <help> maximum number of non linear iterations </help>','  </parameter>');
-				fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="reltol" type="',class(self.reltol),'" default="',convert2str(self.reltol),'">','     <section name="steadystate" />','     <help> relative tolerance criterion [K] </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter" type="',class(self.maxiter),'" default="',convert2str(self.maxiter),'">','     <section name="thermal" />','     <help> stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization) </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_lock" type="',class(self.penalty_lock),'" default="',convert2str(self.penalty_lock),'">','     <section name="thermal" />','     <help> threshold to declare convergence of thermal solution (default is 0)  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_factor" type="',class(self.penalty_factor),'" default="',convert2str(self.penalty_factor),'">','     <section name="thermal" />','     <help> scaling exponent (default is 3) </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isenthalpy" type="',class(self.isenthalpy),'" default="',convert2str(self.isenthalpy),'">','     <section name="thermal" />','     <help> use an enthalpy formulation to include temperate ice (default is 0) </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isdynamicbasalspc" type="',class(self.isdynamicbasalspc),'" default="',convert2str(self.isdynamicbasalspc),'">','     <section name="thermal" />','     <help> enable dynamic setting of basal forcing. recommended for enthalpy formulation (default is 0)  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="requested_outputs" type="',class(self.requested_outputs),'" default="',convert2str(self.requested_outputs),'">','     <section name="thermal" />','     <help> additional outputs requested </help>','  </parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');    
-        
-        end % }}}
+			% thermal solution parameters
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Thermal solution parameters">','<section name="thermal" />');                    
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="spctemperature" type="',class(self.spctemperature),'" default="',convert2str(self.spctemperature),'">','     <section name="thermal" />','     <help> temperature constraints (NaN means no constraint) [K] </help>','  </parameter>');
+
+			% penalty_threshold drop-down (0, 1, or 2)
+			fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="penalty_threshold" type="alternative"  optional="false">','     <section name="thermal" />','     <help> 0: no, 1: artificial_diffusivity, 2: SUPG </help>');
+			fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
+			fprintf(fid,'%s\n','       <option value="1" type="string" default="false"> </option>');
+			fprintf(fid,'%s\n%s\n','       <option value="2" type="string" default="false"> </option>','</parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="stabilization" type="',class(self.stabilization),'" default="',convert2str(self.stabilization),'">','     <section name="thermal" />','     <help> maximum number of non linear iterations </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="reltol" type="',class(self.reltol),'" default="',convert2str(self.reltol),'">','     <section name="steadystate" />','     <help> relative tolerance criterion [K] </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter" type="',class(self.maxiter),'" default="',convert2str(self.maxiter),'">','     <section name="thermal" />','     <help> stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization) </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_lock" type="',class(self.penalty_lock),'" default="',convert2str(self.penalty_lock),'">','     <section name="thermal" />','     <help> threshold to declare convergence of thermal solution (default is 0)  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_factor" type="',class(self.penalty_factor),'" default="',convert2str(self.penalty_factor),'">','     <section name="thermal" />','     <help> scaling exponent (default is 3) </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isenthalpy" type="',class(self.isenthalpy),'" default="',convert2str(self.isenthalpy),'">','     <section name="thermal" />','     <help> use an enthalpy formulation to include temperate ice (default is 0) </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isdynamicbasalspc" type="',class(self.isdynamicbasalspc),'" default="',convert2str(self.isdynamicbasalspc),'">','     <section name="thermal" />','     <help> enable dynamic setting of basal forcing. recommended for enthalpy formulation (default is 0)  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="requested_outputs" type="',class(self.requested_outputs),'" default="',convert2str(self.requested_outputs),'">','     <section name="thermal" />','     <help> additional outputs requested </help>','  </parameter>');
+			fprintf(fid,'%s\n%s\n','</frame>');    
+
+		end % }}}
+		function self = extrude(self,md) % {{{
+			self.spctemperature=project3d(md,'vector',self.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN);
+			if (length(md.initialization.temperature)==md.mesh.numberofvertices),
+				self.spctemperature=NaN(md.mesh.numberofvertices,1);
+				pos=find(md.mesh.vertexonsurface);
+				self.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
+			end
+		end % }}}
 		function self = thermal(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/initialization.py
===================================================================
--- ../trunk-jpl/src/m/classes/initialization.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/initialization.py	(revision 19048)
@@ -1,4 +1,5 @@
 import numpy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
 from checkfield import checkfield
@@ -47,6 +48,22 @@
 
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.vx=project3d(md,'vector',self.vx,'type','node')
+		self.vy=project3d(md,'vector',self.vy,'type','node')
+		self.vz=project3d(md,'vector',self.vz,'type','node')
+		self.vel=project3d(md,'vector',self.vel,'type','node')
+		self.temperature=project3d(md,'vector',self.temperature,'type','node')
+		self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node')
+		self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node')
+		self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1)
+		self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1)
+		self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1)
+
+		#Lithostatic pressure by default
+		self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1))
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		return self
 	#}}}
Index: ../trunk-jpl/src/m/classes/mask.py
===================================================================
--- ../trunk-jpl/src/m/classes/mask.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/mask.py	(revision 19048)
@@ -1,5 +1,6 @@
 import numpy
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from checkfield import checkfield
 from WriteData import WriteData
@@ -27,6 +28,11 @@
 		string="%s\n%s"%(string,fielddisplay(self,"ice_levelset","presence of ice if < 0, icefront position if = 0, no ice if > 0"))
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.ice_levelset=project3d(md,'vector',self.ice_levelset,'type','node')
+		self.groundedice_levelset=project3d(md,'vector',self.groundedice_levelset,'type','node')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		return self
 	#}}}
Index: ../trunk-jpl/src/m/classes/calvinglevermann.py
===================================================================
--- ../trunk-jpl/src/m/classes/calvinglevermann.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/calvinglevermann.py	(revision 19048)
@@ -29,6 +29,11 @@
 
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.coeff=project3d(md,'vector',self.coeff,'type','node')
+		self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 
 		#stabilization = 2 by default
Index: ../trunk-jpl/src/m/classes/qmu.m
===================================================================
--- ../trunk-jpl/src/m/classes/qmu.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/qmu.m	(revision 19048)
@@ -23,35 +23,38 @@
 		vertex_weight               = NaN;
 	end
 	methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '<!-- qmu -->\n');            
-                    
-            % qmu solution parameters
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="qmu parameters">','<section name="qmu" />');                    
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isdakota" type="',class(self.isdakota),'" default="',convert2str(self.isdakota),'">','     <section name="qmu" />','     <help> is qmu analysis activated? </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="variables" type="',class(self.variables),'" default="',convert2str(self.variables),'">','     <section name="qmu" />','     <help> (arrays of each variable class) </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="responses" type="',class(self.responses),'" default="',convert2str(self.responses),'">','     <section name="qmu" />','     <help> (arrays of each response class) </help>','  </parameter>');
-                
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="numberofresponses" type="',class(self.numberofresponses),'" default="',convert2str(self.numberofresponses),'">','     <section name="qmu" />','     <help> number of responses </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="params" type="',class(self.params),'" default="',convert2str(self.params),'">','     <section name="qmu" />','     <help> (array of method-independent parameters)  </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="results" type="',class(self.results),'" default="',convert2str(self.results),'">','     <section name="qmu" />','     <help> (information from dakota files) </help>','  </parameter>');
-                
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="partition" type="',class(self.partition),'" default="',convert2str(self.partition),'">','     <section name="qmu" />','     <help> user provided mesh partitioning, defaults to metis if not specified </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="numberofpartitions" type="',class(self.numberofpartitions),'" default="',convert2str(self.numberofpartitions),'">','     <section name="qmu" />','     <help> number of partitions for semi-discrete qmu  </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="variabledescriptors" type="',class(self.variabledescriptors),'" default="',convert2str(self.variabledescriptors),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
-            
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="responsedescriptors" type="',class(self.responsedescriptors),'" default="',convert2str(self.responsedescriptors),'">','     <section name="qmu" />','     <help> use an enthalpy formulation to include temperate ice (default is 0) </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="method" type="',class(self.method),'" default="',convert2str(self.method),'">','     <section name="qmu" />','     <help> array of dakota_method class </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_profile_directory" type="',class(self.mass_flux_profile_directory),'" default="',convert2str(self.mass_flux_profile_directory),'">','     <section name="qmu" />','     <help> directory for mass flux profiles </help>','  </parameter>');
-                
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_profiles" type="',class(self.mass_flux_profiles),'" default="',convert2str(self.mass_flux_profiles),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_segments" type="',class(self.mass_flux_segments),'" default="',convert2str(self.mass_flux_segments),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="adjacency" type="',class(self.adjacency),'" default="',convert2str(self.adjacency),'">','     <section name="qmu" />','     <help> additional outputs requested </help>','  </parameter>');
-                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vertex_weight" type="',class(self.vertex_weight),'" default="',convert2str(self.vertex_weight),'">','     <section name="qmu" />','     <help> weight applied to each mesh vertex </help>','  </parameter>');
-            
-            fprintf(fid,'%s\n%s\n','</frame>');    
-        
-        end % }}}        
+		function createxml(self,fid) % {{{
+			fprintf(fid, '<!-- qmu -->\n');            
+
+			% qmu solution parameters
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="qmu parameters">','<section name="qmu" />');                    
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isdakota" type="',class(self.isdakota),'" default="',convert2str(self.isdakota),'">','     <section name="qmu" />','     <help> is qmu analysis activated? </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="variables" type="',class(self.variables),'" default="',convert2str(self.variables),'">','     <section name="qmu" />','     <help> (arrays of each variable class) </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="responses" type="',class(self.responses),'" default="',convert2str(self.responses),'">','     <section name="qmu" />','     <help> (arrays of each response class) </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="numberofresponses" type="',class(self.numberofresponses),'" default="',convert2str(self.numberofresponses),'">','     <section name="qmu" />','     <help> number of responses </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="params" type="',class(self.params),'" default="',convert2str(self.params),'">','     <section name="qmu" />','     <help> (array of method-independent parameters)  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="results" type="',class(self.results),'" default="',convert2str(self.results),'">','     <section name="qmu" />','     <help> (information from dakota files) </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="partition" type="',class(self.partition),'" default="',convert2str(self.partition),'">','     <section name="qmu" />','     <help> user provided mesh partitioning, defaults to metis if not specified </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="numberofpartitions" type="',class(self.numberofpartitions),'" default="',convert2str(self.numberofpartitions),'">','     <section name="qmu" />','     <help> number of partitions for semi-discrete qmu  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="variabledescriptors" type="',class(self.variabledescriptors),'" default="',convert2str(self.variabledescriptors),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="responsedescriptors" type="',class(self.responsedescriptors),'" default="',convert2str(self.responsedescriptors),'">','     <section name="qmu" />','     <help> use an enthalpy formulation to include temperate ice (default is 0) </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="method" type="',class(self.method),'" default="',convert2str(self.method),'">','     <section name="qmu" />','     <help> array of dakota_method class </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_profile_directory" type="',class(self.mass_flux_profile_directory),'" default="',convert2str(self.mass_flux_profile_directory),'">','     <section name="qmu" />','     <help> directory for mass flux profiles </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_profiles" type="',class(self.mass_flux_profiles),'" default="',convert2str(self.mass_flux_profiles),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mass_flux_segments" type="',class(self.mass_flux_segments),'" default="',convert2str(self.mass_flux_segments),'">','     <section name="qmu" />','     <help>  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="adjacency" type="',class(self.adjacency),'" default="',convert2str(self.adjacency),'">','     <section name="qmu" />','     <help> additional outputs requested </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vertex_weight" type="',class(self.vertex_weight),'" default="',convert2str(self.vertex_weight),'">','     <section name="qmu" />','     <help> weight applied to each mesh vertex </help>','  </parameter>');
+
+			fprintf(fid,'%s\n%s\n','</frame>');    
+
+		end % }}}        
+		function self = extrude(self,md) % {{{
+			self.partition=project3d(md,'vector',self.partition','type','node');
+		end % }}}
 		function self = qmu(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/matice.m
===================================================================
--- ../trunk-jpl/src/m/classes/matice.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/matice.m	(revision 19048)
@@ -70,6 +70,10 @@
 
 
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node');
+			self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element');
+		end % }}}
 		function self = matice(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/masstransport.py
===================================================================
--- ../trunk-jpl/src/m/classes/masstransport.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/masstransport.py	(revision 19048)
@@ -1,4 +1,5 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from StringToEnum import StringToEnum
 from checkfield import checkfield
@@ -37,6 +38,10 @@
 
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.spcthickness=project3d(md,'vector',self.spcthickness,'type','node')
+		return self
+	#}}}
 	def defaultoutputs(self,md): # {{{
 
 		return ['Thickness','Surface','Base']
Index: ../trunk-jpl/src/m/classes/calving.m
===================================================================
--- ../trunk-jpl/src/m/classes/calving.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/calving.m	(revision 19048)
@@ -28,6 +28,10 @@
 					error('constructor not supported');
 			end
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.calvingrate=project3d(md,'vector',self.calvingrate,'type','node');
+			self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
+		end % }}}
 		function self = setdefaultparameters(self) % {{{
 
 			%stabilization = 2 by default
Index: ../trunk-jpl/src/m/classes/inversionvalidation.m
===================================================================
--- ../trunk-jpl/src/m/classes/inversionvalidation.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/inversionvalidation.m	(revision 19048)
@@ -21,6 +21,15 @@
 		surface_obs                 = NaN
 	end
 	methods
+		function self = extrude(self,md) % {{{
+			self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
+			self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
+			self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
+			self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
+			if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
+			if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
+			if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
+		end % }}}
 		function self = inversionvalidation(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/geometry.py
===================================================================
--- ../trunk-jpl/src/m/classes/geometry.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/geometry.py	(revision 19048)
@@ -1,3 +1,4 @@
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
 from checkfield import checkfield
@@ -31,6 +32,14 @@
 		string="%s\n%s"%(string,fielddisplay(self,'bed','bed elevation [m]'))
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.surface=project3d(md,'vector',self.surface,'type','node')
+		self.thickness=project3d(md,'vector',self.thickness,'type','node')
+		self.hydrostatic_ratio=project3d(md,'vector',self.hydrostatic_ratio,'type','node')
+		self.base=project3d(md,'vector',self.base,'type','node')
+		self.bed=project3d(md,'vector',self.bed,'type','node')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		return self
 	#}}}
Index: ../trunk-jpl/src/m/classes/stressbalance.m
===================================================================
--- ../trunk-jpl/src/m/classes/stressbalance.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/stressbalance.m	(revision 19048)
@@ -69,6 +69,14 @@
 
 			fprintf(fid,'%s\n%s\n','</frame>'); 
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.spcvx=project3d(md,'vector',self.spcvx,'type','node');
+			self.spcvy=project3d(md,'vector',self.spcvy,'type','node');
+			self.spcvz=project3d(md,'vector',self.spcvz,'type','node');
+			self.referential=project3d(md,'vector',self.referential,'type','node');
+			self.loadingforce=project3d(md,'vector',self.loadingforce,'type','node');
+
+		end % }}}
 		function self = stressbalance(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/hydrologyshreve.py
===================================================================
--- ../trunk-jpl/src/m/classes/hydrologyshreve.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/hydrologyshreve.py	(revision 19048)
@@ -26,6 +26,9 @@
 		string="%s\n%s"%(string,fielddisplay(self,'stabilization','artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.'))
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		
 		#Type of stabilization to use 0:nothing 1:artificial_diffusivity
Index: ../trunk-jpl/src/m/classes/hydrologyshreve.m
===================================================================
--- ../trunk-jpl/src/m/classes/hydrologyshreve.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/hydrologyshreve.m	(revision 19048)
@@ -9,17 +9,19 @@
 		stabilization  = 0;
 	end
 	methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- Hydrology -->');
-            
-            % Convergence criteria            
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Hydrologyshreve solution parameters">','<section name="hydrologyshreve" />');                    
-            
+		function createxml(self,fid) % {{{
+			fprintf(fid, '\n\n');
+			fprintf(fid, '%s\n', '<!-- Hydrology -->');
+
+			% Convergence criteria            
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Hydrologyshreve solution parameters">','<section name="hydrologyshreve" />');                    
+
 			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="spcwatercolumn" type="',      class(self.spcwatercolumn),'" default="',        convert2str(self.spcwatercolumn),'">',	'     <section name="hydrologyshreve" />','     <help> water thickness constraints (NaN means no constraint) [m] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="stabilization" type="',         class(self.stabilization),'" default="',           convert2str(self.stabilization),'">',	'     <section name="hydrologyshreve" />','     <help> artificial diffusivity (default is 1). can be more than 1 to increase diffusivity. </help>','</parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');
-        end % }}}
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="stabilization" type="',         class(self.stabilization),'" default="',           convert2str(self.stabilization),'">',	'     <section name="hydrologyshreve" />','     <help> artificial diffusivity (default is 1). can be more than 1 to increase diffusivity. </help>','</parameter>');
+			fprintf(fid,'%s\n%s\n','</frame>');
+		end % }}}
+		function self = extrude(self,md) % {{{
+		end % }}}
 		function self = hydrologyshreve(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/model.m
===================================================================
--- ../trunk-jpl/src/m/classes/model.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/model.m	(revision 19048)
@@ -749,106 +749,28 @@
 
 			%Ok, now deal with the other fields from the 2d mesh:
 
+			%bedinfo and surface info
+			md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
+			md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
+			md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
+
 			%lat long
 			md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
 			md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');
 
-			%drag coefficient is limited to nodes that are on the bedrock.
-			if isa(md.friction,'friction'),
-				md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1);
-				md.friction.p=project3d(md,'vector',md.friction.p,'type','element');
-				md.friction.q=project3d(md,'vector',md.friction.q,'type','element');
-			elseif isa(md.friction,'frictionhydro'),
-				md.friction.q=project3d(md,'vector',md.friction.q,'type','element');
-				md.friction.C=project3d(md,'vector',md.friction.C,'type','element');
-				md.friction.As=project3d(md,'vector',md.friction.As,'type','element');
-				md.friction.effective_pressure=project3d(md,'vector',md.friction.effective_pressure,'type','node','layer',1);
-			elseif isa(md.friction,'frictionwaterlayer'),
-				md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1);
-				md.friction.p=project3d(md,'vector',md.friction.p,'type','element');
-				md.friction.q=project3d(md,'vector',md.friction.q,'type','element');
-				md.friction.water_layer=project3d(md,'vector',md.friction.water_layer,'type','node','layer',1);
-			elseif isa(md.friction,'frictionweertman'),
-				md.friction.C=project3d(md,'vector',md.friction.C,'type','node','layer',1);
-				md.friction.m=project3d(md,'vector',md.friction.m,'type','element');
-	    end
-	 
-			%observations
-			md.inversion.vx_obs=project3d(md,'vector',md.inversion.vx_obs,'type','node');
-			md.inversion.vy_obs=project3d(md,'vector',md.inversion.vy_obs,'type','node');
-			md.inversion.vel_obs=project3d(md,'vector',md.inversion.vel_obs,'type','node');
-			md.inversion.thickness_obs=project3d(md,'vector',md.inversion.thickness_obs,'type','node');
+			md.geometry=extrude(md.geometry,md);
+			md.friction  = extrude(md.friction,md);
+			md.inversion = extrude(md.inversion,md);
 			md.surfaceforcings = extrude(md.surfaceforcings,md);
-			md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node');
+			md.initialization = extrude(md.initialization,md);
 
-			%results
-			if ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end;
-			if ~isnan(md.initialization.vy),md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node');end;
-			if ~isnan(md.initialization.vz),md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node');end;
-			if ~isnan(md.initialization.vel),md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node');end;
-			if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end;
-			if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end;
-      if ~isnan(md.initialization.watercolumn),md.initialization.watercolumn=project3d(md,'vector',md.initialization.watercolumn,'type','node','layer',1);end;
-      if ~isnan(md.initialization.sediment_head),md.initialization.sediment_head=project3d(md,'vector',md.initialization.sediment_head,'type','node','layer',1);end;
-      if ~isnan(md.initialization.epl_head),md.initialization.epl_head=project3d(md,'vector',md.initialization.epl_head,'type','node','layer',1);end;
-      if ~isnan(md.initialization.epl_thickness),md.initialization.epl_thickness=project3d(md,'vector',md.initialization.epl_thickness,'type','node','layer',1);end;
+			md.flowequation.extrude(md);
+			md.stressbalance=extrude(md.stressbalance,md);
+			md.thermal.extrude(md);
+			md.masstransport.extrude(md);
+			md.calving=extrude(md.calving,md);
+			md.hydrology = extrude(md.hydrology,md);
 
-			%bedinfo and surface info
-			md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
-			md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
-
-			%elementstype
-			if ~isnan(md.flowequation.element_equation)
-				oldelements_type=md.flowequation.element_equation;
-				md.flowequation.element_equation=zeros(number_el3d,1);
-				md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element');
-			end
-
-			%verticestype
-			if ~isnan(md.flowequation.vertex_equation)
-				oldvertices_type=md.flowequation.vertex_equation;
-				md.flowequation.vertex_equation=zeros(number_nodes3d,1);
-				md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node');
-			end
-			md.flowequation.borderSSA=project3d(md,'vector',md.flowequation.borderSSA,'type','node');
-			md.flowequation.borderHO=project3d(md,'vector',md.flowequation.borderHO,'type','node');
-			md.flowequation.borderFS=project3d(md,'vector',md.flowequation.borderFS,'type','node');
-
-			%boundary conditions
-			md.stressbalance.spcvx=project3d(md,'vector',md.stressbalance.spcvx,'type','node');
-			md.stressbalance.spcvy=project3d(md,'vector',md.stressbalance.spcvy,'type','node');
-			md.stressbalance.spcvz=project3d(md,'vector',md.stressbalance.spcvz,'type','node');
-			md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN);
-			if (length(md.initialization.temperature)==md.mesh.numberofvertices),
-				md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1);
-				if isprop(md.mesh,'vertexonsurface'),
-					pos=find(md.mesh.vertexonsurface);
-					md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
-				end
-			end
-			md.masstransport.spcthickness=project3d(md,'vector',md.masstransport.spcthickness,'type','node');
-			md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node');
-			md.damage.spcdamage=project3d(md,'vector',md.damage.spcdamage,'type','node');
-			md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node');
-			md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node');
-
-			% Calving variables
-			if isa(md.calving,'calving'),md.calving.calvingrate=project3d(md,'vector',md.calving.calvingrate,'type','node');end;
-			if ~isnan(md.calving.meltingrate),md.calving.meltingrate=project3d(md,'vector',md.calving.meltingrate,'type','node');end;
-			if isa(md.calving,'calvinglevermann'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end;
-			if isa(md.calving,'calvingpi'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end;
-
-			% Hydrologydc variables
-			if isa(md.hydrology,'hydrologydc');
-				md.hydrology.spcsediment_head=project3d(md,'vector',md.hydrology.spcsediment_head,'type','node','layer',1);
-				md.hydrology.mask_eplactive_node=project3d(md,'vector',md.hydrology.mask_eplactive_node,'type','node','layer',1);
-				md.hydrology.sediment_transmitivity=project3d(md,'vector',md.hydrology.sediment_transmitivity,'type','node','layer',1);
-				md.hydrology.basal_moulin_input=project3d(md,'vector',md.hydrology.basal_moulin_input,'type','node','layer',1);
-				if(md.hydrology.isefficientlayer==1);
-					md.hydrology.spcepl_head=project3d(md,'vector',md.hydrology.spcepl_head,'type','node','layer',1);
-		    end
-	    end
-
 			%connectivity
 			if ~isnan(md.mesh.elementconnectivity)
 				md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1);
@@ -860,41 +782,12 @@
 				md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0;
 			end
 
-			%materials
-			md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node');
-			md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element');
-				
-			%damage
-			md.damage.D=project3d(md,'vector',md.damage.D,'type','node');
+			md.materials=extrude(md.materials,md);
+			md.damage=extrude(md.damage,md);
+			md.mask=extrude(md.mask,md);
+			md.qmu=extrude(md.qmu,md);
+			md.basalforcings=extrude(md.basalforcings,md);
 
-			%parameters
-			md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node');
-			md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node');
-			md.gia.mantle_viscosity=project3d(md,'vector',md.gia.mantle_viscosity,'type','node');
-			md.gia.lithosphere_thickness=project3d(md,'vector',md.gia.lithosphere_thickness,'type','node');
-			md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node');
-			md.geometry.base=project3d(md,'vector',md.geometry.base,'type','node');
-			md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node');
-			md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
-			md.mask.groundedice_levelset=project3d(md,'vector',md.mask.groundedice_levelset,'type','node');
-			md.mask.ice_levelset=project3d(md,'vector',md.mask.ice_levelset,'type','node');
-			if ~isnan(md.inversion.cost_functions_coefficients),md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end;
-			if ~isnan(md.inversion.min_parameters),md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node');end;
-			if ~isnan(md.inversion.max_parameters),md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node');end;
-			if ~isnan(md.qmu.partition),md.qmu.partition=project3d(md,'vector',md.qmu.partition','type','node');end
-
-			%Put lithostatic pressure if there is an existing pressure
-			if ~isnan(md.initialization.pressure),
-				md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
-	    end
-
-			%special for thermal modeling:
-			md.basalforcings.groundedice_melting_rate=project3d(md,'vector',md.basalforcings.groundedice_melting_rate,'type','node','layer',1); 
-			md.basalforcings.floatingice_melting_rate=project3d(md,'vector',md.basalforcings.floatingice_melting_rate,'type','node','layer',1); 
-			if ~isnan(md.basalforcings.geothermalflux)
-				md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux
-			end
-
 			%increase connectivity if less than 25:
 			if md.mesh.average_vertex_connectivity<=25,
 				md.mesh.average_vertex_connectivity=100;
Index: ../trunk-jpl/src/m/classes/inversion.py
===================================================================
--- ../trunk-jpl/src/m/classes/inversion.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/inversion.py	(revision 19048)
@@ -1,4 +1,5 @@
 import numpy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
 from StringToEnum import StringToEnum
@@ -71,6 +72,19 @@
 		string="%s\n%s"%(string,'   503: ThicknessAbsGradient')
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node')
+		self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node')
+		self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node')
+		self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node')
+		if not numpy.any(numpy.isnan(self.cost_functions_coefficients)):
+			self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node')
+		if not numpy.any(numpy.isnan(self.min_parameters)):
+			self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node')
+		if not numpy.any(numpy.isnan(self.max_parameters)):
+			self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		
 		#default is incomplete adjoint for now
Index: ../trunk-jpl/src/m/classes/flowequation.py
===================================================================
--- ../trunk-jpl/src/m/classes/flowequation.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/flowequation.py	(revision 19048)
@@ -1,5 +1,6 @@
 import numpy
 import copy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
 from checkfield import checkfield
@@ -57,6 +58,14 @@
 		string="%s\n%s"%(string,fielddisplay(self,'borderFS',"vertices on FS' border (for tiling)"))
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.element_equation=project3d(md,'vector',self.element_equation,'type','element')
+		self.vertex_equation=project3d(md,'vector',self.vertex_equation,'type','node')
+		self.borderSSA=project3d(md,'vector',self.borderSSA,'type','node')
+		self.borderHO=project3d(md,'vector',self.borderHO,'type','node')
+		self.borderFS=project3d(md,'vector',self.borderFS,'type','node')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 
 		#P1 for SSA
Index: ../trunk-jpl/src/m/classes/basalforcings.py
===================================================================
--- ../trunk-jpl/src/m/classes/basalforcings.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/basalforcings.py	(revision 19048)
@@ -1,4 +1,5 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from checkfield import checkfield
 from WriteData import WriteData
@@ -29,6 +30,12 @@
 		string="%s\n%s"%(string,fielddisplay(self,"geothermalflux","geothermal heat flux [W/m^2]"))
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1)
+		self.floatingice_melting_rate=project3d(md,'vector',self.floatingice_melting_rate,'type','node','layer',1)
+		self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','node','layer',1)    #bedrock only gets geothermal flux
+		return self
+	#}}}
 	def initialize(self,md): # {{{
 
 		if numpy.all(numpy.isnan(self.groundedice_melting_rate)):
Index: ../trunk-jpl/src/m/classes/calvingpi.m
===================================================================
--- ../trunk-jpl/src/m/classes/calvingpi.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/calvingpi.m	(revision 19048)
@@ -28,6 +28,10 @@
 					error('constructor not supported');
 			end
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.coeff=project3d(md,'vector',self.coeff,'type','node');
+			self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
+		end % }}}
 		function self = setdefaultparameters(self) % {{{
 
 			%stabilization = 2 by default
Index: ../trunk-jpl/src/m/classes/flowequation.m
===================================================================
--- ../trunk-jpl/src/m/classes/flowequation.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/flowequation.m	(revision 19048)
@@ -106,6 +106,13 @@
 
 			fprintf(fid,'%s\n%s\n','</frame>');
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.element_equation=project3d(md,'vector',self.element_equation,'type','element');
+			self.vertex_equation=project3d(md,'vector',self.vertex_equation,'type','node');
+			self.borderSSA=project3d(md,'vector',self.borderSSA,'type','node');
+			self.borderHO=project3d(md,'vector',self.borderHO,'type','node');
+			self.borderFS=project3d(md,'vector',self.borderFS,'type','node');
+		end % }}}
 		function self = flowequation(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/basalforcings.m
===================================================================
--- ../trunk-jpl/src/m/classes/basalforcings.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/basalforcings.m	(revision 19048)
@@ -10,13 +10,18 @@
 		geothermalflux            = NaN;
 	end
 	methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- basalforcings -->');
-			 fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="melting_rate" type="',            class(self.melting_rate),'" default="',              num2str(self.melting_rate),'">',              '     <section name="basalforcings" />','     <help> basal melting rate (positive if melting) [m/yr] </help>','</parameter>');
-             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="geothermalflux" type="',          class(self.geothermalflux),'" default="',            num2str(self.geothermalflux),'">',            '     <section name="basalforcings" />','     <help> geothermal heat flux [W/m^2] </help>','</parameter>');
-             
-        end % }}}
+		function createxml(self,fid) % {{{
+			fprintf(fid, '\n\n');
+			fprintf(fid, '%s\n', '<!-- basalforcings -->');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="melting_rate" type="',            class(self.melting_rate),'" default="',              num2str(self.melting_rate),'">',              '     <section name="basalforcings" />','     <help> basal melting rate (positive if melting) [m/yr] </help>','</parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="geothermalflux" type="',          class(self.geothermalflux),'" default="',            num2str(self.geothermalflux),'">',            '     <section name="basalforcings" />','     <help> geothermal heat flux [W/m^2] </help>','</parameter>');
+
+		end % }}}
+		function self = extrude(self,md) % {{{
+			self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1); 
+			self.floatingice_melting_rate=project3d(md,'vector',self.floatingice_melting_rate,'type','node','layer',1); 
+			self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux
+		end % }}}
 		function self = basalforcings(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/damage.m
===================================================================
--- ../trunk-jpl/src/m/classes/damage.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/damage.m	(revision 19048)
@@ -84,7 +84,6 @@
 			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',        '<parameter key ="penalty_factor" type="',      class(self.penalty_factor),'" default="',            num2str(self.penalty_factor),'">',            '     <section name="damage" />','     <help> scaling exponent (default is 3) </help>','</parameter>');
 
 		end % }}}
-        
 		function self = damage(varargin) % {{{
 			switch nargin
 				case 0
@@ -103,6 +102,10 @@
 					error('constructor not supported');
 			end
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.D=project3d(md,'vector',self.D,'type','node');
+			self.spcdamage=project3d(md,'vector',self.spcdamage,'type','node');
+		end % }}}
 		function self = setdefaultparameters(self) % {{{
 
 			%damage parameters: 
Index: ../trunk-jpl/src/m/classes/taoinversion.m
===================================================================
--- ../trunk-jpl/src/m/classes/taoinversion.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/taoinversion.m	(revision 19048)
@@ -28,6 +28,15 @@
 		surface_obs               = NaN
 	end
 	methods
+		function self = extrude(self,md) % {{{
+			self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
+			self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
+			self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
+			self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
+			if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
+			if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
+			if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
+		end % }}}
 		function self = taoinversion(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/hydrologydc.m
===================================================================
--- ../trunk-jpl/src/m/classes/hydrologydc.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/hydrologydc.m	(revision 19048)
@@ -36,6 +36,15 @@
 		eplflip_lock             = 0;
   end
 	methods
+		function self = extrude(self,md) % {{{
+			self.spcsediment_head=project3d(md,'vector',self.spcsediment_head,'type','node','layer',1);
+			self.mask_eplactive_node=project3d(md,'vector',self.mask_eplactive_node,'type','node','layer',1);
+			self.sediment_transmitivity=project3d(md,'vector',self.sediment_transmitivity,'type','node','layer',1);
+			self.basal_moulin_input=project3d(md,'vector',self.basal_moulin_input,'type','node','layer',1);
+			if(self.isefficientlayer==1);
+				self.spcepl_head=project3d(md,'vector',self.spcepl_head,'type','node','layer',1);
+			end
+		end % }}}
 		% {{{ function self = hydrologydc(varargin) 
 		function self = hydrologydc(varargin) 
 			switch nargin
@@ -43,7 +52,7 @@
 					self=setdefaultparameters(self);
 				otherwise
 					error('constructor not supported');
-	    end
+			end
 		end 
 		% }}}
 		function self = initialize(self,md) % {{{
@@ -56,7 +65,7 @@
 		% {{{ function self = setdefaultparameters(self) 
 		function self = setdefaultparameters(self) 
 
-		%Parameters from de Fleurian 2014
+			%Parameters from de Fleurian 2014
 			self.water_compressibility    = 5.04e-10;
 			self.isefficientlayer         = 1;
 			self.penalty_factor           = 3;
@@ -85,7 +94,7 @@
 		% }}}
 		% {{{ function md = checkconsistency(self,md,solution,analyses) 
 		function md = checkconsistency(self,md,solution,analyses) 
-		%Early return
+			%Early return
 			if ~ismember(HydrologyDCInefficientAnalysisEnum(),analyses) & ~ismember(HydrologyDCEfficientAnalysisEnum(),analyses),
 				return;
 			end
@@ -100,10 +109,10 @@
 			md = checkfield(md,'fieldname','hydrology.transfer_flag','numel',[1],'values',[0 1]);
 			if self.sedimentlimit_flag==1,
 				md = checkfield(md,'fieldname','hydrology.sedimentlimit','>',0,'numel',1);
-	    end
+			end
 			if self.transfer_flag==1,
 				md = checkfield(md,'fieldname','hydrology.leakage_factor','>',0,'numel',1);
-	    end
+			end
 			md = checkfield(md,'fieldname','hydrology.basal_moulin_input','NaN',1,'timeseries',1);
 
 			md = checkfield(md,'fieldname','hydrology.spcsediment_head','timeseries',1);
@@ -126,7 +135,7 @@
 				if (self.epl_colapse_thickness>self.epl_initial_thickness),
 					md = checkmessage(md,'Colapsing thickness for EPL larger than initial thickness');
 				end 
-	    end
+			end
 		end 
 		% }}}
 		% {{{ function disp(self)
@@ -146,13 +155,13 @@
 			disp(sprintf('%55s  3: normal stress',' '));
 			if self.sedimentlimit_flag==1,
 				fielddisplay(self,'sedimentlimit','user defined upper limit for the inefficient layer [m]');
-	    end
+			end
 			fielddisplay(self,'transfer_flag',['what kind of transfer method is applied between the layers']);
 			disp(sprintf('%55s  0: no transfer',' '));
 			disp(sprintf('%55s  1: constant leakage factor: %s',' ','leakage_factor'));
 			if self.transfer_flag==1,
 				fielddisplay(self,'leakage_factor','user defined leakage factor [m]');
-	    end
+			end
 			fielddisplay(self,'basal_moulin_input','water flux at a given point [m3 s-1]');
 			disp(sprintf('   - for the sediment layer'));
 			fielddisplay(self,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]');
@@ -173,7 +182,7 @@
 				fielddisplay(self,'epl_max_thickness','epl maximal thickness [m]');
 				fielddisplay(self,'epl_conductivity','epl conductivity [m^2/s]');
 				fielddisplay(self,'eplflip_lock','lock the epl activation to avoid fli-floping (default is 0, no stabilization)');
-	    end
+			end
 
 		end 
 		% }}}
@@ -190,10 +199,10 @@
 			WriteData(fid,'object',self,'fieldname','transfer_flag','format','Integer');
 			if self.sedimentlimit_flag==1,
 				WriteData(fid,'object',self,'fieldname','sedimentlimit','format','Double');
-	    end
+			end
 			if self.transfer_flag==1,
 				WriteData(fid,'object',self,'fieldname','leakage_factor','format','Double');
-	    end
+			end
 			WriteData(fid,'object',self,'fieldname','basal_moulin_input','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1)
 
 			WriteData(fid,'object',self,'fieldname','spcsediment_head','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
@@ -213,8 +222,8 @@
 				WriteData(fid,'object',self,'fieldname','epl_max_thickness','format','Double');
 				WriteData(fid,'object',self,'fieldname','epl_conductivity','format','Double');
 				WriteData(fid,'object',self,'fieldname','eplflip_lock','format','Integer');
-	    end
+			end
 		end 
-% }}}
-  end
+		% }}}
+	end
 end
Index: ../trunk-jpl/src/m/classes/calvinglevermann.m
===================================================================
--- ../trunk-jpl/src/m/classes/calvinglevermann.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/calvinglevermann.m	(revision 19048)
@@ -28,6 +28,10 @@
 					error('constructor not supported');
 			end
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.coeff=project3d(md,'vector',self.coeff,'type','node');
+			self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
+		end % }}}
 		function self = setdefaultparameters(self) % {{{
 
 			%stabilization = 2 by default
Index: ../trunk-jpl/src/m/classes/mask.m
===================================================================
--- ../trunk-jpl/src/m/classes/mask.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/mask.m	(revision 19048)
@@ -21,12 +21,16 @@
 		end% }}}
 	end
 	methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- mask -->');
+		function createxml(self,fid) % {{{
+			fprintf(fid, '\n\n');
+			fprintf(fid, '%s\n', '<!-- mask -->');
 			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="groundedice_levelset" type="',class(self.groundedice_levelset),'" default="',self.groundedice_levelset,'">','     <section name="mask" />','     <help> is ice grounded ? grounded ice if > 0, grounding line position if = 0, floating ice if &amp;lt; 0 </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="ice_levelset" type="',class(self.ice_levelset),'" default="',self.ice_levelset,'">','     <section name="mask" />','     <help> presence of ice if > 0, icefront position if = 0, no ice if &amp;lt; 0 </help>','</parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="ice_levelset" type="',class(self.ice_levelset),'" default="',self.ice_levelset,'">','     <section name="mask" />','     <help> presence of ice if > 0, icefront position if = 0, no ice if &amp;lt; 0 </help>','</parameter>');
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.groundedice_levelset=project3d(md,'vector',self.groundedice_levelset,'type','node');
+			self.ice_levelset=project3d(md,'vector',self.ice_levelset,'type','node');
+		end % }}}
 		function self = mask(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/matdamageice.m
===================================================================
--- ../trunk-jpl/src/m/classes/matdamageice.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/matdamageice.m	(revision 19048)
@@ -63,6 +63,10 @@
 
 
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node');
+			self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element');
+		end % }}}
 		function self = matdamageice(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/damage.py
===================================================================
--- ../trunk-jpl/src/m/classes/damage.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/damage.py	(revision 19048)
@@ -1,4 +1,5 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from StringToEnum import StringToEnum
 from checkfield import checkfield
@@ -73,6 +74,11 @@
 
 		return s
 	# }}}
+	def extrude(self,md): # {{{
+		self.D=project3d(md,'vector',self.D,'type','node')
+		self.spcdamage=project3d(md,'vector',self.spcdamage,'type','node')
+		return self
+	#}}}
 	def setdefaultparameters(self):    # {{{
 
 		#damage parameters: 
Index: ../trunk-jpl/src/m/classes/inversion.m
===================================================================
--- ../trunk-jpl/src/m/classes/inversion.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/inversion.m	(revision 19048)
@@ -84,6 +84,15 @@
 					error('constructor not supported');
 			end
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
+			self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
+			self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
+			self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
+			if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
+			if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
+			if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
+		end % }}}
 		function self = setdefaultparameters(self) % {{{
 
 			%default is incomplete adjoint for now
Index: ../trunk-jpl/src/m/classes/geometry.m
===================================================================
--- ../trunk-jpl/src/m/classes/geometry.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/geometry.m	(revision 19048)
@@ -32,18 +32,24 @@
 		end% }}}
 	end
 	methods
-        %function createxml(self,fid)
-        function createxml(self, fid)% {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- geometry -->');
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Geometry parameters">','<section name="geometry" />'); 
+		function createxml(self, fid)% {{{
+			fprintf(fid, '\n\n');
+			fprintf(fid, '%s\n', '<!-- geometry -->');
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Geometry parameters">','<section name="geometry" />'); 
 			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="surface" type="','path','" optional="','false','">','     <section name="geometry" />','     <help> surface elevation [m] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="thickness" type="','path','" optional="','false','">','     <section name="geometry" />','     <help> ice thickness [m] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="bed" type="','path','" optional="','false','">','     <section name="geometry" />','     <help> bed elevation [m] </help>','</parameter>');
-            %fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="bathymetry" type="',class(self.bathymetry),'" default="',convert2str(self.bathymetry),'">','     <section name="geometry" />','     <help> bathymetry elevation [m] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="hydrostatic_ratio" type="',class(self.hydrostatic_ratio),'" default="',convert2str(self.hydrostatic_ratio),'">','     <section name="geometry" />','     <help> coefficient for ice shelves; thickness correction: hydrostatic_ratio H_obs+ (1-hydrostatic_ratio) H_hydro </help>','</parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');
-        end % }}}
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="thickness" type="','path','" optional="','false','">','     <section name="geometry" />','     <help> ice thickness [m] </help>','</parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="bed" type="','path','" optional="','false','">','     <section name="geometry" />','     <help> bed elevation [m] </help>','</parameter>');
+			%fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="bathymetry" type="',class(self.bathymetry),'" default="',convert2str(self.bathymetry),'">','     <section name="geometry" />','     <help> bathymetry elevation [m] </help>','</parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="hydrostatic_ratio" type="',class(self.hydrostatic_ratio),'" default="',convert2str(self.hydrostatic_ratio),'">','     <section name="geometry" />','     <help> coefficient for ice shelves; thickness correction: hydrostatic_ratio H_obs+ (1-hydrostatic_ratio) H_hydro </help>','</parameter>');
+			fprintf(fid,'%s\n%s\n','</frame>');
+		end % }}}
+		function self = extrude(self,md) % {{{
+			self.surface=project3d(md,'vector',self.surface,'type','node');
+			self.thickness=project3d(md,'vector',self.thickness,'type','node');
+			self.hydrostatic_ratio=project3d(md,'vector',self.hydrostatic_ratio,'type','node');
+			self.base=project3d(md,'vector',self.base,'type','node');
+			self.bed=project3d(md,'vector',self.bed,'type','node');
+		end % }}}
 		function self = geometry(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/gia.m
===================================================================
--- ../trunk-jpl/src/m/classes/gia.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/gia.m	(revision 19048)
@@ -10,21 +10,25 @@
 		cross_section_shape           = 0;
 	end
 	methods
-         function createxml(self,fid) % {{{
-            fprintf(fid, '<!-- gia -->\n');            
-                    
-            % gia solution parameters
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="gia parameters">','<section name="gia" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mantle_viscosity" type="',class(self.mantle_viscosity),'" default="',convert2str(self.mantle_viscosity),'">','     <section name="gia" />','     <help> mantle viscosity[Pa s] </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="lithosphere_thickness" type="',class(self.lithosphere_thickness),'" default="',convert2str(self.lithosphere_thickness),'">','     <section name="gia" />','     <help> lithosphere thickness (km) </help>','  </parameter>');
-            %cross_section_shape drop-down (1 o r 2)
-            fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="cross_section_shape" type="alternative" optional="false">','     <section name="gia" />','     <help> 1: square-edged (default). 2: elliptical.  See iedge in GiaDeflectionCore </help>');
-            fprintf(fid,'%s\n','       <option value="1" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="2" type="string" default="false"> </option>','</parameter>');
-            
-            fprintf(fid,'%s\n%s\n','</frame>');    
-        
-        end % }}}
+		function createxml(self,fid) % {{{
+			fprintf(fid, '<!-- gia -->\n');            
+
+			% gia solution parameters
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="gia parameters">','<section name="gia" />');                    
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="mantle_viscosity" type="',class(self.mantle_viscosity),'" default="',convert2str(self.mantle_viscosity),'">','     <section name="gia" />','     <help> mantle viscosity[Pa s] </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="lithosphere_thickness" type="',class(self.lithosphere_thickness),'" default="',convert2str(self.lithosphere_thickness),'">','     <section name="gia" />','     <help> lithosphere thickness (km) </help>','  </parameter>');
+			%cross_section_shape drop-down (1 o r 2)
+			fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="cross_section_shape" type="alternative" optional="false">','     <section name="gia" />','     <help> 1: square-edged (default). 2: elliptical.  See iedge in GiaDeflectionCore </help>');
+			fprintf(fid,'%s\n','       <option value="1" type="string" default="true"> </option>');
+			fprintf(fid,'%s\n%s\n','       <option value="2" type="string" default="false"> </option>','</parameter>');
+
+			fprintf(fid,'%s\n%s\n','</frame>');    
+
+		end % }}}
+		function self = extrude(self,md) % {{{
+			self.mantle_viscosity=project3d(md,'vector',self.mantle_viscosity,'type','node');
+			self.lithosphere_thickness=project3d(md,'vector',self.lithosphere_thickness,'type','node');
+		end % }}}
 		function self = gia(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/friction.py
===================================================================
--- ../trunk-jpl/src/m/classes/friction.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/friction.py	(revision 19048)
@@ -1,4 +1,5 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from checkfield import checkfield
 from WriteData import WriteData
@@ -28,6 +29,12 @@
 		string="%s\n%s"%(string,fielddisplay(self,"q","q exponent"))
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1)
+		self.p=project3d(md,'vector',self.p,'type','element')
+		self.q=project3d(md,'vector',self.q,'type','element')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		return self
 	#}}}
Index: ../trunk-jpl/src/m/classes/friction.m
===================================================================
--- ../trunk-jpl/src/m/classes/friction.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/friction.m	(revision 19048)
@@ -10,15 +10,20 @@
 		q           = NaN;
 	end
 	methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p -->');
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p">','<section name="friction" />');   
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="coefficient" type="',   	class(self.coefficient),'" default="',     	convert2str(self.coefficient),'">',              '     <section name="friction" />','     <help> friction coefficient [SI] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="p" type="',               class(self.p),'" default="',                 convert2str(self.p),'">',   '     <section name="friction" />','     <help> p exponent </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="q" type="',               class(self.q),'" default="',                 convert2str(self.q),'">',            '     <section name="friction" />','     <help> q exponent </help>','</parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');
-        end % }}}
+		function createxml(self,fid) % {{{
+			fprintf(fid, '\n\n');
+			fprintf(fid, '%s\n', '<!-- Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p -->');
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p">','<section name="friction" />');   
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="coefficient" type="',   	class(self.coefficient),'" default="',     	convert2str(self.coefficient),'">',              '     <section name="friction" />','     <help> friction coefficient [SI] </help>','</parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="p" type="',               class(self.p),'" default="',                 convert2str(self.p),'">',   '     <section name="friction" />','     <help> p exponent </help>','</parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="q" type="',               class(self.q),'" default="',                 convert2str(self.q),'">',            '     <section name="friction" />','     <help> q exponent </help>','</parameter>');
+			fprintf(fid,'%s\n%s\n','</frame>');
+		end % }}}
+		function self = extrude(self,md) % {{{
+			self.coefficient=project3d(md,'vector',self.coefficient,'type','node','layer',1);
+			self.p=project3d(md,'vector',self.p,'type','element');
+			self.q=project3d(md,'vector',self.q,'type','element');
+		end % }}}
 		function self = friction(varargin) % {{{
 			switch nargin
 				case 0
Index: ../trunk-jpl/src/m/classes/frictionweertman.m
===================================================================
--- ../trunk-jpl/src/m/classes/frictionweertman.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/frictionweertman.m	(revision 19048)
@@ -9,15 +9,15 @@
 		m = NaN;
 	end
 	methods
-        function createxml(self,fid) % {{{
-            fprintf(fid, '\n\n');
-            fprintf(fid, '%s\n', '<!-- Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p -->');
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p">','<section name="frictionweertman" />');   
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="coefficient" type="',   	class(self.coefficient),'" default="',     	convert2str(self.coefficient),'">',              '     <section name="frictionweertman" />','     <help> frictionweertman coefficient [SI] </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="p" type="',               class(self.p),'" default="',                 convert2str(self.p),'">',   '     <section name="frictionweertman" />','     <help> p exponent </help>','</parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="q" type="',               class(self.q),'" default="',                 convert2str(self.q),'">',            '     <section name="frictionweertman" />','     <help> q exponent </help>','</parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');
-        end % }}}
+		function createxml(self,fid) % {{{
+			fprintf(fid, '\n\n');
+			fprintf(fid, '%s\n', '<!-- Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p -->');
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Friction: Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p">','<section name="frictionweertman" />');   
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="coefficient" type="',   	class(self.coefficient),'" default="',     	convert2str(self.coefficient),'">',              '     <section name="frictionweertman" />','     <help> frictionweertman coefficient [SI] </help>','</parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n',    '<parameter key ="p" type="',               class(self.p),'" default="',                 convert2str(self.p),'">',   '     <section name="frictionweertman" />','     <help> p exponent </help>','</parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n',        '<parameter key ="q" type="',               class(self.q),'" default="',                 convert2str(self.q),'">',            '     <section name="frictionweertman" />','     <help> q exponent </help>','</parameter>');
+			fprintf(fid,'%s\n%s\n','</frame>');
+		end % }}}
 		function self = frictionweertman(varargin) % {{{
 			switch nargin
 				case 0
@@ -26,6 +26,10 @@
 					error('constructor not supported');
 			end
 		end % }}}
+		function self = extrude(self,md) % {{{
+			md.friction.C=project3d(md,'vector',md.friction.C,'type','node','layer',1);
+			md.friction.m=project3d(md,'vector',md.friction.m,'type','element');
+		end % }}}
 		function self = setdefaultparameters(self) % {{{
 
 		end % }}}
Index: ../trunk-jpl/src/m/classes/hydrologydc.py
===================================================================
--- ../trunk-jpl/src/m/classes/hydrologydc.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/hydrologydc.py	(revision 19048)
@@ -1,4 +1,5 @@
 import numpy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
 from checkfield import checkfield
@@ -46,7 +47,6 @@
 		#set defaults
 		self.setdefaultparameters()
 	#}}}
-
 	def __repr__(self): # {{{
 		string='   hydrology Dual Porous Continuum Equivalent parameters:'
 		string='   - general parameters'
@@ -94,6 +94,16 @@
 			string="%s\n%s"%(string,fielddisplay(self,'eplflip_lock','lock epl activity to avoid flip-floping (default is 0, no stabilization)'))
 		return string
 #}}}
+	def extrude(self,md): # {{{
+		self.spcsediment_head=project3d(md,'vector',self.spcsediment_head,'type','node','layer',1)
+		self.spcepl_head=project3d(md,'vector',self.spcepl_head,'type','node','layer',1)
+		self.mask_eplactive_node=project3d(md,'vector',self.mask_eplactive_node,'type','node','layer',1)
+		self.sediment_transmitivity=project3d(md,'vector',self.sediment_transmitivity,'type','node','layer',1)
+		self.basal_moulin_input=project3d(md,'vector',self.basal_moulin_input,'type','node','layer',1)
+		if self.isefficientlayer==1 :
+			self.spcepl_head=project3d(md,'vector',self.spcepl_head,'type','node','layer',1)
+		return self
+	#}}}
 	def setdefaultparameters(self): #{{{ 
 
 		#Parameters from de Fleurian 2014
@@ -124,7 +134,6 @@
 		
 		return self
 	# }}}
-
 	def initialize(self,md): # {{{
 		if numpy.all(numpy.isnan(self.basal_moulin_input)):
 			self.basal_moulin_input=numpy.zeros((md.mesh.numberofvertices,1))
Index: ../trunk-jpl/src/m/classes/frictionhydro.m
===================================================================
--- ../trunk-jpl/src/m/classes/frictionhydro.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/frictionhydro.m	(revision 19048)
@@ -30,6 +30,12 @@
 			md = checkfield(md,'fieldname','friction.C','NaN',1,'size',[md.mesh.numberofelements 1]);
 			md = checkfield(md,'fieldname','friction.As','NaN',1,'size',[md.mesh.numberofelements 1]);
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.q=project3d(md,'vector',self.q,'type','element');
+			self.C=project3d(md,'vector',self.C,'type','element');
+			self.As=project3d(md,'vector',self.As,'type','element');
+			self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1);
+		end % }}}
 		function disp(self) % {{{
 			disp(sprintf('Effective Pressure based friction law described in Gagliardini 2007'));
 			fielddisplay(self,'q','friction law exponent q>=1');
Index: ../trunk-jpl/src/m/classes/model.py
===================================================================
--- ../trunk-jpl/src/m/classes/model.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/model.py	(revision 19048)
@@ -594,103 +594,30 @@
 
 		#Ok, now deal with the other fields from the 2d mesh:
 
+		#bedinfo and surface info
+		md.mesh.vertexonbase=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
+		md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
+		md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node')
+
 		#lat long
 		md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node')
 		md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node')
 
-		#drag coefficient is limited to nodes that are on the bedrock.
-		md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1)
-
-		#p and q (same deal, except for element that are on the bedrock: )
-		md.friction.p=project3d(md,'vector',md.friction.p,'type','element')
-		md.friction.q=project3d(md,'vector',md.friction.q,'type','element')
-
-		#observations
-		md.inversion.vx_obs=project3d(md,'vector',md.inversion.vx_obs,'type','node')
-		md.inversion.vy_obs=project3d(md,'vector',md.inversion.vy_obs,'type','node')
-		md.inversion.vel_obs=project3d(md,'vector',md.inversion.vel_obs,'type','node')
-		md.inversion.thickness_obs=project3d(md,'vector',md.inversion.thickness_obs,'type','node')
+		md.geometry.extrude(md)
+		md.friction.extrude(md)
+		md.inversion.extrude(md)
 		md.surfaceforcings.extrude(md)
-		md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node')
+		md.initialization.extrude(md)
+		md.flowequation.extrude(md)
 
-		#results
-		if not numpy.any(numpy.isnan(md.initialization.vx)):
-			md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node')
-		if not numpy.any(numpy.isnan(md.initialization.vy)):
-			md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node')
-		if not numpy.any(numpy.isnan(md.initialization.vz)):
-			md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node')
-		if not numpy.any(numpy.isnan(md.initialization.vel)):
-			md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node')
-		if not numpy.any(numpy.isnan(md.initialization.temperature)):
-			md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node')
-		if not numpy.any(numpy.isnan(md.initialization.waterfraction)):
-			md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node')
-		if not numpy.any(numpy.isnan(md.initialization.watercolumn)):
-			md.initialization.watercolumn=project3d(md,'vector',md.initialization.watercolumn,'type','node')
-                if not numpy.any(numpy.isnan(md.initialization.sediment_head)):
-                        md.initialization.sediment_head=project3d(md,'vector',md.initialization.sediment_head,'type','node','layer',1)
-                if not numpy.any(numpy.isnan(md.initialization.epl_head)):
-                        md.initialization.epl_head=project3d(md,'vector',md.initialization.epl_head,'type','node','layer',1)
-                if not numpy.any(numpy.isnan(md.initialization.epl_thickness)):
-                        md.initialization.epl_thickness=project3d(md,'vector',md.initialization.epl_thickness,'type','node','layer',1)
+		md.stressbalance.extrude(md)
+		md.thermal.extrude(md)
+		md.masstransport.extrude(md)
 
-		#bedinfo and surface info
-		md.mesh.vertexonbase=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
-		md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
-
-		#elementstype
-		if not numpy.any(numpy.isnan(md.flowequation.element_equation)):
-			oldelements_type=md.flowequation.element_equation
-			md.flowequation.element_equation=numpy.zeros(number_el3d,int)
-			md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element')
-
-		#verticestype
-		if not numpy.any(numpy.isnan(md.flowequation.vertex_equation)):
-			oldvertices_type=md.flowequation.vertex_equation
-			md.flowequation.vertex_equation=numpy.zeros(number_nodes3d,int)
-			md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node')
-
-		md.flowequation.borderSSA=project3d(md,'vector',md.flowequation.borderSSA,'type','node')
-		md.flowequation.borderHO=project3d(md,'vector',md.flowequation.borderHO,'type','node')
-		md.flowequation.borderFS=project3d(md,'vector',md.flowequation.borderFS,'type','node')
-
-		#boundary conditions
-		md.stressbalance.spcvx=project3d(md,'vector',md.stressbalance.spcvx,'type','node')
-		md.stressbalance.spcvy=project3d(md,'vector',md.stressbalance.spcvy,'type','node')
-		md.stressbalance.spcvz=project3d(md,'vector',md.stressbalance.spcvz,'type','node')
-		md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',numpy.nan)
-		if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
-			md.thermal.spctemperature=numpy.nan*numpy.ones((md.mesh.numberofvertices,1))
-			if hasattr(md.mesh,'vertexonsurface'):
-				pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
-				md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
-		md.masstransport.spcthickness=project3d(md,'vector',md.masstransport.spcthickness,'type','node')
-		md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node')
-		md.damage.spcdamage=project3d(md,'vector',md.damage.spcdamage,'type','node')
-		md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node')
-		md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node')
-
 		# Calving variables
-		if isinstance(md.calving,calving):
-			md.calving.calvingrate=project3d(md,'vector',md.calving.calvingrate,'type','node')
-		if isinstance(md.calving,calvinglevermann):
-			md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node')
-		#if isinstance(md.calving,calvingpi):
-		#	md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node')
-		if not numpy.any(numpy.isnan(md.calving.meltingrate)):
-			md.calving.meltingrate=project3d(md,'vector',md.calving.meltingrate,'type','node')
+		md.hydrology.extrude(md)
+		md.calving.extrude(md)
 
-		# Hydrologydc variables
-		if hasattr(md.hydrology,'hydrologydc'):
-			md.hydrology.spcsediment_head=project3d(md,'vector',md.hydrology.spcsediment_head,'type','node','layer',1)
-			md.hydrology.spcepl_head=project3d(md,'vector',md.hydrology.spcepl_head,'type','node','layer',1)
-			md.hydrology.mask_eplactive_node=project3d(md,'vector',md.hydrology.mask_eplactive_node,'type','node','layer',1)
-			md.hydrology.sediment_transmitivity=project3d(md,'vector',md.hydrology.sediment_transmitivity,'type','node','layer',1)
-			md.hydrology.basal_moulin_input=project3d(md,'vector',md.hydrology.basal_moulin_input,'type','node','layer',1)
-			if md.hydrology.isefficientlayer==1 :
-				md.hydrology.spcepl_head=project3d(md,'vector',md.hydrology.spcepl_head,'type','node','layer',1)
-
 		#connectivity
 		md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1))
 		md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1
@@ -700,43 +627,13 @@
 						=md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d
 				md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity<0)]=0
 
-		#materials
-		md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node')
-		md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element')
+		md.materials.extrude(md)
+		md.damage.extrude(md)
+		md.gia.extrude(md)
+		md.mask.extrude(md)
+		md.qmu.extrude(md)
+		md.basalforcings.extrude(md)
 
-		#damage
-		md.damage.D=project3d(md,'vector',md.damage.D,'type','node')
-
-		#parameters
-		md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node')
-		md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node')
-		md.gia.mantle_viscosity=project3d(md,'vector',md.gia.mantle_viscosity,'type','node')
-		md.gia.lithosphere_thickness=project3d(md,'vector',md.gia.lithosphere_thickness,'type','node')
-		md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node')
-		md.geometry.base=project3d(md,'vector',md.geometry.base,'type','node')
-		md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node')
-		md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node')
-		md.mask.ice_levelset=project3d(md,'vector',md.mask.ice_levelset,'type','node')
-		md.mask.groundedice_levelset=project3d(md,'vector',md.mask.groundedice_levelset,'type','node')
-		if not numpy.any(numpy.isnan(md.inversion.cost_functions_coefficients)):
-			md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node')
-		if not numpy.any(numpy.isnan(md.inversion.min_parameters)):
-			md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node')
-		if not numpy.any(numpy.isnan(md.inversion.max_parameters)):
-			md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node')
-		if not numpy.any(numpy.isnan(md.qmu.partition)):
-			md.qmu.partition=project3d(md,'vector',numpy.transpose(md.qmu.partition),'type','node')
-
-		#Put lithostatic pressure if there is an existing pressure
-		if not numpy.any(numpy.isnan(md.initialization.pressure)):
-			md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1))
-
-		#special for thermal modeling:
-		md.basalforcings.groundedice_melting_rate=project3d(md,'vector',md.basalforcings.groundedice_melting_rate,'type','node','layer',1)
-		md.basalforcings.floatingice_melting_rate=project3d(md,'vector',md.basalforcings.floatingice_melting_rate,'type','node','layer',1)
-		if not numpy.any(numpy.isnan(md.basalforcings.geothermalflux)):
-			md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1)    #bedrock only gets geothermal flux
-
 		#increase connectivity if less than 25:
 		if md.mesh.average_vertex_connectivity<=25:
 			md.mesh.average_vertex_connectivity=100
Index: ../trunk-jpl/src/m/classes/stressbalance.py
===================================================================
--- ../trunk-jpl/src/m/classes/stressbalance.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/stressbalance.py	(revision 19048)
@@ -1,6 +1,7 @@
 import numpy
 import sys
 import copy
+from project3d import project3d
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
 from checkfield import checkfield
@@ -74,6 +75,15 @@
 
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.spcvx=project3d(md,'vector',self.spcvx,'type','node')
+		self.spcvy=project3d(md,'vector',self.spcvy,'type','node')
+		self.spcvz=project3d(md,'vector',self.spcvz,'type','node')
+		self.referential=project3d(md,'vector',self.referential,'type','node')
+		self.loadingforce=project3d(md,'vector',self.loadingforce,'type','node')
+
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		#maximum of non-linear iterations.
 		self.maxiter=100
Index: ../trunk-jpl/src/m/classes/qmu.py
===================================================================
--- ../trunk-jpl/src/m/classes/qmu.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/qmu.py	(revision 19048)
@@ -1,4 +1,5 @@
 import numpy
+from project3d import project3d
 from collections import OrderedDict
 from fielddisplay import fielddisplay
 from EnumDefinitions import *
@@ -36,38 +37,13 @@
 		self.setdefaultparameters()
 
 		#}}}
-	def setdefaultparameters(self): # {{{
-		return self
-	#}}}
-	def checkconsistency(self,md,solution,analyses):    # {{{
-
-		#Early return
-		if not md.qmu.isdakota:
-			return
-
-		if not md.qmu.params.evaluation_concurrency==1:
-			md.checkmessage("concurrency should be set to 1 when running dakota in library mode")
-		if md.qmu.partition:
-			if not numpy.size(md.qmu.partition)==md.mesh.numberofvertices:
-				md.checkmessage("user supplied partition for qmu analysis should have size md.mesh.numberofvertices x 1")
-			if not min(md.qmu.partition)==0:
-				md.checkmessage("partition vector not indexed from 0 on")
-			if max(md.qmu.partition)>=md.qmu.numberofpartitions:
-				md.checkmessage("for qmu analysis, partitioning vector cannot go over npart, number of partition areas")
-
-		if not m.strcmpi(md.cluster.name,'none'):
-			if not md.settings.waitonlock:
-				md.checkmessage("waitonlock should be activated when running qmu in parallel mode!")
-
-		return md
-	# }}}
 	def __repr__(self):    # {{{
 		s ='   qmu parameters:\n'
 
 		s+="%s\n" % fielddisplay(self,'isdakota','is qmu analysis activated?')
 		for i,variable in enumerate(self.variables.iteritems()):
 			s+="         variables%s:  (arrays of each variable class)\n" % \
-				string_dim(self.variables,i)
+					string_dim(self.variables,i)
 			fnames=vars(variable)
 			maxlen=0
 			for fname in fnames:
@@ -75,11 +51,11 @@
 
 			for fname in fnames:
 				s+="'            %-*s:    [%ix%i]    '%s'\n" % \
-					(maxlen+1,fname,size(getattr(variable,fname)),type(getattr(variable,fname)))
+						(maxlen+1,fname,size(getattr(variable,fname)),type(getattr(variable,fname)))
 
 		for i,response in enumerate(self.responses.iteritems()):
 			s+="         responses%s:  (arrays of each response class)\n" % \
-				string_dim(self.responses,i)
+					string_dim(self.responses,i)
 			fnames=vars(response)
 			maxlen=0
 			for fname in fnames:
@@ -87,18 +63,18 @@
 
 			for fname in fnames:
 				s+="            %-*s:    [%ix%i]    '%s'\n" % \
-					(maxlen+1,fname,size(getattr(response,fname)),type(getattr(response,fname)))
+						(maxlen+1,fname,size(getattr(response,fname)),type(getattr(response,fname)))
 
 		s+="%s\n" % fielddisplay(self,'numberofresponses','number of responses') 
 
 		for i,method in enumerate(self.method.iteritems()):
 			if isinstance(method,'dakota_method'):
 				s+="            method%s :    '%s'\n" % \
-					(string_dim(method,i),method.method)
+						(string_dim(method,i),method.method)
 
 		for i,param in enumerate(self.params.iteritems()):
 			s+="         params%s:  (array of method-independent parameters)\n" % \
-				string_dim(self.params,i)
+					string_dim(self.params,i)
 			fnames=vars(param)
 			maxlen=0
 			for fname in fnames:
@@ -106,11 +82,11 @@
 
 			for fname in fnames:
 				s+="            %-*s: %s\n" % \
-					(maxlen+1,fname,any2str(getattr(param,fname)))
+						(maxlen+1,fname,any2str(getattr(param,fname)))
 
 		for i,result in enumerate(self.results.iteritems()):
 			s+="         results%s:  (information from dakota files)\n" % \
-				string_dim(self.results,i)
+					string_dim(self.results,i)
 			fnames=vars(result)
 			maxlen=0
 			for fname in fnames:
@@ -118,7 +94,7 @@
 
 			for fname in fnames:
 				s+="            %-*s:    [%ix%i]    '%s'\n" % \
-					(maxlen+1,fname,size(getattr(result,fname)),type(getattr(result,fname)))
+						(maxlen+1,fname,size(getattr(result,fname)),type(getattr(result,fname)))
 
 		s+="%s\n" % fielddisplay(self,'partition','user provided mesh partitioning, defaults to metis if not specified') 
 		s+="%s\n" % fielddisplay(self,'numberofpartitions','number of partitions for semi-discrete qmu') 
@@ -133,6 +109,35 @@
 
 		return s
 	# }}}
+	def extrude(self,md): # {{{
+		self.partition=project3d(md,'vector',numpy.transpose(self.partition),'type','node')
+		return self
+	#}}}
+	def setdefaultparameters(self): # {{{
+		return self
+	#}}}
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		#Early return
+		if not md.qmu.isdakota:
+			return
+
+		if not md.qmu.params.evaluation_concurrency==1:
+			md.checkmessage("concurrency should be set to 1 when running dakota in library mode")
+		if md.qmu.partition:
+			if not numpy.size(md.qmu.partition)==md.mesh.numberofvertices:
+				md.checkmessage("user supplied partition for qmu analysis should have size md.mesh.numberofvertices x 1")
+			if not min(md.qmu.partition)==0:
+				md.checkmessage("partition vector not indexed from 0 on")
+			if max(md.qmu.partition)>=md.qmu.numberofpartitions:
+				md.checkmessage("for qmu analysis, partitioning vector cannot go over npart, number of partition areas")
+
+		if not m.strcmpi(md.cluster.name,'none'):
+			if not md.settings.waitonlock:
+				md.checkmessage("waitonlock should be activated when running qmu in parallel mode!")
+
+		return md
+	# }}}
 	def marshall(self,md,fid):    # {{{
 		WriteData(fid,'object',self,'fieldname','isdakota','format','Boolean')
 		if not self.isdakota:
Index: ../trunk-jpl/src/m/classes/matice.py
===================================================================
--- ../trunk-jpl/src/m/classes/matice.py	(revision 19047)
+++ ../trunk-jpl/src/m/classes/matice.py	(revision 19048)
@@ -1,4 +1,5 @@
 from fielddisplay import fielddisplay
+from project3d import project3d
 from EnumDefinitions import *
 from StringToEnum import StringToEnum
 from checkfield import checkfield
@@ -62,6 +63,11 @@
 
 		return string
 		#}}}
+	def extrude(self,md): # {{{
+		self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node')
+		self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element')
+		return self
+	#}}}
 	def setdefaultparameters(self): # {{{
 		#ice density (kg/m^3)
 		self.rho_ice=917.
Index: ../trunk-jpl/src/m/classes/masstransport.m
===================================================================
--- ../trunk-jpl/src/m/classes/masstransport.m	(revision 19047)
+++ ../trunk-jpl/src/m/classes/masstransport.m	(revision 19048)
@@ -58,6 +58,9 @@
 					error('constructor not supported');
 			end
 		end % }}}
+		function self = extrude(self,md) % {{{
+			self.spcthickness=project3d(md,'vector',self.spcthickness,'type','node');
+		end % }}}
 		function list = defaultoutputs(self,md) % {{{
 
 			list = {'Thickness','Surface','Base'};
