diff --git a/nutils/numeric.py b/nutils/numeric.py
index da55e0698..a9b72f0d0 100644
--- a/nutils/numeric.py
+++ b/nutils/numeric.py
@@ -145,6 +145,54 @@ def meshgrid(*args, dtype=None):
   assert n == 0
   return grid
 
+def _simplex_grid_helper(bounds):
+  if bounds.ndim != 1 or len(bounds) == 0:
+    raise ValueError
+  nd = len(bounds)
+  spacing = [numpy.sqrt((1+i/2)/(1+i)) for i in range(nd)] # simplex height orthogonal to lower dimension
+  grid = meshgrid(*[numpy.arange(bound, step=step) for step, bound in zip(spacing, bounds)])
+  out_of_bounds = []
+  for idim in range(nd-1):
+    stripes = grid[(idim,)+(slice(None),)*(idim+1)+(slice(1,None,2),)]
+    stripes += spacing[idim] * (idim+1) / (idim+2)
+    if stripes.size and stripes.flat[-1] >= bounds[idim]:
+      out_of_bounds.append(idim)
+  if out_of_bounds:
+    select = numpy.ones(grid.shape[1:], dtype=bool)
+    for idim in out_of_bounds:
+      select[(slice(None),)*(idim)+(-1,)+(slice(1,None,2),)] = False
+    points = grid[:,select].T
+  else:
+    points = grid.reshape(nd, -1).T
+  d = numpy.subtract(bounds, points.max(axis=0))
+  assert (d > 0).all()
+  points += d / 2
+  return points
+
+def simplex_grid(shape, spacing):
+  '''Multi-dimensional generator for equilateral simplex grids.
+
+  Simplex_grid generates a point cloud within an n-dimensional orthotope, which
+  ranges from zero to a specified shape. The point coordinates are spaced in
+  such a way that the nearest neighbours are at distance `spacing`, thus
+  forming vertices of regular simplices. The returned array is two-dimensional,
+  with the first axis being the spatial dimension (matching `shape`) and the
+  second a stacking of the generated points.
+
+  Parameters
+  ----------
+  shape : :class:`tuple`
+    list or tuple of dimensions of the orthotope to be filled.
+  spacing : :class:`float`
+    minimum spacing in the generated point cloud.
+
+  Returns
+  -------
+  :class:`numpy.ndarray`
+  '''
+
+  return _simplex_grid_helper(numpy.divide(shape, spacing)) * spacing
+
 def takediag(A, axis=-2, rmaxis=-1):
   axis = normdim(A.ndim, axis)
   rmaxis = normdim(A.ndim, rmaxis)
diff --git a/tests/test_numeric.py b/tests/test_numeric.py
index d18999c06..b24f725a7 100644
--- a/tests/test_numeric.py
+++ b/tests/test_numeric.py
@@ -226,6 +226,31 @@ def test_dtype(self):
     self.assertEqual(m.dtype, float)
     self.assertAllEqual(m, numpy.ones((1,)))
 
+class simplex_grid(TestCase):
+
+  def simplex_grid(self, shape, spacing):
+    coords = numeric.simplex_grid(shape, spacing)
+    self.assertEqual(coords.ndim, 2)
+    self.assertEqual(coords.shape[1], len(shape))
+    self.assertTrue((coords > 0).all())
+    self.assertTrue((coords < shape).all())
+    mindist = min(numpy.linalg.norm(c1 - c2) for i, c1 in enumerate(coords) for c2 in coords[:i])
+    self.assertAlmostEqual(mindist, spacing)
+    return coords
+
+  def test_1d(self):
+    coords = self.simplex_grid([2], .8)
+    self.assertEqual(len(coords), 3)
+    self.assertAllAlmostEqual(coords[:,0], [.2,1,1.8])
+
+  def test_2d(self):
+    coords = self.simplex_grid([2,3], .8)
+    self.assertEqual(len(coords), 13)
+
+  def test_3d(self):
+    coords = self.simplex_grid([2,3,4], .8)
+    self.assertEqual(len(coords), 82)
+
 class overlapping(TestCase):
 
   def test_pairwise(self):