diff --git a/src/Meshes.jl b/src/Meshes.jl
index e7b067e91..02b3e98b6 100644
--- a/src/Meshes.jl
+++ b/src/Meshes.jl
@@ -79,6 +79,9 @@ include("domains.jl")
 # utilities
 include("utils.jl")
 
+# iteration
+include("iteration.jl")
+
 # domain indices
 include("indices.jl")
 
@@ -330,6 +333,9 @@ export
   # indices
   indices,
 
+  # iteration
+  eachvertex,
+
   # partitions
   Partition,
   indices,
diff --git a/src/domains/meshes.jl b/src/domains/meshes.jl
index f98cad96b..e1f68fe9f 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)
diff --git a/src/geometries/multigeom.jl b/src/geometries/multigeom.jl
index 451e6f747..86777fb50 100644
--- a/src/geometries/multigeom.jl
+++ b/src/geometries/multigeom.jl
@@ -20,7 +20,12 @@ struct Multi{M<:Manifold,C<:CRS,G<:Geometry{M,C}} <: Geometry{M,C}
 end
 
 # constructor with iterator of geometries
-Multi(geoms) = Multi(collect(geoms))
+
+Multi(geoms) = begin
+  all(gT <: Geometry{M,C} where {M<:Manifold,C<:CRS} for gT in typeof.(geoms)) ||
+  throw("All elements used must be subtypes of the same Geometry{<:Manifold, <:CRS}")
+  Multi(collect(geoms))
+end
 
 # type aliases for convenience
 const MultiPoint{M<:Manifold,C<:CRS} = Multi{M,C,<:Point{M,C}}
@@ -48,9 +53,9 @@ Base.isapprox(m₁::Multi, m₂::Multi; atol=atol(lentype(m₁)), kwargs...) =
 # POLYTOPE
 # ---------
 
-vertex(m::MultiPolytope, ind) = vertices(m)[ind]
+vertex(m::MultiPolytope, ind) = first(Base.Iterators.drop(eachvertex(m), ind - 1)) # nth(itr, n)
 
-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)
 
diff --git a/src/iteration.jl b/src/iteration.jl
new file mode 100644
index 000000000..0b6d22582
--- /dev/null
+++ b/src/iteration.jl
@@ -0,0 +1,32 @@
+struct VertexItr{T}
+  el::T
+end
+
+eachvertex(e::T) where {T} = error("Vertex iterator not implemented for type $T")
+
+eachvertex(e::Union{Mesh,Polytope,MultiPolytope}) = VertexItr(e)
+
+_v_iterate(el::Polytope, i) =
+  (@inline; (i - 1) % UInt < nvertices(el) % UInt ? (@inbounds vertex(el, i), i + 1) : nothing)
+
+_v_iterate(el::Mesh, i) = (@inline; (i - 1) % UInt < nvertices(el) % UInt ? (@inbounds vertex(el, i), i + 1) : nothing)
+
+Base.iterate(itr::VertexItr{<:Mesh}, i=1) = _v_iterate(itr.el, i)
+
+Base.iterate(itr::VertexItr{<:Polytope}, i=1) = _v_iterate(itr.el, i)
+
+Base.iterate(itr::VertexItr{<:Meshes.MultiPolytope}, state=(1, 1)) = begin
+  ig, ivg = state
+  ig > length(itr.el.geoms) && return nothing
+  is = _v_iterate(itr.el.geoms[ig], ivg)
+  is === nothing && return Base.iterate(itr, (ig + 1, 1))
+  v, ivg = is
+  return (v, (ig, ivg))
+end
+
+Base.length(itr::VertexItr{<:Mesh}) = nvertices(itr.el)
+Base.length(itr::VertexItr{<:Polytope}) = nvertices(itr.el)
+Base.length(itr::VertexItr{<:Meshes.MultiPolytope}) = sum(nvertices, itr.el.geoms)
+
+Base.IteratorSize(::VertexItr) = Base.HasLength()
+Base.eltype(::VertexItr) = Point
diff --git a/test/meshes.jl b/test/meshes.jl
index ce1070f8b..fcb44fa39 100644
--- a/test/meshes.jl
+++ b/test/meshes.jl
@@ -127,6 +127,14 @@
   @test minimum(sub) == Point(Polar(T(1), T(2)))
   @test maximum(sub) == Point(Polar(T(4), T(7)))
 
+  # vertex iteration
+  const rg = RegularGrid(Point(1, 1), Point(2, 2), (1, 1))
+  @test collect(eachvertex(rg)) == Point.([(1, 1), (2, 1), (1, 2), (2, 2)])
+  @test @allocated(begin
+    for _ in eachvertex(rg)
+    end
+  end) == 0
+
   # type stability
   grid = RegularGrid((10, 20), Point(Polar(T(0), T(0))), T.((1, 1)))
   @inferred vertex(grid, (1, 1))
@@ -789,6 +797,7 @@ end
   @test crs(mesh) <: Cartesian{NoDatum}
   @test Meshes.lentype(mesh) == ℳ
   @test vertices(mesh) == points
+  @test collect(eachvertex(mesh)) == points
   @test collect(faces(mesh, 2)) == triangles
   @test collect(elements(mesh)) == triangles
   @test nelements(mesh) == 4
@@ -800,6 +809,12 @@ end
   @test area(mesh) ≈ T(1) * u"m^2"
   @test extrema(mesh) == (cart(0, 0), cart(1, 1))
 
