import numpy as npy

def project2d(md3d,value,layer):
	'''
	returns the value of a field for a given layer of the mesh
	
	project 'value' vector taken at layer 'layer' from extruded 2d-3d
	mesh onto 2d mesh used to do the extrusion. This routine is used to
	compare values between a 2d-3d mesh at a certain layer, and the
	equivalent value (if it exists), on the original 2d mesh.	

	Note: for consistency with usage in ISSM, layers are indexed starting
	from one.  This routine accounts for the fact that python indexes from
	zero.  In other words, when using this function consider the bottom
	layer as layer=1.  

   Usage:
      projection_value=project2d(md3d,value,layer)

   Example:
      vel2=project2d(md3d,md3d.vel,2)
	'''

	if md3d.mesh.domaintype().lower() != '3d':
		raise StandardError("model passed to project2d function should be 3D")

	if layer<1 or layer>md3d.mesh.numberoflayers:
		raise ValueError("layer must be between 0 and %i" % md3d.mesh.numberoflayers)
	
	# coerce to array in case float is passed
	if type(value)!=npy.ndarray:
		print 'coercing array'
		value=npy.array(value)

	vec2d=False
	if value.ndim==2 and value.shape[1]==1: 
		value=value.reshape(-1,)
		vec2d=True

	if value.size==1:
		projection_value=value[(layer-1)*md3d.mesh.numberofelements2d:layer*md3d.mesh.numberofelements2d]
	elif value.shape[0]==md3d.mesh.numberofvertices:
		#print 'indices: ', (layer-1)*md3d.mesh.numberofvertices2d, layer*md3d.mesh.numberofvertices2d
		projection_value=value[(layer-1)*md3d.mesh.numberofvertices2d:layer*md3d.mesh.numberofvertices2d]
	elif value.shape[0]==md3d.mesh.numberofvertices+1:
		projection_value=[value[(layer-1)*md3d.mesh.numberofvertices2d:layer*md3d.mesh.numberofvertices2d], value[-1]]
	else:
		projection_value=value[(layer-1)*md3d.mesh.numberofelements2d:layer*md3d.mesh.numberofelements2d]

	if vec2d:
		projection_value=projection_value.reshape(-1,1)

	return projection_value
