Improve efficiency of cut-cell MeshData

[jlchan]

Chasing allocations

[codecov]

Codecov Report

Attention: Patch coverage is 94.11765% with 1 lines in your changes are missing coverage. Please review.

Project coverage is 96.34%. Comparing base (2c56f42) to head (3e3ff7d).
Report is 21 commits behind head on dev.

❗ Current head 3e3ff7d differs from pull request most recent head 2a857d9. Consider uploading reports for the commit 2a857d9 to get more accurate results

Files	Patch %	Lines
src/cut_cell_meshes.jl	90.90%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##              dev     #167      +/-   ##
==========================================
- Coverage   97.12%   96.34%   -0.78%     
==========================================
  Files          25       26       +1     
  Lines        3167     3724     +557     
==========================================
+ Hits         3076     3588     +512     
- Misses         91      136      +45     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[jlchan]

A script to check accuracy of the resulting hybridized SBP operators:

using StartUpDG
using LinearAlgebra

cells_per_dimension = 2
circle = PresetGeometries.Circle(R=0.6, x0=0.0, y0=0.0)

N = 3
rd = RefElemData(Quad(), N) #, quad_rule_face=gauss_quad(0,0,N)) 

@time md = MeshData(rd, (circle, ), cells_per_dimension; 
                    quad_rule_boundary = gauss_quad(0, 0, N+1),
                    precompute_operators=true)

@profview_allocs MeshData(rd, (circle, ), cells_per_dimension; precompute_operators=false) 
@profview MeshData(rd, (circle, ), cells_per_dimension; precompute_operators=false) 

# target_boundary_degree = 2 * N^2 + (N - 1)
# rq_boundary, wq_boundary = gauss_quad(0, 0, ceil(Int, (N^2+(N-1)-1)/2))
# Vtarget = vandermonde(Line(), target_boundary_degree, rq_boundary)          

# e, f = 1, 1
# (; cut_face_node_indices_by_elem_by_face) = mt.cut_cell_data
# ids = cut_face_node_indices_by_elem_by_face[e][f]
# Vtarget_nxy = [Diagonal(md.nxJ.cut[ids]) * Vtarget Diagonal(md.nyJ.cut[ids]) * Vtarget]
# w_pruned, inds = StartUpDG.caratheodory_pruning_qr(Vtarget_nxy, wq_boundary)                    

mt = md.mesh_type

(; volume_interpolation_matrices, 
   face_interpolation_matrices, 
   differentiation_matrices, 
   mass_matrices) = mt.cut_cell_operators
(; wJf) = mt.cut_cell_data
wf = wJf ./ md.Jf
(; x, nxJ, nyJ) = md

(; cut_face_node_indices_by_elem_by_face) = mt.cut_cell_data
Qxyh, hybridized_project_interp_matrices, 
hybridized_projection_matrices, hybridized_interp_matrices = 
    hybridized_SBP_operators(md)

# for e in eachindex(Qxyh)
#     Dx, Dy = differentiation_matrices[e]
#     M = mass_matrices[e]
#     Qxh, Qyh = Qxyh[e]
#     Vh = hybridized_interp_matrices[e]

#     weak_sbp_error = norm(sum(Qxh, dims=2))
#     accuracy_error = norm(Vh' * Qxh * Vh - M * Dx) + norm(Vh' * Qyh * Vh - M * Dy)
#     @show weak_sbp_error, accuracy_error
# end

(; x, y, xq, yq) = md
u = @. sin(pi * x) * sin(pi * y)
dudx_exact = @. pi * cos(pi * xq) * sin(pi * yq)

dudx = similar(md.xq)
dudx.cartesian .= rd.Vq * ((md.rxJ.cartesian .* (rd.Dr * u.cartesian) + 
                            md.sxJ.cartesian .* (rd.Ds * u.cartesian)) ./ md.J.cartesian)
for e in axes(u.cut, 2)
    Dx, Dy = differentiation_matrices[e]
    Qxh, Qyh = Qxyh[e]
    Vh = hybridized_interp_matrices[e]
    Ph = hybridized_projection_matrices[e]
    Vq = volume_interpolation_matrices[e]
    dudx.cut[:, e] .= Vq * Ph * Qxh * Vh * u.cut[:, e]
end

@show sqrt(sum(md.wJq .* (dudx - dudx_exact).^2))
# scatter(md.xyz..., zcolor = dudx - dudx_exact, leg=false, ratio=1)