Isoparametric cut-cell meshes

plot_isoparametric_cut_mesh.jl

@@ -0,0 +1,230 @@
+using NodesAndModes: face_basis
+using Plots
+using StartUpDG
+
+N = 4
+rd = RefElemData(Quad(), N, quad_rule_face = gauss_quad(0, 0, 2 * N))
+
+cells_per_dimension = 2
+circle = PresetGeometries.Circle(R=0.66, x0=0, y0=0)
+md = MeshData(rd, (circle, ), cells_per_dimension; precompute_operators=true)
+
+mt = md.mesh_type
+(; cutcells) = mt.cut_cell_data
+
+using Triangulate, StaticArrays
+function triangulate_points(coordinates::AbstractMatrix)
+    triin=Triangulate.TriangulateIO()
+    triin.pointlist = coordinates
+    triout, _ = triangulate("Q", triin)
+    VX, VY = (triout.pointlist[i,:] for i = 1:size(triout.pointlist,1))
+    EToV = permutedims(triout.trianglelist)
+    return (VX, VY), EToV
+end
+
+N_phys_frame_geo = max(2 * rd.N, (rd.N-1) * (rd.N + 2)) # (N-1) * (N + 2)
+
+rd_line = RefElemData(Line(), N=rd.N, quad_rule_vol = gauss_quad(0, 0, N))
+
+rd_tri = RefElemData(Tri(), N=rd.N, quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo))
+fv = NodesAndModes.face_vertices(Tri())
+tri_face_coords = (rd_tri.r[vec(rd_tri.Fmask)], rd_tri.s[vec(rd_tri.Fmask)])
+
+# preallocate face interp node arrays
+r1D = rd_line.r
+tri_face_coords_x = zeros(length(r1D), 3)
+tri_face_coords_y = zeros(length(r1D), 3)
+
+# this operator performs a least squares fit, and is equivalent to isoparametric warp and blend
+warp_face_points_to_interp = 
+    face_basis(Tri(), rd_tri.N, rd_tri.rst...) / face_basis(Tri(), rd_tri.N, tri_face_coords...) 
+
+using StartUpDG: map_to_interval
+
+function compute_geometric_determinant_J(x, y, Dr, Ds)
+    xr, xs = Dr * x, Ds * x
+    yr, ys = Dr * y, Ds * y
+    J = @. -xs * yr + xr * ys
+    return J
+end
+
+x_cutcells, y_cutcells, xq_cutcells, yq_cutcells, Jq_cutcells, wJq_cutcells = 
+    ntuple(_ -> Matrix{eltype(rd_tri.wq)}[], 6)
+
+for cutcell in cutcells
+    # create subtriangulation. note that `cutcell.stop_pts[end] == cutcell.stop_pts[1]`
+    vxy = zeros(2, length(cutcell.stop_pts[1:end-1]))
+    for (i, pt) in enumerate(cutcell.stop_pts[1:end-1])
+        vxy[1, i], vxy[2, i] = cutcell(pt)    
+    end
+    (VX, VY), EToV = triangulate_points(vxy)
+
+    xq, yq, Jq, wJq = ntuple(_ -> zeros(length(rd_tri.wq), size(EToV, 1)), 6)
+
+    # loop over each triangle, map 1D interpolation points to faces
+    for e in axes(EToV, 1)
+        ids = view(EToV, e, :)
+        for (face_index, face_vertices) in enumerate(SVector{2}.(fv))
+            vertices_on_face = sort(ids[face_vertices])
+
+            # map each point to a physical element. 
+            for i in eachindex(r1D)
+                # This assumes a PathIntersections.jl ordering of curve points.
+                # If the vertex indices are far apart, it's the last face/boundary curve
+                if (x->abs(x[2]-x[1]))(vertices_on_face) == length(VX) - 1 
+                    s = map_to_interval(r1D[i], cutcell.stop_pts[end-1:end]...)
+                    point = cutcell(s)
+
+                # if vertex indices are consecutive, it's a boundary face    
+                elseif (x->x[2]-x[1])(vertices_on_face) == 1 
+
+                    curve_id = minimum(ids[face_vertices])
+                    s = map_to_interval(r1D[i], cutcell.stop_pts[curve_id:curve_id+1]...)
+                    point = cutcell(s)
+
+                else # it's an internal face
+                    point = SVector{2}.(map_to_interval(r1D[i], VX[ids[face_vertices]]...),
+                                        map_to_interval(r1D[i], VY[ids[face_vertices]]...))
+                end
+                tri_face_coords_x[i, face_index] = point[1]
+                tri_face_coords_y[i, face_index] = point[2]
+            end
+        end
+
+        # this performs a least squares fit interpolation in the face basis, but is equivalent
+        # to isoparametric warp and blend if the face node locations are continuous.
+        tri_warped_coords_x = warp_face_points_to_interp * vec(tri_face_coords_x) 
+        tri_warped_coords_y = warp_face_points_to_interp * vec(tri_face_coords_y)
+
+        Jq_e = abs.(compute_geometric_determinant_J(tri_warped_coords_x, tri_warped_coords_y, 
+                                                    rd_tri.Vq * rd_tri.Dr, rd_tri.Vq * rd_tri.Ds))                                                            
+
+        view(xq, :, e) .= rd_tri.Vq * tri_warped_coords_x
+        view(yq, :, e) .= rd_tri.Vq * tri_warped_coords_y
+        view(Jq, :, e) .= Jq_e
+        @. wJq[:,e] = rd_tri.wq * Jq_e
+    end
+    push!(xq_cutcells, xq)
+    push!(yq_cutcells, yq)
+    push!(Jq_cutcells, Jq)
+    push!(wJq_cutcells, wJq)
+end
+
+# Caratheodory pruning
+function basic_removal(V, w_in)
+
+    if length(w_in) <= size(V, 2)
+        return w_in, eachindex(w_in)
+    end
+    w = copy(w_in)
+    M, N = size(V)
+    inds = collect(1:M)
+    m = M-N
+    Q, _ = qr(V)
+    Q = copy(Q)
+    for _ in 1:m
+        kvec = Q[:,end]
+
+        # for subtracting the kernel vector
+        idp = findall(@. kvec > 0)
+        alphap, k0p = findmin(w[inds[idp]] ./ kvec[idp])
+        k0p = idp[k0p]
+
+        # for adding the kernel vector
+        idn = findall(@. kvec < 0);
+        alphan, k0n = findmax(w[inds[idn]] ./ kvec[idn])
+        k0n = idn[k0n];
+
+        alpha, k0 = abs(alphan) < abs(alphap) ? (alphan, k0n) : (alphap, k0p)
+        w[inds] = w[inds] - alpha * kvec
+        deleteat!(inds, k0)
+        Q, _ = qr(V[inds, :])
+        Q = copy(Q)
+    end
+    return w, inds
+end
+
+xq_pruned, yq_pruned, wJq_pruned = ntuple(_ -> Vector{Float64}[], 3)
+for e in eachindex(xq_cutcells)
+    V2N = vandermonde(mt.physical_frame_elements[e], 2 * rd.N, vec(xq_cutcells[e]), vec(yq_cutcells[e]))
+    w = copy(vec(wJq_cutcells[e]))
+    w_pruned, inds = basic_removal(V2N, w)
+
+    V = vandermonde(mt.physical_frame_elements[e], rd.N, vec(xq_cutcells[e]), vec(yq_cutcells[e]))
+    # @show size(V[inds,:])
+    # @show length(w), length(inds)
+    @show norm(V' * diagm(w) * V - V' * diagm(w_pruned) * V)
+
+    push!(yq_pruned, vec(yq_cutcells[e])[inds])
+    push!(xq_pruned, vec(xq_cutcells[e])[inds])
+    push!(wJq_pruned, vec(wJq_cutcells[e])[inds])
+end
+
+plot()
+for e in eachindex(xq_cutcells, yq_cutcells)
+    scatter!(vec(xq_cutcells[e]), vec(yq_cutcells[e]), label="Reference quadrature"); 
+    scatter!(xq_pruned[e], yq_pruned[e], markersize=8, marker=:circle, 
+             z_order=:back, label="Caratheodory pruning", leg=false)
+end
+display(plot!(leg=false))
+
+e = 1
+path = md.mesh_type.cut_cell_data.cutcells[e]
+xy = path.(LinRange(0,1,2000))
+plot(getindex.(xy, 1),getindex.(xy, 2), linewidth=2, color=:black, label="")
+
+scatter!(vec(xq_cutcells[e]), vec(yq_cutcells[e]), label="Reference quadrature"); 
+scatter!(xq_pruned[e], yq_pruned[e], markersize=8, marker=:circle, 
+         z_order=:back, label="Caratheodory pruning", axis=([], false), ratio=1)
+
+
+
+# # compute normals
+# cutcell = cutcells[1]
+# plot()
+# xf, yf, nxJ, nyJ = ntuple(_ -> zeros(size(rd_line.Vq, 1), length(cutcell.stop_pts)-1), 4)
+# for f in 1:length(cutcell.stop_pts)-1
+#     points = map(s -> cutcell(map_to_interval(s, cutcell.stop_pts[f:f+1]...)), r1D)
+#     x = getindex.(points, 1)
+#     y = getindex.(points, 2)
+
+#     # compute tangent vector
+#     (; Vq, Dr) = rd_line
+#     dxdr = Vq * Dr * x
+#     dydr = Vq * Dr * y
+
+#     tangent_vector = SVector.(dxdr, dydr)
+#     scaling = (cutcell.stop_pts[f+1] - cutcell.stop_pts[f]) / 2
+#     Jf = norm.(tangent_vector) .* scaling
+#     raw_normal = SVector.(-dydr, dxdr)
+#     scaled_normal = (raw_normal) / norm(raw_normal) .* Jf
+
+#     @. nxJ[:, f] = getindex(scaled_normal, 1)
+#     @. nyJ[:, f] = getindex(scaled_normal, 2)
+
+#     # interp face coordinates to face quad nodes
+#     xf[:, f] .= Vq * x
+#     yf[:, f] .= Vq * y
+
+#     scatter!(Vq * x, Vq * y)
+#     quiver!(Vq * x, Vq * y, quiver=(getindex.(scaled_normal, 1), getindex.(scaled_normal, 2)))
+# end
+# plot!(leg=false, ratio=1)
+
+# # test weak SBP property
+# (; x, y) = md
+# elem = md.mesh_type.physical_frame_elements[1]
+# VDM = vandermonde(elem, rd.N, x.cut[:, 1], y.cut[:, 1])
+# Vq, Vxq, Vyq = map(A -> A / VDM, basis(elem, rd.N, xq_pruned[1], yq_pruned[1]))
+# M = Vq' * diagm(wJq_pruned[1]) * Vq
+# Qx = Vq' * diagm(wJq_pruned[1]) * Vxq
+# Vf = vandermonde(elem, rd.N, xf, yf) / VDM
+
+# Dx, Dy = md.mesh_type.cut_cell_operators.differentiation_matrices[1]
+# Qx, Qy = M * Dx, M * Dy
+
+# Bx = Diagonal(vec(Diagonal(rd_line.wq) * nxJ))
+# # Bx = Diagonal(vec(Diagonal(rd_line.wq) * reshape(md.nxJ.cut[md.mesh_type.cut_face_nodes[1]], length(rd_line.wq), :)))
+
+# e = ones(size(Vf, 2))
+# [e' * Qx; e' * Vf' * Bx * Vf]'

