diff --git a/ext/PlotlyBaseExt.jl b/ext/PlotlyBaseExt.jl index 672a81b..0cf2ee9 100644 --- a/ext/PlotlyBaseExt.jl +++ b/ext/PlotlyBaseExt.jl @@ -13,6 +13,7 @@ function CompScienceMeshes.patch(Γ::CompScienceMeshes.AbstractMesh, fcr=nothing x = v[:,1]; y = v[:,2]; z = v[:,3] i = c[:,1].-1; j = c[:,2].-1; k = c[:,3].-1 + # @show typeof(v) if fcr != nothing m, M = extrema(fcr) @@ -119,4 +120,19 @@ function PlotlyBase.scatter3d(surf::CompScienceMeshes.AbstractMesh{3,3}; kwargs. PlotlyBase.scatter3d(CompScienceMeshes.skeleton(surf, 1); kwargs...) end +function PlotlyBase.scatter3d(nodes::CompScienceMeshes.AbstractMesh{3,1}; kwargs...) + T = CompScienceMeshes.coordtype(nodes) + x = T[] + y = T[] + z = T[] + for node in nodes + chrt = CompScienceMeshes.chart(nodes, node) + v1, = chrt.vertices + push!(x, v1[1]) + push!(y, v1[2]) + push!(z, v1[3]) + end + return PlotlyBase.scatter3d(x=x, y=y, z=z, mode="markers"; kwargs...) +end + end \ No newline at end of file diff --git a/src/CompScienceMeshes.jl b/src/CompScienceMeshes.jl index c267744..e931a9d 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 @@ -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") @@ -115,6 +117,7 @@ include("rectangle.jl") include("charts.jl") include("charts/quadrilateral.jl") include("charts/hexahedron.jl") +include("charts/curvilinear_chart.jl") include("subd_chart.jl") include("sphere.jl") include("overlap.jl") @@ -123,6 +126,13 @@ 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_neighborhood.jl") + include("quadpoints.jl") include("submesh.jl") diff --git a/src/baryref.jl b/src/baryref.jl index a1021b1..15a6e85 100644 --- a/src/baryref.jl +++ b/src/baryref.jl @@ -1,5 +1,39 @@ abstract type AbstractRefinement{U,D1,T} <: AbstractMesh{U,D1,T} end +function ==(m1::M, m2::M) where {M<:AbstractRefinement} m1.mesh == m2.mesh && m1.parent == m2.parent end +vertextype(br::AbstractRefinement) = vertextype(br.mesh) +vertices(br::AbstractRefinement) = vertices(br.mesh) +vertices(m::AbstractRefinement, cell) = vertices(m.mesh,cell) +numvertices(br::AbstractRefinement) = numvertices(br.mesh) +cells(br::AbstractRefinement) = cells(br.mesh) +indices(br::AbstractRefinement, p) = indices(br.mesh, p) +numcells(br::AbstractRefinement) = numcells(br.mesh) +skeleton(br::AbstractRefinement, dim::Int) = skeleton(br.mesh, dim) +chart(br::AbstractRefinement, cell) = chart(br.mesh, cell) +vertexarray(m::AbstractRefinement) = vertexarray(m.mesh) +cellarray(m::AbstractRefinement) = cellarray(m.mesh) +function celltype(m::AbstractRefinement{U,D1}) where {U,D1} celltype(m.mesh) end +function celltype(m::AbstractRefinement{U,D1}, ::Type{Val{M}}) where {U,D1,M} celltype(m.mesh, Val{M}) end +function indextype(m::AbstractRefinement{U,D1}) where {U,D1} indextype(m.mesh) end +function indextype(m::AbstractRefinement{U,D1}, ::Type{Val{M}}) where {U,D1,M} indextype(m.mesh, Val{M}) end + + +struct RefinementParent{U,D1,T,M,N} <: AbstractMesh{U,D1,T} + mesh::M + refinement::N +end + +Base.:(==)(m1, m2::RefinementParent) = (m1 == m2.mesh) +Base.:(==)(m1::RefinementParent, m2) = m1.mesh == m2 + +RefinementParent{U,D1,T}(mesh::M, parent::N) where {U,D1,T,M,N} = RefinementParent{U,D1,T,M,N}(mesh, parent) + +parent(m::AbstractRefinement{U,D1,T}) where {U,D1,T} = RefinementParent{U,D1,T}(m.parent, m) +Base.length(m::RefinementParent) = length(m.mesh) +vertices(m::RefinementParent, p) = vertices(m.mesh, p) +cells(m::RefinementParent) = cells(m.mesh) +function chart(m::RefinementParent, ptr) chart(m.mesh, ptr) end + """ True if m1 is a direct refinement of m2. """ @@ -17,48 +51,13 @@ struct BarycentricRefinement{U,D1,T,M,MP} <: AbstractRefinement{U,D1,T} mesh::M parent::MP end -function ==(m1::BarycentricRefinement{U,D1,T,M,MP}, m2::BarycentricRefinement{U,D1,T,M,MP}) where {U,D1,T,M,MP} - return m1.mesh == m2.mesh && m1.parent == m2.parent -end -vertextype(br::BarycentricRefinement) = vertextype(br.mesh) -vertices(br::BarycentricRefinement) = vertices(br.mesh) -vertices(m::BarycentricRefinement, cell) = vertices(m.mesh,cell) -numvertices(br::BarycentricRefinement) = numvertices(br.mesh) -cells(br::BarycentricRefinement) = cells(br.mesh) -indices(br::BarycentricRefinement, p) = indices(br.mesh, p) -numcells(br::BarycentricRefinement) = numcells(br.mesh) -skeleton(br::BarycentricRefinement, dim::Int) = skeleton(br.mesh, dim) -chart(br::BarycentricRefinement, cell) = chart(br.mesh, cell) -vertexarray(m::BarycentricRefinement) = vertexarray(m.mesh) -cellarray(m::BarycentricRefinement) = cellarray(m.mesh) - -function celltype(m::BarycentricRefinement{U,D1}) where {U,D1} SimplexGraph{D1} end -function celltype(m::BarycentricRefinement{U,D1}, ::Type{Val{M}}) where {U,D1,M} SimplexGraph{M+1} end -function indextype(m::BarycentricRefinement{U,D1}) where {U,D1} SVector{D1,Int} end -function indextype(m::BarycentricRefinement{U,D1}, ::Type{Val{M}}) where {U,D1,M} SVector{M+1,Int} end function parent(m::BarycentricRefinement, p) numchildren = factorial(dimension(m)+1) return mod1(p, numchildren) end -struct BarycentricRefinementParent{U,D1,T,M,N} <: AbstractMesh{U,D1,T} - mesh::M - refinement::N -end - -Base.:(==)(m1, m2::BarycentricRefinementParent) = (m1 == m2.mesh) -Base.:(==)(m1::BarycentricRefinementParent, m2) = m1.mesh == m2 - -BarycentricRefinementParent{U,D1,T}(mesh::M, parent::N) where {U,D1,T,M,N} = BarycentricRefinementParent{U,D1,T,M,N}(mesh, parent) - -parent(m::BarycentricRefinement{U,D1,T}) where {U,D1,T} = BarycentricRefinementParent{U,D1,T}(m.parent, m.mesh) -Base.length(m::BarycentricRefinementParent) = length(m.mesh) -vertices(m::BarycentricRefinementParent, p) = vertices(m.mesh, p) -cells(m::BarycentricRefinementParent) = cells(m.mesh) - - -function children(m::BarycentricRefinementParent, p) +function children(m::RefinementParent{<:Any,<:Any,<:Any,<:Any,<:BarycentricRefinement}, p) numchildren = factorial(dimension(m)+1) [(p-1)*numchildren+i for i in 1:numchildren] end @@ -310,7 +309,62 @@ parent(mesh::BarycentricRefinement{U,3} where {U}, cell) = div((cell-1),6)+1 parent(mesh::BarycentricRefinement{U,2} where {U}, cell) = div((cell-1),2)+1 +struct BisectingRefinement{U,D1,T,M,MP} <: AbstractRefinement{U,D1,T} + order::Int + mesh::M + parent::MP +end + +function parent(m::BisectingRefinement, p) + order = m.order + numchildren = (div(order * (order+1), 2) + div((order-1) * order,2)) + return mod1(p, numchildren) +end + +function children(m::RefinementParent{<:Any,<:Any,<:Any,<:Any,<:BisectingRefinement}, p) + order = m.refinement.order + numchildren = (div(order * (order+1), 2) + div((order-1) * order,2)) + [(p-1)*numchildren+i for i in 1:numchildren] +end + +function injections(m::RefinementParent{U,D1,T,<:Any,<:BisectingRefinement}, p) where {U,D1,T} + + order = m.refinement.order + numchildren = (div(order * (order+1), 2) + div((order-1) * order,2)) + + h = 1 / order + # set1 = [ + # simplex(point(i*h, j*h), point((i+1)*h, j*h), point(i*h, (j+1)*h)) + # for j in 0:order-1 for i in 0:order-1-j + # ] + # set2 = [ + # simplex(point((i+1)*h, j*h), point((i+1)*h, (j+1)*h), point(i*h, (j+1)*h)) + # for j in 0:order-2 for i in 0:order-2-j + # ] + # set = [set1; set2] + + + S = Simplex{2,D1-1,2-D1+1,D1,T} + fcs = Vector{S}() + for m in 0:order-1 + for n in 0:order-1-m + push!(fcs, simplex(point(T,m*h,n*h), point(T,(m+1)*h, n*h), point(T, m*h, (n+1)*h))) + if m + n <= order - 2 + push!(fcs, simplex(point(T,(m+1)*h, n*h), point(T,(m+1)*h, (n+1)*h), point(T, m*h, (n+1)*h))) + end + end + end + + # @show length(fcs) + # @show numchildren + # @assert length(fcs) == numchildren + return fcs +end + """ + bisecting_refinement(mesh::AbstractMesh{U,3}; order=2) where U + + bisecting_refinement(mesh) -> refinement Construct a refinement of `mesh` by inserting a new vertex on every existing edge. @@ -318,64 +372,116 @@ Every face is subdived in four small faces. Only defined for 2D meshes. """ -function bisecting_refinement(mesh::Mesh{U,3}) where U +function bisecting_refinement(mesh::AbstractMesh{U,3,T}; order=2) where {U,T} - D1 = 3 + MP = typeof(mesh) edges = skeleton(mesh,1) - faces = skeleton(mesh,2) + # faces = skeleton(mesh,2) + faces = mesh NV, NE, NF = numvertices(mesh), numcells(edges), numcells(faces) - nv = NV + NE - verts = Array{vertextype(mesh)}(undef,nv) + verts = deepcopy(vertices(mesh)) - # copy over the existing vertices - for V in 1 : numvertices(mesh) - verts[V] = mesh.vertices[V] + for (i, E) in enumerate(edges) + ch = chart(edges, E) + for n in 1:order-1 + p = neighborhood(ch, SVector(1 - n/order)) + push!(verts, cartesian(p)) + end end - # add a vertex in each edge centroid - # for (E, Edge) in enumerate(cells(edges)) - for (i,E) in enumerate(edges) - verts[NV+i] = cartesian(center(chart(edges, E))) + idx = Dict{Tuple{Int,Int}, Int}() + k = 1 + for n in 1:order-2 + for m in 1:order-1-n + idx[(n,m)] = k + k += 1 + end + end + + for F in faces + ch = chart(faces, F) + for n in 1:order-2 + for m in 1:order-1-n + p = neighborhood(ch, SVector(n/order,m/order)) + push!(verts, cartesian(p)) + end + end end # Build a matrix that given a coarse face gives # access to the coarse edges making up its boundary - D = copy(transpose(connectivity(edges, faces, identity))) + D = connectivity(faces, edges, identity) rows, vals = rowvals(D), nonzeros(D) # add four faces in each coarse face - nf = 4NF - # fcs = zeros(celltype(mesh), nf) - fcs = Vector{celltype(mesh)}(undef, nf) + # nf = (div(order * (order+1), 2) + div((order-1) * order,2)) * NF + S = SimplexGraph{3} + fcs = Vector{S}() for F in 1 : numcells(faces) - # retrieve the indices of the fine vertices in the centroids - # of the edges in the boundary of coarse face F - es = zeros(Int,3) + # v = zeros(Int, order+1, order+1) + v = fill(-200, order+1, order+1) + + # retrieve the indices of the fine vertices in the nodes of face F + v[order+1,1] = cells(faces)[F][1] + v[1,order+1] = cells(faces)[F][2] + v[1,1] = cells(faces)[F][3] + + # retrieve the indices of the fine vertices along the edges of face F + # e = zeros(Int,3, order-1) for i in nzrange(D, F) E = rows[i] - j = abs(vals[i]) # local index of edge E in face F - @assert 1 <= j <= 3 - e = NV + E - es[j] = e + j, σ = abs(vals[i]), sign(vals[i]) # j = +/-1, +/-2, +/-3 + if j == 1 + for n in 1:order-1 + v[1, order-n+1] = (σ >0) ? + NV + (E-1) * (order-1) + n : + NV + (E-1) * (order-1) + (order-n) + end + elseif j == 2 + for n in 1:order-1 + v[n+1, 1] = (σ >0) ? + NV + (E-1) * (order-1) + n : + NV + (E-1) * (order-1) + (order-n) + end + elseif j == 3 + for n in 1:order-1 + v[order-n+1, n+1] = (σ >0) ? + NV + (E-1) * (order-1) + n : + NV + (E-1) * (order-1) + (order-n) + end + end end + # @show e - # build the four triangles - r1 = mesh.faces[F][1] - r2 = mesh.faces[F][2] - r3 = mesh.faces[F][3] + # retrieve the indices of the fine vertices in the interior of face F + for n in 1:order-2 + for m in 1:order-1-n + k = idx[(n,m)] + v[n+1,m+1] = NV + (NE * (order-1)) + (F-1) * div((order-1)*(order-2),2) + k + end + end - fcs[4(F-1)+1] = SimplexGraph(r1,es[3],es[2]) - fcs[4(F-1)+2] = SimplexGraph(r2,es[1],es[3]) - fcs[4(F-1)+3] = SimplexGraph(r3,es[2],es[1]) - fcs[4(F-1)+4] = SimplexGraph(es[1],es[2],es[3]) + # @show v + + # Build the fine faces in the coarse face F + for n in 1:order + for m in 1:order-n+1 + push!(fcs, S(SVector(v[n,m], v[n+1,m], v[n,m+1]))) + if m + n <= order + push!(fcs, S(SVector(v[n+1,m], v[n+1,m+1], v[n,m+1]))) + end + end + end end - Mesh(verts, fcs) + fine = Mesh(verts, fcs) + MR = typeof(fine) + BisectingRefinement{U,3,T,MR,MP}(order, fine, mesh) end """ diff --git a/src/charts.jl b/src/charts.jl index beb761d..8d3d79e 100644 --- a/src/charts.jl +++ b/src/charts.jl @@ -1,12 +1,15 @@ #export Simplex +abstract type AbstractSimplex{U,D} +end +function nodes(chart) vertices(chart) 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}} @@ -28,7 +31,7 @@ function permute_vertices(s::Simplex{3,2,1,3,T},permutation::Union{Vector{P},SVe Simplex(simp.vertices,simp.tangents,s.normals,simp.volume) end -""" +"""F flip_normal(simplex) Flips the normal of the simplex. Only on triangles embedded in 3D space @@ -274,7 +277,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 +485,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/charts/curvilinear_chart.jl b/src/charts/curvilinear_chart.jl new file mode 100644 index 0000000..870d00e --- /dev/null +++ b/src/charts/curvilinear_chart.jl @@ -0,0 +1,434 @@ +#### 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{<:Any,1}) SVector(s.vertices[begin], s.vertices[begin + 1]) end +function nodes(s::CurvilinearSimplex{<:Any,2}) SVector(s.vertices[begin], s.vertices[begin + 1], s.vertices[begin + 2]) end + +# convenience +function vertices(s::CurvilinearSimplex) s.vertices end + +""" + 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 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`. +""" +function center(ch::CurvilinearSimplex{U,0}) where {U} ch.vertices[1] end +function center(ch::CurvilinearSimplex{U,1}) where {U} neighborhood(ch, SVector(0.5)) end +function center(ch::CurvilinearSimplex{U,2}) where {U} neighborhood(ch, SVector(0.5, 0.5)) end + +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 + +function faces(ch::CurvilinearSimplex{U,2,C,N,T,O}, ::Type{Val{1}}) where {U,C,N,T,O} + + D = 2 + @assert N == binomial(D+O,O) + P = vertextype(ch) + + Nout = binomial(1+O,O) + Cout = U - 1 + + I = ch.vertices + F = SVector{Nout, P} + E = CurvilinearSimplex{U,1,Cout,Nout,T,O} + + ( + E(F(I[2], I[3], I[3+1*(O-1)+1:3+2*(O-1)]...)), + E(F(I[3], I[1], I[3+2*(O-1)+1:3+3*(O-1)]...)), + E(F(I[1], I[2], I[3+0*(O-1)+1:3+1*(O-1)]...)), + ) +end \ No newline at end of file diff --git a/src/fileio/gmsh.jl b/src/fileio/gmsh.jl index 877c46d..5fb500f 100644 --- a/src/fileio/gmsh.jl +++ b/src/fileio/gmsh.jl @@ -22,6 +22,8 @@ function load_gmsh_mesh(meshfile; gmsh.initialize() + @show meshfile + try gmsh.open(meshfile) catch e @@ -35,6 +37,8 @@ function load_gmsh_mesh(meshfile; true, # includeBoundary (default false) true # returnParametricCoord (default true) ) + # @show node_tags + # @show node_coords if physical===nothing element_types, element_tags, element_node_tags = gmsh.model.mesh.getElements() @@ -85,10 +89,16 @@ function load_gmsh_mesh(meshfile; vertices = vertices2d end + # @show elements + # TODO: once we support hexahedrons, we need to catch that one too 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) diff --git a/src/graphs/hosimplex.jl b/src/graphs/hosimplex.jl new file mode 100644 index 0000000..442dc55 --- /dev/null +++ b/src/graphs/hosimplex.jl @@ -0,0 +1,68 @@ +# 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 + + +function faces(g::HOSimplexGraph{3,P,Q}, ::Type{Val{2}}) where {P,Q} + ( + g, + ) +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], I[3+1*(P-1)+1:3+2*(P-1)]...), + F(I[3], I[1], I[3+2*(P-1)+1:3+3*(P-1)]...), + F(I[1], I[2], I[3+0*(P-1)+1:3+1*(P -1)]...), + ) +end + +# nodes of a triangle +function faces(g::HOSimplexGraph{3,P,Q}, ::Type{Val{0}}) where {P,Q} + + N=3 + + Nout = 1 + Pout = P + + F = HOSimplexGraph{Nout,Pout} + I = g.indices + + + ( + F(I[1]), + F(I[2]), + F(I[3]), + ) +end + +function indices(g::HOSimplexGraph{N,P,Q}) where {N,P,Q} + SVector(ntuple(i -> g.indices[i], Val(N))) +end + + +function relorientation( + g::HOSimplexGraph{N1,P1,Q1}, + f::HOSimplexGraph{N2,P2,Q2}) where {N1,P1,Q1,N2,P2,Q2} + + relorientation(indices(g), indices(f)) +end \ No newline at end of file diff --git a/src/mesh.jl b/src/mesh.jl index 65e0460..1e40321 100644 --- a/src/mesh.jl +++ b/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) @@ -362,13 +370,15 @@ function boundary(mesh) nzr = nzrange(conn,e) length(nzr) != 1 && continue relop = vals[nzr[1]] - inds = indices(edges, edge) + # inds = indices(edges, edge) + inds = cells(edges)[e].indices bnd_edges[i] = (relop > 0) ? I(inds) : flip(I(inds)) i += 1 end resize!(bnd_edges, i-1) - bnd = Mesh(vertices(mesh), bnd_edges) + # bnd = Mesh(vertices(mesh), bnd_edges) + return similar_mesh(mesh, vertices(mesh), bnd_edges) end @@ -401,11 +411,14 @@ end Computed an V×M array `vertextocells` where V is the number of vertices and M is the maximum number of cells adjacent to any given vertex such that `vertextocells[v,i]` is the index in the cells of `mesh` of the `i`th -cell adjacent to teh `v`-th vertex. `numneighbors[v]` contains the number -of cells adjacent to the `v`-th vertex. +cell adjacent to the `v`-th vertex. Rows are indexed directly by vertex ID, +including vertices that are not referenced by any cell. Unused entries are +padded with `-1`. `numneighbors[v]` contains the number of cells adjacent +to the `v`-th vertex. This method allows e.g. for the efficient computation of the connectivity -matrix of the mesh. +matrix of the mesh. Use [`vertextocell`](@ref) instead when many vertices +are unreferenced and compact rows are preferred. """ function vertextocellmap(mesh) @@ -435,6 +448,20 @@ function vertextocellmap(mesh) end +""" + vertextocell(mesh) -> vertextocells, numneighbors, LG, GL + +Returns a compact vertex-to-cell adjacency map. `vertextocells` has one row +for each vertex referenced by a cell, and `numneighbors[row]` gives the +number of adjacent cells in that row. Entries beyond that count are padded +with `-1`. Rows are ordered by the first encounter of each vertex while +iterating over the cells. + +The rows are not indexed by vertex ID: `LG[row]` gives the vertex ID for a +row, and `GL[vertex_id]` gives the row for a referenced vertex. Unlike +[`vertextocellmap`](@ref), this does not allocate rows for unreferenced +vertices, so it can use less memory when many vertices are unreferenced. +""" function vertextocell(mesh) row = Dict{Int,Tuple{Int,Int}}() @@ -525,7 +552,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 @@ -534,13 +561,12 @@ function skeleton(mesh, dim::Int; sort=:spacefillingcurve) # sort the simplices on a SFC simplices = cells(sk) - # ctrs = [sum(cellvertices(sk,c))/(dim+1) for c in sk] - ctrs = [cartesian(center(chart(sk, cell))) for cell in simplices] + ctrs = [cartesian(center(chart(sk, cellptr))) for cellptr in sk] if length(sk) > 0 simplices = simplices[sort_sfc(ctrs)] end - Mesh(vertices(mesh), simplices) + similar_mesh(mesh, vertices(mesh), simplices) end @@ -573,7 +599,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 +610,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 +621,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) @@ -646,18 +672,21 @@ positive or negative depending on the relative orientation of face `k` in cell For `op=sign`, the matrix returned is the classic connectivity matrix, i.e. the graph version of the exterior derivative. """ -function connectivity(kcells::AbstractMesh, mcells::AbstractMesh, op = sign) +function connectivity(kmesh::AbstractMesh, mmesh::AbstractMesh, op = sign) - K = celltype(kcells) - M = celltype(mcells) + K = celltype(kmesh) + M = celltype(mmesh) - vtok, _, lgk, _ = vertextocell(kcells) - vtom, _, _, glm = vertextocell(mcells) + vtok, _, lgk, _ = vertextocell(kmesh) + vtom, _, _, glm = vertextocell(mmesh) npos = -1 - dimk = numcells(kcells) - dimm = numcells(mcells) + dimk = length(kmesh) + dimm = length(mmesh) + + kcells = cells(kmesh) + mcells = cells(mmesh) Rows = Int[] Cols = Int[] @@ -669,12 +698,15 @@ function connectivity(kcells::AbstractMesh, mcells::AbstractMesh, op = sign) for q in axes(vtok,2) i = vtok[vk,q] i == npos && break - kcell = K(indices(kcells,i)) + kcell = kcells[i] for s in axes(vtom,2) j = vtom[vm,s] j == npos && break - mcell = M(indices(mcells, j)) + mcell = mcells[j] val = op(relorientation(kcell, mcell)) + # @show indices(kcell) + # @show indices(mcell) + # @show val iszero(val) && continue push!(Rows, j) push!(Cols, i) 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..770d176 --- /dev/null +++ b/src/meshes/curvilinear_mesh.jl @@ -0,0 +1,153 @@ +# 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,D1,T} + vertices::Vector{SVector{U,T}} # coordinates in R^U + # faces::Vector{SVector{N,Int}} # connectivity (local order = Gmsh order) + faces::Vector{HOSimplexGraph{D1,O,N}} # connectivity (local order = Gmsh order) + + +end + +function CurvilinearMesh( + vertices::AbstractVector{SVector{U,T}}, + # faces::AbstractVector{SVector{N,Int}}, + faces::AbstractVector{HOSimplexGraph{D1,O,N}}, + manifolddimension, + order::Integer +) where {U,N,T,D1,O} + order > 0 || throw(ArgumentError("order must be positive.")) + order == O || throw(ArgumentError("order inconsistent with cell type.")) + manifolddimension + 1 == D1 || throw(ArgumentError("manifold dimension inconsistent with cell type.")) + CurvilinearMesh{U,D1,N,T,O}(vertices, faces) +end + +function CurvilinearMesh( + vertices::AbstractVector{SVector{U,T}}, + faces::AbstractVector{SVector{N,Int}}, + # faces::AbstractVector{HOSimplexGraph{D1,O,N}}, + manifolddimension, + order::Integer +) where {U,N,T} + order > 0 || throw(ArgumentError("order must be positive.")) + # order == O || throw(ArgumentError("order inconsistent with cell type.")) + # manifolddimension + 1 == D1 || throw(ArgumentError("manifold dimension inconsistent with cell type.")) + D1 = manifolddimension + 1 + O = order + faces = [HOSimplexGraph{D1,O,N}(f) for f in faces] + CurvilinearMesh{U,D1,N,T,O}(vertices, faces) +end + +# This constructor acts like `similar` for arrays. It allows to override the type +# parameters of the mesh, without referring to the abstract parent type. +function similar_mesh( + mesh::CurvilinearMesh{U,D1,N,T,O}, + vertices::AbstractVector{SVector{U,T}}, + faces::AbstractVector{HOSimplexGraph{D1new,O,Nnew}}, +) where {U,N,T,D1,O,D1new,Nnew} + CurvilinearMesh{U,D1new,Nnew,T,O}(vertices, faces) +end + + +function CurvilinearMesh( + vertices::AbstractVector{SVector{U,T}}, + faces::AbstractVector{HOSimplexGraph{D1,O,N}} +) where {U,N,T,D1,O} + CurvilinearMesh(vertices, faces, D1-1, O) +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 vertexarray(m::CurvilinearMesh) + [v[i] for v in m.vertices, i in 1:universedimension(m)] +end +function cellarray(m::CurvilinearMesh) + [ k[i] for k in m.faces, i in 1:dimension(m)+1 ] +end + +function celltype(m::CurvilinearMesh{U,D1,N,T,O}) where {U,D1,N,T,O} + # SimplexGraph{N-O+1} + HOSimplexGraph{D1,O,N} +end + +function celltype(m::CurvilinearMesh{U,D1,N,T,O}, ::Type{Val{M}}) where {U,D1,N,T,O,M} + # SimplexGraph{M+1} + HOSimplexGraph{M+1,O,binomial(M+O,O)} +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,1,N,T,O}, cellptr) where {U,N,T,O} + cell = m.faces[cellptr] + return SVector(cell.indices[begin]) +end + +function indices(m::CurvilinearMesh{U,2,N,T,O}, cellptr) where {U,N,T,O} + cell = m.faces[cellptr] + return SVector(cell.indices[begin], cell.indices[begin+1]) +end + +function indices(m::CurvilinearMesh{U,3,N,T,O}, cellptr) where {U,N,T,O} + cell = m.faces[cellptr] + return SVector(cell.indices[begin], cell.indices[begin+1], cell.indices[begin+2]) +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)) +function cells(m::CurvilinearMesh) m.faces end +function cell(m::CurvilinearMesh, i::Int) m.faces[i] end + +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 + +function chart( + m::CompScienceMeshes.CurvilinearMesh{U,D1,N,T,O}, + faceid::Int +) where {U,D1,N,T,O} + + vindices = m.faces[faceid].indices + X = m.vertices[vindices] + + CompScienceMeshes.CurvilinearSimplex{U,D1-1,U-1,N,T,O}(X) +end \ No newline at end of file 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/meshes/quadrilateralmesh.jl b/src/meshes/quadrilateralmesh.jl index 31d59ee..8565d87 100644 --- a/src/meshes/quadrilateralmesh.jl +++ b/src/meshes/quadrilateralmesh.jl @@ -11,6 +11,10 @@ function QuadMesh(V::Vector{SVector{3,T}}, F::Vector{SVector{4,Int}}) where {T} return QuadMesh{T}(V, [QuadrilateralGraph(f) for f in F]) end +function similar_mesh(mesh::QuadMesh{T}, vertices::Vector{SVector{U,T}}, faces::Vector{SimplexGraph{2}}) where {U,T} + Mesh{U,2,T}(vertices, faces) +end + function indices(m::QuadMesh, cellptr) m.faces[cellptr].indices end function vertextype(m::QuadMesh) eltype(m.vertices) end function dimension(m::QuadMesh) 2 end diff --git a/src/primitives/gmsh.jl b/src/primitives/gmsh.jl index e6c66ad..ddb1e3e 100644 --- a/src/primitives/gmsh.jl +++ b/src/primitives/gmsh.jl @@ -4,7 +4,7 @@ using GmshTools meshgeo(geofile; physical=nothing, dim=2, tempname=tempname(), kwargs...) """ -function meshgeo(geofile; physical=nothing, dim=2, tempname=tempname(), kwargs...) +function meshgeo(geofile; physical=nothing, dim=2, order=1, tempname=tempname(), kwargs...) io = open(geofile) str = read(io, String) @@ -31,6 +31,7 @@ function meshgeo(geofile; physical=nothing, dim=2, tempname=tempname(), kwargs.. gmsh.option.setNumber("Mesh.MshFileVersion",2) gmsh.open(temp_geo) gmsh.model.mesh.generate(dim) + order > 1 && gmsh.model.mesh.setOrder(order) gmsh.write(temp_msh) gmsh.finalize() @@ -38,9 +39,11 @@ function meshgeo(geofile; physical=nothing, dim=2, tempname=tempname(), kwargs.. # run(`gmsh $temp_geo -2 -format msh2 -o $temp_msh`) if dim == 2 - m = read_gmsh_mesh(temp_msh, physical=physical) + # m = read_gmsh_mesh(temp_msh, physical=physical) + m = load_gmsh_mesh(temp_msh; element=:triangle, physical, order) elseif dim == 3 - m = read_gmsh3d_mesh(temp_msh, physical=physical) + # m = read_gmsh3d_mesh(temp_msh, physical=physical) + m = load_gmsh_mesh(temp_msh; element=:tetrahedron, physical, order) else error("gmsh files of dimension $(dim) cannot be read.") end diff --git a/src/primitives/linemeshes/linemeshes.jl b/src/primitives/linemeshes/linemeshes.jl index 902527e..0609f60 100644 --- a/src/primitives/linemeshes/linemeshes.jl +++ b/src/primitives/linemeshes/linemeshes.jl @@ -116,6 +116,64 @@ function meshrect(a::T, b::T, delta::T) where T<:Real 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 + """ meshcurve(curve, delta::T) where T<:Real 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/src/submesh.jl b/src/submesh.jl index a1cf80f..7ac1add 100644 --- a/src/submesh.jl +++ b/src/submesh.jl @@ -55,6 +55,11 @@ function celltype(m::SubMesh{U,D1}, ::Type{Val{M}}) where {U,D1,M} SimplexGraph{ function indextype(m::SubMesh{U,D1}) where {U,D1} SVector{D1,Int} end function indextype(m::SubMesh{U,D1}, ::Type{Val{M}}) where {U,D1,M} SVector{M+1,Int} end +function similar_mesh(mesh::SubMesh{U,D1,T,S}, vertices, faces::Vector{SimplexGraph{D1new}} +) where {U,D1,T,S,D1new} + similar_mesh(mesh.supermesh, vertices, faces) +end + function SubMesh(supermesh, sub2sup) diff --git a/test/charts/curvilinear_chart.jl b/test/charts/curvilinear_chart.jl new file mode 100644 index 0000000..76b2010 --- /dev/null +++ b/test/charts/curvilinear_chart.jl @@ -0,0 +1,50 @@ +@testitem "chart and graph consistency [p=2]" begin + using SparseArrays + + const CSM = CompScienceMeshes + + verts = [ + point(1.0, 0.0, 0.0), + point(0.0, 1.0, 0.0), + point(0.0, 0.0, 0.0), + point(0.0, 0.5000000000013304, 0.0), + point(0.4999999999986718, 0.0, 0.0), + point(0.5000000000013293, 0.4999999999986707, 0.0), + + point(1.0, 1.0, 0.0), + point(1.0, 0.5, 0.0), + point(0.5, 1.0, 0.0), + + ] + + I = CSM.HOSimplexGraph{3,2} + faces = [ + I(1, 2, 3, 6, 4, 5), + I(1, 7, 2, 8, 9, 6), + ] + + mesh = CSM.CurvilinearMesh(verts, faces, 2, 2) + edges = CSM.skeleton_fast(mesh, 1) + nodes = CSM.skeleton_fast(mesh, 0) + + L = connectivity(nodes, edges, abs) + S = connectivity(mesh, edges, abs) + + ch1 = CSM.chart(mesh, 1) + cl1 = CSM.cell(mesh, 1) + + nnzs = nonzeros(S) + rows = rowvals(S) + for k = nzrange(S, 1)[1] + local_idx = nnzs[k] + global_idx = rows[k] + + fc1 = CSM.faces(ch1, Val{1})[local_idx] + fc2 = CSM.chart(edges, global_idx) + + x1 = CSM.center(fc1) |> cartesian + x2 = CSM.center(fc2) |> cartesian + + @test x1 ≈ x2 + end +end \ No newline at end of file diff --git a/test/meshes/curvilinear_mesh.jl b/test/meshes/curvilinear_mesh.jl new file mode 100644 index 0000000..ceb7a26 --- /dev/null +++ b/test/meshes/curvilinear_mesh.jl @@ -0,0 +1,220 @@ +@testitem "curvy mesh: skeleton [p=2]" begin + const CSM = CompScienceMeshes + + verts = [ + point(1.0, 0.0, 0.0), + point(0.0, 1.0, 0.0), + point(0.0, 0.0, 0.0), + point(0.0, 0.5000000000013304, 0.0), + point(0.4999999999986718, 0.0, 0.0), + point(0.5000000000013293, 0.4999999999986707, 0.0), + + point(1.0, 1.0, 0.0), + point(1.0, 0.5, 0.0), + point(0.5, 1.0, 0.0), + + ] + + I = CSM.HOSimplexGraph{3,2} + faces = [ + I(1, 2, 3, 6, 4, 5), + I(1, 7, 2, 8, 9, 6), + ] + + mesh = CSM.CurvilinearMesh(verts, faces, 2, 2) + + edges = CSM.skeleton_fast(mesh, 1) + @test length(edges) == 5 + + nodes = CSM.skeleton_fast(mesh, 0) + @test length(nodes) == 4 +end + + +@testitem "curvy mesh: connectivity [p=2]" begin + const CSM = CompScienceMeshes + + verts = [ + point(1.0, 0.0, 0.0), + point(0.0, 1.0, 0.0), + point(0.0, 0.0, 0.0), + point(0.0, 0.5000000000013304, 0.0), + point(0.4999999999986718, 0.0, 0.0), + point(0.5000000000013293, 0.4999999999986707, 0.0), + + point(1.0, 1.0, 0.0), + point(1.0, 0.5, 0.0), + point(0.5, 1.0, 0.0), + + ] + + I = CSM.HOSimplexGraph{3,2} + faces = [ + I(1, 2, 3, 6, 4, 5), + I(1, 7, 2, 8, 9, 6), + ] + + mesh = CSM.CurvilinearMesh(verts, faces, 2, 2) + edges = CSM.skeleton_fast(mesh, 1) + nodes = CSM.skeleton_fast(mesh, 0) + + L = connectivity(nodes, edges, identity) + S = connectivity(mesh, edges, identity) + + @test L == [ + 0 1 -2 0; + -2 0 1 0; + 1 -2 0 0; + 0 -2 0 1; + 1 0 0 -2; + ] + + @test S == [ + 1 0; + 2 0; + 3 -2; + 0 1; + 0 3; + ] +end + + +@testitem "curvy mesh boundary [p=2]" begin + const CSM = CompScienceMeshes + + verts = [ + point(1.0, 0.0, 0.0), + point(0.0, 1.0, 0.0), + point(0.0, 0.0, 0.0), + point(0.0, 0.5000000000013304, 0.0), + point(0.4999999999986718, 0.0, 0.0), + point(0.5000000000013293, 0.4999999999986707, 0.0), + + point(1.0, 1.0, 0.0), + point(1.0, 0.5, 0.0), + point(0.5, 1.0, 0.0), + + ] + + I = CSM.HOSimplexGraph{3,2} + faces = [ + I(1, 2, 3, 6, 4, 5), + I(1, 7, 2, 8, 9, 6), + ] + + mesh = CSM.CurvilinearMesh(verts, faces, 2, 2) + boundary = CSM.boundary(mesh) + + @test dimension(boundary) == 1 + @test length(boundary) == 4 +end + + +@testitem "curvy mesh bisecting refinement [rectangle, p=2]" begin + const CSM = CompScienceMeshes + + verts = [ + point(1.0, 0.0, 0.0), + point(0.0, 1.0, 0.0), + point(0.0, 0.0, 0.0), + point(0.0, 0.5000000000013304, 0.0), + point(0.4999999999986718, 0.0, 0.0), + point(0.5000000000013293, 0.4999999999986707, 0.0), + + point(1.0, 1.0, 0.0), + point(1.0, 0.5, 0.0), + point(0.5, 1.0, 0.0), + + ] + + I = CSM.HOSimplexGraph{3,2} + faces = [ + I(1, 2, 3, 6, 4, 5), + I(1, 7, 2, 8, 9, 6), + ] + + mesh = CSM.CurvilinearMesh(verts, faces, 2, 2) + @show typeof(mesh) + @test mesh isa CSM.AbstractMesh{3,3} + fine = CSM.bisecting_refinement(mesh; order=3) + @show typeof(fine) + + @test length(skeleton(mesh,0)) - length(skeleton(mesh,1)) + length(mesh) == 1 + @test length(skeleton(fine,0)) - length(skeleton(fine,1)) + length(fine) == 1 + + @test dimension(fine) == 2 + @test length(fine) == 18 +end + + +@testitem "curvy mesh bisecting refinement [sphere, p=2]" begin + const CSM = CompScienceMeshes + + mesh = gmshsphere(1.0, 0.5; order=2) + @show typeof(mesh) + @test mesh isa CSM.AbstractMesh{3,3} + fine = CSM.bisecting_refinement(mesh; order=3) + + @test length(skeleton(mesh,0)) - length(skeleton(mesh,1)) + length(mesh) == 2 + @test length(skeleton(fine,0)) - length(skeleton(fine,1)) + length(fine) == 2 +end + + +@testitem "curvy mesh bisecting refinement injections [rectangle, p=2]" begin + const CSM = CompScienceMeshes + + verts = [ + point(1.0, 0.0, 0.0), + point(0.0, 1.0, 0.0), + point(0.0, 0.0, 0.0), + point(0.0, 0.5000000000013304, 0.0), + point(0.4999999999986718, 0.0, 0.0), + point(0.5000000000013293, 0.4999999999986707, 0.0), + + point(1.0, 1.0, 0.0), + point(1.0, 0.5, 0.0), + point(0.5, 1.0, 0.0), + + ] + + I = CSM.HOSimplexGraph{3,2} + faces = [ + I(1, 2, 3, 6, 4, 5), + I(1, 7, 2, 8, 9, 6), + ] + + mesh = CSM.CurvilinearMesh(verts, faces, 2, 2) + fine = CSM.bisecting_refinement(mesh; order=3) + parent = CSM.parent(fine) + + @test length(skeleton(mesh,0)) - length(skeleton(mesh,1)) + length(mesh) == 1 + @test length(skeleton(fine,0)) - length(skeleton(fine,1)) + length(fine) == 1 + + @show length(skeleton(fine,1)) == 33 + + # @test dimension(fine) == 2 + # @test length(fine) == 18 + + # c1 = cells(parent)[1] + for t in parent + chart1 = CSM.chart(parent, t) + for (injct, child) in zip(CSM.injections(parent, t), CSM.children(parent, t)) + + # injct1 = CSM.injections(parent, 1)[1] + # child1 = chart(fine, CSM.children(parent, 1)[1]) + + uv = cartesian(CSM.center(injct)) + @show uv + p1 = neighborhood(chart1, uv) + # p2 = CSM.center(child) + p3 = CSM.center(chart(fine, child)) + + x1 = cartesian(p1) + # @show x2 = cartesian(p2) + x3 = cartesian(p3) + + # @test isapprox(x1, x2, atol=1e-12) + @test isapprox(x1, x3, atol=1e-12) + end + end +end \ No newline at end of file diff --git a/test/primitives/gmsh.jl b/test/primitives/gmsh.jl new file mode 100644 index 0000000..6368777 --- /dev/null +++ b/test/primitives/gmsh.jl @@ -0,0 +1,18 @@ +@testitem "meshgeo [p=1,dim=2]" begin + + fn = joinpath(pkgdir(CompScienceMeshes), "test", "assets", "cuboid.geo") + + m = meshgeo(fn; element=:triangle) + @test m isa CompScienceMeshes.Mesh + @test dimension(m) == 2 +end + + +@testitem "meshgeo [p=2,dim=2]" begin + + fn = joinpath(pkgdir(CompScienceMeshes), "test", "assets", "cuboid.geo") + + m = meshgeo(fn; element=:triangle, order=2) + @test m isa CompScienceMeshes.CurvilinearMesh{3,3,6,Float64,2} + @test dimension(m) == 2 +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index a4411b7..65a7c62 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -56,6 +56,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..aa6c5ae --- /dev/null +++ b/test/test_curvilinear_elements.jl @@ -0,0 +1,344 @@ +using Test +using CompScienceMeshes +using StaticArrays + +const CSM = CompScienceMeshes + +## Line order 1 +v1 = SVector(1.0, 0.0) +v2 = SVector(0.0, 0.0) + +verts1 = [v1, v2] + +faces1 = [CSM.SimplexGraph(1,2)] +# faces1c = [SVector(1,2)] +faces1c = [CSM.HOSimplexGraph{2,1}(1,2)] + +m1 = CompScienceMeshes.Mesh(verts1, faces1) +m2 = CompScienceMeshes.CurvilinearMesh(verts1, faces1c, 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 = [CSM.SimplexGraph(1,2,3)] +# faces2c = [SVector(1,2,3)] +faces2c = [CSM.HOSimplexGraph{2,2}(1,2,3)] + +m2 = CompScienceMeshes.CurvilinearMesh(verts2, faces2c, 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] +faces_parabola = [CSM.HOSimplexGraph{2,2}(1,2,3)] + +m_parabola = CompScienceMeshes.CurvilinearMesh(verts_parabola, faces_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 = [CSM.SimplexGraph(1, 2, 3)] +faces1c = [CSM.HOSimplexGraph{3,1}(1, 2, 3)] + +m1 = CompScienceMeshes.Mesh(verts1, faces1) +m2 = CompScienceMeshes.CurvilinearMesh(verts1, faces1c, 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 = [CSM.SimplexGraph(1,2)] +# faces1c = [SVector(1,2)] + +verts2 = [v1, v2, v3, v4, v5] +faces2 = [CSM.HOSimplexGraph{2,4}(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 = [CSM.SimplexGraph(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)] +faces2 = [CSM.HOSimplexGraph{3,2}(1, 2, 3, 6, 4, 5)] + +m1 = CompScienceMeshes.Mesh(verts1, faces1) +m2 = CompScienceMeshes.CurvilinearMesh(verts2, faces2, 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 = [CSM.SimplexGraph(1, 2, 3)] + +verts2 = [ + SVector(1.0, 0.0, 0.0), # Nodes + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 0.0), + SVector(0.0, 0.8571428571429518, 0.0), # Edge based vertices [edge 1] + 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), # Edge based vertices [edge 2] + 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), # Edge based vertices [edge 3] + 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), # Internal vertices + 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 = [CSM.HOSimplexGraph{3,7}(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 = [CSM.HOSimplexGraph{2,4}(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