Index: /issm/trunk-jpl/src/m/partition/AreaAverageOntoPartition.py
===================================================================
--- /issm/trunk-jpl/src/m/partition/AreaAverageOntoPartition.py	(revision 23270)
+++ /issm/trunk-jpl/src/m/partition/AreaAverageOntoPartition.py	(revision 23270)
@@ -0,0 +1,59 @@
+import numpy as np
+
+def AreaAverageOntoPartition(md,vector,layer=None):
+	'''AREAAVERAGEONTOPARTITION 
+   compute partition values for a certain vector expressed on the vertices of the mesh.
+   Use area weighted average.
+
+   Usage:
+      average=AreaAverageOntoPartition(md,vector)
+      average=AreaAverageOntoPartition(md,vector,layer) #if in 3D, chose which layer is partitioned
+'''
+	#some checks
+	if len(np.shape(md.mesh)) == 3:
+		if layer == None:
+			raise RuntimeError('AreaAverageOntoPartition: layer should be provided onto which Area Averaging occurs')
+		
+		#save 3D model
+		md3d = md
+
+		md.mesh.elements = md.mesh.elements2d
+		md.mesh.x = md.mesh.x2d
+		md.mesh.y = md.mesh.y2d
+		md.mesh.numberofvertices = md.mesh.numberofvertices2d
+		md.mesh.numberofelements = md.mesh.numberofelements2d
+		md.qmu.vertex_weight = []
+		md.mesh.vertexconnectivity = []
+
+		#run connectivity routine
+		md = adjacency(md)
+
+		#finally, project vector: 
+		vector = project2d(md3d,vector,layer)
+		md.qmu.partition = project2d(md3d,md3d.qmu.partition,layer)
+
+
+	#ok, first check that part is Matlab indexed
+	#part = md.qmu.partition
+
+	#some check: 
+	if md.qmu.numberofpartitions != max(md.qmu.partition):
+		raise RuntimeError('AreaAverageOntoPartition error message: ''npart'' should be equal to max(md.qmu.partition)')
+
+
+	#initialize output
+	partvector = np.zeros((max(md.qmu.partition),))
+
+	#start weight average
+	weightedvector = vector*md.qmu.vertex_weight
+	for i in range(max(md.qmu.partition)):
+		pos=np.where(md.qmu.partition==i)
+		partvector[i]=sum(weightedvector[pos])/sum(md.qmu.vertex_weight[pos])
+
+
+	#in 3D, restore 3D model:
+	if len(np.shape(md.mesh)) == 3:
+		md = md3d
+
+	return partvector
+
