What is the best way to deal with multiple `mesh_type`s?

[jlchan]

We will be able to support cut-cell, hybrid, and non-conforming meshes using mesh_type in MeshData, but how should we handle multiple mesh_types? If we want a hybrid non-conforming mesh, how should this be represented in a single mesh_type field? A combined type?

I think we will either need to

• factor out parts of the MeshData initialization which depend on each mesh_type in the "combined" mesh type, or
• just specialize on the "combined" mesh type.

I think the latter may just be easier - reusing code sounds like it would require a fair bit of refactoring (and depend closely on the mesh types in consideration.)

[jmark]

Hi @jlchan ! I am playing with you hybrid mesh implementation in StartupDG.jl package. I have a question regarding following code:

using StartUpDG

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

#   1  3_______4______5
#      |   4   |  6  / 
#      |1     2|7  5
#      |   3   |  /  
#  -1  1 ----- 2        1

VX = [-1; 0; -1; 0; 1]
VY = [-1; -1; 1; 1; 1]
EToV = [[1 2 3 4], [2 4 5]]
md = MeshData(VX, VY, EToV, rds)

Now I want the geometric terms. According to your documentation md.rstxyzJ should be a 2x2 matrix. But not in this case. I do not know how to interpretmd.rstxyzJ.

Also, when I try to access md.rxJ I get the following error:

MethodError: no method matching getindex(::Tuple{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(Tri = ViewAxis(1:10, ShapedAxis((10, 1), NamedTuple())), Quad = ViewAxis(11:26, ShapedAxis((16, 1), NamedTuple())))}}}, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(Tri = ViewAxis(1:10, ShapedAxis((10, 1), NamedTuple())), Quad = ViewAxis(11:26, ShapedAxis((16, 1), NamedTuple())))}}}}, ::Int64, ::Int64)
Closest candidates are:
  getindex(::Tuple, ::Int64) at tuple.jl:29
  getindex(::Tuple, ::Integer) at tuple.jl:30
  getindex(::Tuple, ::Real) at deprecated.jl:70
  ...

Stacktrace:
 [1] getproperty(x::MeshData{2, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(Tri = ViewAxis(1:10, ShapedAxis((10, 1), NamedTuple())), Quad = ViewAxis(11:26, ShapedAxis((16, 1), NamedTuple())))}}}, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(Tri = ViewAxis(1:12, ShapedAxis((12, 1), NamedTuple())), Quad = ViewAxis(13:28, ShapedAxis((16, 1), NamedTuple())))}}}, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(Tri = ViewAxis(1:12, ShapedAxis((12, 1), NamedTuple())), Quad = ViewAxis(13:28, ShapedAxis((16, 1), NamedTuple())))}}}, Vector{Int64}, Vector{Matrix{Int64}}, Vector{Int64}, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(Tri = ViewAxis(1:12, ShapedAxis((12, 1), NamedTuple())), Quad = ViewAxis(13:28, ShapedAxis((16, 1), NamedTuple())))}}}, Tuple{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(Tri = ViewAxis(1:10, ShapedAxis((10, 1), NamedTuple())), Quad = ViewAxis(11:26, ShapedAxis((16, 1), NamedTuple())))}}}, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(Tri = ViewAxis(1:10, ShapedAxis((10, 1), NamedTuple())), Quad = ViewAxis(11:26, ShapedAxis((16, 1), NamedTuple())))}}}}, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(Tri = ViewAxis(1:10, ShapedAxis((10, 1), NamedTuple())), Quad = ViewAxis(11:26, ShapedAxis((16, 1), NamedTuple())))}}}, Vector{Int64}, Vector{Int64}}, s::Symbol)
   @ StartUpDG ~/.julia/packages/StartUpDG/kgt4Q/src/MeshData.jl:123
 [2] top-level scope
   @ In[27]:1
 [3] eval
   @ ./boot.jl:368 [inlined]
 [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1428

[jlchan]

Hi @jmark - can you try updating StartUpDG? This bug should have been fixed in e0ba25bbcb0550d9e51978a00767144ebba94067 and in v0.14.

[jmark]

I had v0.13.9 installed. It simply did not come to my mind to try an update first. Sorry for that. :') Now it works as expected! Thank you!

[jmark]

@jlchan The next thing I would like to have is to apply mappings on hybrid meshes. The following code does not work of course since there is no implementation yet for hybrid meshes, or is there?

# Apply mapping to `meshdata`.
@unpack x, y = meshdata

x_ = similar(x)
y_ = similar(y)

for i = 1:length(x)
    x_[i],y_[i] = mapping(x[i],y[i])
end

meshdata = MeshData(rds, meshdata, x_, y_)

I tried to implement this but got a little bit stuck since I need to repeatedly transpose between matrices of ComponentVectors and ComponentVectors of matrices.

function MeshData(rds::LittleDict{AbstractElemShape, <:RefElemData}, md::MeshData{Dim}, xyz...) where {Dim}

    tn(x) = typeof(x).name.name

    xyzf = ntuple(i -> ComponentArray(NamedTuple(Pair.(tn.(keys(rds)), map(k -> rds[k].Vf * xyz[i][tn(k)], rds.keys) ))), Dim)
    xyzq = ntuple(i -> ComponentArray(NamedTuple(Pair.(tn.(keys(rds)), map(k -> rds[k].Vq * xyz[i][tn(k)], rds.keys) ))), Dim)

    # Compute geometric factors and surface normals    
    geo = ComponentArray(NamedTuple(Pair.(tn.(keys(rds)), map(k -> geometric_factors(xyz[1][tn(k)],xyz[2][tn(k)], rds[k].Drst...) , rds.keys) )))

    # 4 entries in the geometric term matrix for 2D hybrid meshes
    rstxyzJ = ntuple(i -> ComponentArray(NamedTuple(Pair.(tn.(keys(rds)), map(k -> geo[tn(k)][i], rds.keys) ))), Dim*Dim)
    rstxyzJ = SMatrix{Dim,Dim}(rstxyzJ)

    # How to cleverly elevate this function call to ComponentVectors?        
    geof = compute_normals(rstxyzJ, rd.Vf, rd.nrstJ...)

#     J = last(geo)
#     wJq = diagm(rd.wq) * (rd.Vq * J)

#     # TODO: should we warp VXYZ as well? Or just set it to nothing since it no longer determines geometric terms?
#     setproperties(md, (xyz, xyzq, xyzf, rstxyzJ, J, wJq,
#                        nxyzJ=geof[1:Dim], Jf=last(geof)))
end

Is there an elegant way to do this? Or do we have to reimplement the helper routines for ComponentVectors or do we have to make these

[jlchan]

I think I should be able to add that in this morning.

I might try to rewrite the original curving routine so we can reuse it for hybrid meshes too.

[jlchan]

Added in e5f145439e5818ac5cdc1d0624e71a907baf45cb. Can you check that this works for you?

[jmark]

Works like a charm! Thank you! :)

[jlchan]

Thanks! Would you mind adding a test for this functionality then?

[jmark]

Sure. I'll try.

[jlchan]

Let me know if I can help with it. I think since it's an isoperarametric curved mesh, we should be able to differentiate linear exactly...

[jlchan]

Added tests in 6bde84faa11ce259669ccb7f2d79687562a47ec7. @jmark OK if I close this?

[jmark]

Wow, you are fast! If you wish to close this issue, then go ahead! :)