Index: /issm/trunk-jpl/src/m/interp/averaging.m
===================================================================
--- /issm/trunk-jpl/src/m/interp/averaging.m	(revision 26375)
+++ /issm/trunk-jpl/src/m/interp/averaging.m	(revision 26376)
@@ -73,5 +73,5 @@
 linesize=rep*numberofelements;
 
-%update weights that holds the volume of all the element holding the node i
+%update weights that hold the volume of all the element holding the node i
 weights=sparse(line,ones(linesize,1),repmat(areas,rep,1),numberofnodes,1);
 
Index: /issm/trunk-jpl/src/m/interp/averaging.py
===================================================================
--- /issm/trunk-jpl/src/m/interp/averaging.py	(revision 26375)
+++ /issm/trunk-jpl/src/m/interp/averaging.py	(revision 26376)
@@ -39,9 +39,9 @@
     # Initialization
     if layer == 0:
-        weights = np.zeros(md.mesh.numberofvertices, )
+        weights = np.zeros((md.mesh.numberofvertices, ))
         data = np.asarray(data).flatten()
     else:
-        weights = np.zeros(md.mesh.numberofvertices2d, )
-        data = data[(layer - 1) * md.mesh.numberofvertices2d + 1:layer * md.mesh.numberofvertices2d, :]
+        weights = np.zeros((md.mesh.numberofvertices2d, ))
+        data = data[(layer - 1) * md.mesh.numberofvertices2d + 1:layer * md.mesh.numberofvertices2d - 1, :]
 
     # Load some variables (it is much faster if the variables are loaded from md once for all)
@@ -56,5 +56,5 @@
 
     # Build some variables
-    if md.mesh.dimension() == 3 and layer == 0:
+    if (md.mesh.dimension() == 3) & (layer == 0):
         rep = 6
         areas = GetAreas(index, md.mesh.x, md.mesh.y, md.mesh.z)
@@ -67,16 +67,16 @@
 
     index = index - 1  # Python indexes from zero
-    line = index.flatten()
+    line=index.T.flatten()
     areas = np.vstack(areas).reshape(-1, )
-    summation = 1. / rep * np.ones(rep, )
+    summation = 1. / rep * np.ones((rep,1) )
     linesize = rep * numberofelements
 
-    # Update weights that holds the volume of all the element holding the node i
-    weights = csc_matrix((np.tile(areas, (rep, 1)).reshape(-1, ), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
+    # Update weights that hold the volume of all the element holding the node i
+    weights = csc_matrix((np.tile(areas, (1, rep)).reshape(-1,), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
 
     # Initialization
     if len(data) == numberofelements:
-        average_node = csc_matrix((np.tile(areas * data, (rep, 1)).reshape(-1, ), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
-        average_node = average_node / weights
+        average_node = csc_matrix((np.tile(np.multiply(areas,data), (1, rep)).reshape(-1, ), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
+        average_node = np.divide(average_node,weights)
         average_node = csc_matrix(average_node)
     else:
@@ -85,7 +85,7 @@
     # Loop over iteration
     for i in np.arange(1, iterations + 1):
-        average_el = np.asarray(np.dot(average_node.todense()[index].reshape(numberofelements, rep), np.vstack(summation))).reshape(-1, )
-        average_node = csc_matrix((np.tile(areas * average_el.reshape(-1), (rep, 1)).reshape(-1, ), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
-        average_node = average_node / weights
+        average_el = np.asarray(average_node.todense()[index].reshape(numberofelements, rep)*summation).reshape(-1, )
+        average_node = csc_matrix((np.tile(np.multiply(areas,average_el.reshape(-1)), (1, rep)).reshape(-1, ), (line, np.zeros(linesize, ))), shape=(numberofnodes, 1))
+        average_node = np.divide(average_node,weights)
         average_node = csc_matrix(average_node)