Merge pull request #2

Implement reference shapes

Project.toml

@@ -1,7 +1,10 @@
 name = "Inti"
 uuid = "fb74042b-437e-4c5b-88cf-d4e2beb394d5"
-authors = ["Luiz M. Faria"]
 version = "0.1.0"
 
+[deps]
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
 [compat]
+StaticArrays = "1"
 julia = "1.6"

src/Geometry/Geometry.jl

@@ -0,0 +1,26 @@
+#=
+
+Definition of basic geometrical concepts.
+
+=#
+
+
+"""
+    ambient_dimension(x)
+
+Dimension of the ambient space where `x` lives. For geometrical objects this can
+differ from its [`geometric_dimension`](@ref); for example a triangle in `ℝ³`
+has ambient dimension `3` but geometric dimension `2`, while a curve in `ℝ³` has
+ambient dimension 3 but geometric dimension 1.
+"""
+function ambient_dimension end
+
+"""
+    geometric_dimension(x)
+
+Number of independent coordinates needed to describe `x`. Lines have geometric
+dimension 1, triangles have geometric dimension 2, etc.
+"""
+function geometric_dimension end
+
+include("referenceshapes.jl")

src/Geometry/referenceshapes.jl

@@ -0,0 +1,77 @@
+"""
+    abstract type AbstractReferenceShape
+
+A fixed reference domain/shape. Used mostly for defining more complex shapes as
+transformations mapping an `AbstractReferenceShape` to some region of `ℜᴹ`.
+
+See e.g. [`ReferenceLine`](@ref) or [`ReferenceTriangle`](@ref) for some
+examples of concrete subtypes.
+"""
+abstract type AbstractReferenceShape end
+
+"""
+    struct ReferenceSimplex{N}
+
+Singleton type representing the N-simplex with N+1 vertices
+`(0,...,0),(0,...,0,1),(0,...,0,1,0),(1,0,...,0)`
+"""
+struct ReferenceSimplex{N} <: AbstractReferenceShape end
+geometric_dimension(::ReferenceSimplex{N}) where {N} = N
+ambient_dimension(::ReferenceSimplex{N}) where {N} = N
+function Base.in(x, ::ReferenceSimplex{N}) where {N}
+    for i in 1:N
+        0 ≤ x[i] ≤ 1 - sum(x[1:i-1]) || return false
+    end
+    return true
+end
+vertices(::ReferenceSimplex{N}) where {N} = svector(i -> i == 1 ? zero(SVector{N, Int64}) : standard_basis_vector(i-1, Val(N)), N+1)
+center(::ReferenceSimplex{N}) where {N} = svector(i -> 1 / (N + 1), N )
+
+"""
+    struct ReferenceTriangle
+
+Singleton type representing the triangle with vertices `(0,0),(1,0),(0,1)`
+"""
+const ReferenceTriangle = ReferenceSimplex{2}
+
+"""
+    struct ReferenceTetrahedron
+
+Singleton type representing the tetrahedron with vertices
+`(0,0,0),(0,0,1),(0,1,0),(1,0,0)`
+"""
+const ReferenceTetrahedron = ReferenceSimplex{3}
+
+"""
+    struct ReferenceHyperCube{N} <: AbstractReferenceShape{N}
+
+Singleton type representing the axis-aligned hypercube in `N` dimensions with
+the lower corner at the origin and the upper corner at `(1,1,…,1)`.
+"""
+struct ReferenceHyperCube{N} <: AbstractReferenceShape end
+geometric_dimension(::ReferenceHyperCube{N}) where {N} = N
+ambient_dimension(::ReferenceHyperCube{N}) where {N} = N
+vertices(::ReferenceHyperCube{N}) where {N} = ntuple(i -> SVector(ntuple(j -> (i >> j) & 1, N)), 2^N)
+Base.in(x, ::ReferenceHyperCube{N}) where {N} = all(0 .≤ x .≤ 1)
+center(::ReferenceHyperCube{N}) where {N} = svector(i -> 0.5, N)
+
+"""
+    const ReferenceLine = ReferenceHyperCube{1}
+
+Singleton type representing the `[0,1]` segment.
+"""
+const ReferenceLine = ReferenceHyperCube{1}
+
+"""
+    const ReferenceSquare = ReferenceHyperCube{2}
+
+Singleton type representing the unit square `[0,1]²`.
+"""
+const ReferenceSquare = ReferenceHyperCube{2}
+
+"""
+    const ReferenceCube = ReferenceHyperCube{3}
+
+Singleton type representing the unit cube `[0,1]³`.
+"""
+const ReferenceCube = ReferenceHyperCube{3}

src/Inti.jl

@@ -1,5 +1,11 @@
 module Inti
 
-# Write your package code here.
+const PROJECT_ROOT = pkgdir(Inti)
+
+using StaticArrays
+
+include("utils.jl")
+include("Geometry/Geometry.jl")
+
 
 end

src/utils.jl

