Create Simplex data type. (#744)

* Add Ksimplex data type.

* Rename `Ksimplex` -> `KSimplex`.

* Rename KSimplex -> Simplex.

* Overhaul constructors for Simplex type.

* Define `issimplex` on Simplex type.

* Add some more tests for Simplex.

* Change `materialize` to always splat.

* Rename file `ksimplex.jl` to `simplex.jl`

* Remove mention of IndAdjTopology.

* Some beauty improvements.

* Nit: docstring format.

Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>

* Docstring rewrite.

Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>

* Nit: docstring.

Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>

* Replace `K_` by `N`.

Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com>

* Remove parametric constructors.

Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com>

* Re-add parametric constructor for `K`.

The reason we need this is for `materialize` to work.
Specifically, consider the following code:
```julia
P3 = Point{3}
pts = [P3(0,0,0), P3(-1,-1,0), P3(-1,1,0), P3(1,1,0), P3(1,-1,0)]
simplicies = connect.([(1, 2, 3), (1, 3, 4), (1, 4, 5), (1, 5, 2)],
Simplex)  # or Simplex{2}
splx = materialize(simplicies[1], pts)
```

`Connection` requires `nvertices` to be defined on the Polytope type, so
`connect` constructs the type `Simplex{K}`.
Then, `materialize` tries to construct `PL(pts)` where `PL ==
Simplex{K}`. Therefore we need this constructor to be defined.

* Refactor construction test.

Refer to the last commit for the rationale behind the parametric constructor.

* Add `display` test.

* Fixup whitespace 4 spaces -> 2 spaces.

* Revert "Change `materialize` to always splat."

This reverts commit 74f76ed.

* Support parametric constructor also for vararg, used by materialize.

See also f74d34e.

* Replace `Vararg` by splatted NTuple constructor.

* Change NTuple constructor to Vararg constructor.

* Remove mention of missing constructors.

* Run formatter on all files.

* Apply suggestions from code review

Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>

* Apply suggestions & Update tests

* Add header

* Add more tests

* Apply suggestions

* Update tests

---------

Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com>
Co-authored-by: Elias Carvalho <73039601+eliascarv@users.noreply.github.com>

src/Meshes.jl

@@ -232,6 +232,7 @@ export
 
   # polytopes
   Polytope,
+  Simplex,
   Chain,
   Segment,
   Rope,

src/connectivities.jl

@@ -48,8 +48,8 @@ indices(c::Connectivity) = c.indices
 Connect a list of `indices` from a global vector of [`Point`](@ref)
 into a [`Polytope`](@ref) of type `PL`.
 
-The type `PL` can be a [`Ngon`](@ref) in which case the length of
-the indices is used to identify the actual N-gon type (e.g. Triangle).
+The type `PL` can be a [`Ngon`](@ref) or [`Simplex`](@ref) in which case the length of
+the indices is used to identify the actual polytope type.
 
 Finally, the type `PL` can be ommitted. In this case, the indices are
 assumed to be connected as a [`Ngon`](@ref) or as a [`Segment`](@ref).
@@ -83,6 +83,11 @@ function connect(indices::Tuple, ::Type{Ngon})
   Connectivity{Ngon{N},N}(indices)
 end
 
+function connect(indices::Tuple, ::Type{Simplex})
+  N = length(indices)
+  Connectivity{Simplex{N - 1},N}(indices)
+end
+
 function connect(indices::Tuple)
   N = length(indices)
   N > 2 ? connect(indices, Ngon) : connect(indices, Segment)

src/polytopes.jl

@@ -189,6 +189,12 @@ include("polytopes/tetrahedron.jl")
 include("polytopes/hexahedron.jl")
 include("polytopes/pyramid.jl")
 
+# --------
+# SIMPLEX
+# --------
+
+include("polytopes/simplex.jl")
+
 # -----------------------
 # N-POLYTOPE (FALLBACKS)
 # -----------------------

src/polytopes/simplex.jl

@@ -0,0 +1,39 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    Simplex(p₁, p₂, ..., pₖ, pₖ₊₁)
+
+A K-simplex is a polytope with parametric dimension K and K+1 vertices.
+It generalizes the [`Segment`](@ref), the  [`Triangle`](@ref) and the
+[`Tetrahedron`](@ref) to higher dimensions (K > 3).
+
+See also [`issimplex`](@ref).
+"""
+struct Simplex{K,Dim,T,N} <: Polytope{K,Dim,T}
+  vertices::NTuple{N,Point{Dim,T}}
+end
+
+function Simplex{K}(vertices::NTuple{N,Point{Dim,T}}) where {K,Dim,T,N}
+  N == K + 1 || throw(ArgumentError("number of vertices must be equal to rank K plus one"))
+  K ≤ Dim || throw(ArgumentError("rank K must be less or equal to embedding dimension"))
+  Simplex{K,Dim,T,N}(vertices)
+end
+
+Simplex{K}(vertices::Point{Dim,T}...) where {K,Dim,T} = Simplex{K}(vertices)
+Simplex{K}(vertices::Tuple...) where {K} = Simplex{K}(Point.(vertices))
+
+function Simplex(vertices::NTuple{N,Point{Dim,T}}) where {Dim,T,N}
+  K = N - 1
+  K ≤ Dim || throw(ArgumentError("rank (number of vertices minus one) must be less or equal to embedding dimension"))
+  Simplex{K,Dim,T,N}(vertices)
+end
+
+Simplex(vertices::Point{Dim,T}...) where {Dim,T} = Simplex(vertices)
+Simplex(vertices::Tuple...) = Simplex(Point.(vertices))
+
+nvertices(::Type{<:Simplex{K}}) where {K} = K + 1
+
+Base.isapprox(s₁::Simplex, s₂::Simplex; kwargs...) =
+  all(isapprox(v₁, v₂; kwargs...) for (v₁, v₂) in zip(vertices(s₁), vertices(s₂)))

src/predicates/issimplex.jl

@@ -15,6 +15,8 @@ issimplex(::Type{<:Triangle}) = true
 
 issimplex(::Type{<:Tetrahedron}) = true
 
+issimplex(::Type{<:Simplex}) = true
+
 """
     issimplex(connectivity)
 

test/polytopes.jl

@@ -873,4 +873,53 @@
       └─ Point(0.0, 0.0, 1.0)"""
     end
   end
+
+  @testset "Simplex" begin
+    pts = (P3(1, 2, 3), P3(1, 2, 4))
+    tups = (T.((1, 2, 3)), T.((1, 2, 4)))
+    @test vertices(Simplex(pts)) == pts
+    @test vertices(Simplex(pts...)) == pts
+    @test vertices(Simplex(tups...)) == pts
+
+    @test vertices(Simplex{1}(pts)) == pts
+    @test vertices(Simplex{1}(pts...)) == pts
+    @test vertices(Simplex{1}(tups...)) == pts
+
+    # higher dimensions
+    pts = ntuple(i -> rand(Point{12,T}), 6)
+    @test vertices(Simplex(pts)) == pts
+
+    s = Simplex(P3(1, 2, 3), P3(1, 2, 4), P3(1, 3, 4))
+    @test issimplex(s)
+    @test paramdim(s) == 2
+    @test nvertices(s) == 3
+    @test s ≈ s
+    @test !(s ≈ Simplex(P3(1, 1, 1), P3(2, 2, 2), P3(3, 3, 3)))
+
+    # error: number of vertices must be equal to rank K plus one
+    pts = ntuple(i -> rand(P3), 3)
+    @test_throws ArgumentError Simplex{3}(pts)
+    # error: rank K must be less or equal to embedding dimension
+    pts = ntuple(i -> rand(P3), 5)
+    @test_throws ArgumentError Simplex(pts)
+    @test_throws ArgumentError Simplex{4}(pts)
+
+    s = Simplex(P3(0, 0, 0), P3(1, 0, 0), P3(0, 1, 0), P3(0, 0, 1))
+    @test sprint(show, s) == "Simplex((0.0, 0.0, 0.0), ..., (0.0, 0.0, 1.0))"
+    if T === Float32
+      @test sprint(show, MIME("text/plain"), s) == """
+      Simplex{3,Float32}
+      ├─ Point(0.0f0, 0.0f0, 0.0f0)
+      ├─ Point(1.0f0, 0.0f0, 0.0f0)
+      ├─ Point(0.0f0, 1.0f0, 0.0f0)
+      └─ Point(0.0f0, 0.0f0, 1.0f0)"""
+    else
+      @test sprint(show, MIME("text/plain"), s) == """
+      Simplex{3,Float64}
+      ├─ Point(0.0, 0.0, 0.0)
+      ├─ Point(1.0, 0.0, 0.0)
+      ├─ Point(0.0, 1.0, 0.0)
+      └─ Point(0.0, 0.0, 1.0)"""
+    end
+  end
 end