+  const cmesh = mesh
+  @test @allocated(begin
+    for _ in eachvertex(cmesh)
+    end
+  end) == 0
+
   # test constructors
   coords = [T.((0, 0)), T.((1, 0)), T.((0, 1)), T.((1, 1)), T.((0.5, 0.5))]
   connec = connect.([(1, 2, 5), (2, 4, 5), (4, 3, 5), (3, 1, 5)], Triangle)
@@ -873,6 +888,8 @@ end
   mesh₂ = SimpleMesh(cart.([(1, 0), (1, 1), (0, 1)]), connect.([(1, 2, 3)]))
   mesh = merge(mesh₁, mesh₂)
   @test vertices(mesh) == [vertices(mesh₁); vertices(mesh₂)]
+  @test collect(eachvertex(mesh)) == vertices(mesh)
+  @test @allocated(iterate(eachvertex(mesh))) == 0
   @test collect(elements(topology(mesh))) == connect.([(1, 2, 3), (4, 5, 6)])
 
   # merge operation with 3D geometries
@@ -880,6 +897,7 @@ end
   mesh₂ = SimpleMesh(cart.([(1, 0, 0), (1, 1, 0), (0, 1, 0), (1, 1, 1)]), connect.([(1, 2, 3, 4)], Tetrahedron))
   mesh = merge(mesh₁, mesh₂)
   @test vertices(mesh) == [vertices(mesh₁); vertices(mesh₂)]
+  @test collect(eachvertex(mesh)) == vertices(mesh)
   @test collect(elements(topology(mesh))) == connect.([(1, 2, 3, 4), (5, 6, 7, 8)], Tetrahedron)
 
   # convert any mesh to SimpleMesh
diff --git a/test/multigeoms.jl b/test/multigeoms.jl
index 0a585fc38..c122277c7 100644
--- a/test/multigeoms.jl
+++ b/test/multigeoms.jl
@@ -10,16 +10,24 @@
   @test crs(multi) <: Cartesian{NoDatum}
   @test Meshes.lentype(multi) == ℳ
   @test vertex(multi, 1) == vertex(poly, 1)
+  @test collect(eachvertex(multi)) == [vertices(poly); vertices(poly)]
   @test vertices(multi) == [vertices(poly); vertices(poly)]
   @test nvertices(multi) == nvertices(poly) + nvertices(poly)
   @test boundary(multi) == merge(boundary(poly), boundary(poly))
   @test rings(multi) == [rings(poly); rings(poly)]
 
+  const cmulti = multi
+  @test @allocated(begin
+    for _ in eachvertex(cmulti)
+    end
+  end) == 0
+
   poly1 = PolyArea(cart.([(0, 0), (1, 0), (1, 1), (0, 1)]))
   poly2 = PolyArea(cart.([(1, 1), (2, 1), (2, 2), (1, 2)]))
   multi = Multi([poly1, poly2])
   @test vertices(multi) == [vertices(poly1); vertices(poly2)]
   @test nvertices(multi) == nvertices(poly1) + nvertices(poly2)
+  @test collect(eachvertex(multi)) == [vertices(poly1); vertices(poly2)]
   @test area(multi) == area(poly1) + area(poly2)
   @test perimeter(multi) == perimeter(poly1) + perimeter(poly2)
   @test centroid(multi) == cart(1, 1)
diff --git a/test/polytopes.jl b/test/polytopes.jl
index 8f6a43b26..ca765adcf 100644
--- a/test/polytopes.jl
+++ b/test/polytopes.jl
@@ -82,6 +82,14 @@
   s = Segment(latlon(0, 135), latlon(0, 45))
   @test_broken measure(s) ≈ 3C / 4
 
+  # vertex iteration
+  const s = Segment(Point(1.0, 1.0, 1.0, 1.0), Point(2.0, 2.0, 2.0, 2.0))
+  @test collect(eachvertex(s)) == [Point(1.0, 1.0, 1.0, 1.0), Point(2.0, 2.0, 2.0, 2.0)]
+  @test @allocated(begin
+    for _ in eachindex(s)
+    end
+  end) == 0
+
   # parameterization
   s = Segment(latlon(45, 0), latlon(45, 90))
   @test s(T(0)) == latlon(45, 0)
@@ -257,6 +265,14 @@ end
   @test @elapsed(issimple(r)) < 0.02
   @test @allocated(issimple(r)) < 950000
 
+  # vertex iteration
+  const r = Ring(Point.([(0, 0), (1, 0), (1, 1), (0, 1)]))
+  @test collect(eachvertex(r)) == Point.([(0, 0), (1, 0), (1, 1), (0, 1)])
+  @test @allocated(begin
+    for _ in eachvertex(r)
+    end
+  end) == 0
+
   # CRS propagation
   r = Ring(merc.([(0, 0), (1, 0), (1, 1), (0, 1)]))
   @test crs(centroid(r)) === crs(r)
@@ -302,6 +318,14 @@ end
   @test vertices(Ngon{3}(pts...)) == verts
   @test vertices(Ngon{3}(tups...)) == verts
 
+  # vertex iteration
+  const ngon = Ngon(pts)
+  @test collect(eachvertex(ngon)) == verts
+  @test @allocated(begin
+    for v in eachvertex(ngon)
+    end
+  end) == 0
+
   NGONS = [Triangle, Quadrangle, Pentagon, Hexagon, Heptagon, Octagon, Nonagon, Decagon]
   NVERT = 3:10
   for (i, NGON) in enumerate(NGONS)