diff --git a/src/Meshes.jl b/src/Meshes.jl
index e7b067e91..cf5b60247 100644
--- a/src/Meshes.jl
+++ b/src/Meshes.jl
@@ -221,6 +221,7 @@ export
   vertex,
   vertices,
   nvertices,
+  eachvertex,
   rings,
   segments,
   angles,
@@ -313,6 +314,7 @@ export
   vertex,
   vertices,
   nvertices,
+  eachvertex,
   element,
   elements,
   nelements,
diff --git a/src/domains/meshes.jl b/src/domains/meshes.jl
index f98cad96b..5b3f73381 100644
--- a/src/domains/meshes.jl
+++ b/src/domains/meshes.jl
@@ -23,7 +23,7 @@ function vertex end
 
 Return the vertices of the `mesh`.
 """
-vertices(m::Mesh) = [vertex(m, ind) for ind in 1:nvertices(m)]
+vertices(m::Mesh) = collect(eachvertex(m))
 
 """
     nvertices(mesh)
@@ -32,6 +32,13 @@ Return the number of vertices of the `mesh`.
 """
 nvertices(m::Mesh) = nvertices(topology(m))
 
+"""
+    eachvertex(mesh)
+
+Return an iterator for the vertices of the `mesh`.
+"""
+eachvertex(m::Mesh) = (vertex(m, i) for i in 1:nvertices(m))
+
 """
     faces(mesh, rank)
 
diff --git a/src/geometries/multigeom.jl b/src/geometries/multigeom.jl
index 451e6f747..19696ae0a 100644
--- a/src/geometries/multigeom.jl
+++ b/src/geometries/multigeom.jl
@@ -48,12 +48,14 @@ Base.isapprox(m₁::Multi, m₂::Multi; atol=atol(lentype(m₁)), kwargs...) =
 # POLYTOPE
 # ---------
 
-vertex(m::MultiPolytope, ind) = vertices(m)[ind]
+vertex(m::MultiPolytope, ind) = first(Iterators.drop(eachvertex(m), ind - 1))
 
-vertices(m::MultiPolytope) = [vertex for geom in m.geoms for vertex in vertices(geom)]
+vertices(m::MultiPolytope) = collect(eachvertex(m))
 
 nvertices(m::MultiPolytope) = sum(nvertices, m.geoms)
 
+eachvertex(m::MultiPolytope) = (v for g in m.geoms for v in eachvertex(g))
+
 Base.unique(m::MultiPolytope) = unique!(deepcopy(m))
 
 function Base.unique!(m::MultiPolytope)
diff --git a/src/geometries/polytopes.jl b/src/geometries/polytopes.jl
index 3785a6d39..ac16288f9 100644
--- a/src/geometries/polytopes.jl
+++ b/src/geometries/polytopes.jl
@@ -247,10 +247,17 @@ vertices(p::Polytope) = p.vertices
 """
     nvertices(polytope)
 
-Return the number of vertices in the `polytope`.
+Return the number of vertices of the `polytope`.
 """
 nvertices(p::Polytope) = nvertices(typeof(p))
 
