Index: /issm/trunk-jpl/src/m/interp/SectionValues.m
===================================================================
--- /issm/trunk-jpl/src/m/interp/SectionValues.m	(revision 17732)
+++ /issm/trunk-jpl/src/m/interp/SectionValues.m	(revision 17733)
@@ -92,7 +92,7 @@
 	%vertically extrude mesh
 
-	%Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system
+	%Get base and surface for each 2d point, offset to make sure that it is inside the glacier system
 	offset=10^-3;
-	bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.base,1),X,Y)+offset;
+	base=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.base,1),X,Y)+offset;
 	surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),X,Y)-offset;
 
@@ -111,5 +111,5 @@
 		X3(i:layers:end)=X;
 		Y3(i:layers:end)=Y;
-		Z3(i:layers:end)=bed+(i-1)*(surface-bed)/(layers-1);
+		Z3(i:layers:end)=base+(i-1)*(surface-base)/(layers-1);
 		S3(i:layers:end)=S;
 
Index: /issm/trunk-jpl/src/m/interp/SectionValues.py
===================================================================
--- /issm/trunk-jpl/src/m/interp/SectionValues.py	(revision 17733)
+++ /issm/trunk-jpl/src/m/interp/SectionValues.py	(revision 17733)
@@ -0,0 +1,140 @@
+import os
+from expread import expread
+import numpy as npy
+from project2d import project2d
+#from InterpFromMesh2d import InterpFromMesh2d
+from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
+#from InterpFromMeshToMesh3d import InterpFromMeshToMesh3d
+
+def SectionValues(md,data,infile,resolution):
+	'''
+	compute the value of a field on a section
+	
+	This routine gets the value of a given field of the model on points
+	given in the file infile (Argus type file). Resolution must be a list
+	[horizontal_resolution, vertical_resolution]
+	
+	Usage:
+	[elements,x,y,z,s,data]=SectionValues(md,data,filename,resolution)
+	[elements,x,y,z,s,data]=SectionValues(md,data,profile_structure,resolution)
+	'''
+
+	if os.path.isfile(infile):
+		profile=expread(infile)[0]
+		nods=profile['nods']
+		x=profile['x']
+		y=profile['y']
+	else:
+		raise IOError('file %s not found' % infile)
+
+	#get the specified resolution
+	if len(resolution)!=2:
+		raise ValueError('SectionValues error message: Resolution must be a list [horizontal_resolution, vertical_resolution]')
+	else:
+		res_h=resolution[0]
+
+	if md.mesh.domaintype().lower() == '3d':
+		if isinstance(resolution[1],int) or isinstance(resolution[1],float):
+			res_v=resolution[1]
+		else:
+			raise ValueError('SectionValues error: resolution must be a length-2 list of integers or floats')
+
+	#initialization
+	X=npy.array([]) #X-coordinate
+	Y=npy.array([]) #Y-coordinate
+	S=npy.array([0.])  #curvilinear coordinate
+	
+	for i in xrange(nods-1):
+	
+		x_start=x[i]
+		x_end=x[i+1]
+		y_start=y[i]
+		y_end=y[i+1]
+		s_start=S[-1]
+	
+		length_segment=npy.sqrt((x_end-x_start)**2+(y_end-y_start)**2)
+		portion=npy.ceil(length_segment/res_h)
+	
+		x_segment=npy.zeros(portion)
+		y_segment=npy.zeros(portion)
+		s_segment=npy.zeros(portion)
+
+		for j in xrange(int(portion)):
+			x_segment[j]=x_start+(j-1)*(x_end-x_start)/portion
+			y_segment[j]=y_start+(j-1)*(y_end-y_start)/portion
+			s_segment[j]=s_start+j*length_segment/portion
+	
+		#plug into X and Y
+		X=npy.append(X,x_segment)
+		Y=npy.append(Y,y_segment)
+		S=npy.append(S,s_segment)
+
+	X=npy.append(X,x[nods-1])
+	Y=npy.append(Y,y[nods-1])
+	
+	#Number of nodes:
+	numberofnodes=X.shape[0]
+	
+	#Compute Z
+	Z=npy.zeros(numberofnodes)
+	
+	#New mesh and Data interpolation
+	if md.mesh.domaintype().lower() == '2d':
+	
+		#Interpolation of data on specified points
+		#data_interp=InterpFromMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y)[0]
+		data_interp=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y)[0]
+		#data_interp=griddata(md.mesh.x,md.mesh.y,data,X,Y)
+	
+		#Compute index
+		index=npy.array([range(1,numberofnodes),range(2,numberofnodes+1)]).T
+	
+	else:
+	
+		#vertically extrude mesh
+	
+		#Get base and surface for each 2d point, offset to make sure that it is inside the glacier system
+		offset=1.e-3
+		base=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.base,1),X,Y)[0]+offset
+		base=base.reshape(-1,)
+		surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),X,Y)[0]-offset
+		surface=surface.reshape(-1,)
+	
+		#Some useful parameters
+		layers=int(npy.ceil(npy.mean(md.geometry.thickness)/res_v))
+		print 'layers=', layers
+		nodesperlayer=int(numberofnodes)
+		nodestot=int(nodesperlayer*layers)
+		elementsperlayer=int(nodesperlayer-1)
+		elementstot=int((nodesperlayer-1)*(layers-1))
+	
+		#initialization
+		X3=npy.zeros(nodesperlayer*layers) 
+		Y3=npy.zeros(nodesperlayer*layers) 
+		Z3=npy.zeros(nodesperlayer*layers) 
+		S3=npy.zeros(nodesperlayer*layers) 
+		index3=npy.zeros((elementstot,4))
+	
+		#Get new coordinates in 3d
+		for i in xrange(1,layers+1):
+			X3[i-1::layers]=X
+			Y3[i-1::layers]=Y
+			Z3[i-1::layers]=base+(i-1)*(surface-base)/(layers-1)
+			S3[i-1::layers]=S
+	
+			if i<layers-1:  #Build index3 with quads
+				ids=npy.vstack((npy.arange(i,nodestot-layers,layers),npy.arange(i+1,nodestot-layers,layers),npy.arange(i+layers+1,nodestot,layers),npy.arange(i+layers,nodestot,layers))).T
+				index3[(i-1)*elementsperlayer:i*elementsperlayer,:]=ids
+
+		#Interpolation of data on specified points
+		# TODO data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,npy.nan)
+	
+		#build outputs
+		X=X3 
+		Y=Y3 
+		Z=Z3  
+		S=S3 
+
+		index=index3
+
+		return index,X,Y,Z,S,data_interp