From add80d413c9c97f8d6d63c1794b30f7ea3b31418 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Tue, 30 Apr 2024 21:44:10 -0500 Subject: [PATCH] Add new cut cell `MeshData` (#165) * add some temporary test files * fix docstrings * remove cruft * update plot made nicer plot of Caratheodory pruning for Christina's proposal * update file name Makes it mroe clear this file is meant for plotting * fixing cut cell demo for v1.0+ * adding Caratheodory pruning * cleaning up cutcell demo * improve efficiency slightly * refactoring functions * improve comments * more refactoring * generalize map_to_interval * add new version of "generate_sampling_points" * add dispatch to preserve old version of MeshData * remove cruft * add new routines to create a cut cell MeshData with positive weights * refactoring * fix construction of cut cell face node indices previously assumed same number of nodes on all faces * fix precomputation of operators * clean up test of SBP property * add test of SBP property using new MeshData * add MomentFitting dispatch for old cut cell MeshData d * remove outdated todo * changing wJf to wf internally * add face node index array * make face centroid computation more compact * specialize connect_mesh * add new MeshData based on subtriangulations * add Subtriangulation() as a quadrature type d * add quadrature type to CutCellMesh meshtype * test both Subtriangulation and MomentFitting * add some docs * allow specifying the target cut cell quadrature degree via keyword arg * committing scratch testing files before deleting d * removing some scratch test files * add tests for the weak SBP property * fix an error when num_cartesian_cells = 0 * remove scratchpad files --- docs/src/more_meshes.md | 11 +- src/StartUpDG.jl | 1 + src/cut_cell_meshes.jl | 871 ++++++++++++++++++++++++++++++++++-- src/physical_frame_basis.jl | 59 ++- test/cut_mesh_tests.jl | 209 +++++---- 5 files changed, 1027 insertions(+), 124 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: diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index 11acfb57..ab146804 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 Subtriangulation, MomentFitting include("state_redistribution.jl") export StateRedistribution, apply! diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 3a3a1bc4..c5cd77ac 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -20,24 +20,15 @@ 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 -# 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 print(io,"Cut-cell MeshData of dimension $DIM with $(md.num_elements) elements " * @@ -45,7 +36,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)) @@ -71,7 +62,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) @@ -175,6 +166,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)` @@ -401,7 +424,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 @@ -471,8 +494,6 @@ function connect_mesh(rd, face_centroids, cutcell_data; tol = 1e2 * eps()) 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 @@ -525,35 +546,124 @@ function calculate_cutcells(vx, vy, objects, ds = 1e-3, arc_tol = 1e-10, corner_ end """ - function MeshData(rd, geometry, vxyz...) + 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 + +struct MomentFitting end +struct Subtriangulation end -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...) + MeshData(rd::RefElemData, objects, + vx::AbstractVector, vy::AbstractVector, + quadrature_type::MomentFitting; + 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, + 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; - quad_rule_face = get_1d_quadrature(rd), + cells_per_dimension_x::Int, cells_per_dimension_y::Int, + quadrature_type=Subtriangulation(); 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; quad_rule_face, precompute_operators) + return MeshData(rd, objects, vx, vy, quadrature_type; kwargs...) end function MeshData(rd::RefElemData, objects, - vx::AbstractVector, vy::AbstractVector; + vx::AbstractVector, vy::AbstractVector, + quadrature_type::MomentFitting; quad_rule_face=get_1d_quadrature(rd), precompute_operators=false) @@ -787,7 +897,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) @@ -801,3 +915,694 @@ 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, 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 + + # 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)) + wf[face_ids, e] .= w1D + + # 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)) + wf[face_ids, e] .= w1D + + # 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] + wf[face_ids, e] .= w1D + + # 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] + wf[face_ids, e] .= w1D + + e = e + 1 + end + end + + return xf, yf, nxJ, nyJ, wf +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}}` 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. +""" +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 + + # 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 + +# 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) + face_centroids_cartesian = + zeros(SVector{2, Float64}, num_faces(Quad()), num_cartesian_cells) + + 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_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 + + 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_cut[e][f] = SVector(x, y) + end + end + + return face_centroids_cartesian, face_centroids_cut +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 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. +""" + +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 + 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 = 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) + + 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; + target_degree = cut_quadrature_target_degree) + + xf_cut, yf_cut, nxJ_cut, nyJ_cut, wf_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, wf_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, wf = + map((x, y) -> NamedArrayPartition((; cartesian=x, cut=y)), + (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 # + #################################################### + + # 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 = 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) + + for f in eachindex(FToF) + # 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) + + # 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. 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) + + 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_indices_by_cell + 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_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_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) + 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) + + 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) + +end \ No newline at end of file diff --git a/src/physical_frame_basis.jl b/src/physical_frame_basis.jl index e81870aa..55af85a0 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 @@ -144,4 +146,55 @@ 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 + +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) + +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 diff --git a/test/cut_mesh_tests.jl b/test/cut_mesh_tests.jl index 92ccb15a..c15de575 100644 --- a/test/cut_mesh_tests.jl +++ b/test/cut_mesh_tests.jl @@ -1,98 +1,137 @@ 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) +@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))) - rd = RefElemData(Quad(), N=3) - md = MeshData(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 + 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(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 +end - @test (@capture_out Base.show(stdout, MIME"text/plain"(), md)) == "Cut-cell MeshData of dimension 2 with 16 elements (12 Cartesian, 4 cut)" +####################################### +# test weak SBP property # +####################################### +@testset "Cut mesh weak SBP property" begin - # test constructor with only one "cells_per_dimension" argument - @test_nowarn MeshData(rd, (circle, ), cells_per_dimension_x) + cells_per_dimension = 2 + circle = PresetGeometries.Circle(R=0.66, x0=.1, y0=0) - # check the volume of the domain - A = 4 - pi * .33^2 - @test sum(md.wJq) ≈ A + rd = RefElemData(Quad(), N=3) - # 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), :)) - @test sum(wJf[md.mapB]) ≈ (8 + 2 * pi * .33) + md = MeshData(rd, (circle, ), cells_per_dimension; + precompute_operators=true) - # 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 - 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])) - end + mt = md.mesh_type + (; wJf) = md.mesh_type.cut_cell_data + wf = wJf ./ md.Jf - @test dudx ≈ dudx_exact - @test dudy ≈ dudy_exact + for (e, elem) in enumerate(mt.physical_frame_elements) + xq, yq, wq = md.xq.cut[:,e], md.yq.cut[:,e], md.wJq.cut[:,e] - # test normals are unit and non-zero - @test all(@. md.nx^2 + md.ny^2 ≈ 1) + 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 - # 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))) + Bx = Diagonal(wf.cut[face_ids] .* md.nxJ.cut[face_ids]) + By = Diagonal(wf.cut[face_ids] .* md.nyJ.cut[face_ids]) -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(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() + # 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