--- a/src/Grids/parameter_space.jl	Thu Apr 24 09:05:30 2025 +0200
+++ b/src/Grids/parameter_space.jl	Thu Apr 24 09:23:12 2025 +0200
@@ -165,6 +165,26 @@
     return Simplex(Tuple(convert(T,v) for v ∈ verticies))
 end
 
+
+function Base.in(x, s::Simplex)
+    v₁ = s.verticies[1]
+    V = map(s.verticies) do v
+        v - v₁
+    end
+
+    A = hcat(V[2:end]...) # matrix with edge vectors as columns
+    b = x - v₁
+
+    # Solve Aλ = b
+    λ = A \ b
+
+    # Compute full barycentric coordinates: first is 1 - sum(λ), then λ
+    λ_full = (1 - sum(λ), λ...)  # Tuple of length NV
+
+    all(λᵢ -> zero(λᵢ) ≤ λᵢ ≤ one(λᵢ), λ_full)
+end
+
+
 """
     verticies(s::Simplex)
 

--- a/test/Grids/parameter_space_test.jl	Thu Apr 24 09:05:30 2025 +0200
+++ b/test/Grids/parameter_space_test.jl	Thu Apr 24 09:23:12 2025 +0200
@@ -92,4 +92,32 @@
     @test verticies(unittetrahedron()) == ([0,0,0], [1,0,0], [0,1,0],[0,0,1])
 
     @test unitsimplex(4) isa Simplex{Float64,4}
+
+    @testset "Base.in" begin
+        @testset "2D" begin
+            T₂ = Simplex([0.0, 0.0], [1.0, 0.0], [0.0, 1.0])
+            @test [0.1, 0.1] ∈ T₂
+            @test [0.3, 0.3] ∈ T₂
+            @test [1.0, 0.0] ∈ T₂
+            @test [0.0, 0.0] ∈ T₂
+            @test [0.0, 1.0] ∈ T₂
+            @test [0.5, 0.5] ∈ T₂
+
+            @test [0.6, 0.6] ∉ T₂
+            @test [-0.1, 0.1] ∉ T₂
+        end
+
+        @testset "3D" begin
+            tetrahedron = Simplex([0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0])
+            @test [0.1, 0.1, 0.1] ∈ tetrahedron
+            @test [0.0, 0.0, 0.0] ∈ tetrahedron
+            @test [1.0, 0.0, 0.0] ∈ tetrahedron
+            @test [0.25, 0.25, 0.25] ∈ tetrahedron
+            @test [0.5, 0.5, 0.0] ∈ tetrahedron
+            @test [0.3, 0.3, 0.3] ∈ tetrahedron
+
+            @test [0.5, 0.5, 1.0] ∉ tetrahedron
+            @test [0.3, 0.3, 0.5] ∉ tetrahedron
+        end
+    end
 end