Index: /issm/trunk-jpl/src/m/mesh/ComputeHessian.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/ComputeHessian.py	(revision 27146)
+++ /issm/trunk-jpl/src/m/mesh/ComputeHessian.py	(revision 27147)
@@ -39,9 +39,9 @@
 
     #compute weights that hold the volume of all the element holding the node i
-    weights = m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile(areas, (3, 1)), numberofnodes, 1)
+    weights = m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile(areas, (1, 3)), numberofnodes, 1)
 
     #compute field on nodes if on elements
     if np.size(field, axis=0) == numberofelements:
-        field = m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile(areas * field, (3, 1)), numberofnodes, 1) / weights
+        field = m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile(areas * field, (1, 3)), numberofnodes, 1) / weights
 
     #Compute gradient for each element
@@ -50,6 +50,6 @@
 
     #Compute gradient for each node (average of the elements around)
-    gradx = m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile((areas * grad_elx), (3, 1)), numberofnodes, 1)
-    grady = m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile((areas * grad_ely), (3, 1)), numberofnodes, 1)
+    gradx = m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile((areas * grad_elx), (1, 3)), numberofnodes, 1)
+    grady = m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile((areas * grad_ely), (1, 3)), numberofnodes, 1)
     gradx = gradx / weights
     grady = grady / weights
@@ -60,7 +60,7 @@
     if m.strcmpi(type, 'node'):
         #Compute Hessian on the nodes (average of the elements around)
-        hessian = np.hstack((m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile((areas * hessian[:, 0]), (3, 1)), numberofnodes, 1) / weights,
-                             m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile((areas * hessian[:, 1]), (3, 1)), numberofnodes, 1) / weights,
-                             m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile((areas * hessian[:, 2]), (3, 1)), numberofnodes, 1) / weights))
+        hessian = np.hstack((m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile((areas * hessian[:, 0]), (1, 3)), numberofnodes, 1) / weights,
+                             m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile((areas * hessian[:, 1]), (1, 3)), numberofnodes, 1) / weights,
+                             m.sparse(line, np.ones((linesize, 1), dtype=int), np.tile((areas * hessian[:, 2]), (1, 3)), numberofnodes, 1) / weights))
 
     return hessian