jlchan / StartUpDG.jl

Initializes and sets up reference elements and physical meshes for DG.
MIT License
27 stars 9 forks source link

Add new cut cell `MeshData` #165

Closed jlchan closed 2 months ago

jlchan commented 2 months ago

This constructs a MeshData for cut cell meshes using subtriangulations and Caratheodory pruning to construct quadrature rules with positive weights on cut cells.

jlchan commented 2 months ago

The following is a script to plot the pruned quadrature nodes:

using Plots
using StartUpDG
using PathIntersections

N = 4
quad_rule_face = gauss_lobatto_quad(0, 0, N)
rd = RefElemData(Quad(), N; quad_rule_face)

cells_per_dimension = 2
circle = PresetGeometries.Circle(R=0.66, x0=0, y0=0)
objects = (circle, )

md = MeshData(rd, objects, cells_per_dimension)
(; cutcells) = md.mesh_type.cut_cell_data

# volume quadrature should be exact for degree N(N-1) + 2N-2 polynomials, 
# or N(N-1) + 2(N-1) = (N+2) * (N-1) degrees
N_phys_frame_geo = N^2 + N * (N-1) + 2 * (N-1) 

target_degree = 2 * N
rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, 
                    quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo))

plot()
for e in eachindex(cutcells)
    xq, yq, _ = StartUpDG.subtriangulated_cutcell_quadrature(cutcells[e], rd_tri)
    scatter!(vec(xq), vec(yq), label="Reference quadrature"); 
    scatter!(md.xq.cut[:, e], md.yq.cut[:, e], markersize=8, marker=:circle, 
             z_order=:back, label="Caratheodory pruning", leg=false)
end
display(plot!(leg=false))