Additional mesh types

In addition to more "standard" mesh types, StartUpDG.jl also has experimental support for hybrid and cut-cell meshes. Both are currently restricted to two dimensional domains.

Hybrid meshes

There is initial support for hybrid meshes in StartUpDG.jl. The following is a short example where we interpolate a polynomial and compute its derivative.

rds = RefElemData((Tri(), Quad()), N = 3)

# Simple hybrid mesh for testing
#   1  7______8______9
#      |      | 3  / |
#      |   4  |  / 5 |
#   0  4 ---- 5 ---- 6 
#      |      |      |
#      |   1  |   2  |
#   -1 1 ---- 2 ---- 3
#     -1      0      1
VX = [-1; 0; 1; -1; 0; 1; -1; 0; 1]
VY = [-1; -1; -1; 0; 0; 0; 1; 1; 1]
EToV = [[1 2 4 5], [2 3 5 6], [5 8 9], [4 5 7 8], [9 6 5]]

md = MeshData(VX, VY, EToV, rds)

# test that the local derivatives of a polynomial recover the exact derivative
(; x, y ) = md
u = @. x^3 - x^2 * y + 2 * y^3
dudx = @. 3 * x^2 - 2 * x * y

# compute local derivatives
(; rxJ, sxJ, J ) = md
dudr, duds = similar(md.x), similar(md.x)
dudr.Quad .= rds[Quad()].Dr * u.Quad
duds.Quad .= rds[Quad()].Ds * u.Quad
dudr.Tri .= rds[Tri()].Dr * u.Tri
duds.Tri .= rds[Tri()].Ds * u.Tri

@show norm(@. dudx - (rxJ * dudr + sxJ * duds) / J) # should be O(1e-14)

The main difference in the representation of hybrid meshes compared with standard MeshData objects is the use of NamedArrayPartition arrays as storage for the geometric coordinates. These arrays have "fields" corresponding to the element type, for example

md.x.Tri
md.x.Quad

but can still be indexed as linear arrays.

The mapP field behaves similarly. If we interpolate the values of u for each element type to surface quadrature nodes, we can use mapP to linearly index into the array to find neighbors.

uf = similar(md.xf)
uf.Quad .= rds.Quad.Vf * u.Quad
uf.Tri .= rds.Tri.Vf * u.Tri
uf[md.mapP] # this returns the exterior node values

Cut Meshes

Initial support for cut-cell meshes is available via PathIntersections.jl. By passing in a tuple of curves (defined as parametrized functions of one coordinate, see PathIntersections.jl documentation for more detail), StartUpDG.jl can compute a MeshData for a cut-cell mesh.

circle = PresetGeometries.Circle(R=0.33, x0=0, y0=0)
cells_per_dimension_x, cells_per_dimension_y = 4, 4

rd = RefElemData(Quad(), N=3)
md = MeshData(rd, (circle, ), cells_per_dimension_x, cells_per_dimension_y; precompute_operators=true)

The interpolation points on cut cells md.x.cut are determined from sampled points and a pivoted QR decomposition. The quadrature points on cut cells md.xq.cut are determined similarly. However, the cut-cell quadrature weights md.wJq.cut are not currently positive. The optional keyword argument precompute_operators specifies whether to precompute differentiation, face interpolation, mass, and lifting matrices for each cut cell. If precompute_operators=true, these are stored in md.mesh_type.cut_cell_operators.

As with hybrid meshes, the nodal coordinates md.x, md.y are NamedArrayPartitions with cartesian and cut fields. For example, md.x.cartesian and md.x.cut are the x-coordinates of the Cartesian and cut cells, respectively. Likewise, md.mapP indexes linearly into the array of face coordinates and specifies exterior node indices. For example, we can interpolate a function to face nodes and compute exterior values via the following code:

(; x, y) = md
u = @. x^rd.N - x * y^(rd.N-1) - x^(rd.N-1) * y + y^rd.N # some random function 

# interpolate the solution to face nodes
uf = similar(md.xf) 
uf.cartesian .= rd.Vf * u.cartesian
for e in 1:size(md.x.cut, 2)
    ids = md.mesh_type.cut_face_nodes[e]
    Vf = md.mesh_type.cut_cell_operators.face_interpolation_matrices[e]
    uf.cut[ids] .= Vf * u.cut[:, e]
end

uf[md.mapP] # these are "exterior" values for each entry of `uf`