src/cut_cell_meshes.jl

@@ -29,14 +29,6 @@ struct CutCellMesh{T1, T2, T3, T4, T5}
 end
 
 # TODO: add isoparametric cut cell mesh with positive quadrature points
-# # This mesh type has a polynomial representation of objects, so we don't store the curve info
-# struct IsoparametricCutCellMesh{T1, T2, T3, T4}
-#     physical_frame_elements::T1
-#     cut_face_nodes::T2
-#     cut_cell_operators::T3
-#     cut_cell_data::T4
-# end
-
 
 function Base.show(io::IO, ::MIME"text/plain", md::MeshData{DIM, <:CutCellMesh}) where {DIM}
     @nospecialize md

test_cut_derivative.jl

@@ -0,0 +1,51 @@
+using StartUpDG
+
+cells_per_dimension = 32
+circle = PresetGeometries.Circle(R=0.66, x0=0, y0=0)
+
+rd = RefElemData(Quad(), N=3)
+md = MeshData(rd, (circle, ), cells_per_dimension; precompute_operators=true)
+
+(; differentiation_matrices, lift_matrices, face_interpolation_matrices) = 
+    md.mesh_type.cut_cell_operators
+
+(; x, y) = md
+
+u_exact(x, y) = sin(pi * x) * sin(pi * y)
+dudx_exact(x, y) = pi * cos(pi * x) * sin(pi * y)
+dudy_exact(x, y) = pi * sin(pi * x) * cos(pi * y)
+
+(; N) = rd
+u_exact(x, y) = x^N + y^N
+dudx_exact(x, y) = N * x^(N-1)
+dudy_exact(x, y) = N * y^(N-1)
+
+u = u_exact.(x, y)
+(; physical_frame_elements, cut_face_nodes) = md.mesh_type
+
+uf = similar(md.xf)
+uf.cartesian = rd.Vf * u.cartesian
+for e in eachindex(face_interpolation_matrices)
+    ids = cut_face_nodes[e]
+    Vf = face_interpolation_matrices[e]
+    uf.cut[ids] = Vf * u.cut[:, e]
+end
+
+uP = vec(uf[md.mapP])
+flux = @. 0.5 * (uP - uf) 
+
+dudx, dudy = similar(md.x), similar(md.x)     
+dudx.cartesian .= (md.rxJ.cartesian .* (rd.Dr * u.cartesian)) ./ md.J
+dudy.cartesian .= (md.syJ.cartesian .* (rd.Ds * u.cartesian)) ./ md.J
+for (e, elem) in enumerate(physical_frame_elements)
+    Dx, Dy = differentiation_matrices[e]
+    LIFT = lift_matrices[e]
+    ids = cut_face_nodes[e]
+    dudx.cut[:, e] .= Dx * u.cut[:,e] + LIFT * (flux[ids] .* md.nxJ.cut[ids])
+    dudy.cut[:, e] .= Dy * u.cut[:,e] + LIFT * (flux[ids] .* md.nyJ.cut[ids])
+end
+
+@show norm(dudx - dudx_exact.(x,y), Inf)
+@show norm(dudy - dudy_exact.(x,y), Inf)
+
+# scatter(md.xyz..., dudx - dudx_exact.(x,y), zcolor=dudx - dudx_exact.(x,y), leg=false)

