Support for curvilinear (simplicial) meshes.

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

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")

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,72 +309,179 @@ 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.
 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
 
 """