Move VXYZ and EToV into a VertexMappedMesh mesh type (#70)

* adding VertexMappedMesh to MeshData

* moving VXYZ, EToV into hybrid mesh types

* moving VXYZ, EToV into cut cell mesh type

* moving VXYZ, EToV into noncon mesh type

* remove outdated comment

* remove macro for _getproperty, use function with inlining

* update docs

* update news

* add deprecation

* cleanup comments

* comments

* misc formatting

NEWS.md

@@ -10,15 +10,21 @@ StartUpDG.jl follows the interpretation of [semantic versioning (semver)](https:
 
 #### Changed
 
-* Changes related to element types and `MeshData`:
+* Changes related to element types:
   * upstream change in NodesAndModes.jl: the abstract type `AbstractElemShape` is now parametrized by the dimension, e.g., `Line <: AbstractElemShape{1}`, `Tri <: AbstractElemShape{2}`. 
   * the spatial dimension `Dim` parameter in `MeshData` is now inferred from the element type through `elem <: AbstractElemShape{Dim}`
   * the `PhysicalFrame` type now has leading type parameter `PhysicalFrame{NDIMS} <: AbstractElemShape{NDIMS}`
-* the `Nplot` field has been removed from `RefElemData`
 * `PhysicalFrame`'s fields are now restricted to be type `SVector` only (instead of `Union{SVector, Tuple}`)
+* The `MeshData` fields `VXYZ` and `EToV` have been moved into a `VertexMappedMesh` type. However, they can still be accessed as a property of `MeshData` (e.g., `md.VX` will still work). 
+
+#### Removed 
+
+* the `Nplot` field has been removed from `RefElemData`
 * all usages of `ComponentArrays` have been replaced by `NamedArrayPartition`
+* the deprecated `MeshData(md::MeshData, rd::RefElemData, xyz...)` constructor has been removed
 
 #### Deprecated
 
+* the old constructor for `MeshData` with fields `VXYZ` and `EToV` has been deprecated and will be removed in the next breaking release. 
 * the old constructor for `RefElemData` with field `Nplot` has been deprecated and will be removed in the next breaking release. 
 

docs/src/MeshData.md

@@ -9,7 +9,8 @@
 * `J, Jf`: volume and surface Jacobians evaluated at interpolation points and surface quadrature points, respectively. `J` is a matrix of size ``N_p \times num_elements``, while `Jf` is a matrix of size ``N_f \times num_elements``. 
 * `nxyz::NTuple{Dim, ...}` and `nxyzJ::NTuple{Dim, ...}`: normalized and `Jf` scaled outward normals evaluated at surface quadrature points. Each element of `nxyzJ` is a matrix of size ``N_f\times num_elements``. 
 
-These are the main quantities used to construct a DG solver. 
+These are the main quantities used to construct a DG solver. Information specific to the type of mesh used is
+stored in the `md.mesh_type` field. 
 
 # Setting up `md::MeshData`
 
@@ -51,7 +52,7 @@ julia> md_periodic_x.is_periodic
 
 ## Creating curved meshes
 
-It is common to generate curved meshes by first generating a linear mesh, then moving high order nodes on the linear mesh. This can be done by calling [`MeshData`](@ref) again with new `x,y` coordinates:
+It is common to generate curved meshes by first generating a linear mesh, then moving high order nodes on the linear mesh. This can be done by calling [`MeshData`](@ref) again with new `x, y` coordinates:
 ```julia
 md = MeshData((VX, VY), EToV, rd)
 @unpack x, y = md
@@ -62,11 +63,9 @@ md_curved = MeshData(rd, md, x, y)
 
 More generally, one can create a copy of a `MeshData` with certain fields modified by using `@set` or `setproperties` from `Setfield.jl`.
 
-## Unstructured (and pre-defined) triangular meshes using Triangulate
+## Unstructured and pre-defined triangular meshes using Triangulate
 
-If `Triangulate` is also loaded, then StartUpDG will include additional utilities for creating and visualizing meshes. 
-
-Several pre-defined geometries are included in StartUpDG.jl. A few examples are `SquareDomain`, `RectangularDomainWithHole`, `Scramjet`, and `CircularDomain`. See `triangulate_example_meshes.jl` for a more complete list and field arguments. These can each be called using `triangulate_domain`, for example the following code will create a mesh of a scramjet:
+StartUpDG.jl also includes additional utilities based on Triangulate.jl for creating and visualizing meshes. Several pre-defined geometries are included in StartUpDG.jl. A few examples are `SquareDomain`, `RectangularDomainWithHole`, `Scramjet`, and `CircularDomain`. See `triangulate_example_meshes.jl` for a more complete list and field arguments. These can each be called using `triangulate_domain`, for example the following code will create a mesh of a scramjet:
 ```julia
 meshIO = triangulate_domain(Scramjet())
 (VX, VY), EToV = triangulateIO_to_VXYEToV(meshIO)

src/MeshData.jl

@@ -13,27 +13,22 @@ md = MeshElemData(VXY, EToV, rd)
 @unpack x, y = md
 ```
 """
-Base.@kwdef struct MeshData{Dim, MeshType, VolumeType, FaceType, VolumeQType,
-                            VertexType, EToVType, FToFType, 
-                            VolumeWeightType, VolumeGeofacsType, VolumeJType,
+Base.@kwdef struct MeshData{Dim, MeshType, 
+                            VolumeType, FaceType, VolumeQType, FToFType, 
+                            VolumeWeightType, VolumeGeofacsType, VolumeJType, 
                             ConnectivityTypeM, ConnectivityTypeP, BoundaryMapType}
 
-    # this field defaults to the element shape, but can be 
-    # used to specify cut-cell, hybrid, non-conforming meshes. 
+    # this field defaults to the element shape, but can be used to specify 
+    # different mesh types (e.g., vertex-mapped, cut-cell, hybrid, non-conforming meshes. 
     mesh_type::MeshType 
 
-    # TODO: move VXYZ, EToV into a `VertexMappedMesh` mesh_type?
-    VXYZ::NTuple{Dim, VertexType}   # vertex coordinates
-    EToV::EToVType                  # mesh vertex array     
-
     FToF::FToFType                  # face connectivity
 
     xyz::NTuple{Dim, VolumeType}    # physical points
     xyzf::NTuple{Dim, FaceType}     # face nodes
     xyzq::NTuple{Dim, VolumeQType}  # phys quad points, Jacobian-scaled weights
     wJq::VolumeWeightType
 
-    # TODO: move mapP, mapB into a "conforming mesh type"?
     # arrays of connectivity indices between face nodes
     mapM::ConnectivityTypeM
     mapP::ConnectivityTypeP
@@ -53,24 +48,24 @@ Base.@kwdef struct MeshData{Dim, MeshType, VolumeType, FaceType, VolumeQType,
     is_periodic::NTuple{Dim, Bool}
 end
 
-# MeshData constructor where we do not specify `nxyz`
-function MeshData(mesh_type, VXYZ, EToV, FToF, xyz, xyzf, xyzq, wJq, 
-                  mapM, mapP, mapB, rstxyzJ, J, nxyzJ, Jf, is_periodic) 
+# MeshData constructor where we do not specify `nxyz` and instead compute it from `nrstJ` and `Jf`
+function MeshData(mesh_type, FToF, xyz, xyzf, xyzq, wJq, mapM, mapP, mapB, rstxyzJ, J, nxyzJ, Jf, is_periodic) 
 
     nxyz = map(nJ -> nJ ./ Jf, nxyzJ)        
 
-    return MeshData(mesh_type, VXYZ, EToV, FToF, 
-                    xyz, xyzf, xyzq, wJq, mapM, mapP, mapB,
-                    rstxyzJ, J, nxyz, nxyzJ, Jf, is_periodic)
+    return MeshData(mesh_type, FToF, xyz, xyzf, xyzq, wJq, mapM, mapP, mapB, rstxyzJ, J, nxyz, nxyzJ, Jf, is_periodic)
 end
 
+# TODO: remove in next breaking release after v0.15. `MeshData` constructor where `VXYZ`, `EToV` are specified 
+@deprecate MeshData(element_type, VXYZ, EToV, FToF, xyz, xyzf, xyzq, wJq, mapM, mapP, mapB, rstxyzJ, J, nxyzJ, Jf, periodicity) MeshData(VertexMappedMesh(element_type, VXYZ, EToV), FToF, xyz, xyzf, xyzq, wJq, mapM, mapP, mapB, rstxyzJ, J, nxyzJ, Jf, periodicity)
+
 function ConstructionBase.setproperties(md::MeshData, patch::NamedTuple)
     fields = (haskey(patch, symbol) ? getproperty(patch, symbol) : getproperty(md, symbol) for symbol in fieldnames(typeof(md)))
     return MeshData(fields...)
 end
 
 ConstructionBase.getproperties(md::MeshData) = 
-    (; mesh_type=md.mesh_type, VXYZ=md.VXYZ, EToV=md.EToV, FToF=md.FToF, xyz=md.xyz, xyzf=md.xyzf, xyzq=md.xyzq, wJq=md.wJq,
+    (; mesh_type=md.mesh_type, FToF=md.FToF, xyz=md.xyz, xyzf=md.xyzf, xyzq=md.xyzq, wJq=md.wJq,
        mapM=md.mapM, mapP=md.mapP, mapB=md.mapB, rstxyzJ=md.rstxyzJ, J=md.J, nxyz = md.nxyz, nxyzJ=md.nxyzJ, Jf=md.Jf,
        is_periodic=md.is_periodic)
 
@@ -81,7 +76,7 @@ end
 
 function Base.show(io::IO, ::MIME"text/plain", md::MeshData{DIM, MeshType}) where {DIM, MeshType}
     @nospecialize md
-    print(io,"$MeshType MeshData of dimension $DIM with $(md.num_elements) elements")
+    print(io,"MeshData of dimension $DIM with $(md.num_elements) elements")
 end
 
 function Base.propertynames(x::MeshData{1}, private::Bool = false)
@@ -99,17 +94,9 @@ function Base.propertynames(x::MeshData{3}, private::Bool = false)
             :nx, :ny, :nz, :rxJ, :sxJ, :txJ, :ryJ, :syJ, :tyJ, :rzJ, :szJ, :tzJ)
 end
 
-# convenience routines for unpacking individual tuple entries
-function Base.getproperty(x::MeshData, s::Symbol)
-
-    if s==:VX
-        return getfield(x, :VXYZ)[1]
-    elseif s==:VY
-        return getfield(x, :VXYZ)[2]
-    elseif s==:VZ
-        return getfield(x, :VXYZ)[3]
-
-    elseif s==:x
+# seems like we need to use @inline to ensure type stability
+@inline function meshdata_getproperty(x::MeshData, s::Symbol)
+    if s==:x
         return getfield(x, :xyz)[1]
     elseif s==:y
         return getfield(x, :xyz)[2]
@@ -177,7 +164,44 @@ function Base.getproperty(x::MeshData, s::Symbol)
     end
 end
 
-num_elements(md) = size(getfield(md, :EToV), 1)
+# generic fallback
+Base.getproperty(x::MeshData, s::Symbol) = meshdata_getproperty(x, s)
+
+"""
+    struct VertexMappedMesh
+
+The default `MeshData` mesh type, represents a mesh which is defined purely by 
+vertex locations and element-to-vertex connectivities. For example, these include 
+affine triangular meshes or bilinear quadrilateral or trilinear hexahedral meshes.
+
+# Fields
+element_type :: TE <: AbstractElemShape \\
+VXYZ :: TV \\
+EToV :: TEV
+"""
+struct VertexMappedMesh{TE <: AbstractElemShape, TV, TEV}
+    element_type::TE
+    VXYZ::TV
+    EToV::TEV
+end
+
+# convenience accessor routines for `VertexMappedMesh` types (lets you do `md.VX` instead of `md.mesh_type.VX`)
+function Base.getproperty(x::MeshData{Dim, <:VertexMappedMesh}, s::Symbol) where {Dim}
+
+    if s===:VX
+        return getfield(getfield(x, :mesh_type), :VXYZ)[1]
+    elseif s===:VY
+        return getfield(getfield(x, :mesh_type), :VXYZ)[2]
+    elseif s===:VZ
+        return getfield(getfield(x, :mesh_type), :VXYZ)[3]
+    elseif s===:EToV
+        return getfield(getfield(x, :mesh_type), s)
+    else
+        meshdata_getproperty(x, s)
+    end
+end
+
+num_elements(md::MeshData{Dim, <:VertexMappedMesh}) where {Dim} = size(md.mesh_type.EToV, 1)
 
 """
     MeshData(VXYZ, EToV, rd::RefElemData)
@@ -236,7 +260,7 @@ function MeshData(VX::AbstractVector{Tv}, EToV, rd::RefElemData{1}) where {Tv}
 
     is_periodic = (false,)
 
-    return MeshData(rd.element_type, tuple(VX), EToV, FToF,
+    return MeshData(VertexMappedMesh(rd.element_type, tuple(VX), EToV), FToF,
                     tuple(x), tuple(xf), tuple(xq), wJq,
                     collect(mapM), mapP, mapB,
                     SMatrix{1,1}(tuple(rxJ)), J,
@@ -276,7 +300,7 @@ function MeshData(VX, VY, EToV, rd::RefElemData{2})
     nxJ, nyJ, Jf = compute_normals(rstxyzJ, rd.Vf, rd.nrstJ...)
 
     is_periodic = (false, false)
-    return MeshData(rd.element_type, tuple(VX, VY), EToV, FToF,
+    return MeshData(VertexMappedMesh(rd.element_type, tuple(VX, VY), EToV), FToF,
                     tuple(x, y), tuple(xf, yf), tuple(xq, yq), wJq,
                     mapM, mapP, mapB,
                     SMatrix{2, 2}(tuple(rxJ, ryJ, sxJ, syJ)), J,
@@ -314,13 +338,14 @@ function MeshData(VX, VY, VZ, EToV, rd::RefElemData{3})
     nxJ, nyJ, nzJ, Jf = compute_normals(rstxyzJ, rd.Vf, rd.nrstJ...)
 
     is_periodic = (false, false, false)
-    return MeshData(rd.element_type, tuple(VX, VY, VZ), EToV, FToF,
+    return MeshData(VertexMappedMesh(rd.element_type, tuple(VX, VY, VZ), EToV), FToF,
                     tuple(x, y, z), tuple(xf, yf, zf), tuple(xq, yq, zq), wJq,
                     mapM, mapP, mapB,
                     rstxyzJ, J, tuple(nxJ, nyJ, nzJ), Jf,
                     is_periodic)
 end
 
+# TODO: remove with breaking change v0.15
 @deprecate MeshData(md::MeshData, rd::RefElemData, xyz...) MeshData(rd, md, xyz...)
 
 function recompute_geometry(rd::RefElemData{Dim}, xyz) where {Dim}

src/RefElemData_SBP.jl

@@ -75,24 +75,24 @@ function RefElemData(elementType::Tri, approxType::SBP, N; tol = 100*eps())
     (Qrh, Qsh), _ = hybridized_SBP_operators(rd)
     Nq = length(rd.wq)
     Vh_sbp = [I(Nq); Ef]
-    Qr = Vh_sbp'*Qrh*Vh_sbp
-    Qs = Vh_sbp'*Qsh*Vh_sbp
-    Dr,Ds = (x->diagm(1 ./ rd.wq) * x).((Qr,Qs))
+    Qr = Vh_sbp' * Qrh * Vh_sbp
+    Qs = Vh_sbp' * Qsh * Vh_sbp
+    Dr,Ds = (x -> diagm(1 ./ rd.wq) * x).((Qr, Qs))
 
     rd = @set rd.rst = quad_rule_vol[1:2]   # set nodes = SBP nodes
     rd = @set rd.rstq = quad_rule_vol[1:2]  # set quad nodes = SBP nodes
-    rd = @set rd.Drst = (Dr,Ds)
+    rd = @set rd.Drst = (Dr, Ds)
     rd = @set rd.Fmask = vec(Fmask)
 
-    # TODO: make these more efficient with custom operator?
-    rd = @set rd.Vf = droptol!(sparse(Ef),tol)
-    rd = @set rd.LIFT = Diagonal(rd.wq)\(rd.Vf'*Diagonal(rd.wf)) 
+    # TODO: make these more efficient with custom operators?
+    rd = @set rd.Vf = droptol!(sparse(Ef), tol)
+    rd = @set rd.LIFT = Diagonal(rd.wq) \ (rd.Vf' * Diagonal(rd.wf)) 
 
     # make V1 the interpolation matrix from triangle vertices to SBP nodal points
-    rd = @set rd.V1 = vandermonde(elementType,N,rd.rst...)/rd.VDM * rd.V1
+    rd = @set rd.V1 = vandermonde(elementType, N, rd.rst...) / rd.VDM * rd.V1
 
     # Vp operator = projects SBP nodal vector onto degree N polynomial, then interpolates to plotting points
-    rd = @set rd.Vp = vandermonde(elementType,N,rd.rstp...)/rd.VDM * rd.Pq
+    rd = @set rd.Vp = vandermonde(elementType, N, rd.rstp...) / rd.VDM * rd.Pq
 
     return _convert_RefElemData_fields_to_SBP(rd, approxType)
 end

src/connectivity_functions.jl

@@ -387,14 +387,15 @@ function build_periodic_boundary_maps!(xf, yf, zf,
     return mapPB[:]
 end
 
-function compute_boundary_face_centroids(md::MeshData{NDIMS, <:Union{Wedge, Pyr}}) where {NDIMS}
-    element_type = md.mesh_type
+function compute_boundary_face_centroids(md::MeshData{3, <:VertexMappedMesh{<:Union{<:Wedge, <:Pyr}}})
+    (; element_type) = md.mesh_type
     (; node_ids_by_face) = element_type
 
     face_nodes_per_element = maximum(maximum.(node_ids_by_face))
     boundary_faces = findall(eachindex(md.FToF) .== vec(md.FToF))
 
     # compute face centroids
+    NDIMS = 3
     face_centroids = ntuple(_ -> similar(md.xf, (length(boundary_faces), )), NDIMS)
     sk = 1    
     for e in 1:md.num_elements
@@ -414,7 +415,8 @@ function compute_boundary_face_centroids(md::MeshData{NDIMS, <:Union{Wedge, Pyr}
 end
 
 # specialize on 3D elements with different types of faces
-function make_periodic(md::MeshData{3, <:Union{Wedge, Pyr}}, is_periodic::NTuple{3}; tol=1e-12, kwargs...)
+function make_periodic(md::MeshData{3, <:VertexMappedMesh{<:Union{<:Wedge, <:Pyr}}}, 
+                       is_periodic::NTuple{3}; tol=1e-12, kwargs...)
 
     NDIMS = 3
     (; mapM) = md

src/cut_cell_meshes.jl

@@ -611,9 +611,6 @@ function MeshData(rd::RefElemData, curves, cells_per_dimension_x, cells_per_dime
         view(wJq.cut, :, e) .= wq
     end
 
-    VXYZ = ntuple(_ -> nothing, 2)
-    EToV = nothing # dummy field for cut cells
-
     # default to non-periodic 
     is_periodic = (false, false)   
 
@@ -642,7 +639,7 @@ function MeshData(rd::RefElemData, curves, cells_per_dimension_x, cells_per_dime
     cut_face_node_ids = [face_ids(e) for e in 1:num_cut_cells]
 
     return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_ids, curves, cut_cell_data), 
-                    VXYZ, EToV, FToF, (x, y), (xf, yf), (xq, yq), wJq, 
+                    FToF, (x, y), (xf, yf), (xq, yq), wJq, 
                     mapM, mapP, mapB, rstxyzJ, J, (nxJ, nyJ), Jf, is_periodic)
 
 end

src/hybrid_meshes.jl

@@ -1,11 +1,13 @@
 # mesh_type identifier for hybrid meshes
-struct HybridMesh{T}
+struct HybridMesh{T, TV, TE}
     element_types::T
+    VXYZ::TV
+    EToV::TE
     # This is an inner constructor for HybridMesh. It sorts tuples of element types alphabetically 
     # so that `HybridMesh.element_types` always yields the same ordering of element types.
-    function HybridMesh(element_types::T) where {T <: Tuple}
+    function HybridMesh(element_types::T, VXYZ, EToV) where {T <: Tuple}
         p = sortperm(collect(typeof.(element_types)); by=string)
-        return new{typeof(Tuple(element_types[p]))}(Tuple(element_types[p]))
+        return new{typeof(Tuple(element_types[p])), typeof(VXYZ), typeof(EToV)}(Tuple(element_types[p]), VXYZ, EToV)
     end
 end
 
@@ -170,7 +172,7 @@ function MeshData(VX, VY, EToV_unsorted, rds::LittleDict{<:AbstractElemShape, <:
 
     mapM, mapP, mapB = vec.(build_node_maps(FToF, xyzf))
 
-    return MeshData(HybridMesh(Tuple(element_types)), (VX, VY), EToV, FToF, 
+    return MeshData(HybridMesh(Tuple(element_types), (VX, VY), EToV), FToF, 
                     xyz, xyzf, xyzq, wJq,
                     mapM, mapP, mapB,
                     rstxyzJ, J, nxyzJ, Jf, 

src/nonconforming.jl

@@ -34,7 +34,9 @@ The `mortar_projection_matrix` similarly maps values from 2 mortar faces back to
 original non-conforming face. These can be used to create DG solvers on non-conforming meshes.
 
 """
-struct NonConformingMesh{CF, NCF, I, P}
+struct NonConformingMesh{TV, TE, CF, NCF, I, P}
+    VXYZ::TV
+    EToV::TE
     conforming_faces::CF
     non_conforming_faces::NCF
     mortar_interpolation_matrix::I
@@ -91,6 +93,8 @@ end
 # one non-conforming quad face is split into 2 mortar faces
 num_mortars_per_face(rd::RefElemData{2, Quad}) = 2
 
+num_elements(md::MeshData{Dim, <:NonConformingMesh}) where {Dim} = size(getproperty(md.mesh_type, :EToV), 1)
+
 function MeshData(mesh::NonConformingQuadMeshExample, rd::RefElemData{2, Quad})
 
     (VX, VY) = mesh.VXY
@@ -163,10 +167,10 @@ function MeshData(mesh::NonConformingQuadMeshExample, rd::RefElemData{2, Quad})
 
     is_periodic = (false, false)
 
-    mesh_type = NonConformingMesh(conforming_faces, non_conforming_faces, 
+    mesh_type = NonConformingMesh(tuple(VX, VY), EToV, conforming_faces, non_conforming_faces, 
                                   mortar_interpolation_matrix, mortar_projection_matrix)
 
-    return MeshData(mesh_type, tuple(VX, VY), EToV, FToF,
+    return MeshData(mesh_type, FToF,
                     tuple(x, y), tuple(x_mortar, y_mortar), tuple(xq, yq), wJq,
                     mapM, mapP, mapB,
                     SMatrix{2, 2}(tuple(rxJ, ryJ, sxJ, syJ)), J,