diff --git a/ntrfc/math/vectorcalc.py b/ntrfc/math/vectorcalc.py
index b123d1ed0f94b9d3071b55b1deaa1ca9bebe9c67..85d7b158b86135eed93684f23b841ca1b280c724 100644
--- a/ntrfc/math/vectorcalc.py
+++ b/ntrfc/math/vectorcalc.py
@@ -267,3 +267,11 @@ def line_intersection(point_a1, point_a2,
     x = det_2d(d, xdiff) / div
     y = det_2d(d, ydiff) / div
     return x, y
+
+
+def unitvec(vec):
+    return vec / vecAbs(vec)
+
+
+def unitvec_list(vecs):
+    return vecs / vecAbs_list(vecs)[:,None]
diff --git a/ntrfc/meshquality/nondimensionals.py b/ntrfc/meshquality/nondimensionals.py
index 949f91ecc5ff6f2d6753f1d587c7e1e21f439e88..a53f16ae52e5a5e06fe1a71aecb90f98d79a286e 100644
--- a/ntrfc/meshquality/nondimensionals.py
+++ b/ntrfc/meshquality/nondimensionals.py
@@ -8,11 +8,7 @@ Created on Fri Apr 24 23:33:53 2020
 import numpy as np
 import pyvista as pv
 
-from ntrfc.math.vectorcalc import vecAbs, vecProjection, vecAngle
-
-
-def unitvec(vec):
-    return vec / vecAbs(vec)
+from ntrfc.math.vectorcalc import vecAbs,vecAbs_list, vecProjection, vecAngle, unitvec, unitvec_list
 
 
 def readDataSet(grid, dataName):
@@ -120,32 +116,20 @@ def cellSpans( solutionMesh, calcFrom):
 def getWalluTaus(labelChunk, solutionMesh, surfaceMesh, mu_0, rhofieldname, velfieldname):
     uTaus = []
     surfaceNormals = surfaceMesh.extract_surface().compute_normals()
-    for cellIdx in labelChunk:
-        cellCenter = solutionMesh["cellCenters"][cellIdx]
-
-        pointOnWallId = surfaceNormals.find_closest_point(cellCenter)
-        pointOnWall = surfaceNormals.points[pointOnWallId]
-        # surfaceNormal = surfaceNormals.point_arrays["Normals"][pointOnWallId]#        pointOnWall = cellNormal+cellCenter
-
-        PT = pv.PolyData(pointOnWall)
-        PT = PT.sample(surfaceMesh)
-        PT.set_active_scalars("gradient")
-
-        gradUWall = PT.active_scalars[0]
+    cell_centers = solutionMesh["cellCenters"]
 
-        gradUWall = np.array([[gradUWall[0], gradUWall[1], gradUWall[2]],
-                              [gradUWall[3], gradUWall[4], gradUWall[5]],
-                              [gradUWall[6], gradUWall[7], gradUWall[8]]])
+    relatedpointonwallids = [surfaceNormals.find_closest_point(i) for i in cell_centers]
+    relatedpointonwall = [surfaceNormals.points[i] for i in relatedpointonwallids]
+    relatedpointonwall_samples = pv.PolyData(relatedpointonwall).sample(surfaceMesh)
+    graduwall = relatedpointonwall_samples["gradient"].reshape(relatedpointonwall_samples.number_of_points,3,3)
+    flowdirections = unitvec_list(solutionMesh[velfieldname])
+    gradUyWs = vecAbs_list([np.dot(flowdirections[i], graduwall[i]) for i in range(relatedpointonwall_samples.number_of_points)])
+    rhoWs = relatedpointonwall_samples[rhofieldname]
 
-        flowdirection = unitvec(solutionMesh[velfieldname][cellIdx])
-        gradUyW = vecAbs(np.dot(flowdirection, gradUWall))
-        PT.set_active_scalars(rhofieldname)
-        rhoW = PT.active_scalars[0]
-        tauW = gradUyW * mu_0 * rhoW
-        u_tau = (tauW / rhoW) ** 0.5
-        uTaus.append(u_tau)
+    tauWs = gradUyWs * mu_0 * rhoWs
+    u_taus = (tauWs / rhoWs) ** 0.5
 
-    return uTaus
+    return u_taus
 
 
 def gridSpacing(mu_0, volmesh):
diff --git a/tests/math/test_ntrfc_vectorcalc.py b/tests/math/test_ntrfc_vectorcalc.py
index ec61df2a97dc6e470148364f7b7c66645edc0ce0..b939b1a85f7d38bce366f997cd24b9301429ff0d 100644
--- a/tests/math/test_ntrfc_vectorcalc.py
+++ b/tests/math/test_ntrfc_vectorcalc.py
@@ -1,3 +1,8 @@
+import numpy as np
+
+from ntrfc.math.vectorcalc import unitvec_list
+
+
 def test_absVec():
     import numpy as np
     from ntrfc.math.vectorcalc import vecAbs
@@ -82,3 +87,18 @@ def test_vecAbs_list():
     c_mags = vecAbs_list(c)
     assert c_mags[0] == 1, "inaccuracy in magnitudes!"
     assert c_mags[1] == 9, "inaccuracy in magnitudes!"
+
+
+def test_unitvec_list():
+    a = np.array([1,0,0])
+    b = np.array([0,0,9])
+    c = np.stack([a,b])
+    directions = unitvec_list(c)
+    assert directions[0] == np.array([1,0,0])
+    assert directions[1] == np.array([0,0,1])
+
+
+def test_unitvec():
+    a = np.array([9,0,0])
+    am = unitvec(a)
+    b = np.array([0,0,-9])