From 8e01c266d60a52e57de1c56a6a0cfd3cf118c49c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 22 Feb 2023 13:28:22 -0600 Subject: [PATCH 1/5] 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 2/5] 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 3/5] 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 4/5] 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 5/5] 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)