Mercurial > repos > public > sbplib_julia

changeset 2009:fe217ad46f9b feature/grids/geometry_functions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 01 May 2025 15:09:13 +0200
parents 4300c59bbeff (diff) df2cbcb7a2b1 (current diff)
children cf0fa2967361
files
diffstat 3 files changed, 516 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/src/Grids/Grids.jl	Thu May 01 15:05:11 2025 +0200
+++ b/src/Grids/Grids.jl	Thu May 01 15:09:13 2025 +0200
@@ -82,6 +82,7 @@
 include("equidistant_grid.jl")
 include("zero_dim_grid.jl")
 include("mapped_grid.jl")
+include("geometry.jl")
 
 function __init__()
     if !isdefined(Base.Experimental, :register_error_hint)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Grids/geometry.jl	Thu May 01 15:09:13 2025 +0200
@@ -0,0 +1,241 @@
+struct Line{PT}
+    p::PT
+    tangent::PT
+
+    Line{PT}(p::PT, tangent::PT) where PT = new{PT}(p,tangent)
+end
+
+"""
+    Line(p,t)
+
+A line, as a callable object, starting at `p` with tangent `t`.
+
+# Example
+```julia-repl
+julia> l = Grids.Line([1,1],[2,1])
+Diffinitive.Grids.Line{StaticArraysCore.SVector{2, Int64}}([1, 1], [2, 1])
+
+julia> l(0)
+2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
+ 1
+ 1
+
+julia> l(1)
+2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
+ 3
+ 2
+```
+
+See also: [`LineSegment`](@ref).
+"""
+function Line(p, t)
+    p = SVector{length(p)}(p)
+    t = SVector{length(t)}(t)
+    p, t = promote(p, t)
+
+    return Line{typeof(p)}(p,t)
+end
+
+(c::Line)(s) = c.p + s*c.tangent
+
+Grids.jacobian(l::Line, t) = l.tangent
+
+struct LineSegment{PT}
+    a::PT
+    b::PT
+
+    LineSegment{PT}(p::PT, tangent::PT) where PT = new{PT}(p,tangent)
+end
+
+"""
+    LineSegment(a,b)
+
+A line segment, as a callable object, from `a` to `b`.
+
+# Example
+```julia-repl
+julia> l = Grids.LineSegment([1,1],[2,1])
+Diffinitive.Grids.LineSegment{StaticArraysCore.SVector{2, Int64}}([1, 1], [2, 1])
+
+julia> l(0)
+2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
+ 1
+ 1
+
+julia> l(0.5)
+2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
+ 1.5
+ 1.0
+
+julia> l(1)
+2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
+ 2
+ 1
+```
+
+See also: [`Line`](@ref).
+"""
+function LineSegment(a, b)
+    a = SVector{length(a)}(a)
+    b = SVector{length(b)}(b)
+    a, b = promote(a, b)
+
+    return LineSegment{typeof(a)}(a,b)
+end
+
+(c::LineSegment)(s) = (1-s)*c.a + s*c.b
+
+Grids.jacobian(c::LineSegment, s) = c.b - c.a
+
+"""
+    linesegments(ps...)
+
+An array of line segments between the points `ps[1]`, `ps[2]`, and so on.
+
+See also: [`polygon_edges`](@ref).
+"""
+function linesegments(ps...)
+    return [LineSegment(ps[i], ps[i+1]) for i ∈ 1:length(ps)-1]
+end
+
+
+"""
+    polygon_edges(ps...)
+
+An array of line segments between the points `ps[1]`, `ps[2]`, and so on
+including the segment between `ps[end]` and `ps[1]`.
+
+See also: [`linesegments`](@ref).
+"""
+function polygon_edges(ps...)
+    n = length(ps)
+    return [LineSegment(ps[i], ps[mod1(i+1,n)]) for i ∈ eachindex(ps)]
+end
+
+struct Circle{PT,T}
+    c::PT
+    r::T
+
+    Circle{PT,T}(c,r) where {PT,T} = new{PT,T}(c,r)
+end
+
+"""
+    Circle(c,r)
+
+A circle with center `c` and radius `r` paramatrized with the angle to the x-axis.
+
+# Example
+```julia-repl
+julia> c = Grids.Circle([1,1], 2)
+Diffinitive.Grids.Circle{StaticArraysCore.SVector{2, Int64}, Int64}([1, 1], 2)
+
+julia> c(0)
+2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
+ 3.0
+ 1.0
+
+julia> c(π/2)
+2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
+ 1.0000000000000002
+ 3.0
+```
+"""
+function Circle(c,r)
+    c = SVector{2}(c)
+    return Circle{typeof(c), typeof(r)}(c,r)
+end
+
+function (C::Circle)(θ)
+    (;c, r) = C
+    c + r*@SVector[cos(θ), sin(θ)]
+end
+
+function Grids.jacobian(C::Circle, θ)
+    (;r) = C
+    r*@SVector[-sin(θ), cos(θ)]
+end
+
+"""
+    TransfiniteInterpolationSurface(c₁, c₂, c₃, c₄)
+
+A surface defined by the transfinite interpolation of the curves `c₁`, `c₂`, `c₃`, and  `c₄`.
+"""
+struct TransfiniteInterpolationSurface{T1,T2,T3,T4}
+    c₁::T1
+    c₂::T2
+    c₃::T3
+    c₄::T4
+end
+
+function (s::TransfiniteInterpolationSurface)(u,v)
+    c₁, c₂, c₃, c₄ = s.c₁, s.c₂, s.c₃, s.c₄
+    P₀₀ = c₁(0)
+    P₁₀ = c₂(0)
+    P₁₁ = c₃(0)
+    P₀₁ = c₄(0)
+    return (1-v)*c₁(u) + u*c₂(v) + v*c₃(1-u) + (1-u)*c₄(1-v) - (
+        (1-u)*(1-v)*P₀₀ + u*(1-v)*P₁₀ + u*v*P₁₁ + (1-u)*v*P₀₁
+    )
+end
+
+function (s::TransfiniteInterpolationSurface)(ξ̄::AbstractArray)
+    s(ξ̄...)
+end
+
+"""
+    check_transfiniteinterpolation(s::TransfiniteInterpolationSurface)
+
+Throw an error if the ends of the curves in the transfinite interpolation do not match.
+"""
+function check_transfiniteinterpolation(s::TransfiniteInterpolationSurface)
+    if check_transfiniteinterpolation(Bool, s)
+        return nothing
+    else
+        throw(ArgumentError("The end of each curve in the transfinite interpolation should be the same as the beginning of the next curve."))
+    end
+end
+
+"""
+    check_transfiniteinterpolation(Bool, s::TransfiniteInterpolationSurface)
+
+Return true if the ends of the curves in the transfinite interpolation match.
+"""
+function check_transfiniteinterpolation(::Type{Bool}, s::TransfiniteInterpolationSurface)
+    if !isapprox(s.c₁(1), s.c₂(0))
+        return false
+    end
+
+    if !isapprox(s.c₂(1), s.c₃(0))
+        return false
+    end
+
+    if !isapprox(s.c₃(1), s.c₄(0))
+        return false
+    end
+
+    if !isapprox(s.c₄(1), s.c₁(0))
+        return false
+    end
+
+    return true
+end
+
+function Grids.jacobian(s::TransfiniteInterpolationSurface, ξ̄)
+    u, v = ξ̄
+
+    c₁, c₂, c₃, c₄ = s.c₁, s.c₂, s.c₃, s.c₄
+    P₀₀ = c₁(0)
+    P₁₀ = c₂(0)
+    P₁₁ = c₃(0)
+    P₀₁ = c₄(0)
+
+    ∂x̄∂ξ₁ = (1-v)*jacobian(c₁,u) + c₂(v) - v*jacobian(c₃,1-u) -c₄(1-v) - (
+        -(1-v)*P₀₀ + (1-v)*P₁₀ + v*P₁₁ - v*P₀₁
+    )
+
+    ∂x̄∂ξ₂ = -c₁(u) + u*jacobian(c₂,v) + c₃(1-u) - (1-u)*jacobian(c₄,1-v) - (
+        -(1-u)*P₀₀ - u*P₁₀ + u*P₁₁ + (1-u)*P₀₁
+    )
+
+    return [∂x̄∂ξ₁ ∂x̄∂ξ₂]
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/Grids/geometry_test.jl	Thu May 01 15:09:13 2025 +0200
@@ -0,0 +1,274 @@
+using Diffinitive.Grids
+using Diffinitive.Grids: Line, LineSegment, linesegments, polygon_edges, Circle, TransfiniteInterpolationSurface, check_transfiniteinterpolation
+using StaticArrays
+
+@testset "Line" begin
+    @testset "Constructors" begin
+        @test Line([1,2],[2,3]) isa Line{SVector{2,Int}}
+        @test Line((1,2),(2,3)) isa Line{SVector{2,Int}}
+        @test Line(@SVector[1,2],[2,3]) isa Line{SVector{2,Int}}
+        @test Line(@SVector[1,2],@SVector[2,3]) isa Line{SVector{2,Int}}
+
+        @test Line([1,2],[2.,3]) isa Line{SVector{2,Float64}}
+        @test Line(@SVector[1,2.],@SVector[2,3]) isa Line{SVector{2,Float64}}
+        @test Line((1,2.),(2,3)) isa Line{SVector{2,Float64}}
+    end
+
+    @testset "Evaluation" begin
+        l = Line([1,2],[2,3])
+
+        @test l(0) == [1,2]
+        @test l(1) == [1,2] + [2,3]
+        @test l(1/2) == [1,2] + [2,3]/2
+    end
+
+    @testset "Grids.jacobian" begin
+        l = Line([1,2],[2,3])
+
+        @test Grids.jacobian(l,0) == [2,3]
+        @test Grids.jacobian(l,1) == [2,3]
+        @test Grids.jacobian(l,1/2) == [2,3]
+    end
+end
+
+@testset "LineSegment" begin
+    @testset "Constructors" begin
+        @test LineSegment([1,2],[2,3]) isa LineSegment{SVector{2,Int}}
+        @test LineSegment((1,2),(2,3)) isa LineSegment{SVector{2,Int}}
+        @test LineSegment(@SVector[1,2],[2,3]) isa LineSegment{SVector{2,Int}}
+        @test LineSegment(@SVector[1,2],@SVector[2,3]) isa LineSegment{SVector{2,Int}}
+
+        @test LineSegment([1,2],[2.,3]) isa LineSegment{SVector{2,Float64}}
+        @test LineSegment(@SVector[1,2.],@SVector[2,3]) isa LineSegment{SVector{2,Float64}}
+        @test LineSegment((1,2.),(2,3)) isa LineSegment{SVector{2,Float64}}
+    end
+
+    @testset "Evaluation" begin
+        l = LineSegment([1,2],[2,3])
+
+        @test l(0) == [1,2]
+        @test l(1) == [2,3]
+        @test l(1/2) == [1,2]/2 + [2,3]/2
+    end
+
+    @testset "Grids.jacobian" begin
+        a, b = [1,2], [2,3]
+        l = LineSegment(a,b)
+        d = b-a
+
+        @test Grids.jacobian(l,0) == d
+        @test Grids.jacobian(l,1) == d
+        @test Grids.jacobian(l,1/2) == d
+    end
+end
+
+@testset "linesegments" begin
+    a,b,c,d = [1,1],[2,2],[3,3],[4,4]
+    @test linesegments(a,b) == [
+        LineSegment(a,b),
+    ]
+
+    @test linesegments(a,b,c) == [
+        LineSegment(a,b),
+        LineSegment(b,c),
+    ]
+
+    @test linesegments(a,b,c,d) == [
+        LineSegment(a,b),
+        LineSegment(b,c),
+        LineSegment(c,d),
+    ]
+end
+
+@testset "polygon_edges" begin
+    a,b,c,d = [1,1],[2,2],[3,3],[4,4]
+    @test polygon_edges(a,b) == [
+        LineSegment(a,b),
+        LineSegment(b,a),
+    ]
+
+    @test polygon_edges(a,b,c) == [
+        LineSegment(a,b),
+        LineSegment(b,c),
+        LineSegment(c,a),
+    ]
+
+    @test polygon_edges(a,b,c,d) == [
+        LineSegment(a,b),
+        LineSegment(b,c),
+        LineSegment(c,d),
+        LineSegment(d,a),
+    ]
+end
+
+@testset "Circle" begin
+    @testset "Constructors" begin
+        @test Circle([1,2], 1) isa Circle{SVector{2,Int},Int}
+        @test Circle([1,2], 1.) isa Circle{SVector{2,Int},Float64}
+        @test Circle([1,2.], 1.) isa Circle{SVector{2,Float64},Float64}
+        @test Circle([1,2.], 1) isa Circle{SVector{2,Float64},Int}
+        @test Circle((1,2.), 1.) isa Circle{SVector{2,Float64},Float64}
+        @test Circle((1,2), 1.) isa Circle{SVector{2,Int},Float64}
+        @test Circle((1.,2), 1) isa Circle{SVector{2,Float64},Int}
+        @test Circle((1,2), 1) isa Circle{SVector{2,Int},Int}
+        @test Circle(@SVector[1,2], 1.) isa Circle{SVector{2,Int},Float64}
+        @test Circle(@SVector[1,2.], 1.) isa Circle{SVector{2,Float64},Float64}
+    end
+
+    @testset "Evaluation" begin
+        c = Circle([0,0], 1)
+        @test c(0) ≈ [1,0]
+        @test c(π/2) ≈ [0,1]
+        @test c(π) ≈ [-1,0]
+        @test c(3π/2) ≈ [0,-1]
+        @test c(π/4) ≈ [1/√(2),1/√(2)]
+
+        c = Circle([0,0], 2)
+        @test c(0) ≈ [2,0]
+        @test c(π/2) ≈ [0,2]
+        @test c(π) ≈ [-2,0]
+        @test c(3π/2) ≈ [0,-2]
+        @test c(π/4) ≈ [√(2),√(2)]
+    end
+
+    @testset "Grids.jacobian" begin
+        c = Circle([0,0], 1)
+        @test Grids.jacobian(c, 0) ≈ [0,1]
+        @test Grids.jacobian(c, π/2) ≈ [-1,0]
+        @test Grids.jacobian(c, π) ≈ [0,-1]
+        @test Grids.jacobian(c, 3π/2) ≈ [1,0]
+        @test Grids.jacobian(c, π/4) ≈ [-1/√(2),1/√(2)]
+
+        c = Circle([0,0], 2)
+        @test Grids.jacobian(c, 0) ≈ [0,2]
+        @test Grids.jacobian(c, π/2) ≈ [-2,0]
+        @test Grids.jacobian(c, π) ≈ [0,-2]
+        @test Grids.jacobian(c, 3π/2) ≈ [2,0]
+        @test Grids.jacobian(c, π/4) ≈ [-√(2),√(2)]
+
+        c = Circle([-1,1], 1)
+        @test Grids.jacobian(c, 0) ≈ [0,1]
+        @test Grids.jacobian(c, π/2) ≈ [-1,0]
+        @test Grids.jacobian(c, π) ≈ [0,-1]
+        @test Grids.jacobian(c, 3π/2) ≈ [1,0]
+        @test Grids.jacobian(c, π/4) ≈ [-1/√(2),1/√(2)]
+
+        c = Circle([-1,1], 2)
+        @test Grids.jacobian(c, 0) ≈ [0,2]
+        @test Grids.jacobian(c, π/2) ≈ [-2,0]
+        @test Grids.jacobian(c, π) ≈ [0,-2]
+        @test Grids.jacobian(c, 3π/2) ≈ [2,0]
+        @test Grids.jacobian(c, π/4) ≈ [-√(2),√(2)]
+    end
+end
+
+@testset "TransfiniteInterpolationSurface" begin
+    @testset "Constructors" begin
+        @test TransfiniteInterpolationSurface(t->[1,2], t->[2,1], t->[0,0], t->[1,1]) isa TransfiniteInterpolationSurface
+
+        cs = polygon_edges([0,0],[1,0],[1,1],[0,1])
+        @test TransfiniteInterpolationSurface(cs...) isa TransfiniteInterpolationSurface
+    end
+
+    @testset "Evaluation" begin
+        cs = polygon_edges([0,0],[1,0],[1,1],[0,1])
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test ti(0,0) == [0,0]
+        @test ti([0,0]) == [0,0]
+        @test ti(1,0) == [1,0]
+        @test ti([1,0]) == [1,0]
+        @test ti(1,1) == [1,1]
+        @test ti([1,1]) == [1,1]
+        @test ti(0,1) == [0,1]
+        @test ti([0,1]) == [0,1]
+
+        @test ti(1/2, 0) == [1/2, 0]
+        @test ti(1/2, 1) == [1/2, 1]
+        @test ti(0,1/2) == [0,1/2]
+        @test ti(1,1/2) == [1,1/2]
+
+
+        a, b, c, d = [1,0],[2,1/4],[2.5,1],[-1/3,1]
+        cs = polygon_edges(a,b,c,d)
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test ti(0,0) == a
+        @test ti(1,0) == b
+        @test ti(1,1) == c
+        @test ti(0,1) == d
+
+        @test ti(1/2, 0) == (a+b)/2
+        @test ti(1/2, 1) == (c+d)/2
+        @test ti(0, 1/2) == (a+d)/2
+        @test ti(1, 1/2) == (b+c)/2
+
+        # TODO: Some test with curved edges?
+    end
+
+    @testset "check_transfiniteinterpolation" begin
+        cs = polygon_edges([0,0],[1,0],[1,1],[0,1])
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test check_transfiniteinterpolation(ti) == nothing
+        @test check_transfiniteinterpolation(Bool, ti) == true
+
+        bad_sides = [
+            LineSegment([0,0],[1/2,0]),
+            LineSegment([1,0],[1,1/2]),
+            LineSegment([1,1],[0,1/2]),
+            LineSegment([0,1],[0,1/2]),
+        ]
+
+        s1 = TransfiniteInterpolationSurface(bad_sides[1],cs[2],cs[3],cs[4])
+        s2 = TransfiniteInterpolationSurface(cs[1],bad_sides[2],cs[3],cs[4])
+        s3 = TransfiniteInterpolationSurface(cs[1],cs[2],bad_sides[3],cs[4])
+        s4 = TransfiniteInterpolationSurface(cs[1],cs[2],cs[3],bad_sides[4])
+
+        @test check_transfiniteinterpolation(Bool, s1) == false
+        @test check_transfiniteinterpolation(Bool, s2) == false
+        @test check_transfiniteinterpolation(Bool, s3) == false
+        @test check_transfiniteinterpolation(Bool, s4) == false
+
+        @test_throws ArgumentError check_transfiniteinterpolation(s1)
+        @test_throws ArgumentError check_transfiniteinterpolation(s2)
+        @test_throws ArgumentError check_transfiniteinterpolation(s3)
+        @test_throws ArgumentError check_transfiniteinterpolation(s4)
+    end
+
+    @testset "Grids.jacobian" begin
+        cs = polygon_edges([0,0],[1,0],[1,1],[0,1])
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test Grids.jacobian(ti, [0,0]) isa SMatrix
+
+        @test Grids.jacobian(ti, [0,0]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1,0]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1,1]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [0,1]) == [1 0; 0 1]
+
+        @test Grids.jacobian(ti, [1/2, 0]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1/2, 1]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [0,1/2]) == [1 0; 0 1]
+        @test Grids.jacobian(ti, [1,1/2]) == [1 0; 0 1]
+
+
+        a, b, c, d = [1,0],[2,1/4],[2.5,1],[-1/3,1]
+        cs = polygon_edges(a,b,c,d)
+        ti = TransfiniteInterpolationSurface(cs...)
+
+        @test Grids.jacobian(ti, [0,0]) == [b-a d-a]
+        @test Grids.jacobian(ti, [1,0]) == [b-a c-b]
+        @test Grids.jacobian(ti, [1,1]) == [c-d c-b]
+        @test Grids.jacobian(ti, [0,1]) == [c-d d-a]
+
+
+        mid(x,y) = (x+y)/2
+        @test Grids.jacobian(ti, [1/2, 0]) ≈ [b-a mid(c,d)-mid(a,b)]
+        @test Grids.jacobian(ti, [1/2, 1]) ≈ [c-d mid(c,d)-mid(a,b)]
+        @test Grids.jacobian(ti, [0, 1/2]) ≈ [mid(b,c)-mid(a,d) d-a]
+        @test Grids.jacobian(ti, [1, 1/2]) ≈ [mid(b,c)-mid(a,d) c-b]
+
+        # TODO: Some test with curved edges?
+    end
+end