Feature/2d curvilinear support 3

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
 
@@ -91,6 +92,7 @@ Pt{N,T} = StaticArrays.SVector{N,T}
 include("defaults.jl")
 include("utils.jl")
 include("utils/localgraphs.jl")
+include("graphs/hosimplex.jl")
 include("utils/sfc_sort.jl")
 include("utils/circumctr.jl")
 # include("combinatorics.jl")
@@ -123,6 +125,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")

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}}
@@ -274,7 +275,8 @@ function barytocart(mani::Simplex, u)
     return r
 end
 
-function cartesian(ch::Simplex, u) barytocart(ch, u) end
+function cartesian(ch::AbstractSimplex, u) barytocart(ch, u) end
+# function cartesian(ch::Simplex, u::SVector{N,T}) where {N,T} barytocart(ch, u) end
 
 
 """
@@ -481,6 +483,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)
 

src/fileio/gmsh.jl

@@ -89,6 +89,10 @@ function load_gmsh_mesh(meshfile;
     if element == :quadrangle
         elements = [QuadrilateralGraph(el) for el in elements]
         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
         elements = [SimplexGraph(el) for el in elements]
         return Mesh(vertices, elements)

src/graphs/hosimplex.jl

@@ -0,0 +1,32 @@
+# N: the number of nodes (aka the dimension plus one)
+# P: the order of the simplex
+# Q: the number of vertices (depends on N and the order of the simplex)
+struct HOSimplexGraph{N,P,Q} <: AbstractGraph
+    indices::SVector{Q,Int}
+end
+
+function HOSimplexGraph{N,P}(idcs::Vararg{Int,Q}) where {N,P,Q}
+    @assert Q == binomial(N+P-1,P)
+    HOSimplexGraph{N,P,Q}(SVector{Q,Int}(idcs))
+end
+
+# edges of a triangle
+function faces(g::HOSimplexGraph{3,P,Q}, ::Type{Val{1}}) where {P,Q}
+
+    N=3
+
+    Nout = N-1
+    Pout = P
+
+    F = HOSimplexGraph{Nout,Pout}
+    I = g.indices
+
+
+    (
+        F(I[2], I[3+P:3+2P-2]..., I[3]),
+        F(I[3], I[3+2P-1:3+3P-3]..., I[1]),
+        F(I[1], I[3+1:3+P-1]..., I[2]),
+    )
+end
+
+

src/mesh.jl

@@ -45,7 +45,15 @@ function mesh(T, mdim, udim=mdim+1)
     Mesh(Pt{udim,T}[], CompScienceMeshes.SimplexGraph{mdim+1}[])
 end
 
+"""
+    similar(mesh::M, vertices, faces) where {M<:AbstractMesh}
 
+Constructs a mesh with the storage type as `mesh`, but with the order and dimension
+parameters inferred from the `faces` argument. The vertices are taken from `vertices`.
+"""
+function similar_mesh(mesh::Mesh{U,D1,T}, vertices, faces::Vector{SimplexGraph{D1new}}) where {U,D1,T,D1new}
+    Mesh{U,D1new,T}(vertices, faces)
+end
 
 """
     vt = vertextype(mesh)
@@ -525,7 +533,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
 
@@ -573,7 +581,7 @@ end
 # end
 
 
-function skeleton_fast(mesh, dimtype::Type)
+function skeleton_fast(mesh::M, dimtype::Type) where {M <: AbstractMesh}
     # C = celltype(mesh)
     # I = indextype(mesh, dimtype)
     F = celltype(mesh, dimtype)
@@ -584,7 +592,7 @@ function skeleton_fast(mesh, dimtype::Type)
         end
     end
     fcs = unique!(fcs)
-    Mesh(vertices(mesh), fcs)
+    similar_mesh(mesh, vertices(mesh), fcs)
 end
 
 function skeleton_fast(mesh, dim::Int)
@@ -595,12 +603,12 @@ function skeleton_fast(mesh, dim::Int)
 end
 
 
-"""
-    skeleton(pred, mesh, dim)
+# """
+#     skeleton(pred, mesh, dim)
 
-Like `skeleton(mesh, dim)`, but only cells for which `pred(cell)`
-returns true are withheld.
-"""
+# Like `skeleton(mesh, dim)`, but only cells for which `pred(cell)`
+# returns true are withheld.
+# """
 # function skeleton(pred, mesh, dim)
 
 #     meshdim = dimension(mesh)