From 8e01c266d60a52e57de1c56a6a0cfd3cf118c49c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 22 Feb 2023 13:28:22 -0600 Subject: [PATCH 01/39] add some temporary test files --- test_cut_derivative.jl | 51 +++++++ test_cutcell_convergence.jl | 74 ++++++++++ test_cutcell_sbp.jl | 240 +++++++++++++++++++++++++++++++++ test_face_basis.jl | 36 +++++ test_isoparametric_cut_mesh.jl | 227 +++++++++++++++++++++++++++++++ 5 files changed, 628 insertions(+) create mode 100644 test_cut_derivative.jl create mode 100644 test_cutcell_convergence.jl create mode 100644 test_cutcell_sbp.jl create mode 100644 test_face_basis.jl create mode 100644 test_isoparametric_cut_mesh.jl diff --git a/test_cut_derivative.jl b/test_cut_derivative.jl new file mode 100644 index 00000000..e5274524 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/test_cutcell_convergence.jl b/test_cutcell_convergence.jl new file mode 100644 index 00000000..3e84d4ff --- /dev/null +++ b/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) \ No newline at end of file diff --git a/test_cutcell_sbp.jl b/test_cutcell_sbp.jl new file mode 100644 index 00000000..038044dd --- /dev/null +++ b/test_cutcell_sbp.jl @@ -0,0 +1,240 @@ +using NodesAndModes: face_basis +using Plots +using StartUpDG + +N = 4 +quad_rule_face = gauss_quad(0, 0, 2*N+1) + +rd = RefElemData(Quad(), N; quad_rule_face) + +# objects = (PresetGeometries.Rectangle(Lx=.4, Ly=2.5, x0=1.2), ) +# vx = [-1, 0, 1.1] +# vy = [-1, 0, 1] +# md = MeshData(rd, objects, vx, vy; precompute_operators=true) + +cells_per_dimension = 2 +circle = PresetGeometries.Circle(R=0.66, x0=0, y0=0) +md = MeshData(rd, (circle, ), cells_per_dimension; quad_rule_face)#, 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 + +# degree N(N-1) + 2N-2 polynomial = N(N-1) + 2(N-1) = (N+2) * (N-1) +N_phys_frame_geo = max(2 * rd.N, (rd.N-1) * (rd.N + 2)) + +rd_line = RefElemData(Line(), N=rd.N, quad_rule_vol = quad_rule_face) + +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) + 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)) + + +# recompute 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) + +# # compute new normals for curved boundaries +# mapB_boundary_cartesian = findall(@. abs(abs(md.xf) - 1) < 1e2 * eps() || abs(abs(md.yf) - 1) < 1e4 * eps()) +# mapB_cut = setdiff(md.mapB, mapB_boundary_cartesian) +# num_face_nodes = length(quad_rule_face[1]) +# xf, yf = map(x -> reshape(x[mapB_cut], num_face_nodes, :), md.xyzf) +# D1D = rd_line.Vq * rd_line.Dr * rd_line.Pq +# dxdr, dydr = D1D * xf, D1D * yf +# nxJ_boundary_cut, nyJ_boundary_cut = dydr, -dxdr + +# 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]) +# xq, yq, wq = xq_pruned[1], yq_pruned[1], wJq_pruned[1] +xq, yq, wq = vec.((xq_cutcells[1], yq_cutcells[1], wJq_cutcells[1])) +Vq, Vxq, Vyq = map(A -> A / VDM, basis(elem, rd.N, xq, yq)) +M = Vq' * diagm(wq) * Vq +Qx, Qy = Vq' * diagm(wq) * Vxq, Vq' * diagm(wq) * Vyq +Vf = vandermonde(elem, rd.N, xf, yf) / VDM + +# Bx = Diagonal(vec(Diagonal(rd_line.wq) * reshape(md.nxJ.cut[md.mesh_type.cut_face_nodes[1]], length(rd_line.wq), :))) +# By = Diagonal(vec(Diagonal(rd_line.wq) * reshape(md.nyJ.cut[md.mesh_type.cut_face_nodes[1]], length(rd_line.wq), :))) +Bx = Diagonal(vec(Diagonal(rd_line.wq) * reshape(nxJ, length(rd_line.wq), :))) +By = Diagonal(vec(Diagonal(rd_line.wq) * reshape(nyJ, length(rd_line.wq), :))) + +e = ones(size(Qx, 2)) +display(sum(Qx - Vf' * Bx * Vf, dims=1)) +# display(e' * ((Qx + Qx') - Vf' * Bx * Vf)) +# display(e' * ((Qy + Qy') - Vf' * By * Vf)) diff --git a/test_face_basis.jl b/test_face_basis.jl new file mode 100644 index 00000000..fd7b29ec --- /dev/null +++ b/test_face_basis.jl @@ -0,0 +1,36 @@ +using NodesAndModes: face_basis +using StartUpDG + +N=7 +rd = RefElemData(Tri(), N) +md_ref = MeshData(uniform_mesh(Tri(), 1)..., rd) +(; x, y) = md_ref +# x = @. x + 0.5 * cos(pi/2 * x) * cos(pi/2 * y) +# y = @. y + 0.5 * cos(pi/2 * x) * cos(pi/2 * y) + +x = @. x + 0.25 * cos(pi/2 * x) * sin(pi/2 * y + 1) +y = @. y + 0.25 * cos(pi/2 * x + 2) * cos(pi/2 * y) + +# x = x + 0.1 * randn(size(x)) +# y = y + 0.1 * randn(size(x)) +(; Dr, Ds) = rd +xr, xs = Dr * x, Ds * x +yr, ys = Dr * y, Ds * y + +# interpolate to quadrature directly +xrq, xsq, yrq, ysq = (x -> rd.Vq * x).((xr, xs, yr, ys)) +Jq = @. -xsq * yrq + xrq * ysq +wJq = Diagonal(rd.wq) * Jq + +md = MeshData(rd, md_ref, x, y) + +@show sum(wJq), sum(md.wJq) + +Vf_face = face_basis(Tri(), rd.N, rd.rstf...) +norm(Vf_face * (Vf_face \ md.xf) - md.xf) + +@show sum(md_ref.wJq), sum(md.wJq) +@show sum(md_ref.wJq .* md_ref.xq), sum(md.wJq .* md.xq) +@show sum(md_ref.wJq .* md_ref.xq.^2), sum(md.wJq .* md.xq.^2) + + diff --git a/test_isoparametric_cut_mesh.jl b/test_isoparametric_cut_mesh.jl new file mode 100644 index 00000000..e28cf1c0 --- /dev/null +++ b/test_isoparametric_cut_mesh.jl @@ -0,0 +1,227 @@ +using NodesAndModes: face_basis +using Plots +using StartUpDG + + +rd = RefElemData(Quad(), N=1, quad_rule_face = gauss_quad(0, 0, 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) + 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)) + + +# 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)) +@show e' * (Qx + Qx') +@show e' * Vf' * Bx * Vf + +# 3×3 Matrix{Float64}: +# 0.507422 -0.103253 0.253711 +# -0.103253 -0.300917 -0.253711 +# 0.253711 -0.253711 6.80375e-18 + From 76a01cf1409edb095c67f91b7633b6f0018f2406 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 20 Apr 2023 09:36:32 -0500 Subject: [PATCH 02/39] fix docstrings --- src/MeshData.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/MeshData.jl b/src/MeshData.jl index 078af19b..9aa864a5 100644 --- a/src/MeshData.jl +++ b/src/MeshData.jl @@ -201,6 +201,8 @@ end num_elements(md::MeshData) = size(md.x, 2) # number of columns in the "x" coordinate array num_elements(md::MeshData{Dim, <:VertexMappedMesh}) where {Dim} = size(md.mesh_type.EToV, 1) +# splat `uniform_mesh` arguments, e.g., enables `MeshData(uniform_mesh(Line(), 1), rd)` +# TODO: wrap `uniform_mesh` in a custom type so we can dispatch more precisely """ MeshData(VXYZ, EToV, rd::RefElemData) @@ -213,14 +215,9 @@ Given new nodal positions `xyz...` (e.g., from mesh curving), recomputes geometr and outputs a new MeshData struct. Only fields modified are the coordinate-dependent terms `xyz`, `xyzf`, `xyzq`, `rstxyzJ`, `J`, `nxyzJ`, `Jf`. """ - -# splat `uniform_mesh` arguments, e.g., enables `MeshData(uniform_mesh(Line(), 1), rd)` -# TODO: wrap `uniform_mesh` in a custom type so we can dispatch more precisely MeshData(mesh::Tuple{<:Tuple, Matrix{Int64}}, other_args...) = MeshData(mesh..., other_args...) - -# splats VXYZ MeshData(VXYZ::T, EToV, other_args...) where {NDIMS, T <: NTuple{NDIMS}} = - MeshData(VXYZ..., EToV, other_args...) + MeshData(VXYZ..., EToV, other_args...) # splats VXYZ function MeshData(VX::AbstractVector{Tv}, EToV, rd::RefElemData{1}) where {Tv} From 6cbffe9834e85c504275bd607c25b592bcef36c0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 25 Apr 2023 14:12:29 -0500 Subject: [PATCH 03/39] remove cruft --- src/cut_cell_meshes.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 34e3f05a..700872bf 100644 --- a/src/cut_cell_meshes.jl +++ b/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 From caa612e52ec60ecd77fc29a3a497f30175ee45c2 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 1 Dec 2023 09:32:28 -0600 Subject: [PATCH 04/39] update plot made nicer plot of Caratheodory pruning for Christina's proposal --- test_isoparametric_cut_mesh.jl | 113 +++++++++++++++++---------------- 1 file changed, 58 insertions(+), 55 deletions(-) diff --git a/test_isoparametric_cut_mesh.jl b/test_isoparametric_cut_mesh.jl index e28cf1c0..6ac4ba05 100644 --- a/test_isoparametric_cut_mesh.jl +++ b/test_isoparametric_cut_mesh.jl @@ -2,8 +2,8 @@ using NodesAndModes: face_basis using Plots using StartUpDG - -rd = RefElemData(Quad(), N=1, quad_rule_face = gauss_quad(0, 0, N)) +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) @@ -161,67 +161,70 @@ for e in eachindex(xq_cutcells) end plot() -for e in eachindex(xq_cutcells) +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)) - -# 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 +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)) + + + +# # 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)) -@show e' * (Qx + Qx') -@show e' * Vf' * Bx * Vf - -# 3×3 Matrix{Float64}: -# 0.507422 -0.103253 0.253711 -# -0.103253 -0.300917 -0.253711 -# 0.253711 -0.253711 6.80375e-18 +# 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]' From d8fe5f59362f3922a12d799bb98845788c45dd0b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 1 Dec 2023 11:22:21 -0600 Subject: [PATCH 05/39] update file name Makes it mroe clear this file is meant for plotting --- ..._isoparametric_cut_mesh.jl => plot_isoparametric_cut_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename test_isoparametric_cut_mesh.jl => plot_isoparametric_cut_mesh.jl (99%) diff --git a/test_isoparametric_cut_mesh.jl b/plot_isoparametric_cut_mesh.jl similarity index 99% rename from test_isoparametric_cut_mesh.jl rename to plot_isoparametric_cut_mesh.jl index 6ac4ba05..9efe437b 100644 --- a/test_isoparametric_cut_mesh.jl +++ b/plot_isoparametric_cut_mesh.jl @@ -175,7 +175,7 @@ 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)) + z_order=:back, label="Caratheodory pruning", axis=([], false), ratio=1) From f30183334009205f11cbd3d593fc73ef0560a18a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Apr 2024 11:10:38 -0500 Subject: [PATCH 06/39] fixing cut cell demo for v1.0+ --- test_cutcell_sbp.jl | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/test_cutcell_sbp.jl b/test_cutcell_sbp.jl index 038044dd..260162f0 100644 --- a/test_cutcell_sbp.jl +++ b/test_cutcell_sbp.jl @@ -1,9 +1,10 @@ using NodesAndModes: face_basis using Plots +using LinearAlgebra using StartUpDG -N = 4 -quad_rule_face = gauss_quad(0, 0, 2*N+1) +N = 3 +quad_rule_face = gauss_quad(0, 0, 2 * N + 1) rd = RefElemData(Quad(), N; quad_rule_face) @@ -76,21 +77,21 @@ for cutcell in cutcells vertices_on_face = sort(ids[face_vertices]) # map each point to a physical element. + # This assumes PathIntersections.jl uses a clockwise ordering of stop curve points. 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]...) + s = map_to_interval(r1D[i], reverse(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]...) + s = map_to_interval(r1D[i], reverse(cutcell.stop_pts[curve_id:curve_id+1])...) point = cutcell(s) - else # it's an internal face + else # it's a non-boundary face, it's a straight line point = SVector{2}.(map_to_interval(r1D[i], VX[ids[face_vertices]]...), map_to_interval(r1D[i], VY[ids[face_vertices]]...)) end @@ -104,8 +105,8 @@ for cutcell in cutcells 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)) + Jq_e = 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 @@ -186,7 +187,7 @@ for f in 1:length(cutcell.stop_pts)-1 x = getindex.(points, 1) y = getindex.(points, 2) - # compute tangent vector + # compute tangent vector using polynomial mapping (; Vq, Dr) = rd_line dxdr = Vq * Dr * x dydr = Vq * Dr * y From ba2756abc44eb19d1f1b26a8646ee4783fc533f0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Apr 2024 14:39:05 -0500 Subject: [PATCH 07/39] adding Caratheodory pruning --- src/physical_frame_basis.jl | 42 +++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/src/physical_frame_basis.jl b/src/physical_frame_basis.jl index e81870aa..add953d8 100644 --- a/src/physical_frame_basis.jl +++ b/src/physical_frame_basis.jl @@ -144,4 +144,46 @@ function map_nodes_to_background_cell(elem::PhysicalFrame{2}, r, s) x = @. 0.5 * (1 + r) * dx + vx[1] y = @. 0.5 * (1 + s) * dy + vy[1] return x, y +end + +""" + caratheodory_pruning_qr(V, w_in) + +This performs Caratheodory pruning using a naive QR-based algorithm. +Returns (w, inds), where `inds` denotes sub-selected indices for a +reduced quadrature rule, and `w` is a vector of reduced positive weights. + +The original Matlab code this was based on was authored by Akil Narayan. +""" +function caratheodory_pruning_qr(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 \ No newline at end of file From cfc059728de8bfcb4147fccf1a729a7d919a3457 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Apr 2024 14:39:10 -0500 Subject: [PATCH 08/39] cleaning up cutcell demo --- test_cutcell_sbp.jl | 85 ++++++++++++++------------------------------- 1 file changed, 27 insertions(+), 58 deletions(-) diff --git a/test_cutcell_sbp.jl b/test_cutcell_sbp.jl index 260162f0..bf4c0565 100644 --- a/test_cutcell_sbp.jl +++ b/test_cutcell_sbp.jl @@ -78,6 +78,7 @@ for cutcell in cutcells # map each point to a physical element. # This assumes PathIntersections.jl uses a clockwise ordering of stop curve points. + # Since StartUpDG uses a CCW ordering, we reverse the order for for i in eachindex(r1D) # 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 @@ -105,8 +106,9 @@ for cutcell in cutcells 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 = compute_geometric_determinant_J(tri_warped_coords_x, tri_warped_coords_y, - rd_tri.Vq * rd_tri.Dr, rd_tri.Vq * rd_tri.Ds) + rxJq_e, sxJq_e, ryJq_e, syJq_e, Jq_e = + StartUpDG.geometric_factors(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 @@ -119,45 +121,11 @@ for cutcell in cutcells 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) + w_pruned, inds = StartUpDG.caratheodory_pruning_qr(V2N, w) V = vandermonde(mt.physical_frame_elements[e], rd.N, vec(xq_cutcells[e]), vec(yq_cutcells[e])) # @show size(V[inds,:]) @@ -169,16 +137,7 @@ for e in eachindex(xq_cutcells) push!(wJq_pruned, vec(wJq_cutcells[e])[inds]) end -plot() -for e in eachindex(xq_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)) - - -# recompute normals +# recompute normals on cut cells cutcell = cutcells[1] plot() xf, yf, nxJ, nyJ = ntuple(_ -> zeros(size(rd_line.Vq, 1), length(cutcell.stop_pts)-1), 4) @@ -188,28 +147,36 @@ for f in 1:length(cutcell.stop_pts)-1 y = getindex.(points, 2) # compute tangent vector using polynomial mapping - (; Vq, Dr) = rd_line - dxdr = Vq * Dr * x - dydr = Vq * Dr * y + dxdr = rd_line.Vq * rd_line.Dr * x + dydr = rd_line.Vq * rd_line.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 + scaled_normal = SVector.(-dydr, dxdr) @. 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 + xf[:, f] .= rd_line.Vq * x + yf[:, f] .= rd_line.Vq * y - scatter!(Vq * x, Vq * y) - quiver!(Vq * x, Vq * y, quiver=(getindex.(scaled_normal, 1), getindex.(scaled_normal, 2))) + scatter!(rd_line.Vq * x, rd_line.Vq * y) + quiver!(rd_line.Vq * x, rd_line.Vq * y, + quiver=(getindex.(scaled_normal, 1), getindex.(scaled_normal, 2))) end plot!(leg=false, ratio=1) +plot() +for e in eachindex(xq_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)) + + + # # compute new normals for curved boundaries # mapB_boundary_cartesian = findall(@. abs(abs(md.xf) - 1) < 1e2 * eps() || abs(abs(md.yf) - 1) < 1e4 * eps()) # mapB_cut = setdiff(md.mapB, mapB_boundary_cartesian) @@ -236,6 +203,8 @@ Bx = Diagonal(vec(Diagonal(rd_line.wq) * reshape(nxJ, length(rd_line.wq), :))) By = Diagonal(vec(Diagonal(rd_line.wq) * reshape(nyJ, length(rd_line.wq), :))) e = ones(size(Qx, 2)) -display(sum(Qx - Vf' * Bx * Vf, dims=1)) +@show norm(sum(Qx - Vf' * Bx * Vf, dims=1)) +@show norm(sum(Qy - Vf' * By * Vf, dims=1)) +# display(sum(Qx - Vf' * Bx * Vf, dims=1)) # display(e' * ((Qx + Qx') - Vf' * Bx * Vf)) # display(e' * ((Qy + Qy') - Vf' * By * Vf)) From db585a854c1da745a930044ee4c1f3de053dc72f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Apr 2024 17:39:07 -0500 Subject: [PATCH 09/39] improve efficiency slightly --- src/physical_frame_basis.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/physical_frame_basis.jl b/src/physical_frame_basis.jl index add953d8..71788c1f 100644 --- a/src/physical_frame_basis.jl +++ b/src/physical_frame_basis.jl @@ -63,11 +63,13 @@ function NodesAndModes.basis(elem::PhysicalFrame{2}, N, x, y) V, Vr, Vs = ntuple(x->zeros(length(r), Np), 3) for j = 0:N P_j = jacobiP(s, 0, 0, j) + dP_j = grad_jacobiP(s, 0, 0, j) for i = 0:N-j P_i = jacobiP(r, 0, 0, i) - V[:, sk] = P_i .* P_j - Vr[:, sk] = grad_jacobiP(r, 0, 0, i) .* P_j * scaling[1] - Vs[:, sk] = P_i .* grad_jacobiP(s, 0, 0, j) * scaling[2] + dP_i = grad_jacobiP(r, 0, 0, i) + @. V[:, sk] = P_i * P_j + @. Vr[:, sk] = dP_i * P_j * scaling[1] + @. Vs[:, sk] = P_i * dP_j * scaling[2] sk += 1 end end From 4bbe94ca3503127829485745cf148afff24d28ec Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Apr 2024 17:39:49 -0500 Subject: [PATCH 10/39] refactoring functions --- src/physical_frame_basis.jl | 9 +++++++++ test_cutcell_sbp.jl | 21 +++++++-------------- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/physical_frame_basis.jl b/src/physical_frame_basis.jl index 71788c1f..55af85a0 100644 --- a/src/physical_frame_basis.jl +++ b/src/physical_frame_basis.jl @@ -148,6 +148,15 @@ function map_nodes_to_background_cell(elem::PhysicalFrame{2}, r, s) return x, y end +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 + """ caratheodory_pruning_qr(V, w_in) diff --git a/test_cutcell_sbp.jl b/test_cutcell_sbp.jl index bf4c0565..1f3acec6 100644 --- a/test_cutcell_sbp.jl +++ b/test_cutcell_sbp.jl @@ -21,14 +21,7 @@ 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 + # degree N(N-1) + 2N-2 polynomial = N(N-1) + 2(N-1) = (N+2) * (N-1) N_phys_frame_geo = max(2 * rd.N, (rd.N-1) * (rd.N + 2)) @@ -60,13 +53,13 @@ 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) + + # vxy = matrix of the form [x_coordinates, y_coordinates] + vxy = hcat(getindex.(cutcell.(cutcell.stop_pts[1:end-1]), 1), + getindex.(cutcell.(cutcell.stop_pts[1:end-1]), 2)) + (VX, VY), EToV = StartUpDG.triangulate_points(permutedims(vxy)) xq, yq, Jq, wJq = ntuple(_ -> zeros(length(rd_tri.wq), size(EToV, 1)), 6) From c6135cccd2c3c5d11821dd65246cdb2f220a383d Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Apr 2024 17:39:57 -0500 Subject: [PATCH 11/39] improve comments --- test_cutcell_sbp.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test_cutcell_sbp.jl b/test_cutcell_sbp.jl index 1f3acec6..bd6ce335 100644 --- a/test_cutcell_sbp.jl +++ b/test_cutcell_sbp.jl @@ -69,7 +69,8 @@ for cutcell in cutcells for (face_index, face_vertices) in enumerate(SVector{2}.(fv)) vertices_on_face = sort(ids[face_vertices]) - # map each point to a physical element. + # map face interpolation points to a physical element. + # This assumes PathIntersections.jl uses a clockwise ordering of stop curve points. # Since StartUpDG uses a CCW ordering, we reverse the order for for i in eachindex(r1D) @@ -94,8 +95,9 @@ for cutcell in cutcells 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. + # this performs a least squares fit interpolation by the face basis. It's + # equivalent to isoparametric warp and blend if the face node locations are + # representable by the face basis (e.g., polynomial). 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) From 201ac32bea385fb73dfa84f620b0bab0b6535608 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 25 Apr 2024 17:40:08 -0500 Subject: [PATCH 12/39] more refactoring --- test_cutcell_sbp.jl | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/test_cutcell_sbp.jl b/test_cutcell_sbp.jl index bd6ce335..5e528d42 100644 --- a/test_cutcell_sbp.jl +++ b/test_cutcell_sbp.jl @@ -134,16 +134,15 @@ end # recompute normals on cut cells 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) + xf_local = getindex.(points, 1) + yf_local = getindex.(points, 2) # compute tangent vector using polynomial mapping - dxdr = rd_line.Vq * rd_line.Dr * x - dydr = rd_line.Vq * rd_line.Dr * y + dxdr = rd_line.Vq * rd_line.Dr * xf_local + dydr = rd_line.Vq * rd_line.Dr * yf_local tangent_vector = SVector.(dxdr, dydr) scaling = (cutcell.stop_pts[f+1] - cutcell.stop_pts[f]) / 2 @@ -153,15 +152,20 @@ for f in 1:length(cutcell.stop_pts)-1 @. nyJ[:, f] = getindex(scaled_normal, 2) # interp face coordinates to face quad nodes - xf[:, f] .= rd_line.Vq * x - yf[:, f] .= rd_line.Vq * y + xf[:, f] .= rd_line.Vq * xf_local + yf[:, f] .= rd_line.Vq * yf_local + +end - scatter!(rd_line.Vq * x, rd_line.Vq * y) - quiver!(rd_line.Vq * x, rd_line.Vq * y, - quiver=(getindex.(scaled_normal, 1), getindex.(scaled_normal, 2))) +# plot normals +plot() +for f in 1:length(cutcell.stop_pts)-1 + scatter!(xf[:,f], yf[:,f]) + quiver!(xf[:,f], yf[:,f], quiver=(nxJ[:,f], nyJ[:,f])) end plot!(leg=false, ratio=1) +# plot pruned cells plot() for e in eachindex(xq_cutcells) scatter!(vec(xq_cutcells[e]), vec(yq_cutcells[e]), label="Reference quadrature"); From 9bde927e5edf3fc70a989a7d99a5532e4a391250 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 28 Apr 2024 20:28:53 -0500 Subject: [PATCH 13/39] generalize map_to_interval --- src/cut_cell_meshes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index b7b68d4d..3c4bb6e9 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -37,7 +37,7 @@ function Base.show(io::IO, ::MIME"text/plain", md::MeshData{DIM, <:CutCellMesh}) end # maps x ∈ [-1,1] to [a,b] -map_to_interval(x, a, b) = a + (b-a) * 0.5 * (1 + x) +map_to_interval(x, a, b) = @. a + (b-a) * 0.5 * (1 + x) function count_cut_faces(cutcells) num_cut_faces = zeros(Int, length(cutcells)) From 8540e8c0e2f432d44f4a5057e11648b57de21a66 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 28 Apr 2024 20:29:10 -0500 Subject: [PATCH 14/39] add new version of "generate_sampling_points" --- src/cut_cell_meshes.jl | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 3c4bb6e9..fadb648c 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -167,6 +167,38 @@ function generate_sampling_points(objects, elem, rd, Np_target; N_sampled = 4 * return x_sampled[ids_in_element], y_sampled[ids_in_element] end + +# new version assuming that elements are Quad +function generate_sampling_points(objects, elem, Np_target::Int, N_sampled::Int) + + r_sampled, s_sampled = equi_nodes(Quad(), N_sampled) # oversampled nodes + + # map sampled points to the background Cartesian cell + x_sampled, y_sampled = map_nodes_to_background_cell(elem, r_sampled, s_sampled) + is_in_domain = fill(true, length(x_sampled)) + for (index, point) in enumerate(zip(x_sampled, y_sampled)) + is_in_domain[index] = !any(map(obj -> PathIntersections.is_contained(obj, point), objects)) + end + + # increase number of background points until we are left with `Np_target` sampling points + while sum(is_in_domain) < Np_target + + N_sampled *= 2 # double degree of sampling + r_sampled, s_sampled = equi_nodes(Quad(), N_sampled) # oversampled nodes + x_sampled, y_sampled = map_nodes_to_background_cell(elem, r_sampled, s_sampled) + + # check if all the points are in all the objects + is_in_domain = fill(true, length(x_sampled)) + for (index, point) in enumerate(zip(x_sampled, y_sampled)) + is_in_domain[index] = !any(map(obj -> PathIntersections.is_contained(obj, point), objects)) + end + end + + ids_in_element = findall(is_in_domain) + + return x_sampled[ids_in_element], y_sampled[ids_in_element] +end + # returns points (xf, yf), scaled normals (nxJ, nyJ), and face Jacobian (Jf) # for a curve returned from PathIntersections. # `out` should hold `xf, yf, nxJ, nyJ, Jf = ntuple(_ -> similar(points, (length(points), num_faces)), 5)` From 76a4a554fafea91a32e9dc6dc485c3faa22c28f5 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 28 Apr 2024 20:29:48 -0500 Subject: [PATCH 15/39] add dispatch to preserve old version of MeshData --- src/cut_cell_meshes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index fadb648c..842d82dd 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -63,7 +63,7 @@ function neighbor_across_face(f, ex, ey) end end -function compute_face_centroids(rd, xf, yf, cutcell_data) +function compute_face_centroids(rd::RefElemData, xf, yf, cutcell_data) (; region_flags, cut_faces_per_cell, cut_face_offsets ) = cutcell_data num_cut_cells = length(cut_faces_per_cell) @@ -425,7 +425,7 @@ Inputs: - `region_flags`, `cutcells` are return arguments from `PathIntersections.define_regions` The keyword argument `tol` is the tolerance for matches between face centroids. """ -function connect_mesh(rd, face_centroids, cutcell_data; tol = 1e2 * eps()) +function connect_mesh(rd::RefElemData, face_centroids, cutcell_data; tol = 1e2 * eps()) (; region_flags, cut_faces_per_cell, cut_face_offsets ) = cutcell_data From 98dfaf7376c52d1474ce3d9254cbd04b0bd3d8a1 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 28 Apr 2024 20:29:54 -0500 Subject: [PATCH 16/39] remove cruft --- src/cut_cell_meshes.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 842d82dd..3e87d14a 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -495,8 +495,6 @@ function connect_mesh(rd::RefElemData, face_centroids, cutcell_data; tol = 1e2 * xy_nbr = SVector(face_centroids_x[j], face_centroids_y[j]) if norm(xy - xy_nbr) < tol * max(1, norm(xy), norm(xy_nbr)) FToF[i] = j - # println("match found for f = $f, e=($ex, $ey), - # enbr=($ex_nbr, $ey_nbr)") end end end From dc9c4f8637bab6dcb257e29e25146594bcd55521 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 28 Apr 2024 20:30:32 -0500 Subject: [PATCH 17/39] add new routines to create a cut cell MeshData with positive weights --- src/cut_cell_meshes.jl | 739 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 739 insertions(+) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 3e87d14a..75ea7c5f 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -546,6 +546,88 @@ function calculate_cutcells(vx, vy, objects, ds = 1e-3, arc_tol = 1e-10, corner_ return region_flags, cutcells end +""" + subtriangulated_cutcell_quadrature(cutcell, rd_tri::RefElemData, + r1D = gauss_lobatto_quad(0,0,rd_tri.N)) + +Constructs a quadrature from subtriangulations. The degree of both the quadrature +rule and isoparametric mapping are specified by `rd_tri`. + +The optional argument `r1D` specifies the nodal points used for constructing a +curved mapping via interpolatory warp and blend. +""" +function subtriangulated_cutcell_quadrature(cutcell, rd_tri, + r1D = NodesAndModes.nodes(Line(), rd_tri.N)) + # vxy = matrix of the form [x_coordinates, y_coordinates] + vxy = hcat(getindex.(cutcell.(cutcell.stop_pts[1:end-1]), 1), + getindex.(cutcell.(cutcell.stop_pts[1:end-1]), 2)) + (VX, VY), EToV = StartUpDG.triangulate_points(permutedims(vxy)) + + # allocate arrays for face interp coordinates + tri_face_coords = (rd_tri.r[vec(rd_tri.Fmask)], rd_tri.s[vec(rd_tri.Fmask)]) + tri_face_coords_x = zeros(length(r1D), 3) + tri_face_coords_y = zeros(length(r1D), 3) + + # allocate output arrays + xq, yq, wJq = ntuple(_ -> zeros(length(rd_tri.wq), size(EToV, 1)), 3) + + # 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}.(rd_tri.fv)) + vertices_on_face = sort(ids[face_vertices]) + + # map face interpolation points to a physical element. + + # This assumes PathIntersections.jl uses a clockwise ordering of stop curve points. + # Since StartUpDG uses a CCW ordering, we reverse the order for + for i in eachindex(r1D) + # 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], reverse(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], reverse(cutcell.stop_pts[curve_id:curve_id+1])...) + point = cutcell(s) + + else # it's a non-boundary face, it's a straight line + 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 operator performs a least squares fit, and is equivalent to isoparametric + # warp and blend. can be moved outside of the cutcell loop for efficiency. + warp_face_points_to_interp = + face_basis(Tri(), rd_tri.N, rd_tri.rst...) / + face_basis(Tri(), rd_tri.N, tri_face_coords...) + + # this performs a least squares fit interpolation by the face basis. It's + # equivalent to isoparametric warp and blend if the face node locations are + # representable by the face basis (e.g., polynomial). + 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 = + StartUpDG.geometric_factors(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 + @. wJq[:,e] = rd_tri.wq * Jq_e + end + + return xq, yq, wJq +end + + """ function MeshData(rd, geometry, vxyz...) @@ -823,3 +905,660 @@ function num_elements(md::MeshData{DIM, <:CutCellMesh}) where {DIM} # the number of elements is given by the number of columns of each component of x return num_cartesian_elements(md) + num_cut_elements(md) end + +""" + construct_cut_volume_quadrature(N, cutcells; target_degree = 2 * N) + +Constructs volume quadrature using subtriangulations of cut cells and +Caratheodory pruning. The resulting quadrature is exact for polynomials +of degree `target_degree`. +""" +function construct_cut_volume_quadrature(N, cutcells, physical_frame_elements; + target_degree = 2 * N) + + # 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 = max(2 * N, (N-1) * (N + 2)) + + rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, + quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) + + Np_target = StartUpDG.Np_cut(target_degree) + xq_pruned, yq_pruned, wJq_pruned = ntuple(_ -> zeros(Np_target, length(cutcells)), 3) + + for (e, cutcell) in enumerate(cutcells) + + xq, yq, wJq = StartUpDG.subtriangulated_cutcell_quadrature(cutcell, rd_tri) + + # if there are not enough quadrature nodes, refine the target polynomial + # degree until there are. + N_refined = N_phys_frame_geo + while length(wJq) < Np_target + N_refined += 1 + rd_tri_refined = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, + quad_rule_vol=NodesAndModes.quad_nodes_tri(N_refined)) + xq, yq, wJq = StartUpDG.subtriangulated_cutcell_quadrature(cutcell, rd_tri_refined) + end + + # perform Caratheodory pruning + Vtarget = vandermonde(physical_frame_elements[e], target_degree, vec(xq), vec(yq)) + w = vec(wJq) + w_pruned, inds = StartUpDG.caratheodory_pruning_qr(Vtarget, w) + + if target_degree >= 2*N + # test exactness of the pruned quadrature rule + V = vandermonde(physical_frame_elements[e], N, vec(xq), vec(yq)) + @assert norm(V' * diagm(w) * V - V' * diagm(w_pruned) * V) < 100 * eps() + # @show cond(V' * diagm(w_pruned) * V) + end + + # @show e, size(xq[inds]), size(xq_pruned) + @. xq_pruned[:, e] = xq[inds] + @. yq_pruned[:, e] = yq[inds] + @. wJq_pruned[:, e] = w_pruned[inds] + end + return xq_pruned, yq_pruned, wJq_pruned +end + +# construct interpolation points using sampling and approximate Fekete +# points (e.g., QR-DEIM). +function construct_cut_interpolation_nodes(N, objects, physical_frame_elements) + + N_sampled = 4 * N + num_cut_elements = length(physical_frame_elements) + x, y = ntuple(_ -> zeros(Np_cut(N), num_cut_elements), 2) + + # Compute interpolation points on cut elements + for e in eachindex(physical_frame_elements) + + physical_frame_element = physical_frame_elements[e] + + x_sampled, y_sampled = + generate_sampling_points(objects, physical_frame_element, 2 * Np_cut(N), N_sampled) + V = vandermonde(physical_frame_element, N, x_sampled, y_sampled) + + # use pivoted QR to find good interpolation points + QRfac = qr(V', ColumnNorm()) + ids = QRfac.p[1:Np_cut(N)] + + # if the condition number of the VDM is really bad, then increase the + # number of sampled points. + condV_original = cond(V[ids, :]) + condV = condV_original + N_sampled = 10 * N + iter, maxiters = 0, 100 + while condV > 1e8 && iter < maxiters + x_sampled, y_sampled = + generate_sampling_points(objects, physical_frame_element, + 2 * Np_cut(N), N_sampled) + V = vandermonde(physical_frame_element, N, x_sampled, y_sampled) + + # use pivoted QR to find good interpolation points + QRfac = qr(V', ColumnNorm()) + ids = QRfac.p[1:Np_cut(N)] + condV = cond(V[ids, :]) + N_sampled += 5 * N + + if condV < 1e8 + @warn "Conditioning of old VDM for element $e is $condV_original. " * + "After recomputing with $(length(x_sampled)) samples: $condV." + end + + iter += 1 + end + if iter >= maxiters + @warn "Adaptive sampling of cut-cell element $e: maximum number of iterations $maxiters exceeded." + end + + view(x, :, e) .= x_sampled[ids] + view(y, :, e) .= y_sampled[ids] + end + return x, y +end + +function construct_cartesian_volume_quadrature(vx, vy, region_flags, quad_rule_volume) + + rq, sq, wq = quad_rule_volume + + num_cartesian_cells = sum(region_flags .== 0) + xq, yq, wJq = ntuple(_ -> zeros(length(wq), num_cartesian_cells), 3) + + # compute quadrature rules for the Cartesian cells + e = 1 + for ex in axes(region_flags, 1), ey in axes(region_flags, 2) + if is_Cartesian(region_flags[ex, ey]) + dx = vx[ex+1] - vx[ex] + dy = vy[ey+1] - vy[ey] + J = dx * dy / sum(wq) + @. xq[:, e] = dx * 0.5 * (1 + rq) + vx[ex] + @. yq[:, e] = dy * 0.5 * (1 + sq) + vy[ey] + @. wJq[:, e] = wq * J + e += 1 + end + end + return xq, yq, wJq +end + +function construct_cartesian_surface_quadrature(vx, vy, region_flags, + quad_rule_face) + + LX, LY = (x -> x[2] - x[1]).(extrema.((vx, vy))) + cells_per_dimension_x, cells_per_dimension_y = size(region_flags) + r1D, w1D = quad_rule_face + + num_cartesian_cells = sum(region_flags .== 0) + Jf = zeros(num_faces(Quad()), num_cartesian_cells) + xf, yf, nxJ, nyJ, wJf = ntuple(_ -> zeros(num_faces(Quad()) * length(r1D), + num_cartesian_cells), 5) + + # the face Jacobian involves scaling between mapped and reference domain + # this is precomputed here since it's needed to compute the normals + + # face 1 and 2: ±s faces + face_ids = [eachindex(r1D); eachindex(r1D) .+ length(r1D)] + Jf[SVector(1, 2), :] .= (LY / (cells_per_dimension_y * sum(w1D))) + + # face 3 and 4: ±r faces + face_ids = [eachindex(r1D) .+ 2 * length(r1D); eachindex(r1D) .+ 3 * length(r1D)] + Jf[SVector(3, 4), :] .= (LX / (cells_per_dimension_x * sum(w1D))) + + # compute face data + e = 1 + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + if is_Cartesian(region_flags[ex, ey]) + # face 1: -r + face_index = 1 + face_ids = eachindex(r1D) .+ (face_index - 1) * length(r1D) + xf[face_ids, e] .= vx[ex] + yf[face_ids, e] .= map_to_interval(r1D, vy[ey], vy[ey+1]) + nxJ[face_ids, e] .= -Jf[face_index, e] + nyJ[face_ids, e] .= zero(eltype(nyJ)) + wJf[face_ids, e] .= w1D * Jf[face_index, e] + + # face 2: +r + face_index = 2 + face_ids = eachindex(r1D) .+ (face_index - 1) * length(r1D) + xf[face_ids, e] .= vx[ex+1] + yf[face_ids, e] .= map_to_interval(r1D, vy[ey], vy[ey+1]) + nxJ[face_ids, e] .= Jf[face_index, e] + nyJ[face_ids, e] .= zero(eltype(nyJ)) + wJf[face_ids, e] .= w1D * Jf[face_index, e] + + # face 3: -s + face_index = 3 + face_ids = eachindex(r1D) .+ (face_index - 1) * length(r1D) + xf[face_ids, e] .= map_to_interval(r1D, vx[ex], vx[ex+1]) + yf[face_ids, e] .= vy[ey] + nxJ[face_ids, e] .= zero(eltype(nxJ)) + nyJ[face_ids, e] .= -Jf[face_index, e] + wJf[face_ids, e] .= w1D * Jf[face_index, e] + + # face 3: +s + face_index = 4 + face_ids = eachindex(r1D) .+ (face_index - 1) * length(r1D) + xf[face_ids, e] .= map_to_interval(r1D, vx[ex], vx[ex+1]) + yf[face_ids, e] .= vy[ey+1] + nxJ[face_ids, e] .= zero(eltype(nxJ)) + nyJ[face_ids, e] .= Jf[face_index, e] + wJf[face_ids, e] .= w1D * Jf[face_index, e] + + e = e + 1 + end + end + + return xf, yf, nxJ, nyJ, wJf +end + +""" + construct_physical_frame_elements(region_flags, cutcells) + +Computes physical frame shifting and scaling parameters from the +vertices of cut cells and the background cell location. +""" +function construct_physical_frame_elements(region_flags, vx, vy, cutcells) + + physical_frame_elements = PhysicalFrame{2}[] # populate this as we iterate through cut cells + e = 1 + for ex in axes(region_flags, 1), ey in axes(region_flags, 2) + if StartUpDG.is_cut(region_flags[ex, ey]) + + # get extremal points (vertices) of the cut cell + cutcell = cutcells[e] + coordinates = cutcell.(cutcell.stop_pts[1:end-1]) + + # store vertex nodes and coordinates of background Cartesian cell + physical_frame_element = + PhysicalFrame(getindex.(coordinates, 1), getindex.(coordinates, 2), + SVector(vx[ex], vx[ex+1]), SVector(vy[ey], vy[ey+1])) + + push!(physical_frame_elements, physical_frame_element) + + e += 1 + end + end + return physical_frame_elements +end + +""" + construct_cut_surface_quadrature(N, cutcells, quad_rule_1D = gauss_quad(0, 0, N)) + +Constructs cut surface quadrature using a degree `N` geometric mapping and a reference +quadrature rule `quad_rule_1D`. Returns `xf, yf, nxJ, nyJ, wf` which are vectors, and +`face_node_indices`, which is a `Vector{Vector{Int}}` to index into each face. + +On boundaries of cut cells, the surface quadrature is taken to be exact for degree +`N^2 + (N-1)` polynomials. This ensures satisfaction of a weak GSBP property. +""" +function construct_cut_surface_quadrature(N, cutcells, quad_rule_1D = gauss_quad(0, 0, N)) + + rd_line = RefElemData(Line(), N, quad_rule_vol = quad_rule_1D) + + # we want to exactly integrate degree N^2 + (N-1) polynomials since our target + # integrand is nxJ * u (nxJ is degree N-1 on the reference face, and u is degree + # N on the physical face, but mapped back to the reference face it's N^2). + N_boundary = ceil(Int, (N^2 + (N-1) - 1) / 2) + rq1D, wq1D = gauss_quad(0, 0, N_boundary) + interp_to_cut_boundary = vandermonde(Line(), rd_line.N, rq1D) / rd_line.VDM + + xf, yf, nxJ, nyJ, wf = ntuple(_ -> Vector{Float64}[], 5) + + face_node_indices = Vector{UnitRange{Int64}}[] + + # recompute normals using isoparametric mapping on cut cells + for cutcell in cutcells + xf_element, yf_element, nxJ_element, nyJ_element, wf_element = + ntuple(_ -> Vector{Float64}[], 5) + + for f in 1:length(cutcell.stop_pts)-1 + + is_boundary_face = + !(cutcell.subcurves[f] isa PathIntersections.PresetGeometries.Line) + + # switch between lower order and higher order face quadrature + # if it's the curved boundary, it may be higher order + interp_to_face_quad_pts = (is_boundary_face) ? interp_to_cut_boundary : rd_line.Vq + + # map 1D interpolation points to curved faces + points = map(s -> cutcell(StartUpDG.map_to_interval(s, cutcell.stop_pts[f:f+1]...)), rd_line.r) + xf_face = getindex.(points, 1) + yf_face = getindex.(points, 2) + + # interp face coordinates to face quad nodes + xfq_face = interp_to_face_quad_pts * xf_face + yfq_face = interp_to_face_quad_pts * yf_face + + # compute tangent vector at quadrature points using + # the polynomial geometric mapping. + dxdr = interp_to_face_quad_pts * rd_line.Dr * xf_face + dydr = interp_to_face_quad_pts * rd_line.Dr * yf_face + scaled_normal = SVector.(-dydr, dxdr) + nxJ_face = @. getindex(scaled_normal, 1) + nyJ_face = @. getindex(scaled_normal, 2) + wf_face = (is_boundary_face) ? wq1D : rd_line.wq + + push!(xf_element, xfq_face) + push!(yf_element, yfq_face) + push!(nxJ_element, nxJ_face) + push!(nyJ_element, nyJ_face) + push!(wf_element, wf_face) + end + + # compute indices of nodes on each face before flattening + num_nodes_per_face = length.(xf_element) + face_node_indices_on_this_element = + map((length, offset) -> (1:length) .+ offset, + num_nodes_per_face, [0; cumsum(num_nodes_per_face[1:end-1])]) + push!(face_node_indices, face_node_indices_on_this_element) + + # flatten the collection of quadrature points for each cutcell + # and push them to the arrays xf, yf, ... + map((x,y) -> push!(x, vcat(y...)), + (xf, yf, nxJ, nyJ, wf), + (xf_element, yf_element, nxJ_element, nyJ_element, wf_element)) + end + return xf, yf, nxJ, nyJ, wf, face_node_indices +end + +# this assumes Cartesian cells are Quads +function compute_face_centroids(vx, vy, region_flags, + cutcells::AbstractVector{<:PathIntersections.PiecewiseCurve}) + + num_cartesian_cells = sum(region_flags .== 0) + num_cut_faces = sum(map(cutcell->length(cutcell.subcurves), cutcells)) + + face_centroids_x = NamedArrayPartition(cartesian=zeros(num_faces(Quad()), num_cartesian_cells), + cut=zeros(num_cut_faces)) + face_centroids_y = similar(face_centroids_x) + + # calculate Cartesian centroids + e = 1 + for ex in axes(region_flags, 1), ey in axes(region_flags, 2) + if is_Cartesian(region_flags[ex, ey]) + # quad faces are -r, +r, -s, +s + + # face 1 + face_centroids_x.cartesian[1, e] = vx[ex] + face_centroids_y.cartesian[1, e] = 0.5 * (vy[ey+1] + vy[ey]) + + # face 2 + face_centroids_x.cartesian[2, e] = vx[ex+1] + face_centroids_y.cartesian[2, e] = 0.5 * (vy[ey+1] + vy[ey]) + + # face 3 + face_centroids_x.cartesian[3, e] = 0.5 * (vx[ex] + vx[ex+1]) + face_centroids_y.cartesian[3, e] = vy[ey] + + # face 4 + face_centroids_x.cartesian[4, e] = 0.5 * (vx[ex] + vx[ex+1]) + face_centroids_y.cartesian[4, e] = vy[ey+1] + e += 1 + end + end + + face_index = 1 + for cutcell in cutcells + for (f, face_parametrization) in enumerate(cutcell.subcurves) + s_midpoint = 0.5 * sum(cutcell.sub_bounds[f]) + x, y = face_parametrization(s_midpoint) + face_centroids_x.cut[face_index] = x + face_centroids_y.cut[face_index] = y + face_index += 1 + end + end + + return face_centroids_x, face_centroids_y +end + +num_faces(cutcell::PathIntersections.PiecewiseCurve) = length(cutcell.subcurves) + +function connect_mesh(face_centroids, region_flags, + cutcells::Vector{<:PathIntersections.PiecewiseCurve}; tol = 1e2 * eps()) + + num_cartesian_cells = sum(region_flags .== 0) + cut_faces_per_cell = num_faces.(cutcells) + cut_face_offsets = [0; cumsum(cut_faces_per_cell)[1:end-1]] + num_cut_faces = sum(cut_faces_per_cell) + num_total_faces = num_cut_faces + num_faces(Quad()) * num_cartesian_cells + + # element_indices[ex, ey] returns the global (flattened) element index into + # the arrays `xf.cartesian[:, e]` or `xf.cut[:, e]` + element_indices = compute_element_indices(region_flags) + + # compute face centroids for making face matches + face_centroids_x, face_centroids_y = face_centroids + + # To determine face-to-face matches, we work with each background Cartesian element + # and search through the 4 neighboring background Cartesian elements for a match in + # the face centroids of the current cell and the face centroids of its neighbors. + # NOTE: this works because the cut cells can only have Cartesian neighbors across + # flat-sided faces. It wouldn't work for meshes with curved interior interfaces. + + FToF = collect(1:num_total_faces) + for ex in axes(region_flags, 1), ey in axes(region_flags, 2) + + e = element_indices[ex, ey] + + # Determine face indices of current cell. The face indices are determined + # from the flattened (e.g., not ex, ey) element ordering. + # NOTE: we search for matches between all faces of `e` and all faces of + # `e_nbr` because the ordering of faces is different for cut elements + # and Cartesian elements. + if is_Cartesian(region_flags[ex, ey]) + face_ids = (1:num_faces(Quad())) .+ (e-1) * num_faces(Quad()) + elseif is_cut(region_flags[ex, ey]) + face_ids = (1:cut_faces_per_cell[e]) .+ cut_face_offsets[e] + + # we offset by the number of Cartesian faces so we can index globally + # into the arrays `face_centroids_x`, `face_centroids_y`. + face_ids = face_ids .+ length(face_centroids_x.cartesian) + end + + if is_inside_domain(ex, ey, region_flags) + + for f in 1:4 # each Cartesian background element has 4 neighbors + + ex_nbr, ey_nbr = neighbor_across_face(f, ex, ey) + if is_inside_domain(ex_nbr, ey_nbr, region_flags) + e_nbr = element_indices[ex_nbr, ey_nbr] + + # determine face indices of neighboring cells + if is_Cartesian(region_flags[ex_nbr, ey_nbr]) + nbr_face_ids = (1:num_faces(Quad())) .+ (e_nbr-1) * num_faces(Quad()) + elseif is_cut(region_flags[ex_nbr, ey_nbr]) + nbr_face_ids = (1:cut_faces_per_cell[e_nbr]) .+ cut_face_offsets[e_nbr] + + # we offset by the number of Cartesian faces so we can index globally + # into the arrays `face_centroids_x`, `face_centroids_y`. + nbr_face_ids = nbr_face_ids .+ length(face_centroids_x.cartesian) + end + + # check for matches in face and neighbor face centroids. + # note: we index into the global `face_centroids_x/y` container + # rather than the `.cut` or `.cartesian subarrays`. + for i in face_ids + xy = SVector(face_centroids_x[i], face_centroids_y[i]) + for j in nbr_face_ids + xy_nbr = SVector(face_centroids_x[j], face_centroids_y[j]) + if norm(xy - xy_nbr) < tol * max(1, norm(xy), norm(xy_nbr)) + FToF[i] = j + end + end + end + + end # if enbr is_inside_domain + end + end # if e is_inside_domain + end + + return FToF +end + + + +function MeshData2(rd::RefElemData, objects, + vx::AbstractVector, vy::AbstractVector; + quad_rule_face=get_1d_quadrature(rd), + precompute_operators=false) + + cells_per_dimension_x = length(vx) - 1 + cells_per_dimension_y = length(vy) - 1 + + LX, LY = (x -> x[2] - x[1]).(extrema.((vx, vy))) + + # compute mesh intersections and cut cell elements. + # `regions` is a matrix of dimensions `(cells_per_dimension_x, cells_per_dimension_y)` with 3 values: + # * 1: cut cell + # * 0: Cartesian cell + # * -1: excluded cells (in the cut-out region) + region_flags, cutcells = calculate_cutcells(vx, vy, objects) + + # pack useful cut cell information together. + num_cartesian_cells = sum(region_flags .== 0) + num_cut_cells = sum(region_flags .> 0) + cut_faces_per_cell = count_cut_faces(cutcells) + cut_face_offsets = [0; cumsum(cut_faces_per_cell)[1:end-1]] + cutcell_data = (; objects, region_flags, cutcells, cut_faces_per_cell, cut_face_offsets) + + N = rd.N + + #################################################### + # Construct cut cells stuff # + #################################################### + + physical_frame_elements = + construct_physical_frame_elements(region_flags, vx, vy, cutcells) + + x_cut, y_cut = + construct_cut_interpolation_nodes(N, objects, physical_frame_elements) + + xq_cut, yq_cut, wJq_cut = + construct_cut_volume_quadrature(N, cutcells, physical_frame_elements) + + xf_cut, yf_cut, nxJ_cut, nyJ_cut, wJf_cut, cut_face_node_indices = + construct_cut_surface_quadrature(N, cutcells, quad_rule_face) + + #################################################### + # Construct Cartesian cells # + #################################################### + + xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wJf_cartesian = + construct_cartesian_surface_quadrature(vx, vy, region_flags, quad_rule_face) + + xq_cartesian, yq_cartesian, wJq_cartesian = + construct_cartesian_volume_quadrature(vx, vy, region_flags, + (rd.rq, rd.sq, rd.wq)) + + # reuse quadrature mapping to construct interpolation nodes + x_cartesian, y_cartesian, _ = + construct_cartesian_volume_quadrature(vx, vy, region_flags, + (rd.r, rd.s, ones(size(rd.r)))) + + #################################################### + # Assemble combined cut/Cartesian arrays # + #################################################### + + x, y = map((x, y) -> NamedArrayPartition((; cartesian=x, cut=y)), + (x_cartesian, y_cartesian), (x_cut, y_cut)) + + xq, yq, wJq = map((x, y) -> NamedArrayPartition((; cartesian=x, cut=y)), + (xq_cartesian, yq_cartesian, wJq_cartesian), + (xq_cut, yq_cut, wJq_cut)) + + xf, yf, nxJ, nyJ, wJf = + map((x, y) -> NamedArrayPartition((; cartesian=x, cut=y)), + (xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wJf_cartesian), + map(x -> vcat(x...), (xf_cut, yf_cut, nxJ_cut, nyJ_cut, wJf_cut))) + + #################################################### + # Mesh connectivity # + #################################################### + + face_centroids = + StartUpDG.compute_face_centroids(vx, vy, region_flags, cutcells) + + FToF = StartUpDG.connect_mesh(face_centroids, region_flags, cutcells) + + # Compute node-to-node mappings + # !!! Warning: this only works if the same quadrature rule is used + # !!! for both cut/cartesian faces! + num_total_faces = length(FToF) + num_points_per_face = length(rd.rf) ÷ num_faces(rd.element_type) + mapM = collect(reshape(1:num_points_per_face * num_total_faces, + num_points_per_face, num_total_faces)) + mapP = copy(mapM) + p = zeros(Int, num_points_per_face) # temp storage for a permutation vector + for f in eachindex(FToF) + idM, idP = view(mapM, :, f), view(mapM, :, FToF[f]) + xyzM = (view(xf, idM), view(yf, idM)) + xyzP = (view(xf, idP), view(yf, idP)) + StartUpDG.match_coordinate_vectors!(p, xyzM, xyzP) + mapP[p, f] .= idP + end + mapB = findall(vec(mapM) .== vec(mapP)) # determine boundary nodes + + # default to non-periodic + is_periodic = (false, false) + + # compute mapping from linear element indices to Cartesian element indices + cartesian_to_linear_element_indices = compute_element_indices(region_flags) + linear_to_cartesian_element_indices = (; cut=zeros(SVector{2, Int}, num_cut_cells), + cartesian=zeros(SVector{2, Int}, num_cartesian_cells)) + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + e = cartesian_to_linear_element_indices[ex, ey] + if is_cut(region_flags[ex, ey]) + linear_to_cartesian_element_indices.cut[e] = SVector(ex, ey) + elseif is_Cartesian(region_flags[ex, ey]) + linear_to_cartesian_element_indices.cartesian[e] = SVector(ex, ey) + end + end + + cut_cell_data = (; cutcells, region_flags, wJf, + cartesian_to_linear_element_indices, + linear_to_cartesian_element_indices, + vxyz=(vx, vy), + cells_per_dimension=(cells_per_dimension_x, + cells_per_dimension_y), + ) # background Cartesian grid info + + # get flattened indices of cut face nodes + face_ids(e) = (1:(num_points_per_face * cut_faces_per_cell[e])) .+ + cut_face_offsets[e] * num_points_per_face + cut_face_node_ids = [face_ids(e) for e in 1:num_cut_cells] + + if precompute_operators == true + + # precompute cut-cell operators and store them in the `md.mesh_type.cut_cell_operators` field + cut_face_nodes = cut_face_node_ids + face_interpolation_matrices = Matrix{eltype(x)}[] + lift_matrices = Matrix{eltype(x)}[] + differentiation_matrices = Tuple{Matrix{eltype(x)}, Matrix{eltype(x)}}[] + mass_matrices = Matrix{eltype(x)}[] + for (e, elem) in enumerate(physical_frame_elements) + + VDM = vandermonde(elem, rd.N, x.cut[:, e], y.cut[:, e]) + Vq, Vrq, Vsq = map(A -> A / VDM, + basis(elem, rd.N, xq.cut[:,e], yq.cut[:, e])) + + M = Vq' * diagm(wJq.cut[:, e]) * Vq + Qr = Vq' * diagm(wJq.cut[:, e]) * Vrq + Qs = Vq' * diagm(wJq.cut[:, e]) * Vsq + Dx_e, Dy_e = M \ Qr, M \ Qs + + Vf = vandermonde(elem, rd.N, xf.cut[cut_face_nodes[e]], + yf.cut[cut_face_nodes[e]]) / VDM + + # don't include jacobian scaling in LIFT matrix (for consistency with the Cartesian mesh) + _, w1D = quad_rule_face + num_cut_faces = length(cut_face_nodes[e]) ÷ length(w1D) + wf = vec(repeat(w1D, 1, num_cut_faces)) + + push!(lift_matrices, M \ (Vf' * diagm(wf))) + push!(face_interpolation_matrices, Vf) + push!(differentiation_matrices, (Dx_e, Dy_e)) + push!(mass_matrices, M) + end + cut_cell_operators = (; differentiation_matrices, face_interpolation_matrices, + mass_matrices, lift_matrices) + + else + + cut_cell_operators = nothing + end + + + ################################################ + # Construct volume geometric terms # + ################################################ + + rxJ_cartesian = LX / (2 * cells_per_dimension_x) + syJ_cartesian = LY / (2 * cells_per_dimension_y) + J_cartesian = (LX / cells_per_dimension_x) * (LY / cells_per_dimension_y) / 4 + + # Note: the volume Jacobian for cut elements is 1 since the "reference element" is the + # cut element itself. Similarly, geometric terms should be 1 since `basis` computes + # physical derivatives accounting for element scaling + + # Note: we use FillArrays.Fill instead of FillArrays.Ones and FillArrays.Zeros because + # we store `rxJ, sxJ, ryJ, syJ` in a single SMatrix, which assumes one homogeneous + # type for all the entries. Since FillArrays.Ones/Zeros are distinct types, using them + # results in a type instability when accessing entries of md.rstxyzJ + rxJ = NamedArrayPartition(cartesian=Fill(rxJ_cartesian, rd.Np, num_cartesian_cells), + cut=Fill(1.0, Np_cut(rd.N), num_cut_cells)) + syJ = NamedArrayPartition(cartesian=Fill(syJ_cartesian, rd.Np, num_cartesian_cells), + cut=Fill(1.0, Np_cut(rd.N), num_cut_cells)) + sxJ, ryJ = ntuple(_ -> NamedArrayPartition(cartesian=Fill(0.0, rd.Np, num_cartesian_cells), + cut=Fill(0.0, Np_cut(rd.N), num_cut_cells)), 2) + J = NamedArrayPartition(cartesian = Fill(J_cartesian, rd.Np, num_cartesian_cells), + cut = Fill(1.0, Np_cut(rd.N), num_cut_cells)) + + # pack geometric terms together + rstxyzJ = SMatrix{2, 2, typeof(rxJ), 4}(rxJ, sxJ, ryJ, syJ) + + Jf = @. sqrt(nxJ^2 + nyJ^2) + + return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_ids, + objects, cut_cell_operators, cut_cell_data), + FToF, (x, y), (xf, yf), (xq, yq), wJq, + mapM, mapP, mapB, rstxyzJ, J, (nxJ, nyJ), Jf, is_periodic) + +end \ No newline at end of file From 2428cf8d69547379994fe6661c2b85d47ab1b9b7 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 28 Apr 2024 20:30:35 -0500 Subject: [PATCH 18/39] refactoring --- test_cutcell_sbp.jl | 236 ++++++++++++-------------------------------- 1 file changed, 62 insertions(+), 174 deletions(-) diff --git a/test_cutcell_sbp.jl b/test_cutcell_sbp.jl index 5e528d42..fa125aac 100644 --- a/test_cutcell_sbp.jl +++ b/test_cutcell_sbp.jl @@ -1,209 +1,97 @@ using NodesAndModes: face_basis using Plots using LinearAlgebra +using StaticArrays using StartUpDG +using PathIntersections N = 3 -quad_rule_face = gauss_quad(0, 0, 2 * N + 1) - +quad_rule_face = gauss_lobatto_quad(0, 0, N) rd = RefElemData(Quad(), N; quad_rule_face) -# objects = (PresetGeometries.Rectangle(Lx=.4, Ly=2.5, x0=1.2), ) -# vx = [-1, 0, 1.1] -# vy = [-1, 0, 1] -# md = MeshData(rd, objects, vx, vy; precompute_operators=true) - -cells_per_dimension = 2 +cells_per_dimension = 4 circle = PresetGeometries.Circle(R=0.66, x0=0, y0=0) -md = MeshData(rd, (circle, ), cells_per_dimension; quad_rule_face)#, precompute_operators=true) +objects = (circle, ) -mt = md.mesh_type -(; cutcells) = mt.cut_cell_data +# md = MeshData(rd, objects, cells_per_dimension; precompute_operators=true) -using Triangulate, StaticArrays +# TODO: remove eventually +cells_per_dimension_x = cells_per_dimension +cells_per_dimension_y = cells_per_dimension +vx = LinRange(-1, 1, cells_per_dimension_x + 1) +vy = LinRange(-1, 1, cells_per_dimension_y + 1) +N = rd.N -# degree N(N-1) + 2N-2 polynomial = N(N-1) + 2(N-1) = (N+2) * (N-1) -N_phys_frame_geo = max(2 * rd.N, (rd.N-1) * (rd.N + 2)) +#################################################### +# Construct Cartesian cells # +#################################################### -rd_line = RefElemData(Line(), N=rd.N, quad_rule_vol = quad_rule_face) +xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wJf_cartesian = + StartUpDG.construct_cartesian_surface_quadrature(vx, vy, region_flags, quad_rule_face) -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)]) +xq_cartesian, yq_cartesian, wJq_cartesian = + StartUpDG.construct_cartesian_volume_quadrature(vx, vy, region_flags, + (rd.rq, rd.sq, rd.wq)) -# 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) +#################################################### +# Construct cut cells stuff # +#################################################### -# 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...) +region_flags, cutcells = StartUpDG.calculate_cutcells(vx, vy, objects) -using StartUpDG: map_to_interval +physical_frame_elements = + StartUpDG.construct_physical_frame_elements(region_flags, vx, vy, cutcells) -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 +xf_cut, yf_cut, nxJ_cut, nyJ_cut, wJf_cut, cut_face_node_indices = + StartUpDG.construct_cut_surface_quadrature(N, cutcells, quad_rule_face) -x_cutcells, y_cutcells, xq_cutcells, yq_cutcells, Jq_cutcells, wJq_cutcells = - ntuple(_ -> Matrix{eltype(rd_tri.wq)}[], 6) - - -for cutcell in cutcells - - # vxy = matrix of the form [x_coordinates, y_coordinates] - vxy = hcat(getindex.(cutcell.(cutcell.stop_pts[1:end-1]), 1), - getindex.(cutcell.(cutcell.stop_pts[1:end-1]), 2)) - (VX, VY), EToV = StartUpDG.triangulate_points(permutedims(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 face interpolation points to a physical element. - - # This assumes PathIntersections.jl uses a clockwise ordering of stop curve points. - # Since StartUpDG uses a CCW ordering, we reverse the order for - for i in eachindex(r1D) - # 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], reverse(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], reverse(cutcell.stop_pts[curve_id:curve_id+1])...) - point = cutcell(s) - - else # it's a non-boundary face, it's a straight line - 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 by the face basis. It's - # equivalent to isoparametric warp and blend if the face node locations are - # representable by the face basis (e.g., polynomial). - 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) - - rxJq_e, sxJq_e, ryJq_e, syJq_e, Jq_e = - StartUpDG.geometric_factors(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 +xq_pruned, yq_pruned, wJq_pruned = + StartUpDG.construct_cut_volume_quadrature(N, cutcells, physical_frame_elements) -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 = StartUpDG.caratheodory_pruning_qr(V2N, w) +#################################################### +# Mesh connectivity # +#################################################### - 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) +face_centroids = + StartUpDG.compute_face_centroids(vx, vy, region_flags, cutcells) - 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 +FToF = StartUpDG.connect_mesh(face_centroids, region_flags, cutcells) -# recompute normals on cut cells -cutcell = cutcells[1] -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) - xf_local = getindex.(points, 1) - yf_local = getindex.(points, 2) +####################################### +# test weak SBP property # +####################################### - # compute tangent vector using polynomial mapping - dxdr = rd_line.Vq * rd_line.Dr * xf_local - dydr = rd_line.Vq * rd_line.Dr * yf_local +elem = physical_frame_elements[1] +xq, yq, wq = xq_pruned[:,1], yq_pruned[:,1], wJq_pruned[:,1] - tangent_vector = SVector.(dxdr, dydr) - scaling = (cutcell.stop_pts[f+1] - cutcell.stop_pts[f]) / 2 - scaled_normal = SVector.(-dydr, dxdr) +Vq, Vxq, Vyq = basis(elem, rd.N, xq, yq) +Vf = vandermonde(elem, rd.N, xf_cut[1], yf_cut[1]) +Qx, Qy = Vq' * diagm(wq) * Vxq, Vq' * diagm(wq) * Vyq - @. nxJ[:, f] = getindex(scaled_normal, 1) - @. nyJ[:, f] = getindex(scaled_normal, 2) +Bx = Diagonal(wJf_cut[1] .* nxJ_cut[1]) +By = Diagonal(wJf_cut[1] .* nyJ_cut[1]) - # interp face coordinates to face quad nodes - xf[:, f] .= rd_line.Vq * xf_local - yf[:, f] .= rd_line.Vq * yf_local +# represent constant vector in basis +e = Vq \ ones(size(Vq, 1)) +@show norm(e' * Qx - e' * Vf' * Bx * Vf) +@show norm(e' * Qy - e' * Vf' * By * Vf) -end +# plot pruned cells -# plot normals -plot() -for f in 1:length(cutcell.stop_pts)-1 - scatter!(xf[:,f], yf[:,f]) - quiver!(xf[:,f], yf[:,f], quiver=(nxJ[:,f], nyJ[:,f])) -end -plot!(leg=false, ratio=1) +# 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 = max(2 * N, (N-1) * (N + 2)) +target_degree = 2 * N +rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, + quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) +Np_target = StartUpDG.Np_cut(target_degree) -# plot pruned cells plot() -for e in eachindex(xq_cutcells) - scatter!(vec(xq_cutcells[e]), vec(yq_cutcells[e]), label="Reference quadrature"); - scatter!(xq_pruned[e], yq_pruned[e], markersize=8, marker=:circle, +for e in eachindex(cutcells) + xq, yq, _ = StartUpDG.subtriangulated_cutcell_quadrature(cutcells[e], rd_tri) + scatter!(vec(xq), vec(yq), 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)) - - - -# # compute new normals for curved boundaries -# mapB_boundary_cartesian = findall(@. abs(abs(md.xf) - 1) < 1e2 * eps() || abs(abs(md.yf) - 1) < 1e4 * eps()) -# mapB_cut = setdiff(md.mapB, mapB_boundary_cartesian) -# num_face_nodes = length(quad_rule_face[1]) -# xf, yf = map(x -> reshape(x[mapB_cut], num_face_nodes, :), md.xyzf) -# D1D = rd_line.Vq * rd_line.Dr * rd_line.Pq -# dxdr, dydr = D1D * xf, D1D * yf -# nxJ_boundary_cut, nyJ_boundary_cut = dydr, -dxdr - -# 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]) -# xq, yq, wq = xq_pruned[1], yq_pruned[1], wJq_pruned[1] -xq, yq, wq = vec.((xq_cutcells[1], yq_cutcells[1], wJq_cutcells[1])) -Vq, Vxq, Vyq = map(A -> A / VDM, basis(elem, rd.N, xq, yq)) -M = Vq' * diagm(wq) * Vq -Qx, Qy = Vq' * diagm(wq) * Vxq, Vq' * diagm(wq) * Vyq -Vf = vandermonde(elem, rd.N, xf, yf) / VDM - -# Bx = Diagonal(vec(Diagonal(rd_line.wq) * reshape(md.nxJ.cut[md.mesh_type.cut_face_nodes[1]], length(rd_line.wq), :))) -# By = Diagonal(vec(Diagonal(rd_line.wq) * reshape(md.nyJ.cut[md.mesh_type.cut_face_nodes[1]], length(rd_line.wq), :))) -Bx = Diagonal(vec(Diagonal(rd_line.wq) * reshape(nxJ, length(rd_line.wq), :))) -By = Diagonal(vec(Diagonal(rd_line.wq) * reshape(nyJ, length(rd_line.wq), :))) - -e = ones(size(Qx, 2)) -@show norm(sum(Qx - Vf' * Bx * Vf, dims=1)) -@show norm(sum(Qy - Vf' * By * Vf, dims=1)) -# display(sum(Qx - Vf' * Bx * Vf, dims=1)) -# display(e' * ((Qx + Qx') - Vf' * Bx * Vf)) -# display(e' * ((Qy + Qy') - Vf' * By * Vf)) From 238347247842bf826b477855b03ad89185818c66 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 28 Apr 2024 21:32:31 -0500 Subject: [PATCH 19/39] fix construction of cut cell face node indices previously assumed same number of nodes on all faces --- src/cut_cell_meshes.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 75ea7c5f..f4f26753 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -1480,10 +1480,13 @@ function MeshData2(rd::RefElemData, objects, cells_per_dimension_y), ) # background Cartesian grid info - # get flattened indices of cut face nodes - face_ids(e) = (1:(num_points_per_face * cut_faces_per_cell[e])) .+ - cut_face_offsets[e] * num_points_per_face - cut_face_node_ids = [face_ids(e) for e in 1:num_cut_cells] + # get flattened indices of cut face nodes. + # note that `cut_face_node_indices = [[face_1_indices, face_2_indices, ...], ...]` + flattened_face_node_indices = + map(x -> UnitRange(x[1][1], x[end][end]), cut_face_node_indices) + face_index_offsets = [0; cumsum(length.(flattened_face_node_indices)[1:end-1])] + cut_face_node_ids = + map((x, y) -> x .+ y, flattened_face_node_indices, face_index_offsets) if precompute_operators == true From aff16046cf3ef36990f6381cf31451d43f7442c3 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 28 Apr 2024 21:32:42 -0500 Subject: [PATCH 20/39] fix precomputation of operators --- src/cut_cell_meshes.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index f4f26753..fc975a4c 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -1429,6 +1429,8 @@ function MeshData2(rd::RefElemData, objects, (xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wJf_cartesian), map(x -> vcat(x...), (xf_cut, yf_cut, nxJ_cut, nyJ_cut, wJf_cut))) + Jf = @. sqrt(nxJ^2 + nyJ^2) + #################################################### # Mesh connectivity # #################################################### @@ -1510,10 +1512,9 @@ function MeshData2(rd::RefElemData, objects, Vf = vandermonde(elem, rd.N, xf.cut[cut_face_nodes[e]], yf.cut[cut_face_nodes[e]]) / VDM - # don't include jacobian scaling in LIFT matrix (for consistency with the Cartesian mesh) - _, w1D = quad_rule_face - num_cut_faces = length(cut_face_nodes[e]) ÷ length(w1D) - wf = vec(repeat(w1D, 1, num_cut_faces)) + # don't include jacobian scaling in LIFT matrix (for consistency + # with the Cartesian mesh) + wf = wJf.cut[cut_face_nodes[e]] ./ Jf.cut[cut_face_nodes[e]] push!(lift_matrices, M \ (Vf' * diagm(wf))) push!(face_interpolation_matrices, Vf) @@ -1557,8 +1558,6 @@ function MeshData2(rd::RefElemData, objects, # pack geometric terms together rstxyzJ = SMatrix{2, 2, typeof(rxJ), 4}(rxJ, sxJ, ryJ, syJ) - Jf = @. sqrt(nxJ^2 + nyJ^2) - return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_ids, objects, cut_cell_operators, cut_cell_data), FToF, (x, y), (xf, yf), (xq, yq), wJq, From b2b63538cee986cca32a73b445a2b96b4fcc597d Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 28 Apr 2024 21:33:01 -0500 Subject: [PATCH 21/39] clean up test of SBP property --- test_cutcell_sbp.jl | 64 ++++++++++++++++++++------------------------- 1 file changed, 28 insertions(+), 36 deletions(-) diff --git a/test_cutcell_sbp.jl b/test_cutcell_sbp.jl index fa125aac..020ec253 100644 --- a/test_cutcell_sbp.jl +++ b/test_cutcell_sbp.jl @@ -49,49 +49,41 @@ xf_cut, yf_cut, nxJ_cut, nyJ_cut, wJf_cut, cut_face_node_indices = xq_pruned, yq_pruned, wJq_pruned = StartUpDG.construct_cut_volume_quadrature(N, cutcells, physical_frame_elements) -#################################################### -# Mesh connectivity # -#################################################### - -face_centroids = - StartUpDG.compute_face_centroids(vx, vy, region_flags, cutcells) - -FToF = StartUpDG.connect_mesh(face_centroids, region_flags, cutcells) - ####################################### # test weak SBP property # ####################################### -elem = physical_frame_elements[1] -xq, yq, wq = xq_pruned[:,1], yq_pruned[:,1], wJq_pruned[:,1] +e = 5 +elem = physical_frame_elements[e] +xq, yq, wq = xq_pruned[:,e], yq_pruned[:,e], wJq_pruned[:,e] Vq, Vxq, Vyq = basis(elem, rd.N, xq, yq) -Vf = vandermonde(elem, rd.N, xf_cut[1], yf_cut[1]) +Vf = vandermonde(elem, rd.N, xf_cut[e], yf_cut[e]) Qx, Qy = Vq' * diagm(wq) * Vxq, Vq' * diagm(wq) * Vyq -Bx = Diagonal(wJf_cut[1] .* nxJ_cut[1]) -By = Diagonal(wJf_cut[1] .* nyJ_cut[1]) +Bx = Diagonal(wJf_cut[e] .* nxJ_cut[e]) +By = Diagonal(wJf_cut[e] .* nyJ_cut[e]) # represent constant vector in basis -e = Vq \ ones(size(Vq, 1)) -@show norm(e' * Qx - e' * Vf' * Bx * Vf) -@show norm(e' * Qy - e' * Vf' * By * Vf) - -# plot pruned cells - -# 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 = max(2 * N, (N-1) * (N + 2)) -target_degree = 2 * N -rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, - quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) -Np_target = StartUpDG.Np_cut(target_degree) - -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!(xq_pruned[:,e], yq_pruned[:,e], markersize=8, marker=:circle, - z_order=:back, label="Caratheodory pruning", leg=false) -end -display(plot!(leg=false)) +ee = Vq \ ones(size(Vq, 1)) +@show norm(ee' * Qx - ee' * Vf' * Bx * Vf) +@show norm(ee' * Qy - ee' * Vf' * By * Vf) + +# # plot pruned cells + +# # 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 = max(2 * N, (N-1) * (N + 2)) +# target_degree = 2 * N +# rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, +# quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) +# Np_target = StartUpDG.Np_cut(target_degree) + +# 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!(xq_pruned[:,e], yq_pruned[:,e], markersize=8, marker=:circle, +# z_order=:back, label="Caratheodory pruning", leg=false) +# end +# display(plot!(leg=false)) From d6e331152e0f23e4db8a82d6291d5aa060c69154 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Sun, 28 Apr 2024 21:33:15 -0500 Subject: [PATCH 22/39] add test of SBP property using new MeshData --- test_cutcell_MeshData.jl | 64 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 test_cutcell_MeshData.jl diff --git a/test_cutcell_MeshData.jl b/test_cutcell_MeshData.jl new file mode 100644 index 00000000..971ba3a0 --- /dev/null +++ b/test_cutcell_MeshData.jl @@ -0,0 +1,64 @@ +using NodesAndModes: face_basis +using Plots +using LinearAlgebra +using StaticArrays +using StartUpDG +using PathIntersections + +N = 3 +quad_rule_face = gauss_lobatto_quad(0, 0, N) +rd = RefElemData(Quad(), N; quad_rule_face) + +cells_per_dimension = 4 +circle = PresetGeometries.Circle(R=0.66, x0=0, y0=0) +objects = (circle, ) + +# TODO: remove eventually +cells_per_dimension_x = cells_per_dimension +cells_per_dimension_y = cells_per_dimension +vx = LinRange(-1, 1, cells_per_dimension_x + 1) +vy = LinRange(-1, 1, cells_per_dimension_y + 1) + +####################################### +# test weak SBP property # +####################################### + +md = StartUpDG.MeshData2(rd, objects, vx, vy; precompute_operators=true) +mt = md.mesh_type + +for (e, elem) in enumerate(mt.physical_frame_elements) + xq, yq, wq = md.xq.cut[:,e], md.yq.cut[:,e], md.wJq.cut[:,e] + + Vq, Vxq, Vyq = basis(elem, rd.N, xq, yq) + face_ids = mt.cut_face_nodes[e] + Vf = vandermonde(elem, rd.N, md.xf.cut[face_ids], md.yf.cut[face_ids]) + Qx, Qy = Vq' * diagm(wq) * Vxq, Vq' * diagm(wq) * Vyq + + (; wJf) = md.mesh_type.cut_cell_data + Bx = Diagonal(wJf.cut[face_ids] .* md.nxJ.cut[face_ids]) + By = Diagonal(wJf.cut[face_ids] .* md.nyJ.cut[face_ids]) + + # represent constant vector in basis + ee = Vq \ ones(size(Vq, 1)) + @show norm(ee' * Qx - ee' * Vf' * Bx * Vf) + @show norm(ee' * Qy - ee' * Vf' * By * Vf) +end + +# plot pruned cells + +# 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 = max(2 * N, (N-1) * (N + 2)) +target_degree = 2 * N +rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, + quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) +Np_target = StartUpDG.Np_cut(target_degree) + +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!(xq_pruned[:,e], yq_pruned[:,e], markersize=8, marker=:circle, + z_order=:back, label="Caratheodory pruning", leg=false) +end +display(plot!(leg=false)) From a8d2becd5c742ffbcbcc7f74daf6b745c3230e8f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 29 Apr 2024 17:20:54 -0500 Subject: [PATCH 23/39] add MomentFitting dispatch for old cut cell MeshData d --- src/StartUpDG.jl | 1 + src/cut_cell_meshes.jl | 31 ++++++++++++++++--------------- test/cut_mesh_tests.jl | 12 +++++++----- 3 files changed, 24 insertions(+), 20 deletions(-) diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 11acfb57..f85c028a 100644 --- a/src/StartUpDG.jl +++ b/src/StartUpDG.jl @@ -63,6 +63,7 @@ export num_faces, num_vertices, HybridMeshExample include("physical_frame_basis.jl") include("cut_cell_meshes.jl") export PhysicalFrame, equi_nodes +export MomentFitting include("state_redistribution.jl") export StateRedistribution, apply! diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index fc975a4c..b4c900fe 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -627,23 +627,24 @@ function subtriangulated_cutcell_quadrature(cutcell, rd_tri, return xq, yq, wJq end +struct MomentFitting end """ - function MeshData(rd, geometry, vxyz...) - -Creates a cut-cell mesh where the boundary is given by `geometry`, which should be a tuple of functions. -These functions can be generated using PathIntersections.PresetGeometries, for example: -```julia -julia> geometry = (PresetGeometries.Circle(R=0.33, x0=0, y0=0), ) -``` -Here, `coordinates_min`, `coordinates_max` contain `(smallest value of x, smallest value of y)` and -`(largest value of x, largest value of y)`, and `cells_per_dimension_x/y` is the number of Cartesian grid -cells placed along each dimension. + MeshData(quadrature_type::MomentFitting, rd::RefElemData, + objects, vx::AbstractVector, vy::AbstractVector; + quad_rule_face=get_1d_quadrature(rd), + precompute_operators=false) + +Constructor for MeshData utilizing moment fitting. Does not guarantee positive +quadrature weights, and is slower due to the use of adaptive sampling +to construct + +!!! Warning: this may be deprecated or removed in future versions. """ -MeshData(rd::RefElemData, objects, cells_per_dimension; kwargs...) = - MeshData(rd::RefElemData, objects, cells_per_dimension, cells_per_dimension; kwargs...) +MeshData(quadrature_type::MomentFitting, rd::RefElemData, objects, cells_per_dimension; kwargs...) = + MeshData(quadrature_type, rd::RefElemData, objects, cells_per_dimension, cells_per_dimension; kwargs...) -function MeshData(rd::RefElemData, objects, +function MeshData(quadrature_type::MomentFitting, rd::RefElemData, objects, cells_per_dimension_x::Int, cells_per_dimension_y::Int; quad_rule_face = get_1d_quadrature(rd), coordinates_min=(-1.0, -1.0), coordinates_max=(1.0, 1.0), @@ -653,10 +654,10 @@ function MeshData(rd::RefElemData, objects, vx = LinRange(coordinates_min[1], coordinates_max[1], cells_per_dimension_x + 1) vy = LinRange(coordinates_min[2], coordinates_max[2], cells_per_dimension_y + 1) - return MeshData(rd, objects, vx, vy; quad_rule_face, precompute_operators) + return MeshData(quadrature_type, rd, objects, vx, vy; quad_rule_face, precompute_operators) end -function MeshData(rd::RefElemData, objects, +function MeshData(::MomentFitting, rd::RefElemData, objects, vx::AbstractVector, vy::AbstractVector; quad_rule_face=get_1d_quadrature(rd), precompute_operators=false) diff --git a/test/cut_mesh_tests.jl b/test/cut_mesh_tests.jl index 92ccb15a..07922814 100644 --- a/test/cut_mesh_tests.jl +++ b/test/cut_mesh_tests.jl @@ -7,7 +7,9 @@ using StartUpDG: PathIntersections circle = PresetGeometries.Circle(R=0.33, x0=0, y0=0) rd = RefElemData(Quad(), N=3) - md = MeshData(rd, (circle, ), cells_per_dimension_x, cells_per_dimension_y; precompute_operators=true) + md = MeshData(MomentFitting(), rd, (circle, ), + cells_per_dimension_x, cells_per_dimension_y; + precompute_operators=true) @test_throws ErrorException("Face index f = 5 > 4; too large.") StartUpDG.neighbor_across_face(5, nothing, nothing) @@ -16,15 +18,14 @@ using StartUpDG: PathIntersections @test (@capture_out Base.show(stdout, MIME"text/plain"(), md)) == "Cut-cell MeshData of dimension 2 with 16 elements (12 Cartesian, 4 cut)" # test constructor with only one "cells_per_dimension" argument - @test_nowarn MeshData(rd, (circle, ), cells_per_dimension_x) + @test_nowarn MeshData(MomentFitting(), rd, (circle, ), cells_per_dimension_x) # check the volume of the domain A = 4 - pi * .33^2 @test sum(md.wJq) ≈ A # check the length of the boundary of the domain - face_weights = reshape(rd.wf, :, num_faces(rd.element_type))[:, 1] - wJf = vec(Diagonal(face_weights) * reshape(md.Jf, length(face_weights), :)) + wJf = md.mesh_type.cut_cell_data.wJf @test sum(wJf[md.mapB]) ≈ (8 + 2 * pi * .33) # check continuity of a function that's in the global polynomial space @@ -81,7 +82,8 @@ end cells_per_dimension = 4 circle = PresetGeometries.Circle(R=0.6, x0=0, y0=0) rd = RefElemData(Quad(), N=4) - md = MeshData(rd, (circle, ), cells_per_dimension, cells_per_dimension) + md = MeshData(MomentFitting(), rd, (circle, ), + cells_per_dimension, cells_per_dimension) srd = StateRedistribution(rd, md) e = @. 0 * md.x + 1 # constant From 8711b899b9d8c668c67401a1dccb4ca5bbd970ee Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 15:06:41 -0500 Subject: [PATCH 24/39] remove outdated todo --- src/cut_cell_meshes.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index b4c900fe..99da68da 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -28,8 +28,6 @@ struct CutCellMesh{T1, T2, T3, T4, T5} cut_cell_data::T5 end -# TODO: add isoparametric cut cell mesh with positive quadrature points - function Base.show(io::IO, ::MIME"text/plain", md::MeshData{DIM, <:CutCellMesh}) where {DIM} @nospecialize md print(io,"Cut-cell MeshData of dimension $DIM with $(md.num_elements) elements " * From 6585529e57ef8c575a015a54f53a22a2d2fc80d9 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 15:07:43 -0500 Subject: [PATCH 25/39] changing wJf to wf internally --- src/cut_cell_meshes.jl | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 99da68da..753f654e 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -1047,8 +1047,8 @@ function construct_cartesian_surface_quadrature(vx, vy, region_flags, num_cartesian_cells = sum(region_flags .== 0) Jf = zeros(num_faces(Quad()), num_cartesian_cells) - xf, yf, nxJ, nyJ, wJf = ntuple(_ -> zeros(num_faces(Quad()) * length(r1D), - num_cartesian_cells), 5) + xf, yf, nxJ, nyJ, wf = ntuple(_ -> zeros(num_faces(Quad()) * length(r1D), + num_cartesian_cells), 5) # the face Jacobian involves scaling between mapped and reference domain # this is precomputed here since it's needed to compute the normals @@ -1072,7 +1072,7 @@ function construct_cartesian_surface_quadrature(vx, vy, region_flags, yf[face_ids, e] .= map_to_interval(r1D, vy[ey], vy[ey+1]) nxJ[face_ids, e] .= -Jf[face_index, e] nyJ[face_ids, e] .= zero(eltype(nyJ)) - wJf[face_ids, e] .= w1D * Jf[face_index, e] + wf[face_ids, e] .= w1D # face 2: +r face_index = 2 @@ -1081,7 +1081,7 @@ function construct_cartesian_surface_quadrature(vx, vy, region_flags, yf[face_ids, e] .= map_to_interval(r1D, vy[ey], vy[ey+1]) nxJ[face_ids, e] .= Jf[face_index, e] nyJ[face_ids, e] .= zero(eltype(nyJ)) - wJf[face_ids, e] .= w1D * Jf[face_index, e] + wf[face_ids, e] .= w1D # face 3: -s face_index = 3 @@ -1090,7 +1090,7 @@ function construct_cartesian_surface_quadrature(vx, vy, region_flags, yf[face_ids, e] .= vy[ey] nxJ[face_ids, e] .= zero(eltype(nxJ)) nyJ[face_ids, e] .= -Jf[face_index, e] - wJf[face_ids, e] .= w1D * Jf[face_index, e] + wf[face_ids, e] .= w1D # face 3: +s face_index = 4 @@ -1099,13 +1099,13 @@ function construct_cartesian_surface_quadrature(vx, vy, region_flags, yf[face_ids, e] .= vy[ey+1] nxJ[face_ids, e] .= zero(eltype(nxJ)) nyJ[face_ids, e] .= Jf[face_index, e] - wJf[face_ids, e] .= w1D * Jf[face_index, e] + wf[face_ids, e] .= w1D e = e + 1 end end - return xf, yf, nxJ, nyJ, wJf + return xf, yf, nxJ, nyJ, wf end """ @@ -1393,15 +1393,15 @@ function MeshData2(rd::RefElemData, objects, xq_cut, yq_cut, wJq_cut = construct_cut_volume_quadrature(N, cutcells, physical_frame_elements) - xf_cut, yf_cut, nxJ_cut, nyJ_cut, wJf_cut, cut_face_node_indices = - construct_cut_surface_quadrature(N, cutcells, quad_rule_face) + xf_cut, yf_cut, nxJ_cut, nyJ_cut, wf_cut, cut_face_node_indices = + StartUpDG.construct_cut_surface_quadrature(N, cutcells, quad_rule_face) #################################################### # Construct Cartesian cells # #################################################### - xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wJf_cartesian = - construct_cartesian_surface_quadrature(vx, vy, region_flags, quad_rule_face) + xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wf_cartesian = + StartUpDG.construct_cartesian_surface_quadrature(vx, vy, region_flags, quad_rule_face) xq_cartesian, yq_cartesian, wJq_cartesian = construct_cartesian_volume_quadrature(vx, vy, region_flags, @@ -1423,12 +1423,13 @@ function MeshData2(rd::RefElemData, objects, (xq_cartesian, yq_cartesian, wJq_cartesian), (xq_cut, yq_cut, wJq_cut)) - xf, yf, nxJ, nyJ, wJf = + xf, yf, nxJ, nyJ, wf = map((x, y) -> NamedArrayPartition((; cartesian=x, cut=y)), - (xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wJf_cartesian), - map(x -> vcat(x...), (xf_cut, yf_cut, nxJ_cut, nyJ_cut, wJf_cut))) + (xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wf_cartesian), + map(x -> vcat(x...), (xf_cut, yf_cut, nxJ_cut, nyJ_cut, wf_cut))) Jf = @. sqrt(nxJ^2 + nyJ^2) + wJf = @. wf * Jf #################################################### # Mesh connectivity # From 60a2a7cc5ddff52d255637e2e4df3f9248d92cea Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 15:08:21 -0500 Subject: [PATCH 26/39] add face node index array --- src/cut_cell_meshes.jl | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 753f654e..6fce5215 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -1143,7 +1143,8 @@ end Constructs cut surface quadrature using a degree `N` geometric mapping and a reference quadrature rule `quad_rule_1D`. Returns `xf, yf, nxJ, nyJ, wf` which are vectors, and -`face_node_indices`, which is a `Vector{Vector{Int}}` to index into each face. +`face_node_indices`, which is a `Vector{Vector{Int}}` of global face node indices (which +index into `xf.cut`, `yf.cut`, etc) for each face of each cut element. On boundaries of cut cells, the surface quadrature is taken to be exact for degree `N^2 + (N-1)` polynomials. This ensures satisfaction of a weak GSBP property. @@ -1215,6 +1216,16 @@ function construct_cut_surface_quadrature(N, cutcells, quad_rule_1D = gauss_quad (xf, yf, nxJ, nyJ, wf), (xf_element, yf_element, nxJ_element, nyJ_element, wf_element)) end + + # compute global cut face node indices + num_nodes_per_cut_element = map(x -> maximum(maximum.(x)), face_node_indices) + face_node_offsets = [0; cumsum(num_nodes_per_cut_element[1:end-1])] + for e in eachindex(face_node_indices) + for f in eachindex(face_node_indices[e]) + face_node_indices[e][f] = face_node_indices[e][f] .+ face_node_offsets[e] + end + end + return xf, yf, nxJ, nyJ, wf, face_node_indices end From bde44b01331297e6bb8307087f3b6777e648c162 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 15:09:28 -0500 Subject: [PATCH 27/39] make face centroid computation more compact --- src/cut_cell_meshes.jl | 40 +++++++++++++--------------------------- 1 file changed, 13 insertions(+), 27 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 6fce5215..496232bf 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -1233,50 +1233,36 @@ end function compute_face_centroids(vx, vy, region_flags, cutcells::AbstractVector{<:PathIntersections.PiecewiseCurve}) - num_cartesian_cells = sum(region_flags .== 0) - num_cut_faces = sum(map(cutcell->length(cutcell.subcurves), cutcells)) + num_cartesian_cells = sum(region_flags .== 0) + face_centroids_cartesian = + zeros(SVector{2, Float64}, num_faces(Quad()), num_cartesian_cells) - face_centroids_x = NamedArrayPartition(cartesian=zeros(num_faces(Quad()), num_cartesian_cells), - cut=zeros(num_cut_faces)) - face_centroids_y = similar(face_centroids_x) + num_faces_per_cut_cell = map(cutcell->length(cutcell.subcurves), cutcells) + face_centroids_cut = [zeros(SVector{2, Float64}, num_faces_per_cut_cell[e]) + for e in eachindex(cutcells)] # calculate Cartesian centroids e = 1 for ex in axes(region_flags, 1), ey in axes(region_flags, 2) if is_Cartesian(region_flags[ex, ey]) # quad faces are -r, +r, -s, +s - - # face 1 - face_centroids_x.cartesian[1, e] = vx[ex] - face_centroids_y.cartesian[1, e] = 0.5 * (vy[ey+1] + vy[ey]) - - # face 2 - face_centroids_x.cartesian[2, e] = vx[ex+1] - face_centroids_y.cartesian[2, e] = 0.5 * (vy[ey+1] + vy[ey]) - - # face 3 - face_centroids_x.cartesian[3, e] = 0.5 * (vx[ex] + vx[ex+1]) - face_centroids_y.cartesian[3, e] = vy[ey] - - # face 4 - face_centroids_x.cartesian[4, e] = 0.5 * (vx[ex] + vx[ex+1]) - face_centroids_y.cartesian[4, e] = vy[ey+1] + face_centroids_cartesian[1, e] = SVector(vx[ex], 0.5 * (vy[ey+1] + vy[ey])) + face_centroids_cartesian[2, e] = SVector(vx[ex+1], 0.5 * (vy[ey+1] + vy[ey])) + face_centroids_cartesian[3, e] = SVector(0.5 * (vx[ex] + vx[ex+1]), vy[ey]) + face_centroids_cartesian[4, e] = SVector(0.5 * (vx[ex] + vx[ex+1]), vy[ey+1]) e += 1 end end - face_index = 1 - for cutcell in cutcells + for (e, cutcell) in enumerate(cutcells) for (f, face_parametrization) in enumerate(cutcell.subcurves) s_midpoint = 0.5 * sum(cutcell.sub_bounds[f]) x, y = face_parametrization(s_midpoint) - face_centroids_x.cut[face_index] = x - face_centroids_y.cut[face_index] = y - face_index += 1 + face_centroids_cut[e][f] = SVector(x, y) end end - return face_centroids_x, face_centroids_y + return face_centroids_cartesian, face_centroids_cut end num_faces(cutcell::PathIntersections.PiecewiseCurve) = length(cutcell.subcurves) From 14fc1ee61333d20f6c4f72713aa99923c05af557 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 15:09:55 -0500 Subject: [PATCH 28/39] specialize connect_mesh --- src/cut_cell_meshes.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 496232bf..2baa0b8b 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -1268,7 +1268,8 @@ end num_faces(cutcell::PathIntersections.PiecewiseCurve) = length(cutcell.subcurves) function connect_mesh(face_centroids, region_flags, - cutcells::Vector{<:PathIntersections.PiecewiseCurve}; tol = 1e2 * eps()) + cutcells::Vector{<:PathIntersections.PiecewiseCurve}; + tol = 1e2 * eps()) num_cartesian_cells = sum(region_flags .== 0) cut_faces_per_cell = num_faces.(cutcells) From 69c83aa0f298e897bff4e02eb1c28acff51ed78c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 15:10:21 -0500 Subject: [PATCH 29/39] add new MeshData based on subtriangulations --- src/cut_cell_meshes.jl | 127 +++++++++++++++++++++++++++-------------- 1 file changed, 85 insertions(+), 42 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 2baa0b8b..11ee6fdd 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -1350,9 +1350,35 @@ function connect_mesh(face_centroids, region_flags, return FToF end +""" + function MeshData(rd, geometry, vxyz...) + +Creates a cut-cell mesh where the boundary is given by `geometry`, which should be a tuple of functions. +These functions can be generated using PathIntersections.PresetGeometries, for example: +```julia +julia> geometry = (PresetGeometries.Circle(R=0.33, x0=0, y0=0), ) +``` +Here, `coordinates_min`, `coordinates_max` contain `(smallest value of x, smallest value of y)` and +`(largest value of x, largest value of y)`, and `cells_per_dimension_x/y` is the number of Cartesian grid +cells placed along each dimension. +""" +MeshData(rd::RefElemData, objects, cells_per_dimension; kwargs...) = + MeshData(rd::RefElemData, objects, cells_per_dimension, cells_per_dimension; kwargs...) + +function MeshData(rd::RefElemData, objects, + cells_per_dimension_x::Int, cells_per_dimension_y::Int; + quad_rule_face = get_1d_quadrature(rd), + coordinates_min=(-1.0, -1.0), coordinates_max=(1.0, 1.0), + precompute_operators=false) + # compute intersections of curve with a background Cartesian grid. + vx = LinRange(coordinates_min[1], coordinates_max[1], cells_per_dimension_x + 1) + vy = LinRange(coordinates_min[2], coordinates_max[2], cells_per_dimension_y + 1) -function MeshData2(rd::RefElemData, objects, + return MeshData(rd, objects, vx, vy; quad_rule_face, precompute_operators) +end + +function MeshData(rd::RefElemData, objects, vx::AbstractVector, vy::AbstractVector; quad_rule_face=get_1d_quadrature(rd), precompute_operators=false) @@ -1367,14 +1393,11 @@ function MeshData2(rd::RefElemData, objects, # * 1: cut cell # * 0: Cartesian cell # * -1: excluded cells (in the cut-out region) - region_flags, cutcells = calculate_cutcells(vx, vy, objects) + region_flags, cutcells = StartUpDG.calculate_cutcells(vx, vy, objects) # pack useful cut cell information together. num_cartesian_cells = sum(region_flags .== 0) num_cut_cells = sum(region_flags .> 0) - cut_faces_per_cell = count_cut_faces(cutcells) - cut_face_offsets = [0; cumsum(cut_faces_per_cell)[1:end-1]] - cutcell_data = (; objects, region_flags, cutcells, cut_faces_per_cell, cut_face_offsets) N = rd.N @@ -1383,13 +1406,13 @@ function MeshData2(rd::RefElemData, objects, #################################################### physical_frame_elements = - construct_physical_frame_elements(region_flags, vx, vy, cutcells) + StartUpDG.construct_physical_frame_elements(region_flags, vx, vy, cutcells) x_cut, y_cut = - construct_cut_interpolation_nodes(N, objects, physical_frame_elements) + StartUpDG.construct_cut_interpolation_nodes(N, objects, physical_frame_elements) xq_cut, yq_cut, wJq_cut = - construct_cut_volume_quadrature(N, cutcells, physical_frame_elements) + StartUpDG.construct_cut_volume_quadrature(N, cutcells, physical_frame_elements) xf_cut, yf_cut, nxJ_cut, nyJ_cut, wf_cut, cut_face_node_indices = StartUpDG.construct_cut_surface_quadrature(N, cutcells, quad_rule_face) @@ -1402,12 +1425,12 @@ function MeshData2(rd::RefElemData, objects, StartUpDG.construct_cartesian_surface_quadrature(vx, vy, region_flags, quad_rule_face) xq_cartesian, yq_cartesian, wJq_cartesian = - construct_cartesian_volume_quadrature(vx, vy, region_flags, + StartUpDG.construct_cartesian_volume_quadrature(vx, vy, region_flags, (rd.rq, rd.sq, rd.wq)) # reuse quadrature mapping to construct interpolation nodes x_cartesian, y_cartesian, _ = - construct_cartesian_volume_quadrature(vx, vy, region_flags, + StartUpDG.construct_cartesian_volume_quadrature(vx, vy, region_flags, (rd.r, rd.s, ones(size(rd.r)))) #################################################### @@ -1433,29 +1456,51 @@ function MeshData2(rd::RefElemData, objects, # Mesh connectivity # #################################################### - face_centroids = - StartUpDG.compute_face_centroids(vx, vy, region_flags, cutcells) - - FToF = StartUpDG.connect_mesh(face_centroids, region_flags, cutcells) - - # Compute node-to-node mappings - # !!! Warning: this only works if the same quadrature rule is used - # !!! for both cut/cartesian faces! - num_total_faces = length(FToF) - num_points_per_face = length(rd.rf) ÷ num_faces(rd.element_type) - mapM = collect(reshape(1:num_points_per_face * num_total_faces, - num_points_per_face, num_total_faces)) + # note: face_centroids_cartesian has dims [4, num_elements] + # face_centroids_cut has dims [num_elements][num_faces_on_this_element] + face_centroids_cartesian, face_centroids_cut = + StartUpDG.compute_face_centroids(vx, vy, region_flags, cutcells) + + face_centroids_x = NamedArrayPartition(cartesian = getindex.(face_centroids_cartesian, 1), + cut = getindex.(vcat(face_centroids_cut...), 1)) + face_centroids_y = NamedArrayPartition(cartesian = getindex.(face_centroids_cartesian, 2), + cut = getindex.(vcat(face_centroids_cut...), 2)) + face_centroids = (face_centroids_x, face_centroids_y) + FToF = StartUpDG.connect_mesh(face_centroids, region_flags, cutcells) + + # create arrays of node indices on each face, e.g., + # xf.cut[cut_face_node_indices_by_face[f]] gives nodes on face f + cartesian_nodes_per_face = rd.Nfq ÷ rd.num_faces + cartesian_face_node_indices_by_face = + [(1:cartesian_nodes_per_face) .+ (f - 1) * cartesian_nodes_per_face + for f in 1:num_cartesian_cells * rd.num_faces] + cut_face_node_indices_by_face = vcat(cut_face_node_indices...) + + # create global index vector (note that we offset cut cell indices + # by num_cartesian nodes). + face_node_indices_by_face = vcat(cartesian_face_node_indices_by_face, + map(x -> x .+ (num_cartesian_cells * rd.Nfq), + cut_face_node_indices_by_face)) + + # create node mappings + mapM_cartesian = reshape(eachindex(xf.cartesian), :, num_cartesian_cells) + mapM = NamedArrayPartition(cartesian = collect(mapM_cartesian), + cut = collect(eachindex(xf.cut) .+ length(mapM_cartesian))) mapP = copy(mapM) - p = zeros(Int, num_points_per_face) # temp storage for a permutation vector + for f in eachindex(FToF) - idM, idP = view(mapM, :, f), view(mapM, :, FToF[f]) - xyzM = (view(xf, idM), view(yf, idM)) - xyzP = (view(xf, idP), view(yf, idP)) - StartUpDG.match_coordinate_vectors!(p, xyzM, xyzP) - mapP[p, f] .= idP + # if it's not a boundary face + if f !== FToF[f] + idM = mapM[face_node_indices_by_face[f]] + idP = mapM[face_node_indices_by_face[FToF[f]]] + xyzM = (view(xf, idM), view(yf, idM)) + xyzP = (view(xf, idP), view(yf, idP)) + p = StartUpDG.match_coordinate_vectors(xyzM, xyzP) + mapP[idM[p]] .= idP + end end mapB = findall(vec(mapM) .== vec(mapP)) # determine boundary nodes - + # default to non-periodic is_periodic = (false, false) @@ -1480,18 +1525,15 @@ function MeshData2(rd::RefElemData, objects, cells_per_dimension_y), ) # background Cartesian grid info - # get flattened indices of cut face nodes. - # note that `cut_face_node_indices = [[face_1_indices, face_2_indices, ...], ...]` - flattened_face_node_indices = + # get flattened indices of cut face nodes. note that + # `cut_face_node_indices = [[face_1_indices, face_2_indices, ...], ...]` + cut_face_node_indices_by_cell = map(x -> UnitRange(x[1][1], x[end][end]), cut_face_node_indices) - face_index_offsets = [0; cumsum(length.(flattened_face_node_indices)[1:end-1])] - cut_face_node_ids = - map((x, y) -> x .+ y, flattened_face_node_indices, face_index_offsets) if precompute_operators == true - # precompute cut-cell operators and store them in the `md.mesh_type.cut_cell_operators` field - cut_face_nodes = cut_face_node_ids + # precompute cut-cell operators and store them in the `md.mesh_type.cut_cell_operators` field. + cut_face_nodes = cut_face_node_indices_by_cell face_interpolation_matrices = Matrix{eltype(x)}[] lift_matrices = Matrix{eltype(x)}[] differentiation_matrices = Tuple{Matrix{eltype(x)}, Matrix{eltype(x)}}[] @@ -1507,12 +1549,13 @@ function MeshData2(rd::RefElemData, objects, Qs = Vq' * diagm(wJq.cut[:, e]) * Vsq Dx_e, Dy_e = M \ Qr, M \ Qs - Vf = vandermonde(elem, rd.N, xf.cut[cut_face_nodes[e]], - yf.cut[cut_face_nodes[e]]) / VDM + Vf = vandermonde(elem, rd.N, xf.cut[cut_face_node_indices_by_cell[e]], + yf.cut[cut_face_node_indices_by_cell[e]]) / VDM # don't include jacobian scaling in LIFT matrix (for consistency # with the Cartesian mesh) - wf = wJf.cut[cut_face_nodes[e]] ./ Jf.cut[cut_face_nodes[e]] + wf = wJf.cut[cut_face_node_indices_by_cell[e]] ./ + Jf.cut[cut_face_node_indices_by_cell[e]] push!(lift_matrices, M \ (Vf' * diagm(wf))) push!(face_interpolation_matrices, Vf) @@ -1555,8 +1598,8 @@ function MeshData2(rd::RefElemData, objects, # pack geometric terms together rstxyzJ = SMatrix{2, 2, typeof(rxJ), 4}(rxJ, sxJ, ryJ, syJ) - - return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_ids, + + return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_indices_by_cell, objects, cut_cell_operators, cut_cell_data), FToF, (x, y), (xf, yf), (xq, yq), wJq, mapM, mapP, mapB, rstxyzJ, J, (nxJ, nyJ), Jf, is_periodic) From 88d8b9c97930b7c854459f522d1d85fee19a57b7 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 15:53:41 -0500 Subject: [PATCH 30/39] add Subtriangulation() as a quadrature type d --- src/StartUpDG.jl | 2 +- src/cut_cell_meshes.jl | 54 +++++++++++++++++++----------------------- 2 files changed, 25 insertions(+), 31 deletions(-) diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index f85c028a..ab146804 100644 --- a/src/StartUpDG.jl +++ b/src/StartUpDG.jl @@ -63,7 +63,7 @@ export num_faces, num_vertices, HybridMeshExample include("physical_frame_basis.jl") include("cut_cell_meshes.jl") export PhysicalFrame, equi_nodes -export MomentFitting +export Subtriangulation, MomentFitting include("state_redistribution.jl") export StateRedistribution, apply! diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 11ee6fdd..ae7056e9 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -626,10 +626,12 @@ function subtriangulated_cutcell_quadrature(cutcell, rd_tri, end struct MomentFitting end +struct Subtriangulation end """ - MeshData(quadrature_type::MomentFitting, rd::RefElemData, - objects, vx::AbstractVector, vy::AbstractVector; + MeshData(rd::RefElemData, objects, + vx::AbstractVector, vy::AbstractVector, + quadrature_type::MomentFitting; quad_rule_face=get_1d_quadrature(rd), precompute_operators=false) @@ -639,11 +641,14 @@ to construct !!! Warning: this may be deprecated or removed in future versions. """ -MeshData(quadrature_type::MomentFitting, rd::RefElemData, objects, cells_per_dimension; kwargs...) = - MeshData(quadrature_type, rd::RefElemData, objects, cells_per_dimension, cells_per_dimension; kwargs...) +MeshData(rd::RefElemData, objects, cells_per_dimension, + quadrature_type=Subtriangulation(); kwargs...) = + MeshData(rd, objects, cells_per_dimension, cells_per_dimension, + quadrature_type; kwargs...) -function MeshData(quadrature_type::MomentFitting, rd::RefElemData, objects, - cells_per_dimension_x::Int, cells_per_dimension_y::Int; +function MeshData(rd::RefElemData, objects, + cells_per_dimension_x::Int, cells_per_dimension_y::Int, + quadrature_type=Subtriangulation(); quad_rule_face = get_1d_quadrature(rd), coordinates_min=(-1.0, -1.0), coordinates_max=(1.0, 1.0), precompute_operators=false) @@ -652,11 +657,13 @@ function MeshData(quadrature_type::MomentFitting, rd::RefElemData, objects, vx = LinRange(coordinates_min[1], coordinates_max[1], cells_per_dimension_x + 1) vy = LinRange(coordinates_min[2], coordinates_max[2], cells_per_dimension_y + 1) - return MeshData(quadrature_type, rd, objects, vx, vy; quad_rule_face, precompute_operators) + return MeshData(rd, objects, vx, vy, quadrature_type; + quad_rule_face, precompute_operators) end -function MeshData(::MomentFitting, rd::RefElemData, objects, - vx::AbstractVector, vy::AbstractVector; +function MeshData(rd::RefElemData, objects, + vx::AbstractVector, vy::AbstractVector, + quadrature_type::MomentFitting; quad_rule_face=get_1d_quadrature(rd), precompute_operators=false) @@ -1353,33 +1360,20 @@ end """ function MeshData(rd, geometry, vxyz...) -Creates a cut-cell mesh where the boundary is given by `geometry`, which should be a tuple of functions. -These functions can be generated using PathIntersections.PresetGeometries, for example: +Creates a cut-cell mesh where the boundary is given by `geometry`, which should be a +tuple of functions. These functions can be generated using PathIntersections.PresetGeometries, +for example: ```julia julia> geometry = (PresetGeometries.Circle(R=0.33, x0=0, y0=0), ) ``` -Here, `coordinates_min`, `coordinates_max` contain `(smallest value of x, smallest value of y)` and -`(largest value of x, largest value of y)`, and `cells_per_dimension_x/y` is the number of Cartesian grid -cells placed along each dimension. +Here, `coordinates_min`, `coordinates_max` contain `(smallest value of x, smallest value of y)` +and `(largest value of x, largest value of y)`, and `cells_per_dimension_x/y` is the number of +Cartesian grid cells placed along each dimension. """ -MeshData(rd::RefElemData, objects, cells_per_dimension; kwargs...) = - MeshData(rd::RefElemData, objects, cells_per_dimension, cells_per_dimension; kwargs...) - -function MeshData(rd::RefElemData, objects, - cells_per_dimension_x::Int, cells_per_dimension_y::Int; - quad_rule_face = get_1d_quadrature(rd), - coordinates_min=(-1.0, -1.0), coordinates_max=(1.0, 1.0), - precompute_operators=false) - - # compute intersections of curve with a background Cartesian grid. - vx = LinRange(coordinates_min[1], coordinates_max[1], cells_per_dimension_x + 1) - vy = LinRange(coordinates_min[2], coordinates_max[2], cells_per_dimension_y + 1) - - return MeshData(rd, objects, vx, vy; quad_rule_face, precompute_operators) -end function MeshData(rd::RefElemData, objects, - vx::AbstractVector, vy::AbstractVector; + vx::AbstractVector, vy::AbstractVector, + quadrature_type::Subtriangulation; quad_rule_face=get_1d_quadrature(rd), precompute_operators=false) From e29ad8089f47b135413bef175e660b550f1f547e Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 15:53:59 -0500 Subject: [PATCH 31/39] add quadrature type to CutCellMesh meshtype --- src/cut_cell_meshes.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index ae7056e9..d1771d07 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -20,12 +20,13 @@ face interpolation, and lifting operators). The field `cut_cell_data` contains additional data from PathIntersections. """ -struct CutCellMesh{T1, T2, T3, T4, T5} +struct CutCellMesh{T1, T2, T3, T4, T5, T6} physical_frame_elements::T1 cut_face_nodes::T2 objects::T3 cut_cell_operators::T4 cut_cell_data::T5 + quadrature_type::T6 end function Base.show(io::IO, ::MIME"text/plain", md::MeshData{DIM, <:CutCellMesh}) where {DIM} @@ -897,7 +898,11 @@ function MeshData(rd::RefElemData, objects, cut_cell_operators = nothing end - return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_ids, objects, cut_cell_operators, cut_cell_data), + # added for consistency + J = NamedArrayPartition(cartesian=Fill(J, size(x.cartesian)), + cut=Fill(1.0, size(x.cut))) + return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_ids, objects, + cut_cell_operators, cut_cell_data, quadrature_type), FToF, (x, y), (xf, yf), (xq, yq), wJq, mapM, mapP, mapB, rstxyzJ, J, (nxJ, nyJ), Jf, is_periodic) @@ -1593,8 +1598,8 @@ function MeshData(rd::RefElemData, objects, # pack geometric terms together rstxyzJ = SMatrix{2, 2, typeof(rxJ), 4}(rxJ, sxJ, ryJ, syJ) - return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_indices_by_cell, - objects, cut_cell_operators, cut_cell_data), + return MeshData(CutCellMesh(physical_frame_elements, cut_face_node_indices_by_cell, objects, + cut_cell_operators, cut_cell_data, quadrature_type), FToF, (x, y), (xf, yf), (xq, yq), wJq, mapM, mapP, mapB, rstxyzJ, J, (nxJ, nyJ), Jf, is_periodic) From 8e714270eebe647cfa28dcc801e78a3f33468c1f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 15:54:13 -0500 Subject: [PATCH 32/39] test both Subtriangulation and MomentFitting --- test/cut_mesh_tests.jl | 192 +++++++++++++++++++++-------------------- 1 file changed, 97 insertions(+), 95 deletions(-) diff --git a/test/cut_mesh_tests.jl b/test/cut_mesh_tests.jl index 07922814..626a3301 100644 --- a/test/cut_mesh_tests.jl +++ b/test/cut_mesh_tests.jl @@ -1,100 +1,102 @@ using StartUpDG: PathIntersections -@testset "Cut meshes" begin - - cells_per_dimension = 4 - cells_per_dimension_x, cells_per_dimension_y = cells_per_dimension, cells_per_dimension - circle = PresetGeometries.Circle(R=0.33, x0=0, y0=0) - - rd = RefElemData(Quad(), N=3) - md = MeshData(MomentFitting(), rd, (circle, ), - cells_per_dimension_x, cells_per_dimension_y; - precompute_operators=true) - - @test_throws ErrorException("Face index f = 5 > 4; too large.") StartUpDG.neighbor_across_face(5, nothing, nothing) - - @test StartUpDG.num_cartesian_elements(md) + StartUpDG.num_cut_elements(md) == md.num_elements - - @test (@capture_out Base.show(stdout, MIME"text/plain"(), md)) == "Cut-cell MeshData of dimension 2 with 16 elements (12 Cartesian, 4 cut)" - - # test constructor with only one "cells_per_dimension" argument - @test_nowarn MeshData(MomentFitting(), rd, (circle, ), cells_per_dimension_x) - - # check the volume of the domain - A = 4 - pi * .33^2 - @test sum(md.wJq) ≈ A - - # check the length of the boundary of the domain - wJf = md.mesh_type.cut_cell_data.wJf - @test sum(wJf[md.mapB]) ≈ (8 + 2 * pi * .33) - - # check continuity of a function that's in the global polynomial space - (; physical_frame_elements) = md.mesh_type - (; x, y) = md - u = @. x^rd.N - x * y^(rd.N-1) - x^(rd.N-1) * y + y^rd.N - 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] +@testset "Cut meshes ($(typeof(quadrature_type)))" for quadrature_type = [Subtriangulation(), MomentFitting()] + @testset "Correctness" begin + + cells_per_dimension = 4 + cells_per_dimension_x, cells_per_dimension_y = cells_per_dimension, cells_per_dimension + circle = PresetGeometries.Circle(R=0.33, x0=0, y0=0) + + rd = RefElemData(Quad(), N=3) + md = MeshData(rd, (circle, ), + cells_per_dimension_x, cells_per_dimension_y, + quadrature_type; + precompute_operators=true) + + @test_throws ErrorException("Face index f = 5 > 4; too large.") StartUpDG.neighbor_across_face(5, nothing, nothing) + + @test StartUpDG.num_cartesian_elements(md) + StartUpDG.num_cut_elements(md) == md.num_elements + + @test (@capture_out Base.show(stdout, MIME"text/plain"(), md)) == "Cut-cell MeshData of dimension 2 with 16 elements (12 Cartesian, 4 cut)" + + # test constructor with only one "cells_per_dimension" argument + @test_nowarn MeshData(rd, (circle, ), cells_per_dimension_x, quadrature_type) + + # check the volume of the domain + A = 4 - pi * .33^2 + @test abs(sum(md.wJq) - A) < 7e-5 # should be 6.417e-5 + + # check the length of the boundary of the domain + wJf = md.mesh_type.cut_cell_data.wJf + @test abs(sum(wJf[md.mapB]) - (8 + 2 * pi * .33)) < 8e-5 # should be 7.66198e-5 + + # check continuity of a function that's in the global polynomial space + (; physical_frame_elements) = md.mesh_type + (; x, y) = md + u = @. x^rd.N - x * y^(rd.N-1) - x^(rd.N-1) * y + y^rd.N + 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 + @test all(uf .≈ vec(uf[md.mapP])) + + dudx_exact = @. rd.N * x^(rd.N-1) - y^(rd.N-1) - (rd.N-1) * x^(rd.N-2) * y + dudy_exact = @. -(rd.N-1) * x * y^(rd.N-2) - x^(rd.N-1) + rd.N * y^(rd.N-1) + (; physical_frame_elements, cut_face_nodes) = md.mesh_type + dudx, dudy = similar(md.x), similar(md.x) + dudx.cartesian .= (md.rxJ.cartesian .* (rd.Dr * u.cartesian)) ./ md.J.cartesian + dudy.cartesian .= (md.syJ.cartesian .* (rd.Ds * u.cartesian)) ./ md.J.cartesian + for (e, elem) in enumerate(physical_frame_elements) + VDM = vandermonde(elem, rd.N, x.cut[:, e], y.cut[:, e]) + Vq, Vxq, Vyq = map(A -> A / VDM, basis(elem, rd.N, md.xq.cut[:,e], md.yq.cut[:, e])) + + M = Vq' * diagm(md.wJq.cut[:, e]) * Vq + Qx = Vq' * diagm(md.wJq.cut[:, e]) * Vxq + Qy = Vq' * diagm(md.wJq.cut[:, e]) * Vyq + Dx, Dy = M \ Qx, M \ Qy + # LIFT = M \ (Vf' * diagm(wJf)) + + # TODO: add interface flux terms into test + dudx.cut[:, e] .= Dx * u.cut[:,e] # (md.rxJ.cut[:,e] .* (Dr * u.cut[:,e])) + dudy.cut[:, e] .= Dy * u.cut[:,e] # (md.rxJ.cut[:,e] .* (Dr * u.cut[:,e])) + end + + @test dudx ≈ dudx_exact + @test dudy ≈ dudy_exact + + # test normals are unit and non-zero + @test all(@. md.nx^2 + md.ny^2 ≈ 1) + + # test creation of equispaced nodes on cut cells + x, y = equi_nodes(physical_frame_elements[1], circle, 10) + # shouldn't have more points than equispaced points on a quad + @test 0 < length(x) <= length(first(equi_nodes(Quad(), 10))) + # no points should be contained in the circle + @test !all(PathIntersections.is_contained.(circle, zip(x, y))) + end - @test all(uf .≈ vec(uf[md.mapP])) - - dudx_exact = @. rd.N * x^(rd.N-1) - y^(rd.N-1) - (rd.N-1) * x^(rd.N-2) * y - dudy_exact = @. -(rd.N-1) * x * y^(rd.N-2) - x^(rd.N-1) + rd.N * y^(rd.N-1) - (; physical_frame_elements, cut_face_nodes) = md.mesh_type - 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) - VDM = vandermonde(elem, rd.N, x.cut[:, e], y.cut[:, e]) - Vq, Vxq, Vyq = map(A -> A / VDM, basis(elem, rd.N, md.xq.cut[:,e], md.yq.cut[:, e])) - - M = Vq' * diagm(md.wJq.cut[:, e]) * Vq - Qx = Vq' * diagm(md.wJq.cut[:, e]) * Vxq - Qy = Vq' * diagm(md.wJq.cut[:, e]) * Vyq - Dx, Dy = M \ Qx, M \ Qy - # LIFT = M \ (Vf' * diagm(wJf)) - - # TODO: add interface flux terms into test - dudx.cut[:, e] .= Dx * u.cut[:,e] # (md.rxJ.cut[:,e] .* (Dr * u.cut[:,e])) - dudy.cut[:, e] .= Dy * u.cut[:,e] # (md.rxJ.cut[:,e] .* (Dr * u.cut[:,e])) + + @testset "State redistribution" begin + # test state redistribution + cells_per_dimension = 4 + circle = PresetGeometries.Circle(R=0.6, x0=0, y0=0) + rd = RefElemData(Quad(), N=4) + md = MeshData(rd, (circle, ), cells_per_dimension, quadrature_type) + + srd = StateRedistribution(rd, md) + e = @. 0 * md.x + 1 # constant + u = @. md.x + md.x^3 .* md.y # degree 4 polynomial + ecopy, ucopy = copy.((e, u)) + + # two ways of applying SRD + apply!(u, srd) + srd(e) # in-place application of SRD functor + + # test exactness + @test norm(e .- ecopy) < 1e3 * eps() + @test norm(u .- ucopy) < 1e3 * eps() end - - @test dudx ≈ dudx_exact - @test dudy ≈ dudy_exact - - # test normals are unit and non-zero - @test all(@. md.nx^2 + md.ny^2 ≈ 1) - - # test creation of equispaced nodes on cut cells - x, y = equi_nodes(physical_frame_elements[1], circle, 10) - # shouldn't have more points than equispaced points on a quad - @test 0 < length(x) <= length(first(equi_nodes(Quad(), 10))) - # no points should be contained in the circle - @test !all(PathIntersections.is_contained.(circle, zip(x, y))) - -end - -@testset "State redistribution" begin - # test state redistribution - cells_per_dimension = 4 - circle = PresetGeometries.Circle(R=0.6, x0=0, y0=0) - rd = RefElemData(Quad(), N=4) - md = MeshData(MomentFitting(), rd, (circle, ), - cells_per_dimension, cells_per_dimension) - - srd = StateRedistribution(rd, md) - e = @. 0 * md.x + 1 # constant - u = @. md.x + md.x^3 .* md.y # degree 4 polynomial - ecopy, ucopy = copy.((e, u)) - - # two ways of applying SRD - apply!(u, srd) - srd(e) # in-place application of SRD functor - - # test exactness - @test norm(e .- ecopy) < 1e3 * eps() - @test norm(u .- ucopy) < 1e3 * eps() end \ No newline at end of file From fc1c211e74b9d40ac6f439884402156b81b61e39 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 16:31:16 -0500 Subject: [PATCH 33/39] add some docs --- docs/src/more_meshes.md | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/docs/src/more_meshes.md b/docs/src/more_meshes.md index 8190b26b..5f1e3bda 100644 --- a/docs/src/more_meshes.md +++ b/docs/src/more_meshes.md @@ -72,10 +72,15 @@ 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) +md = MeshData(rd, (circle, ), cells_per_dimension_x, cells_per_dimension_y, Subtriangulation(); 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 +Here, the final argument `quadrature_type = Subtriangulation()` determines how the quadrature on cut cells is determined. For `Subtriangulation()`, the quadrature on cut cells is constructed from a curved isoparametric subtriangulation of the cut cell. The number of quadrature points on a cut cell is then reduced (while preserving positivity) using Caratheodory pruning. If not specified, the `quadrature_type` argument defaults to `Subtriangulation()`. + +Quadrature rules can also be constructed by specifying `quadrature_type = MomentFitting()`. The quadrature points on cut cells `md.xq.cut` are determined from sampling and a pivoted QR decomposition. This is not recommended, as it can be both slower, and the cut-cell quadrature weights `md.wJq.cut` are not guaranteed to be positive. + +The interpolation points on cut cells `md.x.cut` are determined from sampled points and a pivoted QR decomposition. + +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 `NamedArrayPartition`s 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: From d0c8a63ab638e186bcd0e4bcfafba5203b48c668 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 16:31:45 -0500 Subject: [PATCH 34/39] allow specifying the target cut cell quadrature degree via keyword arg --- src/cut_cell_meshes.jl | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index d1771d07..1706e71e 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -643,23 +643,22 @@ to construct !!! Warning: this may be deprecated or removed in future versions. """ MeshData(rd::RefElemData, objects, cells_per_dimension, - quadrature_type=Subtriangulation(); kwargs...) = - MeshData(rd, objects, cells_per_dimension, cells_per_dimension, - quadrature_type; kwargs...) + quadrature_type=Subtriangulation(); kwargs...) = + MeshData(rd, objects, + cells_per_dimension, cells_per_dimension, + quadrature_type; kwargs...) function MeshData(rd::RefElemData, objects, cells_per_dimension_x::Int, cells_per_dimension_y::Int, quadrature_type=Subtriangulation(); - quad_rule_face = get_1d_quadrature(rd), coordinates_min=(-1.0, -1.0), coordinates_max=(1.0, 1.0), - precompute_operators=false) + kwargs...) # compute intersections of curve with a background Cartesian grid. vx = LinRange(coordinates_min[1], coordinates_max[1], cells_per_dimension_x + 1) vy = LinRange(coordinates_min[2], coordinates_max[2], cells_per_dimension_y + 1) - return MeshData(rd, objects, vx, vy, quadrature_type; - quad_rule_face, precompute_operators) + return MeshData(rd, objects, vx, vy, quadrature_type; kwargs...) end function MeshData(rd::RefElemData, objects, @@ -1380,6 +1379,7 @@ function MeshData(rd::RefElemData, objects, vx::AbstractVector, vy::AbstractVector, quadrature_type::Subtriangulation; quad_rule_face=get_1d_quadrature(rd), + cut_quadrature_target_degree = 2 * rd.N, precompute_operators=false) cells_per_dimension_x = length(vx) - 1 @@ -1405,31 +1405,32 @@ function MeshData(rd::RefElemData, objects, #################################################### physical_frame_elements = - StartUpDG.construct_physical_frame_elements(region_flags, vx, vy, cutcells) + construct_physical_frame_elements(region_flags, vx, vy, cutcells) x_cut, y_cut = - StartUpDG.construct_cut_interpolation_nodes(N, objects, physical_frame_elements) + construct_cut_interpolation_nodes(N, objects, physical_frame_elements) xq_cut, yq_cut, wJq_cut = - StartUpDG.construct_cut_volume_quadrature(N, cutcells, physical_frame_elements) + construct_cut_volume_quadrature(N, cutcells, physical_frame_elements; + target_degree = cut_quadrature_target_degree) xf_cut, yf_cut, nxJ_cut, nyJ_cut, wf_cut, cut_face_node_indices = - StartUpDG.construct_cut_surface_quadrature(N, cutcells, quad_rule_face) + construct_cut_surface_quadrature(N, cutcells, quad_rule_face) #################################################### # Construct Cartesian cells # #################################################### xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wf_cartesian = - StartUpDG.construct_cartesian_surface_quadrature(vx, vy, region_flags, quad_rule_face) + construct_cartesian_surface_quadrature(vx, vy, region_flags, quad_rule_face) xq_cartesian, yq_cartesian, wJq_cartesian = - StartUpDG.construct_cartesian_volume_quadrature(vx, vy, region_flags, + construct_cartesian_volume_quadrature(vx, vy, region_flags, (rd.rq, rd.sq, rd.wq)) # reuse quadrature mapping to construct interpolation nodes x_cartesian, y_cartesian, _ = - StartUpDG.construct_cartesian_volume_quadrature(vx, vy, region_flags, + construct_cartesian_volume_quadrature(vx, vy, region_flags, (rd.r, rd.s, ones(size(rd.r)))) #################################################### From a58c85fd1a02fcb9e23be25bb7ec102b1e35e477 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 17:10:32 -0500 Subject: [PATCH 35/39] committing scratch testing files before deleting d --- test_cutcell_MeshData.jl | 19 +++++++-------- test_cutcell_sbp.jl | 50 ++++++++++++++++++++-------------------- 2 files changed, 33 insertions(+), 36 deletions(-) diff --git a/test_cutcell_MeshData.jl b/test_cutcell_MeshData.jl index 971ba3a0..350010f1 100644 --- a/test_cutcell_MeshData.jl +++ b/test_cutcell_MeshData.jl @@ -6,25 +6,21 @@ using StartUpDG using PathIntersections N = 3 -quad_rule_face = gauss_lobatto_quad(0, 0, N) +quad_rule_face = gauss_quad(0, 0, N+1) rd = RefElemData(Quad(), N; quad_rule_face) cells_per_dimension = 4 circle = PresetGeometries.Circle(R=0.66, x0=0, y0=0) objects = (circle, ) -# TODO: remove eventually -cells_per_dimension_x = cells_per_dimension -cells_per_dimension_y = cells_per_dimension -vx = LinRange(-1, 1, cells_per_dimension_x + 1) -vy = LinRange(-1, 1, cells_per_dimension_y + 1) - ####################################### # test weak SBP property # ####################################### -md = StartUpDG.MeshData2(rd, objects, vx, vy; precompute_operators=true) +md = StartUpDG.MeshData(rd, objects, cells_per_dimension; precompute_operators=true) mt = md.mesh_type +(; wJf) = md.mesh_type.cut_cell_data +wf = wJf ./ md.Jf for (e, elem) in enumerate(mt.physical_frame_elements) xq, yq, wq = md.xq.cut[:,e], md.yq.cut[:,e], md.wJq.cut[:,e] @@ -34,9 +30,8 @@ for (e, elem) in enumerate(mt.physical_frame_elements) Vf = vandermonde(elem, rd.N, md.xf.cut[face_ids], md.yf.cut[face_ids]) Qx, Qy = Vq' * diagm(wq) * Vxq, Vq' * diagm(wq) * Vyq - (; wJf) = md.mesh_type.cut_cell_data - Bx = Diagonal(wJf.cut[face_ids] .* md.nxJ.cut[face_ids]) - By = Diagonal(wJf.cut[face_ids] .* md.nyJ.cut[face_ids]) + Bx = Diagonal(wf.cut[face_ids] .* md.nxJ.cut[face_ids]) + By = Diagonal(wf.cut[face_ids] .* md.nyJ.cut[face_ids]) # represent constant vector in basis ee = Vq \ ones(size(Vq, 1)) @@ -54,6 +49,8 @@ rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) Np_target = StartUpDG.Np_cut(target_degree) +(; cutcells) = md.mesh_type.cut_cell_data + plot() for e in eachindex(cutcells) xq, yq, _ = StartUpDG.subtriangulated_cutcell_quadrature(cutcells[e], rd_tri) diff --git a/test_cutcell_sbp.jl b/test_cutcell_sbp.jl index 020ec253..fc084fa8 100644 --- a/test_cutcell_sbp.jl +++ b/test_cutcell_sbp.jl @@ -13,8 +13,6 @@ cells_per_dimension = 4 circle = PresetGeometries.Circle(R=0.66, x0=0, y0=0) objects = (circle, ) -# md = MeshData(rd, objects, cells_per_dimension; precompute_operators=true) - # TODO: remove eventually cells_per_dimension_x = cells_per_dimension cells_per_dimension_y = cells_per_dimension @@ -23,11 +21,13 @@ vy = LinRange(-1, 1, cells_per_dimension_y + 1) N = rd.N +region_flags, cutcells = StartUpDG.calculate_cutcells(vx, vy, objects) + #################################################### # Construct Cartesian cells # #################################################### -xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wJf_cartesian = +xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wf_cartesian = StartUpDG.construct_cartesian_surface_quadrature(vx, vy, region_flags, quad_rule_face) xq_cartesian, yq_cartesian, wJq_cartesian = @@ -43,7 +43,7 @@ region_flags, cutcells = StartUpDG.calculate_cutcells(vx, vy, objects) physical_frame_elements = StartUpDG.construct_physical_frame_elements(region_flags, vx, vy, cutcells) -xf_cut, yf_cut, nxJ_cut, nyJ_cut, wJf_cut, cut_face_node_indices = +xf_cut, yf_cut, nxJ_cut, nyJ_cut, wf_cut, cut_face_node_indices = StartUpDG.construct_cut_surface_quadrature(N, cutcells, quad_rule_face) xq_pruned, yq_pruned, wJq_pruned = @@ -53,7 +53,7 @@ xq_pruned, yq_pruned, wJq_pruned = # test weak SBP property # ####################################### -e = 5 +e = 9 elem = physical_frame_elements[e] xq, yq, wq = xq_pruned[:,e], yq_pruned[:,e], wJq_pruned[:,e] @@ -61,29 +61,29 @@ Vq, Vxq, Vyq = basis(elem, rd.N, xq, yq) Vf = vandermonde(elem, rd.N, xf_cut[e], yf_cut[e]) Qx, Qy = Vq' * diagm(wq) * Vxq, Vq' * diagm(wq) * Vyq -Bx = Diagonal(wJf_cut[e] .* nxJ_cut[e]) -By = Diagonal(wJf_cut[e] .* nyJ_cut[e]) +Bx = Diagonal(wf_cut[e] .* nxJ_cut[e]) +By = Diagonal(wf_cut[e] .* nyJ_cut[e]) # represent constant vector in basis ee = Vq \ ones(size(Vq, 1)) @show norm(ee' * Qx - ee' * Vf' * Bx * Vf) @show norm(ee' * Qy - ee' * Vf' * By * Vf) -# # plot pruned cells - -# # 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 = max(2 * N, (N-1) * (N + 2)) -# target_degree = 2 * N -# rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, -# quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) -# Np_target = StartUpDG.Np_cut(target_degree) - -# 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!(xq_pruned[:,e], yq_pruned[:,e], markersize=8, marker=:circle, -# z_order=:back, label="Caratheodory pruning", leg=false) -# end -# display(plot!(leg=false)) +# plot pruned cells + +# 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 = max(2 * N, (N-1) * (N + 2)) +target_degree = 2 * N +rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, + quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) +Np_target = StartUpDG.Np_cut(target_degree) + +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!(xq_pruned[:,e], yq_pruned[:,e], markersize=8, marker=:circle, + z_order=:back, label="Caratheodory pruning", leg=false) +end +display(plot!(leg=false)) From 696ab771f5df86fa6fe538404886b669cccd8845 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 17:11:03 -0500 Subject: [PATCH 36/39] removing some scratch test files --- plot_isoparametric_cut_mesh.jl | 230 --------------------------------- test_cut_derivative.jl | 51 -------- test_cutcell_convergence.jl | 74 ----------- test_face_basis.jl | 36 ------ 4 files changed, 391 deletions(-) delete mode 100644 plot_isoparametric_cut_mesh.jl delete mode 100644 test_cut_derivative.jl delete mode 100644 test_cutcell_convergence.jl delete mode 100644 test_face_basis.jl diff --git a/plot_isoparametric_cut_mesh.jl b/plot_isoparametric_cut_mesh.jl deleted file mode 100644 index 9efe437b..00000000 --- a/plot_isoparametric_cut_mesh.jl +++ /dev/null @@ -1,230 +0,0 @@ -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]' diff --git a/test_cut_derivative.jl b/test_cut_derivative.jl deleted file mode 100644 index e5274524..00000000 --- a/test_cut_derivative.jl +++ /dev/null @@ -1,51 +0,0 @@ -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) \ No newline at end of file diff --git a/test_cutcell_convergence.jl b/test_cutcell_convergence.jl deleted file mode 100644 index 3e84d4ff..00000000 --- a/test_cutcell_convergence.jl +++ /dev/null @@ -1,74 +0,0 @@ -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) \ No newline at end of file diff --git a/test_face_basis.jl b/test_face_basis.jl deleted file mode 100644 index fd7b29ec..00000000 --- a/test_face_basis.jl +++ /dev/null @@ -1,36 +0,0 @@ -using NodesAndModes: face_basis -using StartUpDG - -N=7 -rd = RefElemData(Tri(), N) -md_ref = MeshData(uniform_mesh(Tri(), 1)..., rd) -(; x, y) = md_ref -# x = @. x + 0.5 * cos(pi/2 * x) * cos(pi/2 * y) -# y = @. y + 0.5 * cos(pi/2 * x) * cos(pi/2 * y) - -x = @. x + 0.25 * cos(pi/2 * x) * sin(pi/2 * y + 1) -y = @. y + 0.25 * cos(pi/2 * x + 2) * cos(pi/2 * y) - -# x = x + 0.1 * randn(size(x)) -# y = y + 0.1 * randn(size(x)) -(; Dr, Ds) = rd -xr, xs = Dr * x, Ds * x -yr, ys = Dr * y, Ds * y - -# interpolate to quadrature directly -xrq, xsq, yrq, ysq = (x -> rd.Vq * x).((xr, xs, yr, ys)) -Jq = @. -xsq * yrq + xrq * ysq -wJq = Diagonal(rd.wq) * Jq - -md = MeshData(rd, md_ref, x, y) - -@show sum(wJq), sum(md.wJq) - -Vf_face = face_basis(Tri(), rd.N, rd.rstf...) -norm(Vf_face * (Vf_face \ md.xf) - md.xf) - -@show sum(md_ref.wJq), sum(md.wJq) -@show sum(md_ref.wJq .* md_ref.xq), sum(md.wJq .* md.xq) -@show sum(md_ref.wJq .* md_ref.xq.^2), sum(md.wJq .* md.xq.^2) - - From 78aaf6ab799f98b6b67987be1375b04e4d98e6fc Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 17:11:12 -0500 Subject: [PATCH 37/39] add tests for the weak SBP property --- test/cut_mesh_tests.jl | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/test/cut_mesh_tests.jl b/test/cut_mesh_tests.jl index 626a3301..c15de575 100644 --- a/test/cut_mesh_tests.jl +++ b/test/cut_mesh_tests.jl @@ -99,4 +99,39 @@ using StartUpDG: PathIntersections @test norm(e .- ecopy) < 1e3 * eps() @test norm(u .- ucopy) < 1e3 * eps() end +end + +####################################### +# test weak SBP property # +####################################### +@testset "Cut mesh weak SBP property" begin + + cells_per_dimension = 2 + circle = PresetGeometries.Circle(R=0.66, x0=.1, y0=0) + + rd = RefElemData(Quad(), N=3) + + md = MeshData(rd, (circle, ), cells_per_dimension; + precompute_operators=true) + + mt = md.mesh_type + (; wJf) = md.mesh_type.cut_cell_data + wf = wJf ./ md.Jf + + for (e, elem) in enumerate(mt.physical_frame_elements) + xq, yq, wq = md.xq.cut[:,e], md.yq.cut[:,e], md.wJq.cut[:,e] + + Vq, Vxq, Vyq = basis(elem, rd.N, xq, yq) + face_ids = mt.cut_face_nodes[e] + Vf = vandermonde(elem, rd.N, md.xf.cut[face_ids], md.yf.cut[face_ids]) + Qx, Qy = Vq' * diagm(wq) * Vxq, Vq' * diagm(wq) * Vyq + + Bx = Diagonal(wf.cut[face_ids] .* md.nxJ.cut[face_ids]) + By = Diagonal(wf.cut[face_ids] .* md.nyJ.cut[face_ids]) + + # represent constant vector in basis + ee = Vq \ ones(size(Vq, 1)) + @test norm(ee' * Qx - ee' * Vf' * Bx * Vf) < 100 * eps() + @test norm(ee' * Qy - ee' * Vf' * By * Vf) < 100 * eps() + end end \ No newline at end of file From a1317fb19d1451876a7b2663cdf2ccadd3fe7f09 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 17:11:31 -0500 Subject: [PATCH 38/39] fix an error when num_cartesian_cells = 0 --- src/cut_cell_meshes.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 1706e71e..c5cd77ac 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -1483,7 +1483,8 @@ function MeshData(rd::RefElemData, objects, cut_face_node_indices_by_face)) # create node mappings - mapM_cartesian = reshape(eachindex(xf.cartesian), :, num_cartesian_cells) + mapM_cartesian = num_cartesian_cells > 0 ? + reshape(eachindex(xf.cartesian), :, num_cartesian_cells) : Int64[] mapM = NamedArrayPartition(cartesian = collect(mapM_cartesian), cut = collect(eachindex(xf.cut) .+ length(mapM_cartesian))) mapP = copy(mapM) From ca13fd3ec31e939946f85bb887aeab96e3e81d60 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 30 Apr 2024 20:47:43 -0500 Subject: [PATCH 39/39] remove scratchpad files --- test_cutcell_MeshData.jl | 61 --------------------------- test_cutcell_sbp.jl | 89 ---------------------------------------- 2 files changed, 150 deletions(-) delete mode 100644 test_cutcell_MeshData.jl delete mode 100644 test_cutcell_sbp.jl diff --git a/test_cutcell_MeshData.jl b/test_cutcell_MeshData.jl deleted file mode 100644 index 350010f1..00000000 --- a/test_cutcell_MeshData.jl +++ /dev/null @@ -1,61 +0,0 @@ -using NodesAndModes: face_basis -using Plots -using LinearAlgebra -using StaticArrays -using StartUpDG -using PathIntersections - -N = 3 -quad_rule_face = gauss_quad(0, 0, N+1) -rd = RefElemData(Quad(), N; quad_rule_face) - -cells_per_dimension = 4 -circle = PresetGeometries.Circle(R=0.66, x0=0, y0=0) -objects = (circle, ) - -####################################### -# test weak SBP property # -####################################### - -md = StartUpDG.MeshData(rd, objects, cells_per_dimension; precompute_operators=true) -mt = md.mesh_type -(; wJf) = md.mesh_type.cut_cell_data -wf = wJf ./ md.Jf - -for (e, elem) in enumerate(mt.physical_frame_elements) - xq, yq, wq = md.xq.cut[:,e], md.yq.cut[:,e], md.wJq.cut[:,e] - - Vq, Vxq, Vyq = basis(elem, rd.N, xq, yq) - face_ids = mt.cut_face_nodes[e] - Vf = vandermonde(elem, rd.N, md.xf.cut[face_ids], md.yf.cut[face_ids]) - Qx, Qy = Vq' * diagm(wq) * Vxq, Vq' * diagm(wq) * Vyq - - Bx = Diagonal(wf.cut[face_ids] .* md.nxJ.cut[face_ids]) - By = Diagonal(wf.cut[face_ids] .* md.nyJ.cut[face_ids]) - - # represent constant vector in basis - ee = Vq \ ones(size(Vq, 1)) - @show norm(ee' * Qx - ee' * Vf' * Bx * Vf) - @show norm(ee' * Qy - ee' * Vf' * By * Vf) -end - -# plot pruned cells - -# 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 = max(2 * N, (N-1) * (N + 2)) -target_degree = 2 * N -rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, - quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) -Np_target = StartUpDG.Np_cut(target_degree) - -(; cutcells) = md.mesh_type.cut_cell_data - -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!(xq_pruned[:,e], yq_pruned[:,e], markersize=8, marker=:circle, - z_order=:back, label="Caratheodory pruning", leg=false) -end -display(plot!(leg=false)) diff --git a/test_cutcell_sbp.jl b/test_cutcell_sbp.jl deleted file mode 100644 index fc084fa8..00000000 --- a/test_cutcell_sbp.jl +++ /dev/null @@ -1,89 +0,0 @@ -using NodesAndModes: face_basis -using Plots -using LinearAlgebra -using StaticArrays -using StartUpDG -using PathIntersections - -N = 3 -quad_rule_face = gauss_lobatto_quad(0, 0, N) -rd = RefElemData(Quad(), N; quad_rule_face) - -cells_per_dimension = 4 -circle = PresetGeometries.Circle(R=0.66, x0=0, y0=0) -objects = (circle, ) - -# TODO: remove eventually -cells_per_dimension_x = cells_per_dimension -cells_per_dimension_y = cells_per_dimension -vx = LinRange(-1, 1, cells_per_dimension_x + 1) -vy = LinRange(-1, 1, cells_per_dimension_y + 1) - -N = rd.N - -region_flags, cutcells = StartUpDG.calculate_cutcells(vx, vy, objects) - -#################################################### -# Construct Cartesian cells # -#################################################### - -xf_cartesian, yf_cartesian, nxJ_cartesian, nyJ_cartesian, wf_cartesian = - StartUpDG.construct_cartesian_surface_quadrature(vx, vy, region_flags, quad_rule_face) - -xq_cartesian, yq_cartesian, wJq_cartesian = - StartUpDG.construct_cartesian_volume_quadrature(vx, vy, region_flags, - (rd.rq, rd.sq, rd.wq)) - -#################################################### -# Construct cut cells stuff # -#################################################### - -region_flags, cutcells = StartUpDG.calculate_cutcells(vx, vy, objects) - -physical_frame_elements = - StartUpDG.construct_physical_frame_elements(region_flags, vx, vy, cutcells) - -xf_cut, yf_cut, nxJ_cut, nyJ_cut, wf_cut, cut_face_node_indices = - StartUpDG.construct_cut_surface_quadrature(N, cutcells, quad_rule_face) - -xq_pruned, yq_pruned, wJq_pruned = - StartUpDG.construct_cut_volume_quadrature(N, cutcells, physical_frame_elements) - -####################################### -# test weak SBP property # -####################################### - -e = 9 -elem = physical_frame_elements[e] -xq, yq, wq = xq_pruned[:,e], yq_pruned[:,e], wJq_pruned[:,e] - -Vq, Vxq, Vyq = basis(elem, rd.N, xq, yq) -Vf = vandermonde(elem, rd.N, xf_cut[e], yf_cut[e]) -Qx, Qy = Vq' * diagm(wq) * Vxq, Vq' * diagm(wq) * Vyq - -Bx = Diagonal(wf_cut[e] .* nxJ_cut[e]) -By = Diagonal(wf_cut[e] .* nyJ_cut[e]) - -# represent constant vector in basis -ee = Vq \ ones(size(Vq, 1)) -@show norm(ee' * Qx - ee' * Vf' * Bx * Vf) -@show norm(ee' * Qy - ee' * Vf' * By * Vf) - -# plot pruned cells - -# 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 = max(2 * N, (N-1) * (N + 2)) -target_degree = 2 * N -rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, - quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) -Np_target = StartUpDG.Np_cut(target_degree) - -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!(xq_pruned[:,e], yq_pruned[:,e], markersize=8, marker=:circle, - z_order=:back, label="Caratheodory pruning", leg=false) -end -display(plot!(leg=false))