From cb2e2d16e2a138e240b8c44a081dc55d328cc9c0 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Mon, 2 Oct 2023 13:25:38 +0200 Subject: [PATCH 1/5] re-organize and simplify package structure --- src/Geometry/Geometry.jl | 26 ------- src/Geometry/referenceshapes.jl | 77 ------------------- src/Inti.jl | 15 +++- src/utils.jl | 19 +++++ ...hapes_test.jl => reference_shapes_test.jl} | 0 test/referenceshapes_test.jl | 37 --------- test/runtests.jl | 8 +- 7 files changed, 35 insertions(+), 147 deletions(-) delete mode 100644 src/Geometry/Geometry.jl delete mode 100644 src/Geometry/referenceshapes.jl rename test/{Geometry/referenceshapes_test.jl => reference_shapes_test.jl} (100%) delete mode 100644 test/referenceshapes_test.jl diff --git a/src/Geometry/Geometry.jl b/src/Geometry/Geometry.jl deleted file mode 100644 index 580e2bbd..00000000 --- a/src/Geometry/Geometry.jl +++ /dev/null @@ -1,26 +0,0 @@ -#= - -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") diff --git a/src/Geometry/referenceshapes.jl b/src/Geometry/referenceshapes.jl deleted file mode 100644 index 46074341..00000000 --- a/src/Geometry/referenceshapes.jl +++ /dev/null @@ -1,77 +0,0 @@ -""" - 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} diff --git a/src/Inti.jl b/src/Inti.jl index 2e10a22a..ec80a5bd 100644 --- a/src/Inti.jl +++ b/src/Inti.jl @@ -4,8 +4,21 @@ const PROJECT_ROOT = pkgdir(Inti) using StaticArrays +# helper functions include("utils.jl") -include("Geometry/Geometry.jl") + +# basic interpolation and integration +include("reference_shapes.jl") +# include("reference_interpolation.jl") +# include("reference_integration.jl") + +# # geometry and meshes +# include("domain.jl") +# include("mesh.jl") + +# # integral operators +# include("integral_potentials.jl") +# include("integral_operators.jl") end diff --git a/src/utils.jl b/src/utils.jl index 2b5911a8..7a7219cb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -18,3 +18,22 @@ svector(f,n)= ntuple(f,n) |> SVector 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) + +""" + 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) + +NNumber of degrees of freedom necessary to locally represent the geometrical +object. For example, lines have geometric dimension of 1 (whether in `ℝ²` or in +`ℝ³`), while surfaces have geometric dimension of 2. +""" +function geometric_dimension end diff --git a/test/Geometry/referenceshapes_test.jl b/test/reference_shapes_test.jl similarity index 100% rename from test/Geometry/referenceshapes_test.jl rename to test/reference_shapes_test.jl diff --git a/test/referenceshapes_test.jl b/test/referenceshapes_test.jl deleted file mode 100644 index 882355f4..00000000 --- a/test/referenceshapes_test.jl +++ /dev/null @@ -1,37 +0,0 @@ -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 diff --git a/test/runtests.jl b/test/runtests.jl index 9b203116..fc868f58 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,10 +2,6 @@ using Inti using SafeTestsets using Aqua -@safetestset "Code quality" begin - include("aqua_test.jl") -end +@safetestset "Code quality" include("aqua_test.jl") -@safetestset "Geometry" begin - @safetestset "Reference shapes" include("Geometry/referenceshapes_test.jl") -end +@safetestset "Reference shapes" include("reference_shapes_test.jl") From 281ecc14a5af5064e5e08828d304f1fb7433e141 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Mon, 2 Oct 2023 13:32:20 +0200 Subject: [PATCH 2/5] implement `LagrangeElement` interpolants --- src/Inti.jl | 2 +- src/reference_interpolation.jl | 252 +++++++++++++++++++++++++++++++++ 2 files changed, 253 insertions(+), 1 deletion(-) create mode 100644 src/reference_interpolation.jl diff --git a/src/Inti.jl b/src/Inti.jl index ec80a5bd..324e0ca9 100644 --- a/src/Inti.jl +++ b/src/Inti.jl @@ -9,7 +9,7 @@ include("utils.jl") # basic interpolation and integration include("reference_shapes.jl") -# include("reference_interpolation.jl") +include("reference_interpolation.jl") # include("reference_integration.jl") # # geometry and meshes diff --git a/src/reference_interpolation.jl b/src/reference_interpolation.jl new file mode 100644 index 00000000..011aa6c1 --- /dev/null +++ b/src/reference_interpolation.jl @@ -0,0 +1,252 @@ +""" + abstract type AbstractElement{D,T} + +Interpolanting function mapping points on the domain `D<:AbstractReferenceShape` +(of singleton type) to a value of type `T`. + +Instances `el` of `AbstractElement` are expected to implement: +- `el(x̂)`: evaluate the interpolation scheme at the (reference) coordinate `x̂ + ∈ D`. +- `jacobian(el,x̂)` : evaluate the jacobian matrix of the interpolation at the + (reference) coordinate `x ∈ D`. + +!!! note + For performance reasons, both `el(x̂)` and `jacobian(el,x̂)` should + take as input a `StaticVector` and output a static vector or static array. +""" +abstract type AbstractElement{D<:AbstractReferenceShape,T} end + +function (el::AbstractElement)(x) + return abstractmethod(el) +end + +""" + jacobian(f,x) + +Given a (possibly vector-valued) functor `f : 𝐑ᵐ → 𝐅ⁿ`, return the `n × m` +matrix `Aᵢⱼ = ∂fᵢ/∂xⱼ`.Both `x` and `f(x)` are expected to be of `SVector` type. +""" +function jacobian(f, x) + return abstractmethod(f) +end + +domain(::AbstractElement{D,T}) where {D,T} = D() +return_type(::AbstractElement{D,T}) where {D,T} = T +domain_dimension(t::AbstractElement) = domain(t) | center |> length +range_dimension(el::AbstractElement{R,T}) where {R,T} = length(T) + +""" + struct LagrangeElement{D,Np,T} <: AbstractElement{D,T} + +A polynomial `p : D → T` uniquely defined by its `Np` values on the `Np` reference nodes +of `D`. + +The return type `T` should be a vector space (i.e. support addition and +multiplication by scalars). For istance, `T` could be a number or a vector, but +not a `Tuple`. +""" +struct LagrangeElement{D,Np,T} <: AbstractElement{D,T} + vals::SVector{Np,T} +end + +vals(el::LagrangeElement) = el.vals + +""" + reference_nodes(el::LagrangeElement) + +Return the reference nodes on `domain(el)` used for the polynomial +interpolation. The function values on these nodes completely determines the +interpolating polynomial. + +We use the same convention as `gmsh` for defining the reference nodes and their +order (see [node +ordering](https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering) on `gmsh` +documentation). +""" +function reference_nodes(el::LagrangeElement) + return abstractmethod(el) +end + +# infer missig information from type of vals +function LagrangeElement{D}(vals::SVector{Np,T}) where {D,Np,T} + return LagrangeElement{D,Np,T}(vals) +end + +# a more convenient syntax +LagrangeElement{D}(x1,xs...) where {D} = LagrangeElement{D}(SVector(x1,xs...)) + +# construct based on a function +function LagrangeElement{D}(f::Function) where {D} + ref_nodes = reference_nodes(D()) + vals = svector(i -> f(ref_nodes[i]), length(ref_nodes)) + return LagrangeElement{D}(vals) +end + +""" + const LagrangeLine = LagrangeElement{ReferenceLine} +""" +const LagrangeLine = LagrangeElement{ReferenceLine} + +""" + const LagrangeTriangle = LagrangeElement{ReferenceTriangle} +""" +const LagrangeTriangle = LagrangeElement{ReferenceTriangle} + +""" + const LagrangeTetrahedron = LagrangeElement{ReferenceTetrahedron} +""" +const LagrangeTetrahedron = LagrangeElement{ReferenceTetrahedron} + +""" + const LagrangeSquare = LagrangeElement{ReferenceSquare} +""" +const LagrangeSquare = LagrangeElement{ReferenceSquare} + +""" + const LagrangeSquare = LagrangeElement{ReferenceSquare} +""" +const LagrangeCube = LagrangeElement{ReferenceCube} + +#= +Hardcode some basic elements. +TODO: Eventually this could/should be automated. +=# + +# P1 for ReferenceLine +function reference_nodes(::LagrangeLine{2}) + return SVector(SVector(0.0), SVector(1.0)) +end + +function (el::LagrangeLine{2})(u) + v = vals(el) + return v[1] + (v[2] - v[1]) * u[1] +end + +function jacobian(el::LagrangeLine{2}, u) + v = vals(el) + return hcat(v[2] - v[1]) +end + +# P2 for ReferenceLine +function reference_nodes(::LagrangeLine{3}) + return SVector(SVector(0.0), SVector(1.0), SVector(0.5)) +end + +function (el::LagrangeLine{3})(u) + v = vals(el) + return v[1] + (4 * v[3] - 3 * v[1] - v[2]) * u[1] + + 2 * (v[2] + v[1] - 2 * v[3]) * u[1]^2 +end + +function jacobian(el::LagrangeLine{3}, u) + v = vals(el) + return hcat(4 * v[3] - 3 * v[1] - v[2] + 4 * (v[2] + v[1] - 2 * v[3]) * u[1]) +end + +# P1 for ReferenceTriangle +function reference_nodes(::LagrangeTriangle{3}) + return SVector(SVector(0.0, 0.0), SVector(1.0, 0.0), SVector(0.0, 1.0)) +end + +function (el::LagrangeTriangle{3})(u) + v = vals(el) + return v[1] + (v[2] - v[1]) * u[1] + (v[3] - v[1]) * u[2] +end + +function jacobian(el::LagrangeTriangle{3}, u) + v = vals(el) + jac = hcat(v[2] - v[1], v[3] - v[1]) + return jac +end + +# P2 for ReferenceTriangle +function reference_nodes(::LagrangeTriangle{6}) + return SVector(SVector(0.0, 0.0), + SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, 0.0), + SVector(0.5, 0.5), + SVector(0.0, 0.5)) +end + +function (el::LagrangeTriangle{6})(u) + v = vals(el) + return (1 + u[2] * (-3 + 2u[2]) + u[1] * (-3 + 2u[1] + 4u[2])) * v[1] + u[1] * + (-v[2] + u[1] * (2v[2] - 4v[4]) + 4v[4] + u[2] * (-4v[4] + 4v[5] - 4v[6])) + + u[2] * (-v[3] + u[2] * (2v[3] - 4v[6]) + 4v[6]) +end + +function jacobian(el::LagrangeTriangle{6}, u) + v = vals(el) + return hcat((-3 + 4u[1] + 4u[2]) * v[1] - v[2] + + u[1] * (4v[2] - 8v[4]) + + 4v[4] + + u[2] * (-4v[4] + 4v[5] - 4v[6]), + (-3 + 4u[1] + 4u[2]) * v[1] - v[3] + + u[2] * (4v[3] - 8v[6]) + + u[1] * (-4v[4] + 4v[5] - 4v[6]) + + 4v[6]) +end + +# P1 for ReferenceSquare +function reference_nodes(::Type{LagrangeSquare{4}}) + return SVector(SVector(0, 0), + SVector(1, 0), + SVector(1, 1), + SVector(0, 1)) +end + +function (el::LagrangeElement{ReferenceSquare,4})(u) + v = vals(el) + return v[1] + + (v[2] - v[1]) * u[1] + + (v[4] - v[1]) * u[2] + + (v[3] + v[1] - v[2] - v[4]) * u[1] * u[2] +end + +function jacobian(el::LagrangeElement{ReferenceSquare,4}, u) + v = vals(el) + return hcat(((v[2] - v[1]) + (v[3] + v[1] - v[2] - v[4]) * u[2]), + ((v[4] - v[1]) + (v[3] + v[1] - v[2] - v[4]) * u[1])) +end + +# P1 for ReferenceTetrahedron +function reference_nodes(::LagrangeTetrahedron{4}) + return SVector(SVector(0, 0, 0), + SVector(1, 0, 0), + SVector(0, 1, 0), + SVector(0, 0, 1)) +end + +function (el::LagrangeElement{ReferenceTetrahedron,4})(u) + v = vals(el) + return v[1] + (v[2] - v[1]) * u[1] + (v[3] - v[1]) * u[2] + (v[4] - v[1]) * u[3] +end + +function jacobian(el::LagrangeElement{ReferenceTetrahedron,4}, u) + v = vals(el) + return hcat((v[2] - v[1]), (v[3] - v[1]), (v[4] - v[1])) +end + +""" + degree(el::LagrangeElement) + +The polynomial degree of the element. A `LagrangeElement` of degree `K` and +domain `D` belongs to the space [`PolynomialSpace{D,K}`](@ref). +""" +function degree(::Type{LagrangeElement{D,Np}})::Int where {D,Np} + if D == ReferenceLine + return Np - 1 + elseif D == ReferenceTriangle + K = (-3 + sqrt(1 + 8 * Np)) / 2 + return K + elseif D == ReferenceTetrahedron + notimplemented() + elseif D == ReferenceSquare + return sqrt(Np) - 1 + elseif D == ReferenceCube + return Np^(1 / 3) - 1 + else + notimplemented() + end +end From 8f32174e1d08cf9b95c9e8fa07a1a8297a66382a Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Mon, 2 Oct 2023 14:03:12 +0200 Subject: [PATCH 3/5] basic tests for lagrange elements --- src/reference_interpolation.jl | 19 ++++--- src/reference_shapes.jl | 77 ++++++++++++++++++++++++++++ src/utils.jl | 25 +++++++++ test/reference_interpolation_test.jl | 53 +++++++++++++++++++ test/reference_shapes_test.jl | 8 +-- 5 files changed, 170 insertions(+), 12 deletions(-) create mode 100644 src/reference_shapes.jl create mode 100644 test/reference_interpolation_test.jl diff --git a/src/reference_interpolation.jl b/src/reference_interpolation.jl index 011aa6c1..4f83adb1 100644 --- a/src/reference_interpolation.jl +++ b/src/reference_interpolation.jl @@ -17,7 +17,7 @@ Instances `el` of `AbstractElement` are expected to implement: abstract type AbstractElement{D<:AbstractReferenceShape,T} end function (el::AbstractElement)(x) - return abstractmethod(el) + return interface_method(el) end """ @@ -27,13 +27,13 @@ Given a (possibly vector-valued) functor `f : 𝐑ᵐ → 𝐅ⁿ`, return the ` matrix `Aᵢⱼ = ∂fᵢ/∂xⱼ`.Both `x` and `f(x)` are expected to be of `SVector` type. """ function jacobian(f, x) - return abstractmethod(f) + return interface_method(f) end domain(::AbstractElement{D,T}) where {D,T} = D() return_type(::AbstractElement{D,T}) where {D,T} = T -domain_dimension(t::AbstractElement) = domain(t) | center |> length -range_dimension(el::AbstractElement{R,T}) where {R,T} = length(T) +domain_dimension(t::AbstractElement) = domain(t) |> center |> length +range_dimension(el::AbstractElement{R,T}) where {R,T} = domain(t) |> center |> el |> length """ struct LagrangeElement{D,Np,T} <: AbstractElement{D,T} @@ -53,6 +53,7 @@ vals(el::LagrangeElement) = el.vals """ reference_nodes(el::LagrangeElement) + reference_nodes(::Type{<:LagrangeElement}) Return the reference nodes on `domain(el)` used for the polynomial interpolation. The function values on these nodes completely determines the @@ -64,7 +65,7 @@ ordering](https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering) on `gmsh` documentation). """ function reference_nodes(el::LagrangeElement) - return abstractmethod(el) + return interface_method(el) end # infer missig information from type of vals @@ -113,7 +114,7 @@ TODO: Eventually this could/should be automated. =# # P1 for ReferenceLine -function reference_nodes(::LagrangeLine{2}) +function reference_nodes(::Type{LagrangeLine{2}}) return SVector(SVector(0.0), SVector(1.0)) end @@ -128,7 +129,7 @@ function jacobian(el::LagrangeLine{2}, u) end # P2 for ReferenceLine -function reference_nodes(::LagrangeLine{3}) +function reference_nodes(::Type{LagrangeLine{3}}) return SVector(SVector(0.0), SVector(1.0), SVector(0.5)) end @@ -144,7 +145,7 @@ function jacobian(el::LagrangeLine{3}, u) end # P1 for ReferenceTriangle -function reference_nodes(::LagrangeTriangle{3}) +function reference_nodes(::Type{LagrangeTriangle{3}}) return SVector(SVector(0.0, 0.0), SVector(1.0, 0.0), SVector(0.0, 1.0)) end @@ -230,6 +231,7 @@ end """ degree(el::LagrangeElement) + degree(el::Type{<:LagrangeElement}) The polynomial degree of the element. A `LagrangeElement` of degree `K` and domain `D` belongs to the space [`PolynomialSpace{D,K}`](@ref). @@ -250,3 +252,4 @@ function degree(::Type{LagrangeElement{D,Np}})::Int where {D,Np} notimplemented() end end +degree(el::LagrangeElement) = typeof(el) |> degree diff --git a/src/reference_shapes.jl b/src/reference_shapes.jl new file mode 100644 index 00000000..46074341 --- /dev/null +++ b/src/reference_shapes.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} diff --git a/src/utils.jl b/src/utils.jl index 7a7219cb..33a64d9e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -12,6 +12,17 @@ index of the element. """ svector(f,n)= ntuple(f,n) |> SVector +""" + interface_method(x) + +A method of an `abstract type` for which concrete subtypes are expected to +provide an implementation. +""" +function interface_method(T::DataType) + return error("this method needs to be implemented by the concrete subtype $T.") +end +interface_method(x) = interface_method(typeof(x)) + """ standard_basis_vector(k, ::Val{N}) @@ -37,3 +48,17 @@ object. For example, lines have geometric dimension of 1 (whether in `ℝ²` or `ℝ³`), while surfaces have geometric dimension of 2. """ function geometric_dimension end + +""" + return_type(f[,args...]) + +The type returned by `f(args...)`, where `args` is a tuple of types. Falls back +to `Base.promote_op` by default. + +A functors of type `T` with a knonw return type should extend +`return_type(::T,args...)` to avoid relying on `promote_op`. +""" +function return_type(f, args...) + @debug "using `Base.promote_op` to infer return type. Consider defining `return_type(::typeof($f),args...)`." + return Base.promote_op(f, args...) +end diff --git a/test/reference_interpolation_test.jl b/test/reference_interpolation_test.jl new file mode 100644 index 00000000..77764011 --- /dev/null +++ b/test/reference_interpolation_test.jl @@ -0,0 +1,53 @@ +using Test +using StaticArrays +using Inti +using LinearAlgebra + +@testset "Lagrange elements" begin + @testset "LagrangeLine" begin + d = Inti.ReferenceLine() + f = x -> x[1]^2 + x̂ = Inti.reference_nodes(Inti.LagrangeLine{3}) + vals = f.(x̂) + p = Inti.LagrangeLine(vals) + @test p(0) ≈ 0 + @test p(1) ≈ 1 + @test p(0.1) ≈ 0.1^2 + ## line in 3d + vtx = SVector( + SVector(0.0, 0.0, 0.0), + SVector(1.0, 1.0, 1.0) + ) + l = Inti.LagrangeLine(vtx) + @test Inti.domain(l) == Inti.ReferenceLine() + @test l(0.1) ≈ SVector(0.1,0.1,0.1) + ## line in 2d + a = SVector(0.0, 0.0) + b = SVector(1.0, 1.0) + l = Inti.LagrangeLine((a, b)) + @test Inti.domain(l) == Inti.ReferenceLine() + end + @testset "LagrangeTriangle" begin + # triangle in 2d + vtx = SVector( + SVector(0.0, 0.0), + SVector(0.0, 1.0), + SVector(-1.0, 0) + ) + t = Inti.LagrangeTriangle(vtx) + @test Inti.domain_dimension(t) == 2 + @test Inti.range_dimension(t) == 2 + # triangle in 3d + vtx = SVector( + SVector(0.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(-1.0, 0, 0.0) + ) + t = Inti.LagrangeTriangle(vtx) + @test Inti.range_dimension(t) == 3 + @test Inti.domain_dimension(t) == 2 + end + @testset "Tetrahedron" begin + # TODO: add tests + end +end diff --git a/test/reference_shapes_test.jl b/test/reference_shapes_test.jl index 033a7ed5..14519db4 100644 --- a/test/reference_shapes_test.jl +++ b/test/reference_shapes_test.jl @@ -26,7 +26,7 @@ end 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]) + @test Inti.vertices(t) == SVector([0, 0], [1, 0], [0, 1]) end @testset "Tetrahedron" begin t = Inti.ReferenceTetrahedron() @@ -39,7 +39,7 @@ end 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]) + @test Inti.vertices(t) == SVector([0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]) end @testset "NSimplex" begin t = Inti.ReferenceSimplex{4}() @@ -52,5 +52,5 @@ end 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 \ No newline at end of file + @test Inti.vertices(t) == SVector([0, 0, 0, 0], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]) +end From 6025f4c213a6e8d70837e88e0b9d5732fd5f68a2 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Mon, 2 Oct 2023 15:48:09 +0200 Subject: [PATCH 4/5] add reference integration with some basic rules --- src/Inti.jl | 3 +- src/quad_rules_tables.jl | 58 +++++ src/quad_rules_tables_tetrahedron.jl | 326 +++++++++++++++++++++++++++ src/quad_rules_tables_triangle.jl | 215 ++++++++++++++++++ src/reference_integration.jl | 232 +++++++++++++++++++ src/reference_interpolation.jl | 25 +- src/reference_shapes.jl | 12 +- src/utils.jl | 14 ++ test/Manifest.toml | 155 ------------- test/Project.toml | 1 + test/reference_integration_test.jl | 68 ++++++ test/runtests.jl | 4 + 12 files changed, 938 insertions(+), 175 deletions(-) create mode 100644 src/quad_rules_tables.jl create mode 100644 src/quad_rules_tables_tetrahedron.jl create mode 100644 src/quad_rules_tables_triangle.jl create mode 100644 src/reference_integration.jl delete mode 100644 test/Manifest.toml create mode 100644 test/reference_integration_test.jl diff --git a/src/Inti.jl b/src/Inti.jl index 324e0ca9..74e118d3 100644 --- a/src/Inti.jl +++ b/src/Inti.jl @@ -10,7 +10,8 @@ include("utils.jl") # basic interpolation and integration include("reference_shapes.jl") include("reference_interpolation.jl") -# include("reference_integration.jl") +include("quad_rules_tables.jl") +include("reference_integration.jl") # # geometry and meshes # include("domain.jl") diff --git a/src/quad_rules_tables.jl b/src/quad_rules_tables.jl new file mode 100644 index 00000000..06c2bdb0 --- /dev/null +++ b/src/quad_rules_tables.jl @@ -0,0 +1,58 @@ +#= +Tabulated Gaussian quadrature rules from GMSH. +Obtained from `https://gitlab.onelab.info/gmsh/gmsh/-/blob/gmsh_4_7_0/Numeric/GaussQuadratureTri.cpp` +and `https://gitlab.onelab.info/gmsh/gmsh/-/blob/gmsh_4_7_0/Numeric/GaussQuadratureTet.cpp`. +=# + +include("quad_rules_tables_tetrahedron.jl") +include("quad_rules_tables_triangle.jl") + +## +# Dictionaries that contains the quadrature rules +# for various `ReferenceShape`. The dictionary +# key represents the number of quadrature nodes. + +const TRIANGLE_GAUSS_QRULES = Dict(1 => TRIANGLE_G1N1, + 3 => TRIANGLE_G2N3, + 4 => TRIANGLE_G3N4, + 6 => TRIANGLE_G4N6, + 7 => TRIANGLE_G5N7, + 12 => TRIANGLE_G6N12, + 13 => TRIANGLE_G7N13, + 16 => TRIANGLE_G8N16, + 19 => TRIANGLE_G9N19, + 25 => TRIANGLE_G10N25, + 33 => TRIANGLE_G12N33) + +# map a desired quadrature order to the number of nodes +const TRIANGLE_GAUSS_ORDER_TO_NPTS = Dict(1 => 1, + 2 => 3, + 3 => 4, + 4 => 6, + 5 => 7, + 6 => 12, + 7 => 13, + 8 => 16, + 9 => 19, + 10 => 25, + 12 => 33) +const TRIANGLE_GAUSS_NPTS_TO_ORDER = Dict((v, k) for (k, v) in TRIANGLE_GAUSS_ORDER_TO_NPTS) + +const TETAHEDRON_GAUSS_QRULES = Dict(1 => TETAHEDRON_G1N1, + 4 => TETAHEDRON_G2N4, + 5 => TETAHEDRON_G3N5, + 11 => TETAHEDRON_G4N11, + 14 => TETAHEDRON_G5N14, + 24 => TETAHEDRON_G6N24, + 31 => TETAHEDRON_G7N31, + 43 => TETAHEDRON_G8N43) + +# map a desired quadrature order to the number of nodes +const TETRAHEDRON_GAUSS_ORDER_TO_NPTS = Dict(1 => 1, 2 => 4, 3 => 5, 4 => 11, 5 => 14, + 6 => 24, 7 => 31, 8 => 43) +const TETRAHEDRON_GAUSS_NPTS_TO_ORDER = Dict((v, k) + for (k, v) in TETRAHEDRON_GAUSS_ORDER_TO_NPTS) + +## +const GAUSS_QRULES = Dict(ReferenceTriangle => TRIANGLE_GAUSS_QRULES, + ReferenceTetrahedron => TETAHEDRON_GAUSS_QRULES) diff --git a/src/quad_rules_tables_tetrahedron.jl b/src/quad_rules_tables_tetrahedron.jl new file mode 100644 index 00000000..758a045a --- /dev/null +++ b/src/quad_rules_tables_tetrahedron.jl @@ -0,0 +1,326 @@ +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 1 on the tetrahedron *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TETAHEDRON_G1N1 = ((SVector(0.25, 0.25, 0.25), 0.166666666666667),) + +## 0 negative weights, 0 points outside of the tetrahedron, total sum of the +## weights is 1/6 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 2 on the tetrahedron *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TETAHEDRON_G2N4 = ((SVector(0.138196601125, 0.138196601125, 0.138196601125), + 0.0416666666666667), + (SVector(0.585410196625, 0.138196601125, 0.138196601125), + 0.0416666666666667), + (SVector(0.138196601125, 0.585410196625, 0.138196601125), + 0.0416666666666667), + (SVector(0.138196601125, 0.138196601125, 0.585410196625), + 0.0416666666666667)) + +## 0 negative weights, 0 points outside of the tetrahedron, total sum of the +## weights is 1/6 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 3 on the tetrahedron *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TETAHEDRON_G3N5 = ((SVector(0.250000000000, 0.250000000000, 0.250000000000), + -0.133333333333333), + (SVector(0.166666666667, 0.166666666667, 0.166666666667), + +0.075000000000000), + (SVector(0.166666666667, 0.166666666667, 0.500000000000), + +0.075000000000000), + (SVector(0.166666666667, 0.500000000000, 0.166666666667), + +0.075000000000000), + (SVector(0.500000000000, 0.166666666667, 0.166666666667), + +0.075000000000000)) + +## 1 negative weights, 0 points outside of the tetrahedron, total sum of the +## weights is 1/6 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 4 on the tetrahedron *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TETAHEDRON_G4N11 = ((SVector(0.2500000000000, 0.2500000000000, 0.2500000000000), + -0.0131555555555), + (SVector(0.0714285714286, 0.0714285714286, 0.0714285714286), + +0.0076222222222), + (SVector(0.0714285714286, 0.0714285714286, 0.7857142857140), + +0.0076222222222), + (SVector(0.0714285714286, 0.7857142857140, 0.0714285714286), + +0.0076222222222), + (SVector(0.7857142857140, 0.0714285714286, 0.0714285714286), + +0.0076222222222), + (SVector(0.3994035761670, 0.3994035761670, 0.1005964238330), + +0.0248888888888), + (SVector(0.3994035761670, 0.1005964238330, 0.3994035761670), + +0.0248888888888), + (SVector(0.1005964238330, 0.3994035761670, 0.3994035761670), + +0.0248888888888), + (SVector(0.3994035761670, 0.1005964238330, 0.1005964238330), + +0.0248888888888), + (SVector(0.1005964238330, 0.3994035761670, 0.1005964238330), + +0.0248888888888), + (SVector(0.1005964238330, 0.1005964238330, 0.3994035761670), + +0.0248888888888)) + +## 1 negative weights, 0 points outside of the tetrahedron, total sum of the +## weights is 1/6 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 5 on the tetrahedron *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TETAHEDRON_G5N14 = ((SVector(0.0927352503109, 0.0927352503109, 0.0927352503109), + 0.01224884051940), + (SVector(0.7217942490670, 0.0927352503109, 0.0927352503109), + 0.01224884051940), + (SVector(0.0927352503109, 0.7217942490670, 0.0927352503109), + 0.01224884051940), + (SVector(0.0927352503109, 0.0927352503109, 0.7217942490670), + 0.01224884051940), + (SVector(0.3108859192630, 0.3108859192630, 0.3108859192630), + 0.01878132095300), + (SVector(0.0673422422101, 0.3108859192630, 0.3108859192630), + 0.01878132095300), + (SVector(0.3108859192630, 0.0673422422101, 0.3108859192630), + 0.01878132095300), + (SVector(0.3108859192630, 0.3108859192630, 0.0673422422101), + 0.01878132095300), + (SVector(0.4544962958740, 0.4544962958740, 0.0455037041256), + 0.00709100346285), + (SVector(0.4544962958740, 0.0455037041256, 0.4544962958740), + 0.00709100346285), + (SVector(0.0455037041256, 0.4544962958740, 0.4544962958740), + 0.00709100346285), + (SVector(0.4544962958740, 0.0455037041256, 0.0455037041256), + 0.00709100346285), + (SVector(0.0455037041256, 0.4544962958740, 0.0455037041256), + 0.00709100346285), + (SVector(0.0455037041256, 0.0455037041256, 0.4544962958740), + 0.00709100346285)) + +## 0 negative weights, 0 points outside of the tetrahedron, total sum of the +## weights is 1 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 6 on the tetrahedron *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TETAHEDRON_G6N24 = ((SVector(0.2146028712590, 0.2146028712590, 0.2146028712590), + 0.006653791709700), + (SVector(0.3561913862230, 0.2146028712590, 0.2146028712590), + 0.006653791709700), + (SVector(0.2146028712590, 0.3561913862230, 0.2146028712590), + 0.006653791709700), + (SVector(0.2146028712590, 0.2146028712590, 0.3561913862230), + 0.006653791709700), + (SVector(0.0406739585346, 0.0406739585346, 0.0406739585346), + 0.001679535175883), + (SVector(0.8779781243960, 0.0406739585346, 0.0406739585346), + 0.001679535175883), + (SVector(0.0406739585346, 0.8779781243960, 0.0406739585346), + 0.001679535175883), + (SVector(0.0406739585346, 0.0406739585346, 0.8779781243960), + 0.001679535175883), + (SVector(0.3223378901420, 0.3223378901420, 0.3223378901420), + 0.009226196923950), + (SVector(0.0329863295732, 0.3223378901420, 0.3223378901420), + 0.009226196923950), + (SVector(0.3223378901420, 0.0329863295732, 0.3223378901420), + 0.009226196923950), + (SVector(0.3223378901420, 0.3223378901420, 0.0329863295732), + 0.009226196923950), + (SVector(0.0636610018750, 0.0636610018750, 0.2696723314580), + 0.008035714285717), + (SVector(0.0636610018750, 0.2696723314580, 0.0636610018750), + 0.008035714285717), + (SVector(0.0636610018750, 0.0636610018750, 0.6030056647920), + 0.008035714285717), + (SVector(0.0636610018750, 0.6030056647920, 0.0636610018750), + 0.008035714285717), + (SVector(0.0636610018750, 0.2696723314580, 0.6030056647920), + 0.008035714285717), + (SVector(0.0636610018750, 0.6030056647920, 0.2696723314580), + 0.008035714285717), + (SVector(0.2696723314580, 0.0636610018750, 0.0636610018750), + 0.008035714285717), + (SVector(0.2696723314580, 0.0636610018750, 0.6030056647920), + 0.008035714285717), + (SVector(0.2696723314580, 0.6030056647920, 0.0636610018750), + 0.008035714285717), + (SVector(0.6030056647920, 0.0636610018750, 0.2696723314580), + 0.008035714285717), + (SVector(0.6030056647920, 0.0636610018750, 0.0636610018750), + 0.008035714285717), + (SVector(0.6030056647920, 0.2696723314580, 0.0636610018750), + 0.008035714285717)) + +## 0 negative weights, 0 points outside of the tetrahedron, total sum of the +## weights is 1 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 7 on the tetrahedron *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TETAHEDRON_G7N31 = ((SVector(0.50000000000000, 0.50000000000000, 0.00000000000000), + +0.000970017636685), + (SVector(0.50000000000000, 0.00000000000000, 0.50000000000000), + +0.000970017636685), + (SVector(0.00000000000000, 0.50000000000000, 0.50000000000000), + +0.000970017636685), + (SVector(0.00000000000000, 0.00000000000000, 0.50000000000000), + +0.000970017636685), + (SVector(0.00000000000000, 0.50000000000000, 0.00000000000000), + +0.000970017636685), + (SVector(0.50000000000000, 0.00000000000000, 0.00000000000000), + +0.000970017636685), + (SVector(0.25000000000000, 0.25000000000000, 0.25000000000000), + +0.018264223466167), + (SVector(0.07821319233030, 0.07821319233030, 0.07821319233030), + +0.010599941524417), + (SVector(0.07821319233030, 0.07821319233030, 0.76536042300900), + +0.010599941524417), + (SVector(0.07821319233030, 0.76536042300900, 0.07821319233030), + +0.010599941524417), + (SVector(0.76536042300900, 0.07821319233030, 0.07821319233030), + +0.010599941524417), + (SVector(0.12184321666400, 0.12184321666400, 0.12184321666400), + -0.062517740114333), + (SVector(0.12184321666400, 0.12184321666400, 0.63447035000800), + -0.062517740114333), + (SVector(0.12184321666400, 0.63447035000800, 0.12184321666400), + -0.062517740114333), + (SVector(0.63447035000800, 0.12184321666400, 0.12184321666400), + -0.062517740114333), + (SVector(0.33253916444600, 0.33253916444600, 0.33253916444600), + +0.004891425263067), + (SVector(0.33253916444600, 0.33253916444600, 0.00238250666074), + +0.004891425263067), + (SVector(0.33253916444600, 0.00238250666074, 0.33253916444600), + +0.004891425263067), + (SVector(0.00238250666074, 0.33253916444600, 0.33253916444600), + +0.004891425263067), + (SVector(0.10000000000000, 0.10000000000000, 0.20000000000000), + +0.027557319224000), + (SVector(0.10000000000000, 0.20000000000000, 0.10000000000000), + +0.027557319224000), + (SVector(0.10000000000000, 0.10000000000000, 0.60000000000000), + +0.027557319224000), + (SVector(0.10000000000000, 0.60000000000000, 0.10000000000000), + +0.027557319224000), + (SVector(0.10000000000000, 0.20000000000000, 0.60000000000000), + +0.027557319224000), + (SVector(0.10000000000000, 0.60000000000000, 0.20000000000000), + +0.027557319224000), + (SVector(0.20000000000000, 0.10000000000000, 0.10000000000000), + +0.027557319224000), + (SVector(0.20000000000000, 0.10000000000000, 0.60000000000000), + +0.027557319224000), + (SVector(0.20000000000000, 0.60000000000000, 0.10000000000000), + +0.027557319224000), + (SVector(0.60000000000000, 0.10000000000000, 0.20000000000000), + +0.027557319224000), + (SVector(0.60000000000000, 0.10000000000000, 0.10000000000000), + +0.027557319224000), + (SVector(0.60000000000000, 0.20000000000000, 0.10000000000000), + +0.027557319224000)) + +## 4 negative weights, 0 points outside of the tetrahedron + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 8 on the tetrahedron *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TETAHEDRON_G8N43 = ((SVector(0.2500000000000, 0.2500000000000, 0.2500000000000), + -0.020500188658667), + (SVector(0.2068299316110, 0.2068299316110, 0.2068299316110), + +0.014250305822867), + (SVector(0.2068299316110, 0.2068299316110, 0.3795102051680), + +0.014250305822867), + (SVector(0.2068299316110, 0.3795102051680, 0.2068299316110), + +0.014250305822867), + (SVector(0.3795102051680, 0.2068299316110, 0.2068299316110), + +0.014250305822867), + (SVector(0.0821035883105, 0.0821035883105, 0.0821035883105), + +0.001967033313133), + (SVector(0.0821035883105, 0.0821035883105, 0.7536892350680), + +0.001967033313133), + (SVector(0.0821035883105, 0.7536892350680, 0.0821035883105), + +0.001967033313133), + (SVector(0.7536892350680, 0.0821035883105, 0.0821035883105), + +0.001967033313133), + (SVector(0.0057819505052, 0.0057819505052, 0.0057819505052), + +0.000169834109093), + (SVector(0.0057819505052, 0.0057819505052, 0.9826541484840), + +0.000169834109093), + (SVector(0.0057819505052, 0.9826541484840, 0.0057819505052), + +0.000169834109093), + (SVector(0.9826541484840, 0.0057819505052, 0.0057819505052), + +0.000169834109093), + (SVector(0.0505327400189, 0.0505327400189, 0.4494672599810), + +0.004579683824467), + (SVector(0.0505327400189, 0.4494672599810, 0.0505327400189), + +0.004579683824467), + (SVector(0.4494672599810, 0.0505327400189, 0.0505327400189), + +0.004579683824467), + (SVector(0.0505327400189, 0.4494672599810, 0.4494672599810), + +0.004579683824467), + (SVector(0.4494672599810, 0.0505327400189, 0.4494672599810), + +0.004579683824467), + (SVector(0.4494672599810, 0.4494672599810, 0.0505327400189), + +0.004579683824467), + (SVector(0.2290665361170, 0.2290665361170, 0.0356395827885), + +0.005704485808683), + (SVector(0.2290665361170, 0.0356395827885, 0.2290665361170), + +0.005704485808683), + (SVector(0.2290665361170, 0.2290665361170, 0.5062273449780), + +0.005704485808683), + (SVector(0.2290665361170, 0.5062273449780, 0.2290665361170), + +0.005704485808683), + (SVector(0.2290665361170, 0.0356395827885, 0.5062273449780), + +0.005704485808683), + (SVector(0.2290665361170, 0.5062273449780, 0.0356395827885), + +0.005704485808683), + (SVector(0.0356395827885, 0.2290665361170, 0.2290665361170), + +0.005704485808683), + (SVector(0.0356395827885, 0.2290665361170, 0.5062273449780), + +0.005704485808683), + (SVector(0.0356395827885, 0.5062273449780, 0.2290665361170), + +0.005704485808683), + (SVector(0.5062273449780, 0.2290665361170, 0.0356395827885), + +0.005704485808683), + (SVector(0.5062273449780, 0.2290665361170, 0.2290665361170), + +0.005704485808683), + (SVector(0.5062273449780, 0.0356395827885, 0.2290665361170), + +0.005704485808683), + (SVector(0.0366077495532, 0.0366077495532, 0.1904860419350), + +0.002140519141167), + (SVector(0.0366077495532, 0.1904860419350, 0.0366077495532), + +0.002140519141167), + (SVector(0.0366077495532, 0.0366077495532, 0.7362984589590), + +0.002140519141167), + (SVector(0.0366077495532, 0.7362984589590, 0.0366077495532), + +0.002140519141167), + (SVector(0.0366077495532, 0.1904860419350, 0.7362984589590), + +0.002140519141167), + (SVector(0.0366077495532, 0.7362984589590, 0.1904860419350), + +0.002140519141167), + (SVector(0.1904860419350, 0.0366077495532, 0.0366077495532), + +0.002140519141167), + (SVector(0.1904860419350, 0.0366077495532, 0.7362984589590), + +0.002140519141167), + (SVector(0.1904860419350, 0.7362984589590, 0.0366077495532), + +0.002140519141167), + (SVector(0.7362984589590, 0.0366077495532, 0.1904860419350), + +0.002140519141167), + (SVector(0.7362984589590, 0.0366077495532, 0.0366077495532), + +0.002140519141167), + (SVector(0.7362984589590, 0.1904860419350, 0.0366077495532), + +0.002140519141167)) + +## 1 negative weights, 0 points outside of the tetrahedron diff --git a/src/quad_rules_tables_triangle.jl b/src/quad_rules_tables_triangle.jl new file mode 100644 index 00000000..b819a033 --- /dev/null +++ b/src/quad_rules_tables_triangle.jl @@ -0,0 +1,215 @@ +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 1 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G1N1 = ((SVector(0.333333333333333, 0.333333333333333), 0.500000000000000),) +## 0 negative weights, 0 points outside of the triangle, total sum of the +## weights is 0.5 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 2 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G2N3 = ((SVector(0.166666666666667, 0.166666666666667), 0.166666666666667), + (SVector(0.166666666666667, 0.666666666666667), 0.166666666666667), + (SVector(0.666666666666667, 0.166666666666667), 0.166666666666667)) +## 0 negative weights, 0 points outside of the triangle, total sum of the +## weights is 0.5 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 3 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G3N4 = ((SVector(0.333333333333333, 0.333333333333333), -0.281250000000000), + (SVector(0.200000000000000, 0.200000000000000), 0.260416666666667), + (SVector(0.200000000000000, 0.600000000000000), 0.260416666666667), + (SVector(0.600000000000000, 0.200000000000000), 0.260416666666667)) +## 1 negative weight, 0 points outside of the triangle, total sum of the +## weights is 0.5 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 4 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G4N6 = ((SVector(0.445948490915965, 0.445948490915965), 0.111690794839005), + (SVector(0.445948490915965, 0.108103018168070), 0.111690794839005), + (SVector(0.108103018168070, 0.445948490915965), 0.111690794839005), + (SVector(0.091576213509771, 0.091576213509771), 0.054975871827661), + (SVector(0.091576213509771, 0.816847572980459), 0.054975871827661), + (SVector(0.816847572980459, 0.091576213509771), 0.054975871827661)) +## 0 negative weights, 0 points outside of the triangle, total sum of the +## weights is 0.5 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 5 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G5N7 = ((SVector(0.333333333333333, 0.333333333333333), 0.112500000000000), + (SVector(0.470142064105115, 0.470142064105115), 0.066197076394253), + (SVector(0.470142064105115, 0.059715871789770), 0.066197076394253), + (SVector(0.059715871789770, 0.470142064105115), 0.066197076394253), + (SVector(0.101286507323456, 0.101286507323456), 0.062969590272414), + (SVector(0.101286507323456, 0.797426985353087), 0.062969590272414), + (SVector(0.797426985353087, 0.101286507323456), 0.062969590272414)) +## 0 negative weights, 0 points outside of the triangle, total sum of the +## weights is 0.5 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 6 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G6N12 = ((SVector(0.249286745170910, 0.249286745170910), 0.058393137863189), + (SVector(0.249286745170910, 0.501426509658179), 0.058393137863189), + (SVector(0.501426509658179, 0.249286745170910), 0.058393137863189), + (SVector(0.063089014491502, 0.063089014491502), 0.025422453185104), + (SVector(0.063089014491502, 0.873821971016996), 0.025422453185104), + (SVector(0.873821971016996, 0.063089014491502), 0.025422453185104), + (SVector(0.310352451033785, 0.636502499121399), 0.041425537809187), + (SVector(0.636502499121399, 0.053145049844816), 0.041425537809187), + (SVector(0.053145049844816, 0.310352451033785), 0.041425537809187), + (SVector(0.310352451033785, 0.053145049844816), 0.041425537809187), + (SVector(0.636502499121399, 0.310352451033785), 0.041425537809187), + (SVector(0.053145049844816, 0.636502499121399), 0.041425537809187)) +## 0 negative weights, 0 points outside of the triangle, total sum of the +## weights is 0.5 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 7 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G7N13 = ((SVector(0.333333333333333, 0.333333333333333), -0.074785022233841), + (SVector(0.260345966079040, 0.260345966079040), 0.087807628716604), + (SVector(0.260345966079040, 0.479308067841920), 0.087807628716604), + (SVector(0.479308067841920, 0.260345966079040), 0.087807628716604), + (SVector(0.065130102902216, 0.065130102902216), 0.026673617804419), + (SVector(0.065130102902216, 0.869739794195568), 0.026673617804419), + (SVector(0.869739794195568, 0.065130102902216), 0.026673617804419), + (SVector(0.312865496004874, 0.638444188569810), 0.038556880445128), + (SVector(0.638444188569810, 0.048690315425316), 0.038556880445128), + (SVector(0.048690315425316, 0.312865496004874), 0.038556880445128), + (SVector(0.312865496004874, 0.048690315425316), 0.038556880445128), + (SVector(0.638444188569810, 0.312865496004874), 0.038556880445128), + (SVector(0.048690315425316, 0.638444188569810), 0.038556880445128)) +## 1 negative weight, 0 points outside of the triangle, total sum of the +## weights is 0.5 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 8 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G8N16 = ((SVector(0.333333333333333, 0.333333333333333), 0.072157803838894), + (SVector(0.459292588292723, 0.459292588292723), 0.047545817133643), + (SVector(0.459292588292723, 0.081414823414554), 0.047545817133643), + (SVector(0.081414823414554, 0.459292588292723), 0.047545817133643), + (SVector(0.170569307751760, 0.170569307751760), 0.051608685267359), + (SVector(0.170569307751760, 0.658861384496480), 0.051608685267359), + (SVector(0.658861384496480, 0.170569307751760), 0.051608685267359), + (SVector(0.050547228317031, 0.050547228317031), 0.016229248811599), + (SVector(0.050547228317031, 0.898905543365938), 0.016229248811599), + (SVector(0.898905543365938, 0.050547228317031), 0.016229248811599), + (SVector(0.263112829634638, 0.728492392955404), 0.013615157087217), + (SVector(0.728492392955404, 0.008394777409958), 0.013615157087217), + (SVector(0.008394777409958, 0.263112829634638), 0.013615157087217), + (SVector(0.263112829634638, 0.008394777409958), 0.013615157087217), + (SVector(0.728492392955404, 0.263112829634638), 0.013615157087217), + (SVector(0.008394777409958, 0.728492392955404), 0.013615157087217)) +## 0 negative weights, 0 points outside of the triangle, total sum of the +## weights is 0.5 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 9 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G9N19 = ((SVector(0.3333333333330, 0.3333333333333), 0.04856789814140), + (SVector(0.4896825191990, 0.4896825191990), 0.01566735011355), + (SVector(0.4896825191990, 0.0206349616025), 0.01566735011355), + (SVector(0.0206349616025, 0.4896825191990), 0.01566735011355), + (SVector(0.4370895914930, 0.4370895914930), 0.03891377050240), + (SVector(0.4370895914930, 0.1258208170140), 0.03891377050240), + (SVector(0.1258208170140, 0.4370895914930), 0.03891377050240), + (SVector(0.1882035356190, 0.1882035356190), 0.03982386946360), + (SVector(0.1882035356190, 0.6235929287620), 0.03982386946360), + (SVector(0.6235929287620, 0.1882035356190), 0.03982386946360), + (SVector(0.0447295133945, 0.0447295133945), 0.01278883782935), + (SVector(0.0447295133945, 0.9105409732110), 0.01278883782935), + (SVector(0.9105409732110, 0.0447295133945), 0.01278883782935), + (SVector(0.2219629891610, 0.7411985987840), 0.02164176968865), + (SVector(0.7411985987840, 0.0368384120547), 0.02164176968865), + (SVector(0.0368384120547, 0.2219629891610), 0.02164176968865), + (SVector(0.2219629891610, 0.0368384120547), 0.02164176968865), + (SVector(0.7411985987840, 0.2219629891610), 0.02164176968865), + (SVector(0.0368384120547, 0.7411985987840), 0.02164176968865)) +## 0 negative weights, 0 points outside of the triangle, total sum of the +## weights is 0.5 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 10 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G10N25 = ((SVector(0.3333333333330, 0.3333333333330), 0.04540899519140), + (SVector(0.4855776333840, 0.4855776333840), 0.01836297887825), + (SVector(0.4855776333840, 0.0288447332327), 0.01836297887825), + (SVector(0.0288447332327, 0.4855776333840), 0.01836297887825), + (SVector(0.1094815754850, 0.1094815754850), 0.02266052971775), + (SVector(0.1094815754850, 0.7810368490300), 0.02266052971775), + (SVector(0.7810368490300, 0.1094815754850), 0.02266052971775), + (SVector(0.3079398387640, 0.5503529418210), 0.03637895842270), + (SVector(0.5503529418210, 0.1417072194150), 0.03637895842270), + (SVector(0.1417072194150, 0.3079398387640), 0.03637895842270), + (SVector(0.3079398387640, 0.1417072194150), 0.03637895842270), + (SVector(0.5503529418210, 0.3079398387640), 0.03637895842270), + (SVector(0.1417072194150, 0.5503529418210), 0.03637895842270), + (SVector(0.2466725606400, 0.7283239045970), 0.01416362126555), + (SVector(0.7283239045970, 0.0250035347627), 0.01416362126555), + (SVector(0.0250035347627, 0.2466725606400), 0.01416362126555), + (SVector(0.2466725606400, 0.0250035347627), 0.01416362126555), + (SVector(0.7283239045970, 0.2466725606400), 0.01416362126555), + (SVector(0.0250035347627, 0.7283239045970), 0.01416362126555), + (SVector(0.0668032510122, 0.9236559335870), 0.00471083348185), + (SVector(0.9236559335870, 0.0095408154003), 0.00471083348185), + (SVector(0.0095408154003, 0.0668032510122), 0.00471083348185), + (SVector(0.0668032510122, 0.0095408154003), 0.00471083348185), + (SVector(0.9236559335870, 0.0668032510122), 0.00471083348185), + (SVector(0.0095408154003, 0.9236559335870), 0.00471083348185)) +## 0 negative weights, 0 points outside of the triangle, total sum of the +## weights is 0.5 + +## ----------------------------------------------------------------------------- +#*! Quadrature rule for an interpolation of order 12 on the triangle *# +#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *# + +const TRIANGLE_G12N33 = ((SVector(0.4882173897740, 0.4882173897740), 0.01286553322025), + (SVector(0.4882173897740, 0.0235652204524), 0.01286553322025), + (SVector(0.0235652204524, 0.4882173897740), 0.01286553322025), + (SVector(0.4397243922940, 0.4397243922940), 0.02184627226900), + (SVector(0.4397243922940, 0.1205512154110), 0.02184627226900), + (SVector(0.1205512154110, 0.4397243922940), 0.02184627226900), + (SVector(0.2712103850120, 0.2712103850120), 0.03142911210895), + (SVector(0.2712103850120, 0.4575792299760), 0.03142911210895), + (SVector(0.4575792299760, 0.2712103850120), 0.03142911210895), + (SVector(0.1275761455420, 0.1275761455420), 0.01739805646535), + (SVector(0.1275761455420, 0.7448477089170), 0.01739805646535), + (SVector(0.7448477089170, 0.1275761455420), 0.01739805646535), + (SVector(0.0213173504532, 0.0213173504532), 0.00308313052578), + (SVector(0.0213173504532, 0.9573652990940), 0.00308313052578), + (SVector(0.9573652990940, 0.0213173504532), 0.00308313052578), + (SVector(0.2757132696860, 0.6089432357800), 0.02018577888320), + (SVector(0.6089432357800, 0.1153434945350), 0.02018577888320), + (SVector(0.1153434945350, 0.2757132696860), 0.02018577888320), + (SVector(0.2757132696860, 0.1153434945350), 0.02018577888320), + (SVector(0.6089432357800, 0.2757132696860), 0.02018577888320), + (SVector(0.1153434945350, 0.6089432357800), 0.02018577888320), + (SVector(0.2813255809900, 0.6958360867880), 0.01117838660115), + (SVector(0.6958360867880, 0.0228383322223), 0.01117838660115), + (SVector(0.0228383322223, 0.2813255809900), 0.01117838660115), + (SVector(0.2813255809900, 0.0228383322223), 0.01117838660115), + (SVector(0.6958360867880, 0.2813255809900), 0.01117838660115), + (SVector(0.0228383322223, 0.6958360867880), 0.01117838660115), + (SVector(0.1162519159080, 0.8580140335440), 0.00865811555435), + (SVector(0.8580140335440, 0.0257340505483), 0.00865811555435), + (SVector(0.0257340505483, 0.1162519159080), 0.00865811555435), + (SVector(0.1162519159080, 0.0257340505483), 0.00865811555435), + (SVector(0.8580140335440, 0.1162519159080), 0.00865811555435), + (SVector(0.0257340505483, 0.8580140335440), 0.00865811555435)) +## 0 negative weights, 0 points outside of the triangle, total sum of the +## weights is 0.5 diff --git a/src/reference_integration.jl b/src/reference_integration.jl new file mode 100644 index 00000000..3e4c33f1 --- /dev/null +++ b/src/reference_integration.jl @@ -0,0 +1,232 @@ +""" + abstract type ReferenceQuadrature{D} + +A quadrature rule for integrating a function over the domain `D <: ReferenceShape`. + +Calling `x,w = q()` returns the nodes `x`, given as `SVector`s, and weights `w`, +for performing integration over `domain(q)`. +""" +abstract type ReferenceQuadrature{D} end + +""" + domain(q::ReferenceQuadrature) + +The domain of integratino for quadrature rule `q`. +""" +domain(q::ReferenceQuadrature{D}) where {D} = D() + +""" + qcoords(q) + +Return the coordinate of the quadrature nodes associated with `q`. +""" +qcoords(q::ReferenceQuadrature) = q()[1] + +""" + qweights(q) + +Return the quadrature weights associated with `q`. +""" +qweights(q::ReferenceQuadrature) = q()[2] + +function (q::ReferenceQuadrature)() + return interface_method(q) +end + +Base.length(q::ReferenceQuadrature) = length(qweights(q)) + +""" + integrate(f,q::ReferenceQuadrature) + integrate(f,x,w) + +Integrate the function `f` using the quadrature rule `q`. This is simply +`sum(f.(x) .* w)`, where `x` and `w` are the quadrature nodes and weights, +respectively. + +The function `f` should take an `SVector` as input. +""" +function integrate(f, q::ReferenceQuadrature) + x, w = q() + if domain(q) == ReferenceLine() + return integrate(x -> f(x[1]), x, w) + else + return integrate(f, x, w) + end +end + +function integrate(f, x, w) + sum(zip(x, w)) do (x, w) + return f(x) * prod(w) + end +end + +## Define some one-dimensional quadrature rules +""" + struct Fejer{N} + +`N`-point Fejer's first quadrature rule for integrating a function over `[0,1]`. +Exactly integrates all polynomials of degree `≤ N-1`. + +```jldoctest +using Inti + +q = Inti.Fejer(;order=10) + +Inti.integrate(cos,q) ≈ sin(1) - sin(0) + +# output + +true +``` + +""" +struct Fejer{N} <: ReferenceQuadrature{ReferenceLine} end + +Fejer(n::Int) = Fejer{n}() + +Fejer(; order::Int) = Fejer(order + 1) + +# N point fejer quadrature integrates all polynomials up to degree N-1 +""" + order(q::ReferenceQuadrature) + +A quadrature of order `p` integrates all polynomials of degree `≤ p`. +""" +order(::Fejer{N}) where {N} = N - 1 + +@generated function (q::Fejer{N})() where {N} + theta = [(2j - 1) * π / (2 * N) for j in 1:N] + x = -cos.(theta) + w = zero(x) + for j in 1:N + tmp = 0.0 + for l in 1:floor(N / 2) + tmp += 1 / (4 * l^2 - 1) * cos(2 * l * theta[j]) + end + w[j] = 2 / N * (1 - 2 * tmp) + end + xs = svector(i -> SVector(0.5 * (x[i] + 1)), N) + ws = svector(i -> w[i] / 2, N) + return xs, ws +end + +""" + struct Gauss{D,N} <: ReferenceQuadrature{D} + +Tabulated `N`-point symmetric Gauss quadrature rule for integration over `D`. +""" +struct Gauss{D,N} <: ReferenceQuadrature{D} + # gauss quadrature should be constructed using the order, and not the number + # of nodes. This ensures you don't instantiate quadratures which are not + # tabulated. + function Gauss(; domain, order) + domain == :triangle && (domain = ReferenceTriangle()) + domain == :tetrehedron && (domain = ReferenceTetrahedron()) + if domain isa ReferenceTriangle + msg = "quadrature of order $order not available for ReferenceTriangle" + haskey(TRIANGLE_GAUSS_ORDER_TO_NPTS, order) || error(msg) + n = TRIANGLE_GAUSS_ORDER_TO_NPTS[order] + elseif domain isa ReferenceTetrahedron + msg = "quadrature of order $order not available for ReferenceTetrahedron" + haskey(TETRAHEDRON_GAUSS_ORDER_TO_NPTS, order) || error(msg) + n = TETRAHEDRON_GAUSS_ORDER_TO_NPTS[order] + else + error("Tabulated Gauss quadratures only available for `ReferenceTriangle` or `ReferenceTetrahedron`") + end + return new{typeof(domain),n}() + end +end + +function order(q::Gauss{ReferenceTriangle,N}) where {N} + return TRIANGLE_GAUSS_NPTS_TO_ORDER[N] +end + +function order(q::Gauss{ReferenceTetrahedron,N}) where {N} + return TETRAHEDRON_GAUSS_NPTS_TO_ORDER[N] +end + +@generated function (q::Gauss{D,N})() where {D,N} + x, w = _get_gauss_qcoords_and_qweights(D, N) + return :($x, $w) +end + +""" + _get_gauss_and_qweights(R::Type{<:ReferenceShape{D}}, N) where D + +Returns the `N`-point symmetric gaussian qnodes and qweights `(x, w)` for integration over `R`. +""" +function _get_gauss_qcoords_and_qweights(R::Type{<:ReferenceShape}, N) + D = ambient_dimension(R()) + if !haskey(GAUSS_QRULES, R) || !haskey(GAUSS_QRULES[R], N) + error("quadrature rule not found") + end + qrule = GAUSS_QRULES[R][N] + @assert length(qrule) == N + # qnodes + qnodestype = SVector{N,SVector{D,Float64}} + x = qnodestype([q[1] for q in qrule]) + # qweights + qweightstype = SVector{N,Float64} + w = qweightstype([q[2] for q in qrule]) + return x, w +end + +""" + TensorProductQuadrature{N,Q} + +A tensor-product of one-dimension quadrature rules. Integrates over `[0,1]^N`. + +# Examples +```julia +qx = Fejer(10) +qy = TrapezoidalOpen(15) +q = TensorProductQuadrature(qx,qy) +``` +""" +struct TensorProductQuadrature{N,Q} <: ReferenceQuadrature{ReferenceHyperCube{N}} + quads1d::Q +end + +function TensorProductQuadrature(q...) + N = length(q) + Q = typeof(q) + return TensorProductQuadrature{N,Q}(q) +end + +function (q::TensorProductQuadrature{N})() where {N} + nodes1d = ntuple(N) do i + x1d, _ = q.quads1d[i]() + return map(x -> x[1], x1d) # convert the `SVector{1,T}` to just `T` + end + weights1d = map(q -> q()[2], q.quads1d) + nodes_iter = (SVector(x) for x in Iterators.product(nodes1d...)) + weights_iter = (prod(w) for w in Iterators.product(weights1d...)) + return nodes_iter, weights_iter +end + +""" + qrule_for_reference_shape(ref,order) + +Given a `ref`erence shape and a desired quadrature `order`, return +an appropiate quadrature rule. +""" +function qrule_for_reference_shape(ref, order) + if ref === ReferenceLine() || ref === :line + return Fejer(; order) + # return Fejer(; order) + elseif ref === ReferenceSquare() || ref === :square + qx = qrule_for_reference_shape(ReferenceLine(), order) + qy = qx + return TensorProductQuadrature(qx, qy) + elseif ref === ReferenceCube() || ref === :cube + qx = qrule_for_reference_shape(ReferenceLine(), order) + qy = qz = qx + return TensorProductQuadrature(qx, qy, qz) + elseif ref isa ReferenceTriangle || ref === :triangle + return Gauss(; domain=ref, order=order) + elseif ref isa ReferenceTetrahedron || ref === :tetrahedron + return Gauss(; domain=ref, order=order) + else + error("no appropriate quadrature rule found.") + end +end diff --git a/src/reference_interpolation.jl b/src/reference_interpolation.jl index 4f83adb1..685e213c 100644 --- a/src/reference_interpolation.jl +++ b/src/reference_interpolation.jl @@ -1,10 +1,10 @@ """ - abstract type AbstractElement{D,T} + abstract type ReferenceInterpolant{D,T} -Interpolanting function mapping points on the domain `D<:AbstractReferenceShape` +Interpolanting function mapping points on the domain `D<:ReferenceShape` (of singleton type) to a value of type `T`. -Instances `el` of `AbstractElement` are expected to implement: +Instances `el` of `ReferenceInterpolant` are expected to implement: - `el(x̂)`: evaluate the interpolation scheme at the (reference) coordinate `x̂ ∈ D`. - `jacobian(el,x̂)` : evaluate the jacobian matrix of the interpolation at the @@ -14,9 +14,9 @@ Instances `el` of `AbstractElement` are expected to implement: For performance reasons, both `el(x̂)` and `jacobian(el,x̂)` should take as input a `StaticVector` and output a static vector or static array. """ -abstract type AbstractElement{D<:AbstractReferenceShape,T} end +abstract type ReferenceInterpolant{D<:ReferenceShape,T} end -function (el::AbstractElement)(x) +function (el::ReferenceInterpolant)(x) return interface_method(el) end @@ -30,13 +30,13 @@ function jacobian(f, x) return interface_method(f) end -domain(::AbstractElement{D,T}) where {D,T} = D() -return_type(::AbstractElement{D,T}) where {D,T} = T -domain_dimension(t::AbstractElement) = domain(t) |> center |> length -range_dimension(el::AbstractElement{R,T}) where {R,T} = domain(t) |> center |> el |> length +domain(::ReferenceInterpolant{D,T}) where {D,T} = D() +return_type(::ReferenceInterpolant{D,T}) where {D,T} = T +domain_dimension(t::ReferenceInterpolant) = domain(t) |> center |> length +range_dimension(el::ReferenceInterpolant{R,T}) where {R,T} = domain(el) |> center |> el |> length """ - struct LagrangeElement{D,Np,T} <: AbstractElement{D,T} + struct LagrangeElement{D,Np,T} <: ReferenceInterpolant{D,T} A polynomial `p : D → T` uniquely defined by its `Np` values on the `Np` reference nodes of `D`. @@ -45,7 +45,7 @@ The return type `T` should be a vector space (i.e. support addition and multiplication by scalars). For istance, `T` could be a number or a vector, but not a `Tuple`. """ -struct LagrangeElement{D,Np,T} <: AbstractElement{D,T} +struct LagrangeElement{D,Np,T} <: ReferenceInterpolant{D,T} vals::SVector{Np,T} end @@ -233,8 +233,7 @@ end degree(el::LagrangeElement) degree(el::Type{<:LagrangeElement}) -The polynomial degree of the element. A `LagrangeElement` of degree `K` and -domain `D` belongs to the space [`PolynomialSpace{D,K}`](@ref). +The polynomial degree `el`. """ function degree(::Type{LagrangeElement{D,Np}})::Int where {D,Np} if D == ReferenceLine diff --git a/src/reference_shapes.jl b/src/reference_shapes.jl index 46074341..9af931f6 100644 --- a/src/reference_shapes.jl +++ b/src/reference_shapes.jl @@ -1,13 +1,13 @@ """ - abstract type AbstractReferenceShape + abstract type ReferenceShape A fixed reference domain/shape. Used mostly for defining more complex shapes as -transformations mapping an `AbstractReferenceShape` to some region of `ℜᴹ`. +transformations mapping an `ReferenceShape` to some region of `ℜᴹ`. See e.g. [`ReferenceLine`](@ref) or [`ReferenceTriangle`](@ref) for some examples of concrete subtypes. """ -abstract type AbstractReferenceShape end +abstract type ReferenceShape end """ struct ReferenceSimplex{N} @@ -15,7 +15,7 @@ abstract type AbstractReferenceShape end 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 +struct ReferenceSimplex{N} <: ReferenceShape end geometric_dimension(::ReferenceSimplex{N}) where {N} = N ambient_dimension(::ReferenceSimplex{N}) where {N} = N function Base.in(x, ::ReferenceSimplex{N}) where {N} @@ -43,12 +43,12 @@ Singleton type representing the tetrahedron with vertices const ReferenceTetrahedron = ReferenceSimplex{3} """ - struct ReferenceHyperCube{N} <: AbstractReferenceShape{N} + struct ReferenceHyperCube{N} <: ReferenceShape{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 +struct ReferenceHyperCube{N} <: ReferenceShape 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) diff --git a/src/utils.jl b/src/utils.jl index 33a64d9e..a86476b8 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -62,3 +62,17 @@ function return_type(f, args...) @debug "using `Base.promote_op` to infer return type. Consider defining `return_type(::typeof($f),args...)`." return Base.promote_op(f, args...) end + +""" + domain(f) + +Given a function-like object `f: Ω → R`, return `Ω`. +""" +function domain end + +""" + image(f) + +Given a function-like object `f: Ω → R`, return `f(Ω)`. +""" +function image end diff --git a/test/Manifest.toml b/test/Manifest.toml deleted file mode 100644 index a4cccccc..00000000 --- a/test/Manifest.toml +++ /dev/null @@ -1,155 +0,0 @@ -# This file is machine-generated - editing it directly is not advised - -julia_version = "1.9.3" -manifest_format = "2.0" -project_hash = "0be457628dfd5d423c82c1e12d8129fa1cf18c88" - -[[deps.Aqua]] -deps = ["Compat", "Pkg", "Test"] -git-tree-sha1 = "26c25461f946ef74512b3d6fd115d732912ba74a" -uuid = "4c88cf16-eb10-579e-8560-4a9242c79595" -version = "0.7.3" - -[[deps.ArgTools]] -uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" -version = "1.1.1" - -[[deps.Artifacts]] -uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" - -[[deps.Base64]] -uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" - -[[deps.Compat]] -deps = ["UUIDs"] -git-tree-sha1 = "8a62af3e248a8c4bad6b32cbbe663ae02275e32c" -uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.10.0" - - [deps.Compat.extensions] - CompatLinearAlgebraExt = "LinearAlgebra" - - [deps.Compat.weakdeps] - Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" - LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" - -[[deps.Dates]] -deps = ["Printf"] -uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" - -[[deps.Downloads]] -deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] -uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -version = "1.6.0" - -[[deps.FileWatching]] -uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" - -[[deps.InteractiveUtils]] -deps = ["Markdown"] -uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" - -[[deps.LibCURL]] -deps = ["LibCURL_jll", "MozillaCACerts_jll"] -uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" -version = "0.6.3" - -[[deps.LibCURL_jll]] -deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] -uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "7.84.0+0" - -[[deps.LibGit2]] -deps = ["Base64", "NetworkOptions", "Printf", "SHA"] -uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" - -[[deps.LibSSH2_jll]] -deps = ["Artifacts", "Libdl", "MbedTLS_jll"] -uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.10.2+0" - -[[deps.Libdl]] -uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" - -[[deps.Logging]] -uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" - -[[deps.Markdown]] -deps = ["Base64"] -uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" - -[[deps.MbedTLS_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.2+0" - -[[deps.MozillaCACerts_jll]] -uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2022.10.11" - -[[deps.NetworkOptions]] -uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" -version = "1.2.0" - -[[deps.Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] -uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.9.2" - -[[deps.Printf]] -deps = ["Unicode"] -uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" - -[[deps.REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] -uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" - -[[deps.Random]] -deps = ["SHA", "Serialization"] -uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" - -[[deps.SHA]] -uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" -version = "0.7.0" - -[[deps.Serialization]] -uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" - -[[deps.Sockets]] -uuid = "6462fe0b-24de-5631-8697-dd941f90decc" - -[[deps.TOML]] -deps = ["Dates"] -uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" -version = "1.0.3" - -[[deps.Tar]] -deps = ["ArgTools", "SHA"] -uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" -version = "1.10.0" - -[[deps.Test]] -deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] -uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[[deps.UUIDs]] -deps = ["Random", "SHA"] -uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" - -[[deps.Unicode]] -uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" - -[[deps.Zlib_jll]] -deps = ["Libdl"] -uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.13+0" - -[[deps.nghttp2_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.48.0+0" - -[[deps.p7zip_jll]] -deps = ["Artifacts", "Libdl"] -uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.4.0+0" diff --git a/test/Project.toml b/test/Project.toml index 5e3af3c7..8215136d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/reference_integration_test.jl b/test/reference_integration_test.jl new file mode 100644 index 00000000..d130ffca --- /dev/null +++ b/test/reference_integration_test.jl @@ -0,0 +1,68 @@ +using Test +using LinearAlgebra +using Inti + +@testset "Fejer quadrature" begin + N = 5 + q = Inti.Fejer{N}() + x, w = q() + D = Inti.domain(q) + @test D == Inti.ReferenceLine() + @test all(qnode ∈ D for qnode in x) + @test sum(w) ≈ 1 + # integrate all polynomial of degree N-1 exactly + for n in 1:(N - 1) + @test Inti.integrate(x -> x[1]^n, q) ≈ 1 / (n + 1) + end +end + +@testset "Gauss quad on triangle" begin + d = Inti.ReferenceTriangle() + # exact value for x^a*y^b integrate over reference triangle + exa = (a, b) -> factorial(a) * factorial(b) / factorial(a + b + 2) + # check all quadrature implemented + orders = keys(Inti.TRIANGLE_GAUSS_ORDER_TO_NPTS) + for p in orders + q = Inti.Gauss(; domain=d, order=p) + x, w = q() + @test Inti.domain(q) == d + @test all(qnode ∈ d for qnode in x) + for i in 0:p, j in 0:p + i + p > p && continue + @test Inti.integrate(x -> x[1]^i * x[2]^j, q) ≈ exa(i, j) + end + end +end + +@testset "Gauss quad on tetrahedron" begin + d = Inti.ReferenceTetrahedron() + # exact value for x^a*y^b*z^c integrate over reference tetrahedron + exa = (a, b, c) -> factorial(a) * factorial(b) * factorial(c) / factorial(a + b + c + 3) + # check all quadrature implemented + orders = keys(Inti.TETRAHEDRON_GAUSS_ORDER_TO_NPTS) + for p in orders + q = Inti.Gauss(; domain=d, order=p) + x, w = q() + @test Inti.domain(q) == d + @test all(qnode ∈ d for qnode in x) + for i in 0:p, j in 0:p, k in 0:p + i + j + k > p && continue + @test Inti.integrate(x -> x[1]^i * x[2]^j * x[3]^k, q) ≈ exa(i, j, k) + end + end +end + +@testset "Tensor product quad on square" begin + px = 10 + py = 12 + qx = Inti.Fejer(;order=px) + qy = Inti.Fejer(;order=py) + q = Inti.TensorProductQuadrature(qx, qy) + x, w = q() + D = Inti.domain(q) + @test D == Inti.ReferenceSquare() + @test all(qnode ∈ D for qnode in x) + for i in 0:px, j in 0:py + @test Inti.integrate(x -> x[1]^i * x[2]^j, q) ≈ 1 / (i + 1) * 1 / (j + 1) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index fc868f58..857dc067 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,3 +5,7 @@ using Aqua @safetestset "Code quality" include("aqua_test.jl") @safetestset "Reference shapes" include("reference_shapes_test.jl") + +@safetestset "Reference interpolation" include("reference_interpolation_test.jl") + +@safetestset "Reference integration" include("reference_integration_test.jl") From d715c4990433b610e0dee1345977d739ab986545 Mon Sep 17 00:00:00 2001 From: "Luiz M. Faria" Date: Mon, 2 Oct 2023 16:06:08 +0200 Subject: [PATCH 5/5] fix tests --- test/reference_interpolation_test.jl | 2 ++ test/runtests.jl | 2 ++ test/utils_test.jl | 10 ++++++++++ 3 files changed, 14 insertions(+) create mode 100644 test/utils_test.jl diff --git a/test/reference_interpolation_test.jl b/test/reference_interpolation_test.jl index 77764011..2c8a5124 100644 --- a/test/reference_interpolation_test.jl +++ b/test/reference_interpolation_test.jl @@ -10,6 +10,7 @@ using LinearAlgebra x̂ = Inti.reference_nodes(Inti.LagrangeLine{3}) vals = f.(x̂) p = Inti.LagrangeLine(vals) + @test Inti.return_type(p) == Float64 @test p(0) ≈ 0 @test p(1) ≈ 1 @test p(0.1) ≈ 0.1^2 @@ -35,6 +36,7 @@ using LinearAlgebra SVector(-1.0, 0) ) t = Inti.LagrangeTriangle(vtx) + @test Inti.return_type(t) == SVector{2,Float64} @test Inti.domain_dimension(t) == 2 @test Inti.range_dimension(t) == 2 # triangle in 3d diff --git a/test/runtests.jl b/test/runtests.jl index 857dc067..3e267fb1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,8 @@ using Aqua @safetestset "Code quality" include("aqua_test.jl") +@safetestset "Utility functions" include("utils_test.jl") + @safetestset "Reference shapes" include("reference_shapes_test.jl") @safetestset "Reference interpolation" include("reference_interpolation_test.jl") diff --git a/test/utils_test.jl b/test/utils_test.jl new file mode 100644 index 00000000..e74bcb0a --- /dev/null +++ b/test/utils_test.jl @@ -0,0 +1,10 @@ +using Inti +using Test + +f = (x,y) -> x*y +@test Inti.return_type(f,Int,Int) == Int +@test Inti.return_type(f,Int,Float64) == Float64 + +struct MyType end +@test_throws ErrorException Inti.interface_method(MyType) +@test_throws ErrorException Inti.interface_method(MyType())