+"""
+    eachvertex(polytope)
+
+Return an iterator for the vertices of the `polytope`.
+"""
+eachvertex(p::Polytope) = (vertex(p, i) for i in 1:nvertices(p))
+
 """
     unique(polytope)
 
diff --git a/test/meshes.jl b/test/meshes.jl
index ce1070f8b..031354233 100644
--- a/test/meshes.jl
+++ b/test/meshes.jl
@@ -127,6 +127,10 @@
   @test minimum(sub) == Point(Polar(T(1), T(2)))
   @test maximum(sub) == Point(Polar(T(4), T(7)))
 
+  # vertex iteration
+  grid = RegularGrid((10, 10), cart(0, 0), T.((1, 1)))
+  vertextest(grid)
+
   # type stability
   grid = RegularGrid((10, 20), Point(Polar(T(0), T(0))), T.((1, 1)))
   @inferred vertex(grid, (1, 1))
@@ -501,6 +505,12 @@ end
   @test topology(rg) == topology(cg)
   @test vertices(rg) == vertices(cg)
 
+  # vertex iteration
+  x = range(zero(T), stop=one(T), length=6)
+  y = T[0.0, 0.1, 0.3, 0.7, 0.9, 1.0]
+  grid = RectilinearGrid(x, y)
+  vertextest(grid)
+
   # type stability
   x = range(zero(T), stop=one(T), length=6) * u"mm"
   y = T[0.0, 0.1, 0.3, 0.7, 0.9, 1.0] * u"cm"
@@ -682,6 +692,12 @@ end
   @test topology(sg) == topology(rg)
   @test vertices(sg) == vertices(rg)
 
+  # vertex iteration
+  X = repeat(range(zero(T), stop=one(T), length=6), 1, 6)
+  Y = repeat(T[0.0, 0.1, 0.3, 0.7, 0.9, 1.0]', 6, 1)
+  grid = StructuredGrid(X, Y)
+  vertextest(grid)
+
   # type stability
   X = repeat(range(zero(T), stop=one(T), length=6), 1, 6) * u"mm"
   Y = repeat(T[0.0, 0.1, 0.3, 0.7, 0.9, 1.0]', 6, 1) * u"cm"
@@ -929,6 +945,12 @@ end
   @test vertex(mesh, 4) == points[4]
   @test vertex(mesh, 5) == points[5]
 
+  # vertex iteration
+  points = cart.([(0, 0), (1, 0), (0, 1), (1, 1), (0.5, 0.5)])
+  connec = connect.([(1, 2, 5), (2, 4, 5), (4, 3, 5), (3, 1, 5)], Triangle)
+  mesh = SimpleMesh(points, connec)
+  vertextest(mesh)
+
   points = cart.([(0, 0), (1, 0), (0, 1), (1, 1), (0.5, 0.5)])
   connec = connect.([(1, 2, 5), (2, 4, 5), (4, 3, 5), (3, 1, 5)], Triangle)
   mesh = SimpleMesh(points, connec)
@@ -988,6 +1010,14 @@ end
   trans2 = Translate(T(-10), T(-10))
   @test TransformedMesh(TransformedMesh(grid, trans1), trans2) == TransformedMesh(grid, trans1 → trans2)
 
+  # vertex iteration
+  trans = Identity()
+  points = latlon.([(0, 0), (0, 1), (1, 0), (1, 1), (0.5, 0.5)])
+  connec = connect.([(1, 2, 5), (2, 4, 5), (4, 3, 5), (3, 1, 5)], Triangle)
+  mesh = SimpleMesh(points, connec)
+  tmesh = TransformedMesh(mesh, trans)
+  vertextest(tmesh)
+
   # transforms that change the Manifold and/or CRS
   points = latlon.([(0, 0), (0, 1), (1, 0), (1, 1), (0.5, 0.5)])
   connec = connect.([(1, 2, 5), (2, 4, 5), (4, 3, 5), (3, 1, 5)], Triangle)
diff --git a/test/multigeoms.jl b/test/multigeoms.jl
index 0a585fc38..668582690 100644
--- a/test/multigeoms.jl
+++ b/test/multigeoms.jl
@@ -97,4 +97,20 @@
   poly2 = PolyArea(merc.([(1, 1), (2, 1), (2, 2), (1, 2)]))
   multi = Multi([poly1, poly2])
   @test crs(centroid(multi)) === crs(multi)
+
+  # vertex iteration
+  ring1 = Ring(cart.([(0, 0), (1, 0), (1, 1), (0, 1)]))
+  ring2 = Ring(cart.([(0, 0), (2, 0), (2, 2), (0, 2)]))
+  ring3 = Ring(cart.([(0.2, 0.2), (0.4, 0.2), (0.4, 0.4), (0.2, 0.4)]))
+  ring4 = Ring(cart.([(0.6, 0.2), (0.8, 0.2), (0.8, 0.4), (0.6, 0.4)]))
+  poly1 = PolyArea(ring1)
+  poly2 = PolyArea(ring2)
+  poly3 = PolyArea([ring1, ring3])
+  poly4 = PolyArea([ring2, ring4])
+  multi1 = Multi([ring1, ring2, ring3, ring4])
+  multi2 = Multi([poly1, poly2])
+  multi3 = Multi([poly3, poly4])
+  vertextest(multi1)
+  vertextest(multi2)
+  vertextest(multi3, bytes=3100)
 end
diff --git a/test/polytopes.jl b/test/polytopes.jl
index 8f6a43b26..2f3f26ea8 100644
--- a/test/polytopes.jl
+++ b/test/polytopes.jl
@@ -42,6 +42,7 @@
   s = Segment(cart(0, 0), cart(1, 1))
   equaltest(s)
   isapproxtest(s)
+  vertextest(s)
 
   s = Segment(Point(1.0, 1.0, 1.0, 1.0), Point(2.0, 2.0, 2.0, 2.0))
   @test all(Point(x, x, x, x) ∈ s for x in 1:0.01:2)
@@ -118,10 +119,12 @@ end
   c = Rope(cart(0, 0), cart(1, 0), cart(0, 1))
   equaltest(c)
   isapproxtest(c)
+  vertextest(c)
 
   c = Ring(cart(0, 0), cart(1, 0), cart(0, 1))
   equaltest(c)
   isapproxtest(c)
+  vertextest(c)
 
   # circular equality
   c1 = Ring(cart.([(1, 1), (2, 2), (3, 3)]))
@@ -352,6 +355,7 @@ end
   t = Triangle(cart(0, 0), cart(1, 0), cart(0, 1))
   equaltest(t)
   isapproxtest(t)
+  vertextest(t)
 
   t = Triangle(cart(0, 0), cart(1, 0), cart(0, 1))
   @test perimeter(t) ≈ T(1 + 1 + √2) * u"m"
@@ -466,6 +470,7 @@ end
   q = Quadrangle(cart(0, 0), cart(1, 0), cart(1, 1), cart(0, 1))
   equaltest(q)
   isapproxtest(q)
+  vertextest(q)
 
   q = Quadrangle(cart(0, 0), cart(1, 0), cart(1, 1), cart(0, 1))
   @test_throws DomainError((T(1.2), T(1.2)), "q(u, v) is not defined for u, v outside [0, 1]².") q(T(1.2), T(1.2))
@@ -537,6 +542,7 @@ end
   p = PolyArea(cart(0, 0), cart(1, 0), cart(0, 1))
   equaltest(p)
   isapproxtest(p)
+  vertextest(p)
 
   # COMMAND USED TO GENERATE TEST FILES (VARY --seed = 1, 2, ..., 5)
   # rpg --cluster 30 --algo 2opt --format line --seed 1 --output poly1
@@ -721,6 +727,7 @@ end
   t = Tetrahedron(cart(0, 0, 0), cart(1, 0, 0), cart(0, 1, 0), cart(0, 0, 1))
   equaltest(t)
   isapproxtest(t)
+  vertextest(t)
 
   # CRS propagation
   c1 = Cartesian{WGS84Latest}(T(0), T(0), T(0))
@@ -786,6 +793,7 @@ end
   )
   equaltest(h)
   isapproxtest(h)
+  vertextest(t)
 
   h = Hexahedron(
     cart(0, 0, 0),
@@ -904,6 +912,7 @@ end
   @test m[5] isa Triangle
   equaltest(p)
   isapproxtest(p)
+  vertextest(p)
 
   p = Pyramid(cart(0, 0, 0), cart(1, 0, 0), cart(1, 1, 0), cart(0, 1, 0), cart(0, 0, 1))
   @test sprint(show, p) == "Pyramid((x: 0.0 m, y: 0.0 m, z: 0.0 m), ..., (x: 0.0 m, y: 0.0 m, z: 1.0 m))"
@@ -942,6 +951,7 @@ end
   @test m[5] isa Quadrangle
   equaltest(w)
   isapproxtest(w)
+  vertextest(w)
 
   w = Wedge(cart(0, 0, 0), cart(1, 0, 0), cart(0, 1, 0), cart(0, 0, 1), cart(1, 0, 1), cart(0, 1, 1))
   @test sprint(show, w) == "Wedge((x: 0.0 m, y: 0.0 m, z: 0.0 m), ..., (x: 0.0 m, y: 1.0 m, z: 1.0 m))"
diff --git a/test/testutils.jl b/test/testutils.jl
index d12c2f374..9fbcf8f90 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -176,3 +176,17 @@ function _isapproxtest(g::Geometry, ::Val{3})
   @test isapprox(g, Translate(0u"m", τ32, 0u"m")(g32), atol=1.1τ32)
   @test isapprox(g, Translate(0u"m", 0u"m", τ32)(g32), atol=1.1τ32)
 end
+
+function eachvertexalloc(g)
+  iterate(eachvertex(g)) # precompile run
+  @allocated for _ in eachvertex(g)
+  end
+end
+
+function vertextest(g; bytes=0)
+  @test collect(eachvertex(g)) == vertices(g)
+  @test eachvertexalloc(g) ≤ bytes
+  # type stability
+  @test isconcretetype(eltype(vertices(g)))
+  @inferred vertices(g)
+end