test_cutcell_convergence.jl

@@ -0,0 +1,74 @@
+using Plots
+using PathIntersections
+using StartUpDG
+
+function compute_L2_error(N, cells_per_dimension; 
+                          u_exact = (x,y) -> sin(pi * x) * sin(pi * y),
+                          use_srd=false, use_interp=false)
+    rd = RefElemData(Quad(), N, quad_rule_vol=quad_nodes(Quad(), N+1))    
+    objects = (PresetGeometries.Circle(R=0.3),)
+    md = MeshData(rd, objects, cells_per_dimension, cells_per_dimension; precompute_operators=true)
+
+    srd = StateRedistribution(rd, md)
+
+    (; physical_frame_elements, cut_face_nodes) = md.mesh_type
+    Vq_cut, Pq_cut, M_cut = ntuple(_ -> Matrix{Float64}[], 3)
+    for (e, elem) in enumerate(physical_frame_elements)
+
+        VDM = vandermonde(elem, rd.N, md.x.cut[:, e], md.y.cut[:, e]) # TODO: should these be md.x, md.y?
+        Vq, _ = map(A -> A / VDM, basis(elem, rd.N, md.xq.cut[:,e], md.yq.cut[:, e]))
+
+        M = Vq' * diagm(md.wJq.cut[:, e]) * Vq        
+        push!(M_cut, M)
+        push!(Vq_cut, Vq)
+        push!(Pq_cut, M \ (Vq' * diagm(md.wJq.cut[:, e])))
+    end
+
+    if use_interp==true 
+        u = u_exact.(md.xyz...)
+    else # projection
+        uq = u_exact.(md.xyzq...)
+        u = similar(md.x)
+        u.cartesian .= rd.Pq * uq.cartesian
+        for e in eachindex(physical_frame_elements)
+            u.cut[:, e] .= Pq_cut[e] * uq.cut[:, e]
+        end
+    end
+
+    if use_srd == true
+        srd(u)
+    end
+
+    # eval solution at quad points
+    uq = similar(md.xq)
+    uq.cartesian .= rd.Vq * u.cartesian
+    for e in eachindex(physical_frame_elements)
+        uq.cut[:, e] .= Vq_cut[e] * u.cut[:, e]
+    end
+
+    L2err = sqrt(abs(sum(md.wJq .* (uq - u_exact.(md.xyzq...)).^2)))
+    return L2err
+end
+
+N = 3
+num_cells = [4, 8, 16, 32, 64]
+
+L2_error, L2_error_srd, L2_error_interp, L2_error_interp_srd  = ntuple(_ -> Float64[], 4)
+for cells_per_dimension in num_cells
+    @show cells_per_dimension
+    use_interp = true
+    push!(L2_error_interp, compute_L2_error(N, cells_per_dimension; use_interp))
+    push!(L2_error_interp_srd, compute_L2_error(N, cells_per_dimension; use_interp, use_srd=true))
+    use_interp = false
+    push!(L2_error, compute_L2_error(N, cells_per_dimension; use_interp))
+    push!(L2_error_srd, compute_L2_error(N, cells_per_dimension; use_interp, use_srd=true))
+end
+
+h = 2 ./ num_cells
+plot()
+plot!(h, L2_error_interp, marker=:dot, label="Interp")
+plot!(h, L2_error_interp_srd, marker=:dot, label="Interp (SRD)")
+plot!(h, L2_error, marker=:dot, label="Projection")
+plot!(h, L2_error_srd, marker=:dot, label="Projection (SRD)")
+plot!(h, 1e-1*h.^(N+1), linestyle=:dash, label="h^{N+1}")
+plot!(xaxis=:log, yaxis=:log, legend = :topleft)