Merge branch 'JuliaGeometry:master' into JRC_addGreinerHormann

Author: souma4

.github/workflows/CI.yml

@@ -54,7 +54,7 @@ jobs:
         env:
           FLOAT_TYPE: ${{ matrix.float }}
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           file: lcov.info
           token: ${{ secrets.CODECOV_TOKEN }}

Project.toml

@@ -1,7 +1,7 @@
 name = "Meshes"
 uuid = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 authors = ["Júlio Hoffimann and contributors"]
-version = "0.52.0"
+version = "0.52.6"
 
 [deps]
 Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
@@ -18,6 +18,7 @@ ScopedValues = "7e506255-f358-4e82-b7e4-beb19740aa63"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 TransformsBase = "28dd2a49-a57a-4bfb-84ca-1a49db9b96b8"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
@@ -37,6 +38,7 @@ ScopedValues = "1.2"
 SparseArrays = "1.9"
 StaticArrays = "1.0"
 StatsBase = "0.33, 0.34"
+TiledIteration = "0.5"
 TransformsBase = "1.6"
 Unitful = "1.17"
 julia = "1.9"

docs/src/algorithms/neighborsearch.md

@@ -10,6 +10,7 @@ point of reference. This can be performed with search methods:
 
 ```@docs
 NeighborSearchMethod
+BoundedNeighborSearchMethod
 search
 search!
 searchdists

ext/grid.jl

@@ -63,7 +63,7 @@ function vizgridfallback!(plot, M, pdim, edim)
     # decide whether or not to reverse connectivity list
     rfunc = Makie.@lift _reverse($grid)
 
-    verts = Makie.@lift map(asmakie, vertices($grid))
+    verts = Makie.@lift map(asmakie, eachvertex($grid))
     quads = Makie.@lift [GB.QuadFace($rfunc(indices(e))) for e in elements(topology($grid))]
 
     dims = Makie.@lift size($grid)

ext/mesh.jl

@@ -64,7 +64,7 @@ function vizmesh!(plot, ::Type{<:𝔼}, ::Val{2}, ::Val)
     dim = embeddim($mesh)
     nvert = nvertices($mesh)
     nelem = nelements($mesh)
-    verts = vertices($mesh)
+    verts = eachvertex($mesh)
     topo = topology($mesh)
     elems = elements(topo)
 

ext/utils.jl

@@ -34,9 +34,9 @@ asmakie(multis::AbstractVector{<:Multi}) = mapreduce(m -> asmakie.(parent(m)), v
 
 function asmakie(poly::Polygon)
   rs = rings(poly)
-  outer = [asmakie(p) for p in vertices(first(rs))]
+  outer = map(asmakie, eachvertex(rs[1]))
   if hasholes(poly)
-    inners = map(i -> [asmakie(p) for p in vertices(rs[i])], 2:length(rs))
+    inners = map(i -> map(asmakie, eachvertex(rs[i])), 2:length(rs))
     Makie.Polygon(outer, inners)
   else
     Makie.Polygon(outer)

src/Meshes.jl

@@ -20,6 +20,7 @@ using Distances: Haversine, SphericalAngle
 using Distances: evaluate, result_type
 using Rotations: Rotation, QuatRotation, Angle2d
 using Rotations: rotation_between
+using TiledIteration: TileIterator
 using CoordRefSystems: Basic, Projected, Geographic
 using NearestNeighbors: KDTree, BallTree
 using NearestNeighbors: knn, inrange
@@ -220,6 +221,7 @@ export
   vertex,
   vertices,
   nvertices,
+  eachvertex,
   rings,
   segments,
   angles,
@@ -312,6 +314,7 @@ export
   vertex,
   vertices,
   nvertices,
+  eachvertex,
   element,
   elements,
   nelements,

src/boundingboxes.jl

@@ -13,7 +13,7 @@ function boundingbox end
 # FALLBACKS
 # ----------
 
-boundingbox(p::Polytope) = _pboxes(vertices(p))
+boundingbox(p::Polytope) = _pboxes(eachvertex(p))
 
 boundingbox(p::Primitive) = boundingbox(boundary(p))
 
@@ -78,7 +78,7 @@ boundingbox(g::TransformedGrid{<:Any,<:Any,<:OrthoRegularGrid}) = boundingbox(pa
 boundingbox(g::TransformedGrid{<:Any,<:Any,<:OrthoRectilinearGrid}) =
   boundingbox(parent(g)) |> transform(g) |> boundingbox
 
-boundingbox(m::Mesh) = _pboxes(vertices(m))
+boundingbox(m::Mesh) = _pboxes(eachvertex(m))
 
 # ----------------
 # IMPLEMENTATIONS

src/coarsening/regular.jl

@@ -22,26 +22,43 @@ RegularCoarsening(factors::Vararg{Int,N}) where {N} = RegularCoarsening(factors)
 
 function coarsen(grid::OrthoRegularGrid, method::RegularCoarsening)
   factors = fitdims(method.factors, paramdim(grid))
-  RegularGrid(minimum(grid), maximum(grid), dims=size(grid) .÷ factors)
+  dims = _coarsesize(grid, factors)
+  RegularGrid(minimum(grid), maximum(grid), dims=dims)
 end
 
 function coarsen(grid::RectilinearGrid, method::RegularCoarsening)
   factors = fitdims(method.factors, paramdim(grid))
-  dims = vsize(grid)
-  rngs = ntuple(i -> 1:factors[i]:dims[i], paramdim(grid))
+  inds = _coarseinds(grid, factors)
   xyzₛ = xyz(grid)
-  xyzₜ = ntuple(i -> xyzₛ[i][rngs[i]], paramdim(grid))
+  xyzₜ = ntuple(i -> xyzₛ[i][inds[i]], paramdim(grid))
   RectilinearGrid{manifold(grid),crs(grid)}(xyzₜ)
 end
 
 function coarsen(grid::StructuredGrid, method::RegularCoarsening)
   factors = fitdims(method.factors, paramdim(grid))
-  dims = vsize(grid)
-  rngs = ntuple(i -> 1:factors[i]:dims[i], paramdim(grid))
+  inds = _coarseinds(grid, factors)
   XYZₛ = XYZ(grid)
-  XYZₜ = ntuple(i -> XYZₛ[i][rngs...], paramdim(grid))
+  XYZₜ = ntuple(i -> XYZₛ[i][inds...], paramdim(grid))
   StructuredGrid{manifold(grid),crs(grid)}(XYZₜ)
 end
 
 coarsen(grid::TransformedGrid, method::RegularCoarsening) =
   TransformedGrid(coarsen(parent(grid), method), transform(grid))
+
+# -----------------
+# HELPER FUNCTIONS
+# -----------------
+
+function _coarsesize(grid, factors)
+  dims = size(grid)
+  axes = ntuple(i -> 1:dims[i], paramdim(grid))
+  size(TileIterator(axes, factors))
+end
+
+_coarsevsize(grid, factors) = _coarsesize(grid, factors) .+ .!isperiodic(grid)
+
+function _coarseinds(grid, factors)
+  dims = vsize(grid)
+  tdims = _coarsevsize(grid, factors)
+  ntuple(i -> floor.(Int, range(start=1, stop=dims[i], length=tdims[i])), paramdim(grid))
+end

src/discretization.jl

@@ -105,8 +105,7 @@ function discretize(polygon::Polygon, method::BoundaryTriangulationMethod)
 
   # handle degenerate polygons
   if nvertices(cpoly) == 1
-    v = first(vertices(cpoly))
-    points = [v, v, v]
+    points = fill(vertex(cpoly, 1), 3)
     connec = [connect((1, 2, 3))]
     return SimpleMesh(points, connec)
   end
@@ -153,8 +152,8 @@ function discretize(polygon::Polygon, method::BoundaryTriangulationMethod)
 end
 
 function discretizewithin(ring::Ring, method::BoundaryTriangulationMethod)
-  # collect vertices to get rid of static containers
-  points = collect(vertices(ring))
+  # retrieve vertices of ring
+  points = collect(eachvertex(ring))
 
   # discretize within 2D ring with given method
   ring2D = Ring(_proj2D(manifold(ring), points))
@@ -198,7 +197,7 @@ function simplexify(chain::Chain)
   np = nvertices(chain) + isclosed(chain)
   ip = isperiodic(chain)
 
-  points = collect(vertices(chain))
+  points = collect(eachvertex(chain))
   topo = GridTopology((np - 1,), ip)
 
   SimpleMesh(points, topo)

src/discretization/dehn.jl

@@ -22,12 +22,8 @@ with small number of vertices.
 struct DehnTriangulation <: BoundaryTriangulationMethod end
 
 function discretizewithin(ring::Ring{𝔼{2}}, ::DehnTriangulation)
-  # points on resulting mesh
-  points = collect(vertices(ring))
-
-  # Dehn's recursion
+  points = collect(eachvertex(ring))
   connec = dehn1899(points, 1:length(points))
-
   SimpleMesh(points, connec)
 end
 

src/discretization/delaunay.jl

@@ -25,7 +25,7 @@ end
 DelaunayTriangulation(rng=Random.default_rng()) = DelaunayTriangulation(rng)
 
 function discretizewithin(ring::Ring{𝔼{2}}, method::DelaunayTriangulation)
-  points = vertices(ring)
+  points = collect(eachvertex(ring))
   coords = map(p -> ustrip.(to(p)), points)
   bnodes = [1:nvertices(ring); 1]
   triang = triangulate(coords, boundary_nodes=bnodes, rng=method.rng)

src/discretization/fan.jl

@@ -12,7 +12,7 @@ See [https://en.wikipedia.org/wiki/Fan_triangulation]
 struct FanTriangulation <: BoundaryTriangulationMethod end
 
 function discretizewithin(ring::Ring, ::FanTriangulation)
-  points = collect(vertices(ring))
+  points = collect(eachvertex(ring))
   connec = [connect((1, i, i + 1)) for i in 2:(nvertices(ring) - 1)]
   SimpleMesh(points, connec)
 end

src/discretization/held.jl

@@ -45,10 +45,10 @@ function discretizewithin(ring::Ring{𝔼{2}}, method::HeldTriangulation)
   𝒫 = ℛ |> StdCoords()
 
   # points of resulting mesh
-  points = collect(vertices(ℛ))
+  points = collect(eachvertex(ℛ))
 
   # standardized points for algorithm
-  stdpts = collect(vertices(𝒫))
+  stdpts = collect(eachvertex(𝒫))
 
   # keep track of global indices
   I = CircularVector(1:nvertices(𝒫))

src/domains.jl

@@ -36,7 +36,7 @@ Base.isapprox(d1::Domain, d2::Domain; kwargs...) =
 
 Base.getindex(d::Domain, ind::Int) = element(d, ind)
 
-Base.getindex(d::Domain, inds::AbstractVector) = [element(d, ind) for ind in inds]
+Base.getindex(d::Domain, inds::AbstractVector) = [d[ind] for ind in inds]
 
 Base.firstindex(d::Domain) = 1
 

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)
 
@@ -220,8 +227,9 @@ function Base.getindex(g::Grid, I::CartesianIndices)
   odims = size(g)
   cinds = first(I):CartesianIndex(Tuple(last(I)) .+ 1)
   inds = vec(LinearIndices(odims .+ 1)[cinds])
+  points = [vertex(g, ind) for ind in inds]
   periodic = isperiodic(topology(g)) .&& dims .== odims
-  SimpleMesh(vertices(g)[inds], GridTopology(dims, periodic))
+  SimpleMesh(points, GridTopology(dims, periodic))
 end
 
 _asrange(::Int, r::UnitRange{Int}) = r

src/domains/meshes/transformedmesh.jl

@@ -18,7 +18,7 @@ struct TransformedMesh{M<:Manifold,C<:CRS,TP<:Topology,MS<:Mesh,TR<:Transform} <
 end
 
 function TransformedMesh(m::Mesh, t::Transform)
-  p = t(first(vertices(m)))
+  p = t(vertex(m, 1))
   TransformedMesh{manifold(p),crs(p)}(m, t)
 end
 

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)

src/geometries/polytopes.jl

@@ -32,20 +32,21 @@ macro polytope(type, K, N)
   structexpr = if K == 3
     quote
       struct $type{C<:CRS,Mₚ<:Manifold} <: Polytope{$K,𝔼{3},C}
-        vertices::NTuple{$N,Point{Mₚ,C}}
+        vertices::SVector{$N,Point{Mₚ,C}}
       end
     end
   else
     quote
       struct $type{M<:Manifold,C<:CRS} <: Polytope{$K,M,C}
-        vertices::NTuple{$N,Point{M,C}}
+        vertices::SVector{$N,Point{M,C}}
       end
     end
   end
 
   expr = quote
     $Base.@__doc__ $structexpr
 
+    $type(vertices::NTuple{$N,P}) where {P<:Point} = $type(SVector(vertices))
     $type(vertices::Vararg{Tuple,$N}) = $type(Point.(vertices))
     $type(vertices::Vararg{P,$N}) where {P<:Point} = $type(vertices)
   end
@@ -246,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)
 

src/geometries/polytopes/ngon.jl

@@ -21,7 +21,7 @@ are `Triangle` (N=3), `Quadrangle` (N=4), `Pentagon` (N=5), etc.
   `Heptagon`, `Octagon`, `Nonagon`, `Decagon`.
 """
 struct Ngon{N,M<:Manifold,C<:CRS} <: Polygon{M,C}
-  vertices::NTuple{N,Point{M,C}}
+  vertices::SVector{N,Point{M,C}}
   function Ngon{N,M,C}(vertices) where {N,M<:Manifold,C<:CRS}
     if N < 3
       throw(ArgumentError("the number of vertices must be greater than or equal to 3"))
@@ -30,11 +30,13 @@ struct Ngon{N,M<:Manifold,C<:CRS} <: Polygon{M,C}
   end
 end
 
-Ngon{N}(vertices::NTuple{N,Point{M,C}}) where {N,M<:Manifold,C<:CRS} = Ngon{N,M,C}(vertices)
+Ngon{N}(vertices::SVector{N,Point{M,C}}) where {N,M<:Manifold,C<:CRS} = Ngon{N,M,C}(vertices)
+Ngon{N}(vertices::NTuple{N,P}) where {N,P<:Point} = Ngon{N}(SVector(vertices))
 Ngon{N}(vertices::Vararg{P,N}) where {N,P<:Point} = Ngon{N}(vertices)
 Ngon{N}(vertices::Vararg{Tuple,N}) where {N} = Ngon{N}(Point.(vertices))
 
-Ngon(vertices::NTuple{N,Point{M,C}}) where {N,M<:Manifold,C<:CRS} = Ngon{N,M,C}(vertices)
+Ngon(vertices::SVector{N,Point{M,C}}) where {N,M<:Manifold,C<:CRS} = Ngon{N,M,C}(vertices)
+Ngon(vertices::NTuple{N,P}) where {N,P<:Point} = Ngon(SVector(vertices))
 Ngon(vertices::P...) where {P<:Point} = Ngon(vertices)
 Ngon(vertices::Tuple...) = Ngon(Point.(vertices))
 

src/geometries/polytopes/polyarea.jl

@@ -31,7 +31,19 @@ PolyArea(outer...) = PolyArea(collect(outer))
 Base.isapprox(p₁::PolyArea, p₂::PolyArea; atol=atol(lentype(p₁)), kwargs...) =
   length(p₁.rings) == length(p₂.rings) && all(isapprox(r₁, r₂; atol, kwargs...) for (r₁, r₂) in zip(p₁.rings, p₂.rings))
 
-vertices(p::PolyArea) = mapreduce(vertices, vcat, p.rings)
+function vertex(p::PolyArea, ind)
+  offset = 0
+  for r in p.rings
+    nverts = nvertices(r)
+    if ind ≤ offset + nverts
+      return vertex(r, ind - offset)
+    end
+    offset += nverts
+  end
+  throw(BoundsError(p, ind))
+end
+
+vertices(p::PolyArea) = collect(eachvertex(p))
 
 nvertices(p::PolyArea) = mapreduce(nvertices, +, p.rings)