Mercurial > repos > public > sbplib_julia
changeset 2003:524a52f190d7 feature/sbp_operators/laplace_curvilinear
Merge feature/grids/geometry_functions
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Tue, 29 Apr 2025 09:03:05 +0200 |
parents | 496b8d3903c6 (current diff) 4300c59bbeff (diff) |
children | |
files | Manifest.toml src/Grids/Grids.jl |
diffstat | 16 files changed, 617 insertions(+), 74 deletions(-) [+] |
line wrap: on
line diff
--- a/Manifest.toml Tue Feb 11 09:03:04 2025 +0100 +++ b/Manifest.toml Tue Apr 29 09:03:05 2025 +0200 @@ -60,9 +60,9 @@ [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "02c8bd479d26dbeff8a7eb1d77edfc10dacabc01" +git-tree-sha1 = "e3be13f448a43610f978d29b7adf78c76022467a" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.11" +version = "1.9.12" [deps.StaticArrays.extensions] StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/benchmark/Manifest.toml Tue Feb 11 09:03:04 2025 +0100 +++ b/benchmark/Manifest.toml Tue Apr 29 09:03:05 2025 +0200 @@ -61,7 +61,7 @@ deps = ["LinearAlgebra", "StaticArrays", "TOML"] path = ".." uuid = "5a373a26-915f-4769-bcab-bf03835de17b" -version = "0.1.3" +version = "0.1.4" [deps.Diffinitive.extensions] DiffinitiveMakieExt = "Makie" @@ -242,9 +242,9 @@ [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "02c8bd479d26dbeff8a7eb1d77edfc10dacabc01" +git-tree-sha1 = "e3be13f448a43610f978d29b7adf78c76022467a" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.11" +version = "1.9.12" [deps.StaticArrays.extensions] StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/docs/Manifest.toml Tue Feb 11 09:03:04 2025 +0100 +++ b/docs/Manifest.toml Tue Apr 29 09:03:05 2025 +0200 @@ -28,9 +28,9 @@ [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] -git-tree-sha1 = "bce6804e5e6044c6daab27bb533d1295e4a2e759" +git-tree-sha1 = "962834c22b66e32aa10f7611c08c8ca4e20749a9" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.7.6" +version = "0.7.8" [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] @@ -46,7 +46,7 @@ deps = ["LinearAlgebra", "StaticArrays", "TOML"] path = ".." uuid = "5a373a26-915f-4769-bcab-bf03835de17b" -version = "0.1.3" +version = "0.1.4" [deps.Diffinitive.extensions] DiffinitiveMakieExt = "Makie" @@ -67,9 +67,9 @@ [[deps.Documenter]] deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"] -git-tree-sha1 = "d0ea2c044963ed6f37703cead7e29f70cba13d7e" +git-tree-sha1 = "182a9a3fe886587ba230a417f1651a4cbc2b92d4" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "1.8.0" +version = "1.8.1" [[deps.Downloads]] deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] @@ -78,9 +78,9 @@ [[deps.Expat_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "e51db81749b0777b2147fbe7b783ee79045b8e99" +git-tree-sha1 = "d55dffd9ae73ff72f1c0482454dcf2ec6c6c4a63" uuid = "2e619515-83b5-522b-bb60-26c02a35a201" -version = "2.6.4+3" +version = "2.6.5+0" [[deps.FileWatching]] uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" @@ -205,9 +205,9 @@ [[deps.OpenSSL_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "7493f61f55a6cce7325f197443aa80d32554ba10" +git-tree-sha1 = "a9697f1d06cc3eb3fb3ad49cc67f2cfabaac31ea" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "3.0.15+3" +version = "3.0.16+0" [[deps.PCRE2_jll]] deps = ["Artifacts", "Libdl"] @@ -276,9 +276,9 @@ [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "02c8bd479d26dbeff8a7eb1d77edfc10dacabc01" +git-tree-sha1 = "e3be13f448a43610f978d29b7adf78c76022467a" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.11" +version = "1.9.12" [deps.StaticArrays.extensions] StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/docs/make.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/docs/make.jl Tue Apr 29 09:03:05 2025 +0200 @@ -30,7 +30,7 @@ "operator_file_format.md", "grids_and_grid_functions.md", "matrix_and_tensor_representations.md", - "manifolds_charts_atlases.md" + "manifolds_charts_atlases.md", "Submodules" => [ "submodules/grids.md", "submodules/lazy_tensors.md",
--- a/docs/src/manifolds_charts_atlases.md Tue Feb 11 09:03:04 2025 +0100 +++ b/docs/src/manifolds_charts_atlases.md Tue Apr 29 09:03:05 2025 +0200 @@ -1,18 +1,25 @@ # Manifolds, Charts, and Atlases To construct grids on more complicated geometries we use manifolds described -by one or more charts. The charts describe a mapping from some parameter space -to the geometry that we are interested in. If there are more than one chart -for a given geometry this collection of charts and their connection is -described by and atlas. +by one or more charts. The charts describe a mapping from a logical parameter +space to the geometry that we are interested in. If there are more than one +chart for a given geometry this collection of charts and how they are +connected is described by an atlas. + +We consider a mapping from the logical coordidinates ``\xi \in \Xi`` to the +physical coordinates ``x \in \Omega``. A `Chart` describes the mapping by a +`ParameterSpace` respresenting ``\Xi`` and some mapping object that takes +arguments ``\xi \in \Xi`` and returns coordinates ``x\in\Omega``. The mapping +object can either be a function or some other callable object. For the construction of differential and difference operators on a manifold -with a chart the library needs to know the Jacobian of the mapping as a -function of coordinates in the logical parameter space. Internally, -Diffinitive.jl uses a local Jacobian function, `Grids.jacobian(f, ξ)`. For -geometry objects provided by the library this function should have fast and -efficient implementations. If you are creating your own mapping functions you -can implement `Grids.jacobian` for your function or type, for example +with a chart the library needs to know the Jacobian, +``\frac{\partial x}{\partial \xi}``, of the mapping as a function of +coordinates in the logical parameter space. Internally, Diffinitive.jl uses a +local Jacobian function, `Grids.jacobian(f, ξ)`. For geometry objects provided +by the library this function should have fast and efficient implementations. +If you are creating your own mapping functions you must implement +`Grids.jacobian` for your function or type, for example ```julia f(x) = 2x @@ -32,5 +39,3 @@ using ForwardDiff Grids.jacobian(f,x) = ForwardDiff.jacobian(f,x) ``` - -<!-- What more needs to be said here? --/>
--- a/ext/DiffinitivePlotsExt.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/ext/DiffinitivePlotsExt.jl Tue Apr 29 09:03:05 2025 +0200 @@ -57,16 +57,4 @@ return x, y end -# get_axis_limits(plt, :x) - - -# ReicpesPipline/src/user_recipe.jl -# @recipe function f(f::FuncOrFuncs{F}) where {F<:Function} - -# @recipe function f(f::Function, xmin::Number, xmax::Number) - -# _scaled_adapted_grid(f, xscale, yscale, xmin, xmax) - end - -
--- a/src/Grids/Grids.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/src/Grids/Grids.jl Tue Apr 29 09:03:05 2025 +0200 @@ -84,4 +84,17 @@ include("mapped_grid.jl") include("geometry.jl") +function __init__() + if !isdefined(Base.Experimental, :register_error_hint) + return + end + + Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs + if exc.f == Grids.jacobian + print(io, "\nThis possibly means that a function used to define a coordinate mapping is missing a method for `Grids.jacobian`.\n") + print(io, "Provide one by for exmple implementing `Grids.jacobian(::$(typeof(exc.args[1])), x) = ...` or `Grids.jacobian(f, x) = ForwardDiff.jacobian(f,x)`") + end + end +end + end # module
--- a/src/Grids/geometry.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/src/Grids/geometry.jl Tue Apr 29 09:03:05 2025 +0200 @@ -1,32 +1,148 @@ 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{T,PT} +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)(θ) @@ -34,6 +150,16 @@ 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 @@ -56,4 +182,60 @@ s(ξ̄...) end -# TODO: Implement jacobian() for the different mapping helpers +""" + 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
--- a/src/Grids/manifolds.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/src/Grids/manifolds.jl Tue Apr 29 09:03:05 2025 +0200 @@ -9,9 +9,15 @@ end Base.ndims(::Chart{D}) where D = D -(c::Chart)(ξ) = c.mapping(ξ) parameterspace(c::Chart) = c.parameterspace +function (c::Chart)(ξ) + if ξ ∉ parameterspace(c) + throw(DomainError(ξ, "chart was called logical coordinates outside the parameterspace. If this was inteded, use the `mapping` field from the Chart struct instead.")) + end + return c.mapping(ξ) +end + """ jacobian(c::Chart, ξ) @@ -30,8 +36,12 @@ ``` which will both allow calling `jacobian(c,ξ)`. """ -jacobian(c::Chart, ξ) = jacobian(c.mapping, ξ) -# TBD: Can we register a error hint for when jacobian is called with a function that doesn't have a registered jacobian? +function jacobian(c::Chart, ξ) + if ξ ∉ parameterspace(c) + throw(DomainError(ξ, "jacobian was called with logical coordinates outside the parameterspace of the chart. If this was inteded, use the `mapping` field from the Chart struct instead.")) + end + return jacobian(c.mapping, ξ) +end boundary_identifiers(c::Chart) = boundary_identifiers(parameterspace(c)) @@ -54,7 +64,7 @@ """ connections(::Atlas) -Collection of pairs of multiblock boundary identifiers. +Collection of 2-tuples of multiblock boundary identifiers. """ function connections end
--- a/src/Grids/parameter_space.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/src/Grids/parameter_space.jl Tue Apr 29 09:03:05 2025 +0200 @@ -18,6 +18,13 @@ abstract type ParameterSpace{D} end Base.ndims(::ParameterSpace{D}) where D = D +@doc """ + in(x, S::ParameterSpace) + ∈(x, S::ParameterSpace) + +Test if the point `x` is in the parameter space `S`. +""" Base.in(x,::ParameterSpace) + """ Interval{T} <: ParameterSpace{1} @@ -47,6 +54,8 @@ boundary_identifiers(::Interval) = (LowerBoundary(), UpperBoundary()) +Base.in(x, i::Interval) = i.a <= x <= i.b + """ unitinterval(T=Float64) @@ -102,6 +111,11 @@ end end +function Base.in(x, box::HyperBox) + return all(eachindex(x)) do i + box.a[i] <= x[i] <= box.b[i] + end +end """ unitsquare(T=Float64) @@ -150,6 +164,20 @@ 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 + λ = A \ (x - v₁) + + λ_full = (1 - sum(λ), λ...) # Full barycentric coordinates + + return all(λᵢ -> zero(λᵢ) ≤ λᵢ ≤ one(λᵢ), λ_full) +end + """ verticies(s::Simplex)
--- a/test/Grids/equidistant_grid_test.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/test/Grids/equidistant_grid_test.jl Tue Apr 29 09:03:05 2025 +0200 @@ -164,7 +164,6 @@ @test equidistant_grid(HyperBox((0,),(2,)),4) == equidistant_grid(@SVector[0], @SVector[2], 4) end - @testset "equidistant_grid(::Chart)" begin c = Chart(unitsquare()) do (ξ,η) @SVector[2ξ, 3η]
--- a/test/Grids/geometry_test.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/test/Grids/geometry_test.jl Tue Apr 29 09:03:05 2025 +0200 @@ -1,23 +1,274 @@ +using Diffinitive.Grids +using Diffinitive.Grids: Line, LineSegment, linesegments, polygon_edges, Circle, TransfiniteInterpolationSurface, check_transfiniteinterpolation +using StaticArrays + @testset "Line" begin - @test_broken false + @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 - @test_broken false + @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 - @test_broken false + 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 - @test_broken false + 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 - @test_broken false + @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 - @test_broken false + @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
--- a/test/Grids/manifolds_test.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/test/Grids/manifolds_test.jl Tue Apr 29 09:03:05 2025 +0200 @@ -15,14 +15,23 @@ @testset "Chart" begin X(ξ) = 2ξ - Grids.jacobian(::typeof(X), ξ) = @SVector[2,2] + Grids.jacobian(::typeof(X), ξ) = @SMatrix[2 0; 0 2] c = Chart(X, unitsquare()) @test c isa Chart{2} - @test c([3,2]) == [6,4] @test parameterspace(c) == unitsquare() @test ndims(c) == 2 - @test jacobian(c, [3,2]) == [2,2] + @testset "Calling" begin + c = Chart(X, unitsquare()) + @test c([.3,.2]) == [.6,.4] + @test_throws DomainError c([3,2]) + end + + @testset "jacobian(::Chart, ⋅)" begin + c = Chart(X, unitsquare()) + @test_throws DomainError jacobian(c, [3,2]) + @test jacobian(c, [.3,.2]) == [2 0; 0 2] + end @test Set(boundary_identifiers(Chart(X,unitsquare()))) == Set([east(),west(),south(),north()]) end
--- a/test/Grids/parameter_space_test.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/test/Grids/parameter_space_test.jl Tue Apr 29 09:03:05 2025 +0200 @@ -1,6 +1,7 @@ using Test using Diffinitive.Grids +using StaticArrays @testset "ParameterSpace" begin @test ndims(HyperBox([1,1], [2,2])) == 2 @@ -22,6 +23,13 @@ @test limits(unitinterval(Int)) == (0,1) @test boundary_identifiers(unitinterval()) == (LowerBoundary(), UpperBoundary()) + + @test 0 ∈ Interval(0,1) + @test 0. ∈ Interval(0,1) + @test 1. ∈ Interval(0,1) + @test 2 ∉ Interval(0,1) + @test -1 ∉ Interval(0,1) + @test -1. ∉ Interval(0,1) end @testset "HyperBox" begin @@ -59,6 +67,13 @@ CartesianBoundary{3,LowerBoundary}(), CartesianBoundary{3,UpperBoundary}(), ] + + @test @SVector[1.5, 3.5] ∈ HyperBox([1,2], [3,4]) + @test @SVector[1, 2] ∈ HyperBox([1,2], [3,4]) + @test @SVector[3, 4] ∈ HyperBox([1,2], [3,4]) + + @test @SVector[0.5, 3.5] ∉ HyperBox([1,2], [3,4]) + @test @SVector[1.5, 4.5] ∉ HyperBox([1,2], [3,4]) end @testset "Simplex" begin @@ -77,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 + T₃ = 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] ∈ T₃ + @test [0.0, 0.0, 0.0] ∈ T₃ + @test [1.0, 0.0, 0.0] ∈ T₃ + @test [0.25, 0.25, 0.25] ∈ T₃ + @test [0.5, 0.5, 0.0] ∈ T₃ + @test [0.3, 0.3, 0.3] ∈ T₃ + + @test [0.5, 0.5, 1.0] ∉ T₃ + @test [0.3, 0.3, 0.5] ∉ T₃ + end + end end
--- a/test/Manifest.toml Tue Feb 11 09:03:04 2025 +0100 +++ b/test/Manifest.toml Tue Apr 29 09:03:05 2025 +0200 @@ -6,9 +6,9 @@ [[deps.Aqua]] deps = ["Compat", "Pkg", "Test"] -git-tree-sha1 = "3f96fac9ed31d16aba9f2f46fb5cbfc98db6b57c" +git-tree-sha1 = "500941611ff835a025f484f55836f6feea61720a" uuid = "4c88cf16-eb10-579e-8560-4a9242c79595" -version = "0.8.10" +version = "0.8.11" [[deps.ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" @@ -101,14 +101,14 @@ version = "0.2.4" [[deps.JET]] -deps = ["CodeTracking", "InteractiveUtils", "JuliaInterpreter", "LoweredCodeUtils", "MacroTools", "Pkg", "PrecompileTools", "Preferences", "Test"] -git-tree-sha1 = "24bdbf3ef611b69d1f5ef9331e69b6007152989e" +deps = ["CodeTracking", "InteractiveUtils", "JuliaInterpreter", "JuliaSyntax", "LoweredCodeUtils", "MacroTools", "Pkg", "PrecompileTools", "Preferences", "Test"] +git-tree-sha1 = "a453c9b3320dd73f5b05e8882446b6051cb254c4" uuid = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" -version = "0.9.14" +version = "0.9.18" [deps.JET.extensions] JETCthulhuExt = "Cthulhu" - ReviseExt = "Revise" + JETReviseExt = "Revise" [deps.JET.weakdeps] Cthulhu = "f68482b8-f384-11e8-15f7-abe071a5a75f" @@ -128,9 +128,14 @@ [[deps.JuliaInterpreter]] deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"] -git-tree-sha1 = "a729439c18f7112cbbd9fcdc1771ecc7f071df6a" +git-tree-sha1 = "4bf4b400a8234cff0f177da4a160a90296159ce9" uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a" -version = "0.9.39" +version = "0.9.41" + +[[deps.JuliaSyntax]] +git-tree-sha1 = "937da4713526b96ac9a178e2035019d3b78ead4a" +uuid = "70703baa-626e-46a2-a12c-08ffd08c73b4" +version = "0.4.10" [[deps.LRUCache]] git-tree-sha1 = "b3cc6698599b10e652832c2f23db3cab99d51b59" @@ -226,9 +231,9 @@ [[deps.NaNMath]] deps = ["OpenLibm_jll"] -git-tree-sha1 = "fe891aea7ccd23897520db7f16931212454e277e" +git-tree-sha1 = "cc0a5deefdb12ab3a096f00a6d42133af4560d71" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "1.1.1" +version = "1.1.2" [[deps.NetworkOptions]] uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" @@ -351,9 +356,9 @@ [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "02c8bd479d26dbeff8a7eb1d77edfc10dacabc01" +git-tree-sha1 = "e3be13f448a43610f978d29b7adf78c76022467a" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.11" +version = "1.9.12" [deps.StaticArrays.extensions] StaticArraysChainRulesCoreExt = "ChainRulesCore"
--- a/update_manifest.jl Tue Feb 11 09:03:04 2025 +0100 +++ b/update_manifest.jl Tue Apr 29 09:03:05 2025 +0200 @@ -1,12 +1,22 @@ using Pkg +function update_and_log(root, d) + printstyled("Updating '$d'"; color=:cyan, bold=true) + println() + dp = joinpath(root, d) + + update_directory(dp) + println() +end + function update_directory(d) Pkg.activate(d) Pkg.update() - println() end -update_directory(".") -update_directory("benchmark") -update_directory("docs") -update_directory("test") +root = @__DIR__ +update_and_log(root, ".") +update_and_log(root, "benchmark") +update_and_log(root, "docs") +update_and_log(root, "test") +