diff --git a/src/CompScienceMeshes.jl b/src/CompScienceMeshes.jl index 6a8cd74..dd262a9 100644 --- a/src/CompScienceMeshes.jl +++ b/src/CompScienceMeshes.jl @@ -60,11 +60,12 @@ export dimension, universedimension export volume export neighborhood export quadpoints +export meshorder export isinside, isinclosure, overlap # specific to simplicial charts -export simplex, center, vertices +export simplex, center, vertices, nodes export faces, edges, circumcenter export barytocart, carttobary @@ -122,6 +123,10 @@ include("isinside.jl") include("findchart.jl") include("neighborhood.jl") include("subd_neighborhood.jl") + +include("meshes/curvilinear_mesh.jl") +include("meshes/curvilinear_chart.jl") +include("meshes/curvilinear_neighborhood.jl") include("quadpoints.jl") include("submesh.jl") diff --git a/src/charts.jl b/src/charts.jl index b74ead2..f6479f8 100644 --- a/src/charts.jl +++ b/src/charts.jl @@ -1,12 +1,13 @@ #export Simplex - +abstract type AbstractSimplex{U,D} +end # U: the dimension of the universe # D: the dimension of the manifold # N: the number of vertices # T: the type of the coordinates # C: the complimentary dimension (should always be U-D) -struct Simplex{U,D,C,N,T} +struct Simplex{U,D,C,N,T} <: AbstractSimplex{U,D} vertices::SVector{N,SVector{U,T}} tangents::SVector{D,SVector{U,T}} normals::SVector{C,SVector{U,T}} @@ -481,6 +482,8 @@ tangents(splx::Simplex, u) = hcat((splx.tangents)...) # vertices(splx::Simplex) = hcat((splx.vertices)...) function vertices(s::Simplex) s.vertices end +function nodes(s::Simplex) s.vertices end + """ verticeslist(simplex) diff --git a/src/fileio/gmsh.jl b/src/fileio/gmsh.jl index a620384..a183d3c 100644 --- a/src/fileio/gmsh.jl +++ b/src/fileio/gmsh.jl @@ -88,6 +88,8 @@ function load_gmsh_mesh(meshfile; # TODO: once we support hexahedrons, we need to catch that one too if element == :quadrangle return QuadMesh(vertices, elements) + elseif element == :line && order > 1 + return CurvilinearMesh(vertices, elements, order) else return Mesh(vertices, elements) end diff --git a/src/mesh.jl b/src/mesh.jl index 6d89e2e..5646a9d 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -522,7 +522,7 @@ function skeleton(mesh, dim::Int; sort=:spacefillingcurve) meshdim = dimension(mesh) @assert 0 <= dim <= meshdim - if dim == meshdim + if dim == meshdim && (mesh isa Mesh || mesh isa QuadMesh) return mesh end diff --git a/src/meshes/curvilinear_chart.jl b/src/meshes/curvilinear_chart.jl new file mode 100644 index 0000000..148ba14 --- /dev/null +++ b/src/meshes/curvilinear_chart.jl @@ -0,0 +1,307 @@ +#TODO: Extend to 2D and 3D simplices. + +""" + CurvilinearSimplex{U,D,C,N,T} + +A curvilinear simplex (currently specialized to 1D) embedded in `U`-dimensional +space, represented by `N` interpolation vertices and `N` parametric nodes. + +# Fields +- `vertices :: SVector{N, SVector{U,T}}` + Coordinates of the interpolation/control points in physical space. + +- `ζnodes :: SVector{N,T}` + Parametric locations of the nodes on the reference interval `[0,1]` used + for high-order interpolation (e.g. Gmsh’s high-order node ordering). + +# Description +A `CurvilinearSimplex` represents a single high-order 1D finite-/boundary-element. +Its geometry is given by an isoparametric map + +x(ζ) = Σ_r ℓ_r(ζ) * vertices[r], + +where `ℓ_r(·)` are the Lagrange polynomials defined by `ζnodes`. + +The type parameters are: +- `U`: ambient (universe) dimension, +- `D`: topological dimension (here always `1`), +- `C`: codimension (= `U-D`), +- `N`: number of local nodes, +- `T`: coordinate type (e.g. `Float64`). + +This is the geometric building block used for curvilinear 1D meshes. +""" +struct CurvilinearSimplex{U,D,C,N,T} <: AbstractSimplex{U,D} + vertices::SVector{N,SVector{U,T}} + ζnodes::SVector{N,T} +end + +""" + nodes(s::CurvilinearSimplex{U,1}) + +Return the physical coordinates of the two endpoints of a 1D curvilinear element, +identified by the parametric values `ζ = 1` (left) and `ζ = 0` (right) under the +barycentric convention `ζ = λ_left`. + +# Returns +`(x_left, x_right) :: (SVector{U,T}, SVector{U,T})` +""" +@inline function nodes(s::CurvilinearSimplex{U,1}) where {U} + return SVector(s.vertices[begin], s.vertices[begin + 1]) +end + +# convenience +vertices(s::CurvilinearSimplex) = s.vertices + +""" + coordtype(s::CurvilinearSimplex) + coordtype(::Type{CurvilinearSimplex}) + +Return the scalar coordinate type `T` used to represent vertex coordinates. +""" +coordtype(::CurvilinearSimplex{U,D,C,N,T}) where {U,D,C,N,T} = T + +""" + dimension(s::CurvilinearSimplex) + dimension(::Type{CurvilinearSimplex}) + +Return the topological dimension of the simplex (`1` for line elements). +""" +dimension(::CurvilinearSimplex{U,D,C,N,T}) where {U,D,C,N,T} = D + +""" + universedimension(s::CurvilinearSimplex) + universedimension(::Type{CurvilinearSimplex}) + +Return the dimension of the embedding space (e.g. `2` for curves in the plane). +""" +universedimension(::CurvilinearSimplex{U,D,C,N,T}) where {U,D,C,N,T} = U + +""" + vertextype(s::CurvilinearSimplex) + vertextype(::Type{CurvilinearSimplex}) + +Return the static vector type used to represent vertex coordinates. +""" +vertextype(::CurvilinearSimplex{U,D,C,N,T}) where {U,D,C,N,T} = SVector{U,T} + +# type-form +coordtype(::Type{CurvilinearSimplex{U,D,C,N,T}}) where {U,D,C,N,T} = T +dimension(::Type{CurvilinearSimplex{U,D,C,N,T}}) where {U,D,C,N,T} = D +universedimension(::Type{CurvilinearSimplex{U,D,C,N,T}}) where {U,D,C,N,T} = U +vertextype(::Type{CurvilinearSimplex{U,D,C,N,T}}) where {U,D,C,N,T} = SVector{U,T} + +""" + neighborhood(ch::CurvilinearSimplex, bary) + +Construct a `MeshPoint` located at the parametric (barycentric) coordinate `bary` +on the given curvilinear simplex. + +The returned object stores: +- the chart `ch`, +- the parametric point `bary`, +- the corresponding physical coordinates `cartesian`. + +This is the canonical way to create a point for evaluating reference functions +or geometric quantities on curvilinear elements. +""" +function neighborhood(p::C, bary) where {C<:CurvilinearSimplex} + D = dimension(p) + T = coordtype(p) + P = SVector{D,T} + cart = barytocart(p, T.(bary)) + Q = typeof(cart) + U = length(cart) + MeshPointNM{T,C,D,U}(p, P(bary), cart) +end + +barytocart(ch::CurvilinearSimplex{U,1}, u) where {U} = pushforward(ch, u[1]) + +cartesian(ch::AbstractSimplex{U,1}, u::SVector{1,<:Real}) where {U} = barytocart(ch, u) +parametric(ch::AbstractSimplex{U,1}, u::SVector{1,<:Real}) where {U} = u + +""" + pushforward(ch::CurvilinearSimplex{U,1,C,N,T}, ζ) + +Evaluate the isoparametric mapping from the scalar barycentric coordinate +`ζ = λ_left ∈ [0,1]` to physical space: + +x(ζ) = Σ_r ℓ_r(ζ) * ch.vertices[r], + +where `ℓ_r` are the Lagrange polynomials associated with `ch.ζnodes`. + +With this convention, `ζ=1` maps to the left endpoint and `ζ=0` to the right. +""" +@inline function pushforward(ch::CurvilinearSimplex{U,1,C,N,T}, ζ::Real) where {U,C,N,T} + z = T(ζ) + s = zero(SVector{U,T}) + @inbounds for r in 1:N + s += ch.vertices[r] * _ℓ(ch.ζnodes, r, z) + end + s +end + +""" + dpushforward(ch::CurvilinearSimplex, ζ) + +Return the derivative of the isoparametric map `x(ζ)` with respect to `ζ`. + +Mathematically: + +dx/dζ = Σ_r ℓ'_r(ζ) * ch.vertices[r]. + + +This quantity is used for computing Jacobians, tangent vectors, and element +volumes on curvilinear 1D elements. +""" +@inline function dpushforward(ch::CurvilinearSimplex{U,1,C,N,T}, ζ::Real) where {U,C,N,T} + z = T(ζ) + s = zero(SVector{U,T}) + @inbounds for r in 1:N + s += ch.vertices[r] * _dℓ(ch.ζnodes, r, z) + end + s +end + +""" + _ℓ(ζnodes, r, ζ) + +Evaluate the `r`-th 1D Lagrange basis polynomial `ℓ_r(ζ)` associated with the +nodes given in `ζnodes`. + +The polynomial satisfies: +- `ℓ_r(ζnodes[r]) = 1` +- `ℓ_r(ζnodes[j]) = 0` for all `j ≠ r` + +This is the standard explicit Lagrange form: + +ℓ_r(ζ) = ∏_{j ≠ r} (ζ - ζ_j) / (ζ_r - ζ_j). + +# Internal +Used for geometric mapping, Jacobians, and interpolation on curvilinear 1D elements. +""" +@inline function _ℓ(ζnodes::SVector{D,T}, r::Int, ζ::T) where {D,T<:Real} + ζr = ζnodes[r] + num = one(T); den = one(T) + @inbounds for j in 1:D + j == r && continue + num *= (ζ - ζnodes[j]) + den *= (ζr - ζnodes[j]) + end + + return num/den +end + +""" + _dℓ(ζnodes, r, ζ) + +Evaluate the derivative of the `r`-th Lagrange basis polynomial at ζ. +This is the analytical derivative of `_ℓ`. + +Mathematically: + +ℓ'r(ζ) = Σ{k ≠ r} +∏{j ≠ r,k} (ζ - ζ_j) +/ (ζ_r - ζ_k) / ∏{j ≠ r} (ζ_r - ζ_j). + + +This is used to compute the derivative of the physical mapping and hence the +Jacobian and tangents. +""" +@inline function _dℓ(ζnodes::SVector{D,T}, r::Int, ζ::T) where {D,T<:Real} + ζr = ζnodes[r] + s = zero(T) + @inbounds for k in 1:D + k == r && continue + num = one(T); den = (ζr - ζnodes[k]) + @inbounds for j in 1:D + (j == r || j == k) && continue + num *= (ζ - ζnodes[j]) + den *= (ζr - ζnodes[j]) + end + s += num/den + end + s +end + + +""" + gmsh_line_ζnodes(::Val{N}) + +Return the parametric node locations `ζ ∈ [0,1]` for a 1D high-order line element, +arranged to match the element’s local node numbering (compatible with Gmsh element +node indices and the barycentric convention used here). + +By convention in this code, the scalar parameter is the *left* barycentric +coordinate: `ζ = 1` corresponds to the left endpoint and `ζ = 0` to the right +endpoint. + +- For `N = 2` (linear): returns `(1.0, 0.0)` → `(left, right)`. +- For `N > 2`: returns `(1.0, 0.0, ζ₃, …, ζ_N)` where the interior nodes `ζ_k` + lie strictly between 0 and 1. (Their order is chosen to be consistent with + Gmsh’s line-element node numbering and this code’s interpolation routines.) + +This tuple is used by the Lagrange map `x(ζ) = Σ_r ℓ_r(ζ) * vertices[r]`. +""" +gmsh_line_ζnodes(::Val{N}) where {N} = + N == 2 ? SVector{2,Float64}(1.0, 0.0) : + SVector{N,Float64}(1.0, 0.0, (k/(N-1) for k in reverse(1:N-2))...) + + +function chart( + m::CompScienceMeshes.CurvilinearMesh{U,N,T,O}, + faceid::Int +) where {U,N,T,O} + + vindices = m.faces[faceid] + X = m.vertices[vindices] + ζnodes = gmsh_line_ζnodes(Val(N)) # match that order + + CompScienceMeshes.CurvilinearSimplex{U,1,U-1,N,T}(X, ζnodes) +end + +""" + jacobian(ch, ζ) + +Return the scalar Jacobian `|dx/dζ|` of the 1D isoparametric mapping at +parametric coordinate `ζ`. +""" +jacobian(ch::CurvilinearSimplex{U,1}, ζ::Real) where {U} = norm(dpushforward(ch, ζ)) + +""" + center(ch::CurvilinearSimplex) + +Return the physical midpoint of the 1D curvilinear element, obtained by +evaluating the mapping at `ζ = 0.5`. +""" +center(ch::CurvilinearSimplex) = pushforward(ch, 0.5) + +domain(::CurvilinearSimplex{U,1,C,N,T}) where {U,C,N,T} = ReferenceSimplex{1,T,2}() + +""" + volume(ch) + +Compute the physical length of the curvilinear element `ch` by integrating +the Jacobian over `[0,1]`. + +A 3-point Gauss–Legendre quadrature rule is used: + +∫₀¹ |dx/dζ| dζ +≈ Σ_k w_k * |dx/dζ(ζ_k)| / 2 + + +This is the element volume used by numerical integration for BEM/FEM operators. +""" +function volume(ch::CurvilinearSimplex{U,1,C,N,T}) where {U,C,N,T} + + ξ, w = legendre(20, 0.0, 1.0) + #ξ = (-sqrt(T(3)/T(5)), zero(T), sqrt(T(3)/T(5))) + #w = (T(5)/T(9), T(8)/T(9), T(5)/T(9)) + #s = zero(T) + #@inbounds for k in 1:3 + # ζ = (ξ[k] + one(T)) / T(2) + # s += w[k] * jacobian(ch, ζ) / T(2) + #end + #s + return dot(w, jacobian.(Ref(ch), ξ)) +end \ No newline at end of file diff --git a/src/meshes/curvilinear_mesh.jl b/src/meshes/curvilinear_mesh.jl new file mode 100644 index 0000000..e2148e2 --- /dev/null +++ b/src/meshes/curvilinear_mesh.jl @@ -0,0 +1,70 @@ +# Mesh: 1D elements embedded in R^U, with N control points per element (N = order+1) +mutable struct CurvilinearMesh{U,N,T,O} <: AbstractMesh{U,N,T} + vertices::Vector{SVector{U,T}} # coordinates in R^U + faces::Vector{SVector{N,Int}} # connectivity (local order = Gmsh order) + + function CurvilinearMesh( + vertices::Vector{SVector{U,T}}, + faces::Vector{SVector{N,Int}}, + order::Integer + ) where {U,N,T} + order > 0 || throw(ArgumentError("order must be positive.")) + new{U,N,T,Int(order)}(vertices, faces) + end +end + + + +# Constructors from abstract containers +CurvilinearMesh(vertices::AbstractVector{<:SVector{U,T}}, + faces::AbstractVector{<:SVector{N,Int}}, + order::Integer) where {U,N,T} = + CurvilinearMesh(Vector{SVector{U,T}}(vertices), + Vector{SVector{N,Int}}(faces), + order) + +# From matrices (verts: U×nV, faces: N×nF) +function CurvilinearMesh(verts::AbstractMatrix{T}, + faces::AbstractMatrix{<:Integer}, + order::Integer) where {T} + U = size(verts, 1) + N = size(faces, 1) + v = [SVector{U,T}(verts[:, i]) for i in 1:size(verts, 2)] + f = [SVector{N,Int}(faces[:, i]) for i in 1:size(faces, 2)] + CurvilinearMesh(v, f, order) +end + +# Mesh interface +coordtype(::CurvilinearMesh{U,N,T,O}) where {U,N,T,O} = T +vertextype(::CurvilinearMesh{U,N,T,O}) where {U,N,T,O} = SVector{U,T} +universedimension(::CurvilinearMesh{U}) where {U} = U +dimension(::CurvilinearMesh) = 1 + +function celltype(m::CurvilinearMesh{U,N,T,O}) where {U,N,T,O} SimplexGraph{N-O+1} end +function celltype(m::CurvilinearMesh{U,N,T,O}, ::Type{Val{M}}) where {U,N,T,O,M} SimplexGraph{M+1} end + +function indextype(m::CurvilinearMesh{U,N}) where {U,N} SVector{N-O+1,Int} end +function indextype(m::CurvilinearMesh{U,N}, ::Type{Val{M}}) where {U,N,M} SVector{M+1,Int} end + +function indices(m::CurvilinearMesh{U,N,T,O}, cell) where {U,N,T,O} + # Currently, we only support lines + # First come the topological nodes + return SVector(m.faces[cell][begin], m.faces[cell][begin+1]) +end + + +vertices(m::CurvilinearMesh) = m.vertices +faces(m::CurvilinearMesh) = m.faces + +numvertices(m::CurvilinearMesh) = length(m.vertices) +numcells(m::CurvilinearMesh) = length(m.faces) +cells(m::CurvilinearMesh) = Base.OneTo(length(m.faces)) +cell(m::CurvilinearMesh, i::Int) = m.faces[i] + +Base.eltype(::Type{CurvilinearMesh{U,N,T,O}}) where {U,N,T,O} = SVector{N,Int} + +meshorder(::CurvilinearMesh{U,N,T,O}) where {U,N,T,O} = O + +function skeleton_fast(mesh::CurvilinearMesh, dim::Int) + skeleton_fast(mesh, Val{dim}) +end diff --git a/src/meshes/curvilinear_neighborhood.jl b/src/meshes/curvilinear_neighborhood.jl new file mode 100644 index 0000000..0047cfb --- /dev/null +++ b/src/meshes/curvilinear_neighborhood.jl @@ -0,0 +1,19 @@ +""" + tangents(neighborhod, i) -> tangent_i + +Return the i-th tangent vector at the neighborhood. +""" +function tangents(mp::MeshPointNM{T,C,D,U}, i) where {T,C<:CurvilinearSimplex,D,U} + return t = dpushforward(mp.patch, mp.bary[1]) +end + +function normal(mp::MeshPointNM{T,C,D,U}) where {T,C<:CurvilinearSimplex,D,U} + t = tangents(mp, 1) + nt = SVector(-t[2], t[1]) + + return nt ./ norm(nt) +end + +function jacobian(mp::MeshPointNM{T,C,D,U}) where {T,C<:CurvilinearSimplex,D,U} + return jacobian(mp.patch, mp.bary[1]) +end \ No newline at end of file diff --git a/src/primitives/linemeshes/linemeshes.jl b/src/primitives/linemeshes/linemeshes.jl index 3d2372f..184c50f 100644 --- a/src/primitives/linemeshes/linemeshes.jl +++ b/src/primitives/linemeshes/linemeshes.jl @@ -57,3 +57,62 @@ CT = SVector{2,Int} return Mesh(vertices, faces) end + +function gmshcircle(radius::T, delta::T, tempname=tempname(); order::Int=1) where T<:Real + s = """ +lc = $delta; + +Point(1) = {0, 0, 0, lc}; +Point(2) = {$radius, 0, 0, lc}; +Point(3) = {0, $radius, 0, lc}; +Point(4) = {-$radius, 0, 0, lc}; +Point(5) = {0, -$radius, 0, lc}; + +Circle(1) = {2, 1, 3}; +Circle(2) = {3, 1, 4}; +Circle(3) = {4, 1, 5}; +Circle(4) = {5, 1, 2}; + +Physical Curve("boundary") = {1, 2, 3, 4}; +""" + + fn = tempname + io = open(fn, "w") + try + print(io, s) + finally + close(io) + end + + fno = tempname * ".msh" + + gmsh.initialize() + gmsh.option.setNumber("Mesh.MshFileVersion", 2.0) + gmsh.open(fn) + + # Choose element order (1: linear segments, 2: quadratic) + @show order + gmsh.option.setNumber("Mesh.ElementOrder", order) + + # Generate only the 1D mesh (curves) + gmsh.model.mesh.generate(1) + + # If you *don’t* define Physicals, enable the next line instead: + # gmsh.option.setNumber("Mesh.SaveAll", 1) + + gmsh.write(fno) + gmsh.finalize() + + m = load_gmsh_mesh(fno, + element=:line, + udim=2, + vertextype=Float64, + order=order, + sort=false + ) + + #rm(fno) + rm(fn) + + return m +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 886acc6..89a148d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -54,6 +54,8 @@ include("primitives/surfacemeshes/test_rectangle.jl") include("primitives/surfacemeshes/test_sphere.jl") include("primitives/volumemeshes/test_cuboid.jl") +include("test_curvilinear_elements.jl") + using TestItemRunner @run_package_tests end diff --git a/test/test_curvilinear_elements.jl b/test/test_curvilinear_elements.jl new file mode 100644 index 0000000..d49c9f2 --- /dev/null +++ b/test/test_curvilinear_elements.jl @@ -0,0 +1,79 @@ +using Test +using CompScienceMeshes +using StaticArrays + +## Verify consistency with linear mesh + +v1 = SVector(0.0, 0.0) +v2 = SVector(1.0, 0.0) +v3 = SVector(0.25, 0.0) +v4 = SVector(0.50, 0.0) +v5 = SVector(0.75, 0.0) + +verts1 = [v1, v2] +faces1 = [SVector(1,2)] + +verts2 = [v1, v2, v3, v4, v5] +faces2 = [SVector(1,2,3,4,5)] + +m1 = CompScienceMeshes.Mesh(verts1, faces1) +m2 = CompScienceMeshes.CurvilinearMesh(verts2, faces2, 4) + +ch1 = chart(m1, 1) +ch2 = chart(m2, 1) + +mp1 = neighborhood(ch1, 0.0) +mp2 = neighborhood(ch2, 0.0) + +@test tangents(mp1, 1) == SVector(-1.0, 0.0) +@test tangents(mp2, 1) == SVector(-1.0, 0.0) # Sign is flipped + +@test normal(mp1) == SVector(0.0, -1.0) +@test normal(mp2) == SVector(0.0, -1.0) # Sign is flipped + +@test cartesian(mp1) == SVector(1.0, 0.0) +@test cartesian(mp2) == SVector(1.0, 0.0) # We are using barycentric coordinates: 0.0 refers to the end vertex + +@test nodes(ch1)[1] == SVector(0.0, 0.0) +@test nodes(ch1)[2] == SVector(1.0, 0.0) + +@test nodes(ch2)[1] == SVector(0.0, 0.0) +@test nodes(ch2)[2] == SVector(1.0, 0.0) + +@test volume(ch2) ≈ 1.0 + + +## + +# Consider parabola t ↦(t, t²) +v1 = SVector(0.0, 0.0^2) +v2 = SVector(1.0, 1.0^2) +v3 = SVector(0.25, 0.25^2) +v4 = SVector(0.50, 0.50^2) +v5 = SVector(0.75, 0.75^2) + +verts_parabola = [v1, v2, v3, v4, v5] +face_parabola = [SVector(1,2,3,4,5)] + +m_parabola = CompScienceMeshes.CurvilinearMesh(verts_parabola, face_parabola, 4) +ch_parabola = chart(m_parabola, 1) + +mp_1 = neighborhood(ch_parabola, 0.0) +mp_2 = neighborhood(ch_parabola, 0.25) +mp_3 = neighborhood(ch_parabola, 0.50) +mp_4 = neighborhood(ch_parabola, 0.75) +mp_5 = neighborhood(ch_parabola, 1.0) + +@test cartesian(mp_1) == SVector(1.00, 1.0) +@test cartesian(mp_2) == SVector(0.75, 0.5625) +@test cartesian(mp_3) == SVector(0.50, 0.25) +@test cartesian(mp_4) == SVector(0.25, 0.0625) +@test cartesian(mp_5) == SVector(0.00, 0.0) + +@test tangents(mp_1, 1) == SVector(-1.0, -2.0) # Sign is flipped +@test tangents(mp_2, 1) ≈ SVector(-1.0, -1.5) # Sign is flipped +@test tangents(mp_3, 1) ≈ SVector(-1.0, -1.0) # Sign is flipped +@test tangents(mp_4, 1) ≈ SVector(-1.0, -0.5) # Sign is flipped +@test tangents(mp_5, 1) ≈ SVector(-1.0, -0.0) # Sign is flipped + +@test volume(ch_parabola) ≈ 1.478942857544597 # (Reference from Mathematica 14) \ No newline at end of file