diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
index 53e92c7f66..375da8525a 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
@@ -42,6 +42,7 @@ SourceFiles
 #include "point.H"
 #include "primitiveFieldsFwd.H"
 #include "pointHit.H"
+#include "cachedRandom.H"
 #include "Random.H"
 #include "FixedList.H"
 #include "UList.H"
@@ -164,6 +165,10 @@ public:
             //  uniform distribution
             inline Point randomPoint(Random& rndGen) const;
 
+            //- Return a random point in the tetrahedron from a
+            //  uniform distribution
+            inline Point randomPoint(cachedRandom& rndGen) const;
+
             //- Calculate the barycentric coordinates of the given
             //  point, in the same order as a, b, c, d.  Returns the
             //  determinant of the solution.
diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
index 0051951e6b..d770df6a50 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
@@ -240,6 +240,42 @@ inline Point Foam::tetrahedron<Point, PointRef>::randomPoint
 }
 
 
+template<class Point, class PointRef>
+inline Point Foam::tetrahedron<Point, PointRef>::randomPoint
+(
+    cachedRandom& rndGen
+) const
+{
+    // Adapted from
+    // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html
+
+    scalar s = rndGen.sample01<scalar>();
+    scalar t = rndGen.sample01<scalar>();
+    scalar u = rndGen.sample01<scalar>();
+
+    if (s + t > 1.0)
+    {
+        s = 1.0 - s;
+        t = 1.0 - t;
+    }
+
+    if (t + u > 1.0)
+    {
+        scalar tmp = u;
+        u = 1.0 - s - t;
+        t = 1.0 - tmp;
+    }
+    else if (s + t + u > 1.0)
+    {
+        scalar tmp = u;
+        u = s + t + u - 1.0;
+        s = 1.0 - t - tmp;
+    }
+
+    return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_;
+}
+
+
 template<class Point, class PointRef>
 Foam::scalar Foam::tetrahedron<Point, PointRef>::barycentric
 (