diff --git a/src/CompScienceMeshes.jl b/src/CompScienceMeshes.jl
index 6a8cd74..8c959a8 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,14 @@ include("isinside.jl")
 include("findchart.jl")
 include("neighborhood.jl")
 include("subd_neighborhood.jl")
+
+include("meshes/curvilinear_silvester.jl")
+include("meshes/curvilinear_gmsh_line.jl")
+include("meshes/curvilinear_gmsh_triangle.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..201e3b0 100644
--- a/src/fileio/gmsh.jl
+++ b/src/fileio/gmsh.jl
@@ -88,6 +88,10 @@ 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, 1, order)
+    elseif element == :triangle && order > 1
+        return CurvilinearMesh(vertices, elements, 2, 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..c91472f
--- /dev/null
+++ b/src/meshes/curvilinear_chart.jl
@@ -0,0 +1,423 @@
+#### Implementation background, todos etc
+# We implement the standard Lagrange interpolated curvilinear elements
+# For the arbitrary order code, we use the Silvester polynomials as described in
+# [1] R. D. Graglia and A. F. Peterson, Higher-order techniques in computational electromagnetics.
+# TODO:
+# 1) Extend to 3D simplices.
+# 2) Add as many specializations as needed
+# 3) Rethink what MeshPointNM should contain as information
+
+"""
+    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.
+
+The type parameters are:
+- `U`: ambient (universe) dimension,
+- `D`: topological dimension,
+- `C`: codimension (= `U-D`),
+- `N`: number of local nodes,
+- `T`: coordinate type (e.g. `Float64`).
+- `O`: element order.
+
+This is the geometric building block used for curvilinear 1D meshes.
+"""
+struct CurvilinearSimplex{U,D,C,N,T,O} <: AbstractSimplex{U,D}
+    vertices::SVector{N,SVector{U,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})`
+"""
+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::Real) where {U} = barytocart(ch, SVector{1}(u))
+barytocart(ch::CurvilinearSimplex{U,2}, u::Real, v::Real) where {U} = barytocart(ch, SVector{2}(u, v))
+barytocart(ch::CurvilinearSimplex{U,3}, u::Real, v::Real, w::Real) where {U} = barytocart(ch, SVector{3}(u, v, w))
+
+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
+
+function barytocart(ch::CurvilinearSimplex{U,1,C,N,T,O}, bary::AbstractVector) where {U,C,N,T,O}
+    u = bary[1]
+
+    s = zero(SVector{U,T})
+    @inbounds for (r, (i, j)) in enumerate(gmsh_line_index_to_tuple[O])
+        s += ch.vertices[r] * silvester(O, u)[i + 1] * silvester(O, 1 - u)[j + 1]
+    end
+
+    return s
+end
+
+# barytocart specializations for 1D curvilinear elements using Silvester interpolation
+# Orders 1 through 10, expressions simplified and verified such that u=1 yields v1
+
+function barytocart(ch::CurvilinearSimplex{U,1,C,N,T,1}, bary::AbstractVector) where {U,C,N,T}
+    u = bary[1]
+    return ch.vertices[1] * u + ch.vertices[2] * (1 - u)
+end
+
+function barytocart(ch::CurvilinearSimplex{U,1,C,N,T,2}, bary::AbstractVector) where {U,C,N,T}
+    u = bary[1]
+    u2 = u*u
+
+    return ch.vertices[1] * (2u2 - u) +
+           ch.vertices[2] * (2u2 - 3u + 1) +
+           ch.vertices[3] * (-4u2 + 4u)
+end
+
+function barytocart(ch::CurvilinearSimplex{U,2,C,N,T,O}, bary) where {U,C,N,T,O}
+    u = bary[1]
+    v = bary[2]
+
+    s = zero(SVector{U,T})
+    @inbounds for (r, (i, j, k)) in enumerate(gmsh_triangle_index_to_tuple[O])
+        s += ch.vertices[r] * 
+            silvester(O, u)[i + 1] *
+            silvester(O, v)[j + 1] *
+            silvester(O, 1 - u - v)[k + 1]
+    end
+
+    return s
+end
+
+function barytocart(ch::CurvilinearSimplex{U,2,C,N,T,1}, bary) where {U,C,N,T}
+    u, v = bary
+    w = 1 - u - v
+
+    return ch.vertices[1] * u +
+           ch.vertices[2] * v +
+           ch.vertices[3] * w
+end
+
+#=
+function barytocart(ch::CurvilinearSimplex{U,2,C,N,T,2}, bary) where {U,C,N,T}
+    u, v = bary
+    w = 1 - u - v
+
+    u2 = u * u
+    v2 = v * v
+    w2 = w * w
+
+    return ch.vertices[1] * (2u2 - u) +
+           ch.vertices[2] * (2v2 - v) +
+           ch.vertices[3] * (2w2 - w) +
+           ch.vertices[4] * (4u*v) +
+           ch.vertices[5] * (4v*w) +
+           ch.vertices[6] * (4u*w)
+end
+
+function barytocart(ch::CurvilinearSimplex{U,2,C,N,T,3}, bary) where {U,C,N,T}
+    u, v = bary
+    w = 1 - u - v
+
+    u2 = u * u
+    u3 = u2 * u
+    v2 = v * v
+    v3 = v2 * v
+    w2 = w * w
+    w3 = w2 * w
+
+    return ch.vertices[1] * (10u3 - 15u2 + 6u) +
+           ch.vertices[2] * (10v3 - 15v2 + 6v) +
+           ch.vertices[3] * (10w3 - 15w2 + 6w) +
+           ch.vertices[4] * (60u2*v - 60u*v) +
+           ch.vertices[5] * (60u*v2 - 60u*v) +
+           ch.vertices[6] * (60v2*w - 60v*w) +
+           ch.vertices[7] * (60v*w2 - 60v*w) +
+           ch.vertices[8] * (60u*w2 - 60u*w) +
+           ch.vertices[9] * (60u2*w - 60u*w) +
+           ch.vertices[10] * (60u*v*w)
+end
+=#
+
+"""
+    unitarybase(ch::CurvilinearSimplex, bary)
+
+Compute the unitary base vectors as defined in (3.38) in Graglia & Peterson.
+"""
+function unitarybase(ch::CurvilinearSimplex{U,1,C,N,T,O}, bary) where {U,C,N,T,O}
+    u = bary[1]
+
+    s = zero(SVector{U,T})
+
+    Ru, dRu = fast_silvester(O, u)
+    Rv, dRv = fast_silvester(O, 1 - u)
+
+    @inbounds for (r, (i, j)) in enumerate(gmsh_line_index_to_tuple[O])
+
+        # Adapted (2.44) in Graglia & Peterson with ξ₁ = u and ξ₂ = 1 - u
+        s += ch.vertices[r] * (dRu[i + 1] * Rv[j + 1] - Ru[i + 1] * dRv[j + 1])
+    end
+
+    return SVector{1}(s)
+end
+
+function unitarybase(ch::CurvilinearSimplex{U,1,C,N,T,1}, bary) where {U,C,N,T}
+    u = bary[1]
+
+    # Adapted (2.44) in Graglia & Peterson with ξ₁ = u and ξ₂ = 1 - u
+    s = ch.vertices[1] - ch.vertices[2]
+
+    return SVector{1}(s)
+end
+
+function unitarybase2(ch::CurvilinearSimplex{U,1,C,N,T,2}, bary) where {U,C,N,T}
+    u = bary[1]
+
+    # Adapted (2.44) in Graglia & Peterson with ξ₁ = u and ξ₂ = 1 - u
+    s = ch.vertices[1] * (+4u - 1) +
+        ch.vertices[2] * (+4u - 3) +
+        ch.vertices[3] * (-8u + 4)
+
+    return SVector{1}(s)
+end
+
+# TODO: Speed-up gains possible by specialization to explicit polynomials
+# do at least for lowest orders to avoid memory allocations
+function unitarybase(ch::CurvilinearSimplex{U,2,C,N,T,O}, bary) where {U,C,N,T,O}
+    u = bary[1]
+    v = bary[2]
+
+    Ru, dRu = fast_silvester(O, u)
+    Rv, dRv = fast_silvester(O, v)
+    Rw, dRw = fast_silvester(O, 1 - u - v)
+
+    𝓵¹ = zero(SVector{U,T})
+    𝓵² = zero(SVector{U,T})
+    @inbounds for (r, (i, j, k)) in enumerate(gmsh_triangle_index_to_tuple[O])
+
+        𝓵¹ += ch.vertices[r] * (
+                dRu[i + 1] * Rv[j + 1] * Rw[k + 1] -
+                Ru[i + 1] * Rv[j + 1] * dRw[k + 1]
+            )
+
+        𝓵² += ch.vertices[r] * (
+                Ru[i + 1] * dRv[j + 1] * Rw[k + 1] -
+                Ru[i + 1] * Rv[j + 1] * dRw[k + 1]
+            )
+    end
+    
+    return SVector(𝓵¹, 𝓵²)
+end
+
+# TODO: Trying to avoid allocations. Even for lowest order however this is much slower. Why???
+function unitarybase_experimental(ch::CurvilinearSimplex{U,2,C,N,T,1}, bary) where {U,C,N,T}
+    u = bary[1]
+    v = bary[2]
+
+    Ru = @SVector[silv(Val(1), Val(i), u) for i=0:1]
+    Rv = @SVector[silv(Val(1), Val(i), v) for i=0:1]
+    Rw = @SVector[silv(Val(1), Val(i), 1 - u - v) for i=0:1]
+
+    dRu = @SVector[dsilv(Val(1), Val(i), u) for i=0:1]
+    dRv = @SVector[dsilv(Val(1), Val(i), v) for i=0:1]
+    dRw = @SVector[dsilv(Val(1), Val(i), 1 - u - v) for i=0:1]
+
+    𝓵¹ = zero(SVector{U,T})
+    𝓵² = zero(SVector{U,T})
+    @inbounds for (r, (i, j, k)) in enumerate(gmsh_triangle_index_to_tuple[1])
+
+        𝓵¹ += ch.vertices[r] * (
+                dRu[i + 1] * Rv[j+1] * Rw[k+1] -
+                Ru[i + 1] * Rv[j+1] * dRw[k+1]
+            )
+
+        𝓵² += ch.vertices[r] * (
+                Ru[i + 1] * dRv[j+1] * Rw[k+1] -
+                Ru[i + 1] * Rv[j+1] * dRw[k+1]
+            )
+    end
+    
+    return SVector(𝓵¹, 𝓵²)
+end
+
+
+function unitarybase2(ch::CurvilinearSimplex{U,2,C,N,T,7}, bary) where {U,C,N,T}
+    u = bary[1]
+    v = bary[2]
+
+    Ru = @SVector[silv(Val(7), Val(i), u) for i=0:7]
+    Rv = @SVector[silv(Val(7), Val(i), v) for i=0:7]
+    Rw = @SVector[silv(Val(7), Val(i), 1 - u - v) for i=0:7]
+
+    dRu = @SVector[dsilv(Val(7), Val(i), u) for i=0:7]
+    dRv = @SVector[dsilv(Val(7), Val(i), v) for i=0:7]
+    dRw = @SVector[dsilv(Val(7), Val(i), 1 - u - v) for i=0:7]
+
+    𝓵¹ = zero(SVector{U,T})
+    𝓵² = zero(SVector{U,T})
+    @inbounds for (r, (i, j, k)) in enumerate(gmsh_triangle_index_to_tuple[7])
+
+        𝓵¹ += ch.vertices[r] * (
+                dRu[i + 1] * Rv[j+1] * Rw[k+1] -
+                Ru[i + 1] * Rv[j+1] * dRw[k+1]
+            )
+
+        𝓵² += ch.vertices[r] * (
+                Ru[i + 1] * dRv[j+1] * Rw[k+1] -
+                Ru[i + 1] * Rv[j+1] * dRw[k+1]
+            )
+    end
+    
+    return SVector(𝓵¹, 𝓵²)
+end
+
+function chart(
+    m::CompScienceMeshes.CurvilinearMesh{U,D1,N,T,O},
+    faceid::Int
+) where {U,D1,N,T,O}
+
+    vindices = m.faces[faceid]
+    X = m.vertices[vindices]
+
+    CompScienceMeshes.CurvilinearSimplex{U,D1-1,U-1,N,T,O}(X)
+end
+
+function normal(ch::CurvilinearSimplex{U,2}, bary) where {U}
+    t = unitarybase(ch, bary)
+
+    nt = cross(t[1], t[2])
+
+    return nt ./ norm(nt)
+end
+
+
+"""
+    jacobian(ch, ζ)
+
+Return the scalar Jacobian determinant.
+"""
+jacobian(ch::CurvilinearSimplex{U,1}, u) where {U} = norm(unitarybase(ch, u))
+
+function jacobian(ch::CurvilinearSimplex{U,2}, bary) where {U} 
+    t = unitarybase(ch, bary)
+
+    return norm(cross(t[1], t[2]))
+end
+
+"""
+    center(ch::CurvilinearSimplex)
+
+Return the physical midpoint of the 1D curvilinear element, obtained by
+evaluating the mapping at `ζ = 0.5`.
+"""
+center(ch::CurvilinearSimplex{U,1}) where {U} = barytocart(ch, SVector(0.5))
+center(ch::CurvilinearSimplex{U,2}) where {U} = barytocart(ch, SVector(0.5, 0.5))
+
+domain(::CurvilinearSimplex{U,1,C,N,T}) where {U,C,N,T} = ReferenceSimplex{1,T,2}()
+domain(::CurvilinearSimplex{U,2,C,N,T}) where {U,C,N,T} = ReferenceSimplex{2,T,3}()
+
+
+"""
+    volume(ch)
+
+Compute the physical length of the curvilinear element `ch` by integrating
+the Jacobian over `[0,1]`.
+"""
+function volume(ch::CurvilinearSimplex{U,1,C,N,T}) where {U,C,N,T}
+
+    ξ, w = legendre(20, 0.0, 1.0)
+
+    return dot(w, jacobian.(Ref(ch), ξ))
+end
+
+"""
+    volume(ch)
+
+Compute the physical length of the curvilinear element `ch` by integrating
+the Jacobian over `[0,1]`.
+"""
+# TODO: volume should be saved in the simplex type. However, Jacobian requires the existence...
+function volume(ch::CurvilinearSimplex{U,2,C,N,T}) where {U,C,N,T}
+
+    ξ, w = legendre(20, 0.0, 1.0)
+
+    vol = T(0)
+
+    for i in eachindex(ξ)
+        for j in eachindex(ξ)
+            vol += w[i] * w[j] * jacobian(ch, SVector(ξ[i], ξ[j]*(1 - ξ[i]))) * (1 - ξ[i])
+        end
+    end
+
+    return vol
+end
\ No newline at end of file
diff --git a/src/meshes/curvilinear_gmsh_line.jl b/src/meshes/curvilinear_gmsh_line.jl
new file mode 100644
index 0000000..c179fb3
--- /dev/null
+++ b/src/meshes/curvilinear_gmsh_line.jl
@@ -0,0 +1,139 @@
+# This file is created by reverse engineering Gmsh msh file
+# and using ChatGPT for analyzing the output.
+# Up to order 5 all results are manually verified.
+
+const gmsh_line_index_to_tuple_order1 = [
+    # Vertices
+    SVector(1, 0), # V1
+    SVector(0, 1), # V2
+]
+
+const gmsh_line_index_to_tuple_order2 = [
+    # Vertices
+    SVector(2, 0), # V1
+    SVector(0, 2), # V2
+
+    # Edge
+    SVector(1, 1),
+]
+
+const gmsh_line_index_to_tuple_order3 = [
+    # Vertices
+    SVector(3, 0), # V1
+    SVector(0, 3), # V2
+
+    # Edge
+    SVector(2, 1),
+    SVector(1, 2),
+]
+
+const gmsh_line_index_to_tuple_order4 = [
+    # Vertices
+    SVector(4, 0), # V1
+    SVector(0, 4), # V2
+
+    # Edge
+    SVector(3, 1),
+    SVector(2, 2),
+    SVector(1, 3),
+]
+
+const gmsh_line_index_to_tuple_order5 = [
+    # Vertices
+    SVector(5, 0), # V1
+    SVector(0, 5), # V2
+
+    # Edge
+    SVector(4, 1),
+    SVector(3, 2),
+    SVector(2, 3),
+    SVector(1, 4),
+]
+
+const gmsh_line_index_to_tuple_order6 = [
+    # Vertices
+    SVector(6, 0), # V1
+    SVector(0, 6), # V2
+
+    # Edge
+    SVector(5, 1),
+    SVector(4, 2),
+    SVector(3, 3),
+    SVector(2, 4),
+    SVector(1, 5),
+]
+
+const gmsh_line_index_to_tuple_order7 = [
+    # Vertices
+    SVector(7, 0), # V1
+    SVector(0, 7), # V2
+
+    # Edge
+    SVector(6, 1),
+    SVector(5, 2),
+    SVector(4, 3),
+    SVector(3, 4),
+    SVector(2, 5),
+    SVector(1, 6),
+]
+
+const gmsh_line_index_to_tuple_order8 = [
+    # Vertices
+    SVector(8, 0), # V1
+    SVector(0, 8), # V2
+
+    # Edge
+    SVector(7, 1),
+    SVector(6, 2),
+    SVector(5, 3),
+    SVector(4, 4),
+    SVector(3, 5),
+    SVector(2, 6),
+    SVector(1, 7),
+]
+
+const gmsh_line_index_to_tuple_order9 = [
+    # Vertices
+    SVector(9, 0), # V1
+    SVector(0, 9), # V2
+
+    # Edge
+    SVector(8, 1),
+    SVector(7, 2),
+    SVector(6, 3),
+    SVector(5, 4),
+    SVector(4, 5),
+    SVector(3, 6),
+    SVector(2, 7),
+    SVector(1, 8),
+]
+
+const gmsh_line_index_to_tuple_order10 = [
+    # Vertices
+    SVector(10, 0), # V1
+    SVector(0, 10), # V2
+
+    # Edge
+    SVector(9, 1),
+    SVector(8, 2),
+    SVector(7, 3),
+    SVector(6, 4),
+    SVector(5, 5),
+    SVector(4, 6),
+    SVector(3, 7),
+    SVector(2, 8),
+    SVector(1, 9),
+]
+
+const gmsh_line_index_to_tuple = [
+    gmsh_line_index_to_tuple_order1,
+    gmsh_line_index_to_tuple_order2,
+    gmsh_line_index_to_tuple_order3,
+    gmsh_line_index_to_tuple_order4,
+    gmsh_line_index_to_tuple_order5,
+    gmsh_line_index_to_tuple_order6,
+    gmsh_line_index_to_tuple_order7,
+    gmsh_line_index_to_tuple_order8,
+    gmsh_line_index_to_tuple_order9,
+    gmsh_line_index_to_tuple_order10,
+]
\ No newline at end of file
diff --git a/src/meshes/curvilinear_gmsh_triangle.jl b/src/meshes/curvilinear_gmsh_triangle.jl
new file mode 100644
index 0000000..b05c839
--- /dev/null
+++ b/src/meshes/curvilinear_gmsh_triangle.jl
@@ -0,0 +1,411 @@
+# This file is created by reverse engineering Gmsh msh file
+# and using ChatGPT for analyzing the output.
+# Up to order 5 all results are manually verified.
+
+const gmsh_triangle_index_to_tuple_order1 = [
+    # Vertices
+    SVector(1, 0, 0), # V1
+    SVector(0, 1, 0), # V2
+    SVector(0, 0, 1), # V3
+]
+
+const gmsh_triangle_index_to_tuple_order2 = [
+    # Vertices
+    SVector(2, 0, 0), # V1
+    SVector(0, 2, 0), # V2
+    SVector(0, 0, 2), # V3
+
+    # Edge 1: V1 -> V2
+    SVector(1, 1, 0),
+
+    # Edge 2: V2 -> V3
+    SVector(0, 1, 1),
+
+    # Edge 3: V2 -> V3
+    SVector(1, 0, 1),
+]
+
+const gmsh_triangle_index_to_tuple_order3 = [
+    # Vertices
+    SVector(3, 0, 0), # V1
+    SVector(0, 3, 0), # V2
+    SVector(0, 0, 3), # V3
+
+    # Edge 1: V1 -> V2
+    SVector(2, 1, 0),
+    SVector(1, 2, 0),
+
+    # Edge 2: V2 -> V3
+    SVector(0, 2, 1),
+    SVector(0, 1, 2),
+
+    # Edge 3: V2 -> V3
+    SVector(1, 0, 2),
+    SVector(2, 0, 1),
+
+    # Interior
+    SVector(1, 1, 1),
+]
+
+const gmsh_triangle_index_to_tuple_order4 = [
+    # Vertices
+    SVector(4, 0, 0), # V1
+    SVector(0, 4, 0), # V2
+    SVector(0, 0, 4), # V3
+
+    # Edge 1: V1 -> V2
+    SVector(3, 1, 0),
+    SVector(2, 2, 0),
+    SVector(1, 3, 0),
+
+    # Edge 2: V2 -> V3
+    SVector(0, 3, 1),
+    SVector(0, 2, 2),
+    SVector(0, 1, 3),
+
+    # Edge 3: V3 -> V1
+    SVector(1, 0, 3),
+    SVector(2, 0, 2),
+    SVector(3, 0, 1),
+
+    # Interior
+    SVector(2, 1, 1),
+    SVector(1, 2, 1),
+    SVector(1, 1, 2),
+]
+
+const gmsh_triangle_index_to_tuple_order5 = [
+    # Vertices
+    SVector(5, 0, 0), # V1
+    SVector(0, 5, 0), # V2
+    SVector(0, 0, 5), # V3
+
+    # Edge 1: V1 -> V2
+    SVector(4, 1, 0),
+    SVector(3, 2, 0),
+    SVector(2, 3, 0),
+    SVector(1, 4, 0),
+
+    # Edge 2: V2 -> V3
+    SVector(0, 4, 1),
+    SVector(0, 3, 2),
+    SVector(0, 2, 3),
+    SVector(0, 1, 4),
+
+    # Edge 3: V3 -> V1
+    SVector(1, 0, 4),
+    SVector(2, 0, 3),
+    SVector(3, 0, 2),
+    SVector(4, 0, 1),
+
+    # Interior
+    SVector(3, 1, 1),
+    SVector(1, 3, 1),
+    SVector(1, 1, 3),
+    SVector(2, 2, 1),
+    SVector(1, 2, 2),
+    SVector(2, 1, 2),
+]
+
+const gmsh_triangle_index_to_tuple_order6 = [
+    # Vertices
+    SVector(6, 0, 0), # V1
+    SVector(0, 6, 0), # V2
+    SVector(0, 0, 6), # V3
+
+    # Edge 1: V1 -> V2
+    SVector(5, 1, 0),
+    SVector(4, 2, 0),
+    SVector(3, 3, 0),
+    SVector(2, 4, 0),
+    SVector(1, 5, 0),
+
+    # Edge 2: V2 -> V3
+    SVector(0, 5, 1),
+    SVector(0, 4, 2),
+    SVector(0, 3, 3),
+    SVector(0, 2, 4),
+    SVector(0, 1, 5),
+
+    # Edge 3: V3 -> V1
+    SVector(1, 0, 5),
+    SVector(2, 0, 4),
+    SVector(3, 0, 3),
+    SVector(4, 0, 2),
+    SVector(5, 0, 1),
+
+    # Interior nodes
+    SVector(4, 1, 1),
+    SVector(1, 4, 1),
+    SVector(1, 1, 4),
+    SVector(3, 2, 1),
+    SVector(2, 3, 1),
+    SVector(1, 3, 2),
+    SVector(1, 2, 3),
+    SVector(2, 1, 3),
+    SVector(3, 1, 2),
+    SVector(2, 2, 2),
+]
+
+const gmsh_triangle_index_to_tuple_order7 = [
+    # Vertices
+    SVector(7, 0, 0), # V1
+    SVector(0, 7, 0), # V2
+    SVector(0, 0, 7), # V3
+
+    # Edge 1: V1 -> V2
+    SVector(6, 1, 0),
+    SVector(5, 2, 0),
+    SVector(4, 3, 0),
+    SVector(3, 4, 0),
+    SVector(2, 5, 0),
+    SVector(1, 6, 0),
+
+    # Edge 2: V2 -> V3
+    SVector(0, 6, 1),
+    SVector(0, 5, 2),
+    SVector(0, 4, 3),
+    SVector(0, 3, 4),
+    SVector(0, 2, 5),
+    SVector(0, 1, 6),
+
+    # Edge 3: V3 -> V1
+    SVector(1, 0, 6),
+    SVector(2, 0, 5),
+    SVector(3, 0, 4),
+    SVector(4, 0, 3),
+    SVector(5, 0, 2),
+    SVector(6, 0, 1),
+
+    # Interior
+    SVector(5, 1, 1),
+    SVector(1, 5, 1),
+    SVector(1, 1, 5),
+    SVector(4, 2, 1),
+    SVector(3, 3, 1),
+    SVector(2, 4, 1),
+    SVector(1, 4, 2),
+    SVector(1, 3, 3),
+    SVector(1, 2, 4),
+    SVector(2, 1, 4),
+    SVector(3, 1, 3),
+    SVector(4, 1, 2),
+    SVector(3, 2, 2),
+    SVector(2, 3, 2),
+    SVector(2, 2, 3),
+]
+
+const gmsh_triangle_index_to_tuple_order8 = [
+    # Vertices
+    SVector(8, 0, 0), # V1
+    SVector(0, 8, 0), # V2
+    SVector(0, 0, 8), # V3
+
+    # Edge 1: V1 -> V2
+    SVector(7, 1, 0),
+    SVector(6, 2, 0),
+    SVector(5, 3, 0),
+    SVector(4, 4, 0),
+    SVector(3, 5, 0),
+    SVector(2, 6, 0),
+    SVector(1, 7, 0),
+
+    # Edge 2: V2 -> V3
+    SVector(0, 7, 1),
+    SVector(0, 6, 2),
+    SVector(0, 5, 3),
+    SVector(0, 4, 4),
+    SVector(0, 3, 5),
+    SVector(0, 2, 6),
+    SVector(0, 1, 7),
+
+    # Edge 3: V3 -> V1
+    SVector(1, 0, 7),
+    SVector(2, 0, 6),
+    SVector(3, 0, 5),
+    SVector(4, 0, 4),
+    SVector(5, 0, 3),
+    SVector(6, 0, 2),
+    SVector(7, 0, 1),
+
+    # Interior
+    SVector(6, 1, 1),
+    SVector(1, 6, 1),
+    SVector(1, 1, 6),
+    SVector(5, 2, 1),
+    SVector(4, 3, 1),
+    SVector(3, 4, 1),
+    SVector(2, 5, 1),
+    SVector(1, 5, 2),
+    SVector(1, 4, 3),
+    SVector(1, 3, 4),
+    SVector(1, 2, 5),
+    SVector(2, 1, 5),
+    SVector(3, 1, 4),
+    SVector(4, 1, 3),
+    SVector(5, 1, 2),
+    SVector(4, 2, 2),
+    SVector(2, 4, 2),
+    SVector(2, 2, 4),
+    SVector(3, 3, 2),
+    SVector(2, 3, 3),
+    SVector(3, 2, 3),
+]
+
+const gmsh_triangle_index_to_tuple_order9 = [
+    # Vertices
+    SVector(9, 0, 0), # V1
+    SVector(0, 9, 0), # V2
+    SVector(0, 0, 9), # V3
+
+    # Edge 1: V1 -> V2
+    SVector(8, 1, 0),
+    SVector(7, 2, 0),
+    SVector(6, 3, 0),
+    SVector(5, 4, 0),
+    SVector(4, 5, 0),
+    SVector(3, 6, 0),
+    SVector(2, 7, 0),
+    SVector(1, 8, 0),
+
+    # Edge 2: V2 -> V3
+    SVector(0, 8, 1),
+    SVector(0, 7, 2),
+    SVector(0, 6, 3),
+    SVector(0, 5, 4),
+    SVector(0, 4, 5),
+    SVector(0, 3, 6),
+    SVector(0, 2, 7),
+    SVector(0, 1, 8),
+
+    # Edge 3: V3 -> V1
+    SVector(1, 0, 8),
+    SVector(2, 0, 7),
+    SVector(3, 0, 6),
+    SVector(4, 0, 5),
+    SVector(5, 0, 4),
+    SVector(6, 0, 3),
+    SVector(7, 0, 2),
+    SVector(8, 0, 1),
+
+    # Interior
+    SVector(7, 1, 1),
+    SVector(1, 7, 1),
+    SVector(1, 1, 7),
+    SVector(6, 2, 1),
+    SVector(5, 3, 1),
+    SVector(4, 4, 1),
+    SVector(3, 5, 1),
+    SVector(2, 6, 1),
+    SVector(1, 6, 2),
+    SVector(1, 5, 3),
+    SVector(1, 4, 4),
+    SVector(1, 3, 5),
+    SVector(1, 2, 6),
+    SVector(2, 1, 6),
+    SVector(3, 1, 5),
+    SVector(4, 1, 4),
+    SVector(5, 1, 3),
+    SVector(6, 1, 2),
+    SVector(5, 2, 2),
+    SVector(2, 5, 2),
+    SVector(2, 2, 5),
+    SVector(4, 3, 2),
+    SVector(3, 4, 2),
+    SVector(2, 4, 3),
+    SVector(2, 3, 4),
+    SVector(3, 2, 4),
+    SVector(4, 2, 3),
+    SVector(3, 3, 3),
+]
+
+const gmsh_triangle_index_to_tuple_order10 = [
+    # Vertices
+    SVector(10, 0, 0),
+    SVector(0, 10, 0),
+    SVector(0, 0, 10),
+
+    # Edge 1
+    SVector(9, 1, 0),
+    SVector(8, 2, 0),
+    SVector(7, 3, 0),
+    SVector(6, 4, 0),
+    SVector(5, 5, 0),
+    SVector(4, 6, 0),
+    SVector(3, 7, 0),
+    SVector(2, 8, 0),
+    SVector(1, 9, 0),
+
+    # Edge 2
+    SVector(0, 9, 1),
+    SVector(0, 8, 2),
+    SVector(0, 7, 3),
+    SVector(0, 6, 4),
+    SVector(0, 5, 5),
+    SVector(0, 4, 6),
+    SVector(0, 3, 7),
+    SVector(0, 2, 8),
+    SVector(0, 1, 9),
+
+    # Edge 3
+    SVector(1, 0, 9),
+    SVector(2, 0, 8),
+    SVector(3, 0, 7),
+    SVector(4, 0, 6),
+    SVector(5, 0, 5),
+    SVector(6, 0, 4),
+    SVector(7, 0, 3),
+    SVector(8, 0, 2),
+    SVector(9, 0, 1),
+
+    # Interior nodes
+    SVector(8, 1, 1),
+    SVector(1, 8, 1),
+    SVector(1, 1, 8),
+    SVector(7, 2, 1),
+    SVector(6, 3, 1),
+    SVector(5, 4, 1),
+    SVector(4, 5, 1),
+    SVector(3, 6, 1),
+    SVector(2, 7, 1),
+    SVector(1, 7, 2),
+    SVector(1, 6, 3),
+    SVector(1, 5, 4),
+    SVector(1, 4, 5),
+    SVector(1, 3, 6),
+    SVector(1, 2, 7),
+    SVector(2, 1, 7),
+    SVector(3, 1, 6),
+    SVector(4, 1, 5),
+    SVector(5, 1, 4),
+    SVector(6, 1, 3),
+    SVector(7, 1, 2),
+    SVector(6, 2, 2),
+    SVector(2, 6, 2),
+    SVector(2, 2, 6),
+    SVector(5, 3, 2),
+    SVector(4, 4, 2),
+    SVector(3, 5, 2),
+    SVector(2, 5, 3),
+    SVector(2, 4, 4),
+    SVector(2, 3, 5),
+    SVector(3, 2, 5),
+    SVector(4, 2, 4),
+    SVector(5, 2, 3),
+    SVector(4, 3, 3),
+    SVector(3, 4, 3),
+    SVector(3, 3, 4),
+]
+
+const gmsh_triangle_index_to_tuple = [
+    gmsh_triangle_index_to_tuple_order1,
+    gmsh_triangle_index_to_tuple_order2,
+    gmsh_triangle_index_to_tuple_order3,
+    gmsh_triangle_index_to_tuple_order4,
+    gmsh_triangle_index_to_tuple_order5,
+    gmsh_triangle_index_to_tuple_order6,
+    gmsh_triangle_index_to_tuple_order7,
+    gmsh_triangle_index_to_tuple_order8,
+    gmsh_triangle_index_to_tuple_order9,
+    gmsh_triangle_index_to_tuple_order10,
+]
\ 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..8ef1737
--- /dev/null
+++ b/src/meshes/curvilinear_mesh.jl
@@ -0,0 +1,75 @@
+# Mesh: 1D elements embedded in R^U, with N control points per element (N = order+1)
+mutable struct CurvilinearMesh{U,D1,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::AbstractVector{SVector{U,T}},
+        faces::AbstractVector{SVector{N,Int}},
+        manifolddimension,
+        order::Integer
+    ) where {U,N,T}
+        order > 0 || throw(ArgumentError("order must be positive."))
+        new{U,manifolddimension+1,N,T,Int(order)}(vertices, faces)
+    end
+end
+
+                    #=
+
+# Constructors from abstract containers
+CurvilinearMesh(vertices::AbstractVector{<:SVector{U,T}},
+                faces::AbstractVector{<:SVector{N,Int}},
+                manifolddimension::Integer,
+                order::Integer) where {U,N,T} =
+    CurvilinearMesh(Vector{SVector{U,T}}(vertices),
+                    Vector{SVector{N,Int}}(faces),
+                    manifolddimension,
+                    order)
+
+# From matrices (verts: U×nV, faces: N×nF)
+function CurvilinearMesh(verts::AbstractMatrix{T},
+                         faces::AbstractMatrix{<:Integer},
+                         manifolddimension,
+                         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,D1,N,T,O}) where {U,D1,N,T,O} = T
+vertextype(::CurvilinearMesh{U,D1,N,T,O}) where {U,D1,N,T,O} = SVector{U,T}
+universedimension(::CurvilinearMesh{U}) where {U} = U
+dimension(::CurvilinearMesh{U,D1}) where {U,D1} = D1-1
+
+function celltype(m::CurvilinearMesh{U,D1,N,T,O}) where {U,D1,N,T,O} SimplexGraph{N-O+1} end
+function celltype(m::CurvilinearMesh{U,D1,N,T,O}, ::Type{Val{M}}) where {U,D1,N,T,O,M} SimplexGraph{M+1} end
+
+function indextype(m::CurvilinearMesh{U,D1,N,T,O}) where {U,D1,N,T,O} SVector{N-O+1,Int} end
+function indextype(m::CurvilinearMesh{U,D1,N,T,O}, ::Type{Val{M}}) where {U,D1,N,T,O,M} SVector{M+1,Int} end
+
+function indices(m::CurvilinearMesh{U,D1,N,T,O}, cell) where {U,D1,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,D1,N,T,O}}) where {U,D1,N,T,O} = SVector{N,Int}
+
+meshorder(::CurvilinearMesh{U,D1,N,T,O}) where {U,D1,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..3ea0a8e
--- /dev/null
+++ b/src/meshes/curvilinear_neighborhood.jl
@@ -0,0 +1,28 @@
+"""
+    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 unitarybase(mp.patch, mp.bary)[i]
+end
+
+function normal(mp::MeshPointNM{T,C,1,U}) where {T,C<:CurvilinearSimplex,U}
+    t = tangents(mp, 1)
+    nt = SVector(-t[2], t[1])
+
+    return nt ./ norm(nt)
+end
+
+function normal(mp::MeshPointNM{T,C,2,U}) where {T,C<:CurvilinearSimplex,U}
+    t1 = tangents(mp, 1)
+    t2 = tangents(mp, 2)
+
+    nt = cross(t1, t2)
+
+    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)
+end
\ No newline at end of file
diff --git a/src/meshes/curvilinear_silvester.jl b/src/meshes/curvilinear_silvester.jl
new file mode 100644
index 0000000..b47e3de
--- /dev/null
+++ b/src/meshes/curvilinear_silvester.jl
@@ -0,0 +1,223 @@
+"""
+    silvester(p, x)
+Compute the Silvester polynomials of order `p` at `x`.
+"""
+function silvester(p, x)
+    T = typeof(x)
+
+    R = zeros(T, p + 1)
+    R[1] = T(1)
+
+    for i = 0:(p-1)
+        R[i + 1 + 1] = (p*x - i)/(i + 1)*R[i + 1]
+    end
+
+    return R
+end
+
+"""
+    dsilvester(p, x)
+Compute the derivatives of the Silvester polynomials of order `p` at `x`.
+"""
+function dsilvester(p, x)
+    T = typeof(x)
+
+    R = silvester(p, x)
+
+    dR = zeros(T, p + 1)
+
+    # (2.6)
+    dR[1] = T(0)
+
+    for i = 0:(p - 1)
+        dR[i + 1 + 1] = p/(i + 1) * R[i + 1] + (p*x - i)/(i + 1) * dR[i + 1]
+    end
+
+    return dR
+end
+
+"""
+    fast_silvester(p, x)
+Compute the Silvester polynomials of order `p` and their derivatives at `x`.
+"""
+# TODO: Consider using pre-allocated memory.
+function fast_silvester(p, x)
+    T = typeof(x)
+
+    R = zeros(T, p + 1)
+    dR = zeros(T, p + 1)
+
+    R[1] = T(1)
+    dR[1] = T(0)
+
+    for i = 0:(p - 1)
+        R[i + 1 + 1] = (p*x - i)/(i + 1)*R[i + 1]
+        dR[i + 1 + 1] = p/(i + 1) * R[i + 1] + (p*x - i)/(i + 1) * dR[i + 1]
+    end
+
+    return R, dR
+end
+
+#### Silvester polynomials up to order 10 ####
+# Explicitly calculated. However, it seems that fast_silvester(p, x) 
+# approach is faster despite memory allocations.
+# These lines should be removed if we settle at the fast_silver approach.
+@inline silv(::Val{1}, ::Val{0}, ξ::T) where {T} = T(1)
+@inline silv(::Val{1}, ::Val{1}, ξ::T) where {T} = ξ
+
+@inline silv(::Val{2}, ::Val{0}, ξ::T) where {T} = T(1)
+@inline silv(::Val{2}, ::Val{1}, ξ::T) where {T} = 2 * ξ
+@inline silv(::Val{2}, ::Val{2}, ξ::T) where {T} = 2 * ξ^2 - ξ
+
+@inline silv(::Val{3}, ::Val{0}, ξ::T) where {T} = T(1)
+@inline silv(::Val{3}, ::Val{1}, ξ::T) where {T} = 3ξ
+@inline silv(::Val{3}, ::Val{2}, ξ::T) where {T} = (9ξ^2 - 3ξ) / 2
+@inline silv(::Val{3}, ::Val{3}, ξ::T) where {T} = (27ξ^3 - 27ξ^2 + 6ξ) / 6
+
+
+@inline silv(::Val{4}, ::Val{0}, ξ::T) where {T} = T(1)
+@inline silv(::Val{4}, ::Val{1}, ξ::T) where {T} = 4ξ
+@inline silv(::Val{4}, ::Val{2}, ξ::T) where {T} = (16ξ^2 - 4ξ) / 2
+@inline silv(::Val{4}, ::Val{3}, ξ::T) where {T} = (64ξ^3 - 48ξ^2 + 8ξ) / 6
+@inline silv(::Val{4}, ::Val{4}, ξ::T) where {T} = (256ξ^4 - 256ξ^3 + 96ξ^2 - 12ξ) / 24
+
+
+@inline silv(::Val{5}, ::Val{0}, ξ::T) where {T} = T(1)
+@inline silv(::Val{5}, ::Val{1}, ξ::T) where {T} = 5ξ
+@inline silv(::Val{5}, ::Val{2}, ξ::T) where {T} = (25ξ^2 - 5ξ) / 2
+@inline silv(::Val{5}, ::Val{3}, ξ::T) where {T} = (125ξ^3 - 75ξ^2 + 10ξ) / 6
+@inline silv(::Val{5}, ::Val{4}, ξ::T) where {T} = (625ξ^4 - 500ξ^3 + 150ξ^2 - 20ξ) / 24
+@inline silv(::Val{5}, ::Val{5}, ξ::T) where {T} = (3125ξ^5 - 3125ξ^4 + 1250ξ^3 - 250ξ^2 + 20ξ) / 120
+
+
+@inline silv(::Val{6}, ::Val{0}, ξ::T) where {T} = T(1)
+@inline silv(::Val{6}, ::Val{1}, ξ::T) where {T} = 6ξ
+@inline silv(::Val{6}, ::Val{2}, ξ::T) where {T} = (36ξ^2 - 6ξ) / 2
+@inline silv(::Val{6}, ::Val{3}, ξ::T) where {T} = (216ξ^3 - 108ξ^2 + 12ξ) / 6
+@inline silv(::Val{6}, ::Val{4}, ξ::T) where {T} = (1296ξ^4 - 864ξ^3 + 192ξ^2 - 24ξ) / 24
+@inline silv(::Val{6}, ::Val{5}, ξ::T) where {T} = (7776ξ^5 - 6480ξ^4 + 2160ξ^3 - 360ξ^2 + 30ξ) / 120
+@inline silv(::Val{6}, ::Val{6}, ξ::T) where {T} = (46656ξ^6 - 46656ξ^5 + 19440ξ^4 - 4320ξ^3 + 540ξ^2 - 30ξ) / 720
+
+
+@inline silv(::Val{7}, ::Val{0}, ξ::T) where {T} = T(1)
+@inline silv(::Val{7}, ::Val{1}, ξ::T) where {T} = 7ξ
+@inline silv(::Val{7}, ::Val{2}, ξ::T) where {T} = (49ξ^2 - 7ξ) / 2
+@inline silv(::Val{7}, ::Val{3}, ξ::T) where {T} = (343ξ^3 - 147ξ^2 + 14ξ) / 6
+@inline silv(::Val{7}, ::Val{4}, ξ::T) where {T} = (2401ξ^4 - 1372ξ^3 + 294ξ^2 - 28ξ) / 24
+@inline silv(::Val{7}, ::Val{5}, ξ::T) where {T} = (16807ξ^5 - 11760ξ^4 + 3528ξ^3 - 560ξ^2 + 35ξ) / 120
+@inline silv(::Val{7}, ::Val{6}, ξ::T) where {T} = (117649ξ^6 - 100842ξ^5 + 38808ξ^4 - 8400ξ^3 + 1050ξ^2 - 42ξ) / 720
+@inline silv(::Val{7}, ::Val{7}, ξ::T) where {T} = (823543ξ^7 - 823543ξ^6 + 352947ξ^5 - 88200ξ^4 + 13860ξ^3 - 1260ξ^2 + 56ξ) / 5040
+
+
+@inline silv(::Val{8}, ::Val{0}, ξ::T) where {T} = T(1)
+@inline silv(::Val{8}, ::Val{1}, ξ::T) where {T} = 8ξ
+@inline silv(::Val{8}, ::Val{2}, ξ::T) where {T} = (64ξ^2 - 8ξ) / 2
+@inline silv(::Val{8}, ::Val{3}, ξ::T) where {T} = (512ξ^3 - 192ξ^2 + 16ξ) / 6
+@inline silv(::Val{8}, ::Val{4}, ξ::T) where {T} = (4096ξ^4 - 2048ξ^3 + 384ξ^2 - 32ξ) / 24
+@inline silv(::Val{8}, ::Val{5}, ξ::T) where {T} = (32768ξ^5 - 20480ξ^4 + 5120ξ^3 - 800ξ^2 + 48ξ) / 120
+@inline silv(::Val{8}, ::Val{6}, ξ::T) where {T} = (262144ξ^6 - 196608ξ^5 + 73728ξ^4 - 16128ξ^3 + 2016ξ^2 - 72ξ) / 720
+@inline silv(::Val{8}, ::Val{7}, ξ::T) where {T} = (2097152ξ^7 - 1835008ξ^6 + 802816ξ^5 - 206080ξ^4 + 32256ξ^3 - 3024ξ^2 + 84ξ) / 5040
+@inline silv(::Val{8}, ::Val{8}, ξ::T) where {T} = (16777216ξ^8 - 16777216ξ^7 + 7520256ξ^6 - 2007040ξ^5 + 345600ξ^4 - 38400ξ^3 + 2688ξ^2 - 96ξ) / 40320
+
+
+@inline silv(::Val{9}, ::Val{0}, ξ::T) where {T} = T(1)
+@inline silv(::Val{9}, ::Val{1}, ξ::T) where {T} = 9ξ
+@inline silv(::Val{9}, ::Val{2}, ξ::T) where {T} = (81ξ^2 - 9ξ) / 2
+@inline silv(::Val{9}, ::Val{3}, ξ::T) where {T} = (729ξ^3 - 243ξ^2 + 18ξ) / 6
+@inline silv(::Val{9}, ::Val{4}, ξ::T) where {T} = (6561ξ^4 - 2916ξ^3 + 486ξ^2 - 36ξ) / 24
+@inline silv(::Val{9}, ::Val{5}, ξ::T) where {T} = (59049ξ^5 - 32805ξ^4 + 8100ξ^3 - 1215ξ^2 + 54ξ) / 120
+@inline silv(::Val{9}, ::Val{6}, ξ::T) where {T} = (531441ξ^6 - 354294ξ^5 + 118098ξ^4 - 24192ξ^3 + 2916ξ^2 - 90ξ) / 720
+@inline silv(::Val{9}, ::Val{7}, ξ::T) where {T} = (4782969ξ^7 - 3835122ξ^6 + 1536795ξ^5 - 376200ξ^4 + 58806ξ^3 - 5292ξ^2 + 108ξ) / 5040
+@inline silv(::Val{9}, ::Val{8}, ξ::T) where {T} = (43046721ξ^8 - 38742048ξ^7 + 16796160ξ^6 - 4469856ξ^5 + 777600ξ^4 - 86400ξ^3 + 6048ξ^2 - 120ξ) / 40320
+@inline silv(::Val{9}, ::Val{9}, ξ::T) where {T} = (387420489ξ^9 - 387420489ξ^8 + 172186884ξ^7 - 48189030ξ^6 + 8845200ξ^5 - 1108800ξ^4 + 95040ξ^3 - 5040ξ^2 + 144ξ) / 362880
+
+@inline silv(::Val{10}, ::Val{0}, ξ::T) where {T} = T(1)
+@inline silv(::Val{10}, ::Val{1}, ξ::T) where {T} = 10ξ
+@inline silv(::Val{10}, ::Val{2}, ξ::T) where {T} = (100ξ^2 - 10ξ) / 2
+@inline silv(::Val{10}, ::Val{3}, ξ::T) where {T} = (1000ξ^3 - 300ξ^2 + 20ξ) / 6
+@inline silv(::Val{10}, ::Val{4}, ξ::T) where {T} = (10000ξ^4 - 4000ξ^3 + 600ξ^2 - 40ξ) / 24
+@inline silv(::Val{10}, ::Val{5}, ξ::T) where {T} = (100000ξ^5 - 50000ξ^4 + 10000ξ^3 - 1250ξ^2 + 50ξ) / 120
+@inline silv(::Val{10}, ::Val{6}, ξ::T) where {T} = (1000000ξ^6 - 600000ξ^5 + 180000ξ^4 - 33750ξ^3 + 3750ξ^2 - 100ξ) / 720
+@inline silv(::Val{10}, ::Val{7}, ξ::T) where {T} = (10000000ξ^7 - 7000000ξ^6 + 2520000ξ^5 - 588000ξ^4 + 92400ξ^3 - 8820ξ^2 + 140ξ) / 5040
+@inline silv(::Val{10}, ::Val{8}, ξ::T) where {T} = (100000000ξ^8 - 80000000ξ^7 + 33600000ξ^6 - 9408000ξ^5 + 1814400ξ^4 - 241920ξ^3 + 20160ξ^2 - 160ξ) / 40320
+@inline silv(::Val{10}, ::Val{9}, ξ::T) where {T} = (1000000000ξ^9 - 900000000ξ^8 + 405000000ξ^7 - 121500000ξ^6 + 27000000ξ^5 - 4320000ξ^4 + 472500ξ^3 - 31500ξ^2 + 180ξ) / 362880
+@inline silv(::Val{10}, ::Val{10}, ξ::T) where {T} = (10000000000ξ^10 - 10000000000ξ^9 + 4500000000ξ^8 - 1350000000ξ^7 + 300000000ξ^6 - 50000000ξ^5 + 6000000ξ^4 - 480000ξ^3 + 24000ξ^2 - 200ξ) / 3628800
+
+@inline dsilv(::Val{1}, ::Val{0}, ξ::T) where {T} = T(0)
+@inline dsilv(::Val{1}, ::Val{1}, ξ::T) where {T} = T(1)
+
+@inline dsilv(::Val{2}, ::Val{0}, ξ::T) where {T} = T(0)
+@inline dsilv(::Val{2}, ::Val{1}, ξ::T) where {T} = T(2)
+@inline dsilv(::Val{2}, ::Val{2}, ξ::T) where {T} = 4 * ξ - 1
+
+@inline dsilv(::Val{3}, ::Val{0}, ξ::T) where {T} = T(0)
+@inline dsilv(::Val{3}, ::Val{1}, ξ::T) where {T} = T(3)
+@inline dsilv(::Val{3}, ::Val{2}, ξ::T) where {T} = (18ξ - 3) / 2
+@inline dsilv(::Val{3}, ::Val{3}, ξ::T) where {T} = (81ξ^2 - 54ξ + 6) / 6
+
+@inline dsilv(::Val{4}, ::Val{0}, ξ::T) where {T} = T(0)
+@inline dsilv(::Val{4}, ::Val{1}, ξ::T) where {T} = T(4)
+@inline dsilv(::Val{4}, ::Val{2}, ξ::T) where {T} = (32ξ - 4) / 2
+@inline dsilv(::Val{4}, ::Val{3}, ξ::T) where {T} = (192ξ^2 - 96ξ + 8) / 6
+@inline dsilv(::Val{4}, ::Val{4}, ξ::T) where {T} = (1024ξ^3 - 768ξ^2 + 192ξ - 12) / 24
+
+@inline dsilv(::Val{5}, ::Val{0}, ξ::T) where {T} = T(0)
+@inline dsilv(::Val{5}, ::Val{1}, ξ::T) where {T} = T(5)
+@inline dsilv(::Val{5}, ::Val{2}, ξ::T) where {T} = (50ξ - 5) / 2
+@inline dsilv(::Val{5}, ::Val{3}, ξ::T) where {T} = (375ξ^2 - 150ξ + 10) / 6
+@inline dsilv(::Val{5}, ::Val{4}, ξ::T) where {T} = (2500ξ^3 - 1500ξ^2 + 300ξ - 20) / 24
+@inline dsilv(::Val{5}, ::Val{5}, ξ::T) where {T} = (15625ξ^4 - 12500ξ^3 + 3750ξ^2 - 500ξ + 20) / 120
+
+@inline dsilv(::Val{6}, ::Val{0}, ξ::T) where {T} = T(0)
+@inline dsilv(::Val{6}, ::Val{1}, ξ::T) where {T} = T(6)
+@inline dsilv(::Val{6}, ::Val{2}, ξ::T) where {T} = (72ξ - 6) / 2
+@inline dsilv(::Val{6}, ::Val{3}, ξ::T) where {T} = (648ξ^2 - 216ξ + 12) / 6
+@inline dsilv(::Val{6}, ::Val{4}, ξ::T) where {T} = (5184ξ^3 - 2592ξ^2 + 384ξ - 24) / 24
+@inline dsilv(::Val{6}, ::Val{5}, ξ::T) where {T} = (38880ξ^4 - 25920ξ^3 + 6480ξ^2 - 720ξ + 30) / 120
+@inline dsilv(::Val{6}, ::Val{6}, ξ::T) where {T} = (279936ξ^5 - 233280ξ^4 + 77760ξ^3 - 12960ξ^2 + 1080ξ - 30) / 720
+
+
+@inline dsilv(::Val{7}, ::Val{0}, ξ::T) where {T} = T(0)
+@inline dsilv(::Val{7}, ::Val{1}, ξ::T) where {T} = T(7)
+@inline dsilv(::Val{7}, ::Val{2}, ξ::T) where {T} = (98ξ - 7) / 2
+@inline dsilv(::Val{7}, ::Val{3}, ξ::T) where {T} = (1029ξ^2 - 294ξ + 14) / 6
+@inline dsilv(::Val{7}, ::Val{4}, ξ::T) where {T} = (9604ξ^3 - 4116ξ^2 + 588ξ - 28) / 24
+@inline dsilv(::Val{7}, ::Val{5}, ξ::T) where {T} = (84035ξ^4 - 47040ξ^3 + 10584ξ^2 - 1120ξ + 35) / 120
+@inline dsilv(::Val{7}, ::Val{6}, ξ::T) where {T} = (705894ξ^5 - 504210ξ^4 + 155232ξ^3 - 25200ξ^2 + 2100ξ - 42) / 720
+@inline dsilv(::Val{7}, ::Val{7}, ξ::T) where {T} = (5764801ξ^6 - 4941258ξ^5 + 1764735ξ^4 - 352800ξ^3 + 41580ξ^2 - 2520ξ + 56) / 5040
+
+
+@inline dsilv(::Val{8}, ::Val{0}, ξ::T) where {T} = T(0)
+@inline dsilv(::Val{8}, ::Val{1}, ξ::T) where {T} = T(8)
+@inline dsilv(::Val{8}, ::Val{2}, ξ::T) where {T} = (128ξ - 8) / 2
+@inline dsilv(::Val{8}, ::Val{3}, ξ::T) where {T} = (1536ξ^2 - 384ξ + 16) / 6
+@inline dsilv(::Val{8}, ::Val{4}, ξ::T) where {T} = (16384ξ^3 - 6144ξ^2 + 768ξ - 32) / 24
+@inline dsilv(::Val{8}, ::Val{5}, ξ::T) where {T} = (163840ξ^4 - 81920ξ^3 + 15360ξ^2 - 1600ξ + 48) / 120
+@inline dsilv(::Val{8}, ::Val{6}, ξ::T) where {T} = (1572864ξ^5 - 983040ξ^4 + 294912ξ^3 - 48384ξ^2 + 4032ξ - 72) / 720
+@inline dsilv(::Val{8}, ::Val{7}, ξ::T) where {T} = (14680064ξ^6 - 11010048ξ^5 + 4014080ξ^4 - 824320ξ^3 + 96768ξ^2 - 6048ξ + 84) / 5040
+@inline dsilv(::Val{8}, ::Val{8}, ξ::T) where {T} = (134217728ξ^7 - 117440512ξ^6 + 45121536ξ^5 - 10035200ξ^4 + 1382400ξ^3 - 115200ξ^2 + 5376ξ - 96) / 40320
+
+
+@inline dsilv(::Val{9}, ::Val{0}, ξ::T) where {T} = T(0)
+@inline dsilv(::Val{9}, ::Val{1}, ξ::T) where {T} = T(9)
+@inline dsilv(::Val{9}, ::Val{2}, ξ::T) where {T} = (162ξ - 9) / 2
+@inline dsilv(::Val{9}, ::Val{3}, ξ::T) where {T} = (2187ξ^2 - 486ξ + 18) / 6
+@inline dsilv(::Val{9}, ::Val{4}, ξ::T) where {T} = (26244ξ^3 - 8748ξ^2 + 972ξ - 36) / 24
+@inline dsilv(::Val{9}, ::Val{5}, ξ::T) where {T} = (295245ξ^4 - 131220ξ^3 + 24300ξ^2 - 2430ξ + 54) / 120
+@inline dsilv(::Val{9}, ::Val{6}, ξ::T) where {T} = (3188646ξ^5 - 1771470ξ^4 + 472392ξ^3 - 72576ξ^2 + 5832ξ - 90) / 720
+@inline dsilv(::Val{9}, ::Val{7}, ξ::T) where {T} = (34012203ξ^6 - 23010732ξ^5 + 7683975ξ^4 - 1504800ξ^3 + 176418ξ^2 - 10584ξ + 108) / 5040
+@inline dsilv(::Val{9}, ::Val{8}, ξ::T) where {T} = (362797056ξ^7 - 310136448ξ^6 + 117573120ξ^5 - 27974400ξ^4 + 4665600ξ^3 - 518400ξ^2 + 30240ξ - 120) / 40320
+@inline dsilv(::Val{9}, ::Val{9}, ξ::T) where {T} = (3874204890ξ^8 - 3486784401ξ^7 + 1205308188ξ^6 - 289134180ξ^5 + 44226000ξ^4 - 4435200ξ^3 + 285120ξ^2 - 10080ξ + 144) / 362880
+
+
+@inline dsilv(::Val{10}, ::Val{0}, ξ::T) where {T} = T(0)
+@inline dsilv(::Val{10}, ::Val{1}, ξ::T) where {T} = T(10)
+@inline dsilv(::Val{10}, ::Val{2}, ξ::T) where {T} = (200ξ - 10) / 2
+@inline dsilv(::Val{10}, ::Val{3}, ξ::T) where {T} = (3000ξ^2 - 600ξ + 20) / 6
+@inline dsilv(::Val{10}, ::Val{4}, ξ::T) where {T} = (40000ξ^3 - 12000ξ^2 + 1200ξ - 40) / 24
+@inline dsilv(::Val{10}, ::Val{5}, ξ::T) where {T} = (500000ξ^4 - 200000ξ^3 + 30000ξ^2 - 2500ξ + 50) / 120
+@inline dsilv(::Val{10}, ::Val{6}, ξ::T) where {T} = (6000000ξ^5 - 3000000ξ^4 + 720000ξ^3 - 101250ξ^2 + 7500ξ - 100) / 720
+@inline dsilv(::Val{10}, ::Val{7}, ξ::T) where {T} = (70000000ξ^6 - 42000000ξ^5 + 12600000ξ^4 - 2352000ξ^3 + 277200ξ^2 - 17640ξ + 140) / 5040
+@inline dsilv(::Val{10}, ::Val{8}, ξ::T) where {T} = (800000000ξ^7 - 560000000ξ^6 + 201600000ξ^5 - 47040000ξ^4 + 7257600ξ^3 - 725760ξ^2 + 40320ξ - 160) / 40320
+@inline dsilv(::Val{10}, ::Val{9}, ξ::T) where {T} = (9000000000ξ^8 - 7200000000ξ^7 + 2835000000ξ^6 - 729000000ξ^5 + 135000000ξ^4 - 17280000ξ^3 + 1417500ξ^2 - 63000ξ + 180) / 362880
+@inline dsilv(::Val{10}, ::Val{10}, ξ::T) where {T} = (100000000000ξ^9 - 90000000000ξ^8 + 36000000000ξ^7 - 9450000000ξ^6 + 1800000000ξ^5 - 250000000ξ^4 + 24000000ξ^3 - 1440000ξ^2 + 48000ξ - 200) / 3628800
\ No newline at end of file
diff --git a/src/primitives/linemeshes/linemeshes.jl b/src/primitives/linemeshes/linemeshes.jl
index 3d2372f..7188d83 100644
--- a/src/primitives/linemeshes/linemeshes.jl
+++ b/src/primitives/linemeshes/linemeshes.jl
@@ -57,3 +57,61 @@ 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)
+    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/src/primitives/surfacemeshes/mesh_sphere.jl b/src/primitives/surfacemeshes/mesh_sphere.jl
index b01325e..c521109 100644
--- a/src/primitives/surfacemeshes/mesh_sphere.jl
+++ b/src/primitives/surfacemeshes/mesh_sphere.jl
@@ -232,7 +232,7 @@ Create a mesh of a sphere of radius `radius` by parsing a .geo script
 
 The target edge size is `delta`.
 """
-function gmshsphere(radius, delta; tempname=tempname())
+function gmshsphere(radius, delta; tempname=tempname(), order=1)
     s = """
 lc = $delta;
 
@@ -277,11 +277,20 @@ Ruled Surface(4)={4} In Sphere{1};
     gmsh.option.setNumber("General.Terminal", 0)
     gmsh.option.setNumber("Mesh.MshFileVersion",2)
     gmsh.open(fn)
+
+    # Choose element order (1: linear segments, 2: quadratic)
+    gmsh.option.setNumber("Mesh.ElementOrder", order)
+
     gmsh.model.mesh.generate(2)
     gmsh.write(fno)
     gmsh.finalize()
 
-    m = read_gmsh_mesh(fno)
+    m = load_gmsh_mesh(fno,
+        element=:triangle,
+        vertextype=Float64,
+        order=order,
+        sort=false
+    )
 
     rm(fno)
     rm(fn)
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..b85a834
--- /dev/null
+++ b/test/test_curvilinear_elements.jl
@@ -0,0 +1,334 @@
+using Test
+using CompScienceMeshes
+using StaticArrays
+
+## Line order 1
+v1 = SVector(1.0, 0.0)
+v2 = SVector(0.0, 0.0)
+
+verts1 = [v1, v2]
+faces1 = [SVector(1,2)]
+
+m1 = CompScienceMeshes.Mesh(verts1, faces1)
+m2 = CompScienceMeshes.CurvilinearMesh(verts1, faces1, 1, 1)
+
+ch1 = chart(m1, 1)
+ch2 = chart(m2, 1)
+
+@test barytocart(ch1, 1.0) == SVector(1.0, 0.0)
+@test barytocart(ch2, 1.0) == SVector(1.0, 0.0)
+@test barytocart(ch1, 0.0) == SVector(0.0, 0.0)
+@test barytocart(ch2, 0.0) == SVector(0.0, 0.0)
+
+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)
+
+## Line order 2
+
+v3 = SVector(0.5, 0.0)
+verts2 = [v1, v2, v3]
+faces2 = [SVector(1,2,3)]
+
+m2 = CompScienceMeshes.CurvilinearMesh(verts2, faces2, 1, 2)
+
+ch2 = chart(m2, 1)
+
+@test barytocart(ch1, 1.0) == SVector(1.0, 0.0)
+@test barytocart(ch2, 1.0) == SVector(1.0, 0.0)
+@test barytocart(ch1, 0.0) == SVector(0.0, 0.0)
+@test barytocart(ch2, 0.0) == SVector(0.0, 0.0)
+
+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
+
+## Line order 2  (parabola)
+
+v1 = SVector(0.0, 0.0^2)
+v2 = SVector(1.0, 1.0^2)
+v3 = SVector(0.50, 0.50^2)
+
+verts_parabola = [v1, v2, v3]
+face_parabola = [SVector(1,2,3)]
+
+m_parabola = CompScienceMeshes.CurvilinearMesh(verts_parabola, face_parabola, 1, 2)
+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)
+@test tangents(mp_2, 1) ≈ SVector(-1.0, -1.5)
+@test tangents(mp_3, 1) ≈ SVector(-1.0, -1.0)
+@test tangents(mp_4, 1) ≈ SVector(-1.0, -0.5)
+@test tangents(mp_5, 1) ≈ SVector(-1.0, -0.0)
+
+@test volume(ch_parabola) ≈ 1.478942857544597 # (Reference from Mathematica 14)
+
+
+
+## Triangle order 1
+v1 = SVector(1.0, 0.0, 0.0)
+v2 = SVector(0.0, 1.0, 0.0)
+v3 = SVector(0.0, 0.0, 0.0)
+
+verts1 = [v1, v2, v3]
+faces1 = [SVector(1, 2, 3)]
+
+m1 = CompScienceMeshes.Mesh(verts1, faces1)
+m2 = CompScienceMeshes.CurvilinearMesh(verts1, faces1, 2, 1)
+
+ch1 = chart(m1, 1)
+ch2 = chart(m2, 1)
+
+@test barytocart(ch1, SVector(1.0, 0.0)) == SVector(1.0, 0.0, 0.0)
+@test barytocart(ch1, SVector(0.0, 1.0)) == SVector(0.0, 1.0, 0.0)
+@test barytocart(ch1, SVector(0.0, 0.0)) == SVector(0.0, 0.0, 0.0)
+@test barytocart(ch2, SVector(1.0, 0.0)) == SVector(1.0, 0.0, 0.0)
+@test barytocart(ch2, SVector(0.0, 1.0)) == SVector(0.0, 1.0, 0.0)
+@test barytocart(ch2, SVector(0.0, 0.0)) == SVector(0.0, 0.0, 0.0)
+
+CompScienceMeshes.unitarybase(ch2, SVector(0.0, 0.0))
+CompScienceMeshes.unitarybase(ch2, SVector(1.0, 0.0))
+CompScienceMeshes.unitarybase(ch2, SVector(0.0, 1.0))
+CompScienceMeshes.unitarybase(ch2, SVector(0.5, 0.5))
+
+## 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, 1, 4)
+
+ch1 = chart(m1, 1)
+ch2 = chart(m2, 1)
+
+@test barytocart(ch1, 1.0) == SVector(0.0, 0.0)
+@test barytocart(ch2, 1.0) == SVector(0.0, 0.0)
+@test barytocart(ch1, 0.0) == SVector(1.0, 0.0)
+@test barytocart(ch2, 0.0) == SVector(1.0, 0.0)
+
+mp1 = neighborhood(ch1, 0.0)
+mp2 = neighborhood(ch2, 0.0)
+
+@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 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 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
+
+##
+v1 = SVector(1.0, 0.0, 0.0)
+v2 = SVector(0.0, 1.0, 0.0)
+v3 = SVector(0.0, 0.0, 0.0)
+
+verts1 = [v1, v2, v3]
+faces1 = [SVector(1, 2, 3)]
+
+
+verts2 = [
+    SVector(1.0, 0.0, 0.0),                         # 1
+    SVector(0.0, 1.0, 0.0),                         # 2
+    SVector(0.0, 0.0, 0.0),                         # 3
+    SVector(0.0, 0.5000000000013304, 0.0),          # 4
+    SVector(0.4999999999986718, 0.0, 0.0),          # 5
+    SVector(0.5000000000013293, 0.4999999999986707, 0.0) # 6
+]
+
+face2 = [SVector(1, 2, 3, 6, 4, 5)]
+
+m1 = CompScienceMeshes.Mesh(verts1, faces1)
+m2 = CompScienceMeshes.CurvilinearMesh(verts2, face2, 2, 2)
+
+ch1 = chart(m1, 1)
+ch2 = chart(m2, 1)
+
+## Order 7
+
+v1 = SVector(1.0, 0.0, 0.0)
+v2 = SVector(0.0, 1.0, 0.0)
+v3 = SVector(0.0, 0.0, 0.0)
+
+verts1 = [v1, v2, v3]
+faces1 = [SVector(1, 2, 3)]
+
+verts2 = [
+    SVector(1.0, 0.0, 0.0),
+    SVector(0.0, 1.0, 0.0),
+    SVector(0.0, 0.0, 0.0),
+    SVector(0.0, 0.8571428571429518, 0.0),
+    SVector(0.0, 0.7142857142863779, 0.0),
+    SVector(0.0, 0.5714285714298644, 0.0),
+    SVector(0.0, 0.428571428572595, 0.0),
+    SVector(0.0, 0.2857142857150634, 0.0),
+    SVector(0.0, 0.1428571428575318, 0.0),
+    SVector(0.1428571428568644, 0.0, 0.0),
+    SVector(0.2857142857135911, 0.0, 0.0),
+    SVector(0.4285714285703177, 0.0, 0.0),
+    SVector(0.5714285714274409, 0.0, 0.0),
+    SVector(0.7142857142849605, 0.0, 0.0),
+    SVector(0.8571428571424802, 0.0, 0.0),
+    SVector(0.8571428571430437, 0.1428571428569563, 0.0),
+    SVector(0.7142857142863934, 0.2857142857136066, 0.0),
+    SVector(0.5714285714297733, 0.4285714285702267, 0.0),
+    SVector(0.428571428572577, 0.571428571427423, 0.0),
+    SVector(0.2857142857150514, 0.7142857142849486, 0.0),
+    SVector(0.1428571428575258, 0.8571428571424742, 0.0),
+    SVector(0.7142857142855676, 0.1428571428569817, 0.0),
+    SVector(0.1428571428574321, 0.7142857142854938, 0.0),
+    SVector(0.1428571428569261, 0.1428571428574504, 0.0),
+    SVector(0.5714285714289259, 0.2857142857136924, 0.0),
+    SVector(0.4285714285721156, 0.4285714285706285, 0.0),
+    SVector(0.2857142857148444, 0.5714285714280445, 0.0),
+    SVector(0.1428571428573385, 0.5714285714289371, 0.0),
+    SVector(0.1428571428571978, 0.4285714285721468, 0.0),
+    SVector(0.1428571428570315, 0.2857142857148569, 0.0),
+    SVector(0.285714285713734, 0.1428571428573691, 0.0),
+    SVector(0.428571428570712, 0.1428571428572502, 0.0),
+    SVector(0.571428571428096, 0.1428571428570705, 0.0),
+    SVector(0.4285714285714073, 0.2857142857140783, 0.0),
+    SVector(0.2857142857145434, 0.4285714285713746, 0.0),
+    SVector(0.2857142857140778, 0.2857142857145756, 0.0)
+]
+
+face2 = [SVector(1, 2, 3, 16, 17, 18, 19, 20, 21, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
+                22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36)]
+
+m1 = CompScienceMeshes.Mesh(verts1, faces1)
+m2 = CompScienceMeshes.CurvilinearMesh(verts2, face2, 2, 7)
+
+ch1 = chart(m1, 1)
+ch2 = chart(m2, 1)
+
+@test barytocart(ch1, SVector(1.0, 0.0)) == SVector(1.0, 0.0, 0.0)
+@test barytocart(ch1, SVector(0.0, 1.0)) == SVector(0.0, 1.0, 0.0)
+@test barytocart(ch1, SVector(0.0, 0.0)) == SVector(0.0, 0.0, 0.0)
+@test barytocart(ch2, SVector(1.0, 0.0)) == SVector(1.0, 0.0, 0.0)
+@test barytocart(ch2, SVector(0.0, 1.0)) == SVector(0.0, 1.0, 0.0)
+@test barytocart(ch2, SVector(0.0, 0.0)) == SVector(0.0, 0.0, 0.0)
+
+mp1 = neighborhood(ch1, SVector(0.0, 0.0))
+mp2 = neighborhood(ch2, SVector(0.0, 0.0))
+
+@test volume(ch1) == 0.5
+@test volume(ch2) ≈ 0.5
+
+@test cartesian(mp1) == SVector(0.0, 0.0, 0.0)
+@test cartesian(mp2) == SVector(0.0, 0.0, 0.0) 
+
+@test tangents(mp1, 1) == SVector(1.0, 0.0, 0.0)
+@test tangents(mp1, 2) == SVector(0.0, 1.0, 0.0)
+
+CompScienceMeshes.unitarybase(ch2, SVector(0.0, 0.0))
+CompScienceMeshes.unitarybase(ch2, SVector(1.0, 0.0))
+CompScienceMeshes.unitarybase(ch2, SVector(0.0, 1.0))
+
+@test normal(ch2, SVector(0.0, 0.0)) == SVector(0.0, 0.0, 1.0)
+
+@test tangents(mp2, 1) ≈ SVector(1.0, 0.0, 0.0) # Sign is flipped
+@test tangents(mp2, 2) ≈ SVector(0.0, 1.0, 0.0) # Sign is flipped
+
+@test normal(mp1) == SVector(0.0, 0.0, 1.0)
+@test normal(mp2) == SVector(0.0, 0.0, 1.0)
+
+@test jacobian(mp1) == 1.0
+@test jacobian(mp2) ≈ 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, 1, 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)
+
+#=
+##
+using CompScienceMeshes
+using Test
+using StaticArrays
+fn = joinpath(dirname(@__FILE__),"triangle.msh")
+m = CompScienceMeshes.load_gmsh_mesh(fn,  order=5)
+
+m.faces
+
+# Test that we are not missing a node
+for p = 1:10
+    triplets = CompScienceMeshes.gmsh_triangle_index_to_triplet[p]
+    @test length(triplets) == div((p+1)*(p+2), 2)
+    #map = gmsh_triangle_index_to_triplet(p)
+    #for (i, t) in enumerate(triplets)
+    #    @test map[t] == i - 1  # Gmsh uses 0-based indexing
+    #end
+end
+
+function gmsh_triplet_to_index(p::Int)
+    triplets = gmsh_index_to_triplet(p)
+    map = Dict{Tuple{Int,Int,Int}, Int}()
+    for (i, t) in enumerate(triplets)
+        map[t] = i - 1  # Gmsh uses 0-based indexing
+    end
+    return map
+end
+
+=#
\ No newline at end of file