@@ -0,0 +1,20 @@
+#=
+
+Utility functions that have nowhere obvious to go.
+
+=#
+
+"""
+    svector(f,n)
+
+Create an `SVector` of length n, computing each element as f(i), where i is the
+index of the element.
+"""
+svector(f,n)=  ntuple(f,n) |> SVector
+
+"""
+    standard_basis_vector(k, ::Val{N})
+
+Create an `SVector` of length N with a 1 in the kth position and zeros elsewhere.
+"""
+standard_basis_vector(k, n::Val{N}) where {N} = svector(i-> i == k ? 1 : 0, n)

test/Geometry/referenceshapes_test.jl

@@ -0,0 +1,56 @@
+using Test
+using Inti
+using StaticArrays
+
+@testset "Line" begin
+    l = Inti.ReferenceLine()
+    @test Inti.ambient_dimension(l) == 1
+    @test Inti.geometric_dimension(l) == 1
+    x = SVector(0.5)
+    @test x ∈ l
+    x = SVector(1.0)
+    @test x ∈ l
+    x = SVector(1.1)
+    @test !in(x, l)
+    @test Inti.center(l) == SVector(1/2)
+    @test Inti.vertices(l) == (SVector(0), SVector(1))
+end
+@testset "Triangle" begin
+    t = Inti.ReferenceTriangle()
+    @test Inti.ambient_dimension(t) == 2
+    @test Inti.geometric_dimension(t) == 2
+    x = SVector(0.5, 0.5)
+    @test x ∈ t
+    x = SVector(1.0, 0.0)
+    @test x ∈ t
+    x = SVector(1.1, 0.0)
+    @test !in(x, t)
+    @test Inti.center(t) == SVector(1/3, 1/3)
+    @test Inti.vertices(t) == SVector{3, SVector{2, Int64}}([0, 0], [1, 0], [0, 1])
+end
+@testset "Tetrahedron" begin
+    t = Inti.ReferenceTetrahedron()
+    @test Inti.ambient_dimension(t) == 3
+    @test Inti.geometric_dimension(t) == 3
+    x = SVector(0.5, 0.5, 0.0)
+    @test x ∈ t
+    x = SVector(1.0, 0.0, 0.0) # point on edge
+    @test x ∈ t
+    x = SVector(1.1, 0.0, 0.0)
+    @test !in(x, t)
+    @test Inti.center(t) == SVector(1/4, 1/4, 1/4)
+    @test Inti.vertices(t) == SVector{4, SVector{3, Int64}}([0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1])
+end
+@testset "NSimplex" begin
+    t = Inti.ReferenceSimplex{4}()
+    @test Inti.ambient_dimension(t) == 4
+    @test Inti.geometric_dimension(t) == 4
+    x = SVector(0.5, 0.4, 0.0, 0.05)
+    @test x ∈ t
+    x = SVector(1.0, 0.0, 0.0, 0.0)
+    @test x ∈ t
+    x = SVector(1.1, 0.0, 0.0, 0.0)
+    @test !in(x, t)
+    @test Inti.center(t) == SVector(1/5, 1/5, 1/5, 1/5)
+    @test Inti.vertices(t) == SVector{5, SVector{4, Int64}}([0, 0, 0, 0], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1])
+end

test/Project.toml

@@ -1,3 +1,5 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/aqua_test.jl

@@ -0,0 +1,7 @@
+using Inti
+using Test
+using Aqua
+
+@testset "Aqua" begin
+    Aqua.test_all(Inti)
+end

test/referenceshapes_test.jl

@@ -0,0 +1,37 @@
+using Test
+import Inti
+using StaticArrays
+
+@testset "Line" begin
+    l = Inti.ReferenceLine()
+    @test Inti.ambient_dimension(l) == 1
+    @test Inti.geometric_dimension(l) == 1
+    x = SVector(0.5)
+    @test x ∈ l
+    x = SVector(1.0)
+    @test x ∈ l
+    x = SVector(1.1)
+    @test !in(x, l)
+end
+@testset "Triangle" begin
+    t = Inti.ReferenceTriangle()
+    @test Inti.ambient_dimension(t) == 2
+    @test Inti.geometric_dimension(t) == 2
+    x = SVector(0.5, 0.5)
+    @test x ∈ t
+    x = SVector(1.0, 0.0)
+    @test x ∈ t
+    x = SVector(1.1, 0.0)
+    @test !in(x, t)
+end
+@testset "Tetrahedron" begin
+    t = Inti.ReferenceTetrahedron()
+    @test Inti.ambient_dimension(t) == 3
+    @test Inti.geometric_dimension(t) == 3
+    x = SVector(0.5, 0.5, 0.0)
+    @test x ∈ t
+    x = SVector(1.0, 0.0, 0.0) # point on edge
+    @test x ∈ t
+    x = SVector(1.1, 0.0, 0.0)
+    @test !in(x, t)
+end

test/runtests.jl

@@ -1,10 +1,11 @@
 using Inti
-using Test
+using SafeTestsets
 using Aqua
 
-@testset "Inti.jl" begin
-    @testset "Code quality (Aqua.jl)" begin
-        Aqua.test_all(Inti)
-    end
-    # Write your tests here.
+@safetestset "Code quality" begin
+    include("aqua_test.jl")
+end
+
+@safetestset "Geometry" begin
+    @safetestset "Reference shapes" include("Geometry/referenceshapes_test.jl")
 end