From 5c250c9b6371d0bd32cf8cd933e79bf03782c59b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 7 May 2024 15:25:48 -0500 Subject: [PATCH 01/10] move old MomentFitting cut cell code into separate file --- src/StartUpDG.jl | 6 +- src/cut_cell_meshes.jl | 636 +-------------------------------- src/cut_cell_moment_fitting.jl | 632 ++++++++++++++++++++++++++++++++ 3 files changed, 639 insertions(+), 635 deletions(-) create mode 100644 src/cut_cell_moment_fitting.jl diff --git a/src/StartUpDG.jl b/src/StartUpDG.jl index ab146804..45c70bd4 100644 --- a/src/StartUpDG.jl +++ b/src/StartUpDG.jl @@ -63,7 +63,11 @@ export num_faces, num_vertices, HybridMeshExample include("physical_frame_basis.jl") include("cut_cell_meshes.jl") export PhysicalFrame, equi_nodes -export Subtriangulation, MomentFitting +export Subtriangulation + +# ! this will be deprecated in a future release +include("cut_cell_moment_fitting.jl") +export MomentFitting include("state_redistribution.jl") export StateRedistribution, apply! diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index c5cd77ac..f038a9cc 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -38,16 +38,6 @@ end # maps x ∈ [-1,1] to [a,b] 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)) - for e in eachindex(cutcells) - curve = cutcells[e] - stop_points = curve.stop_pts - num_cut_faces[e] = length(stop_points) - 1 - end - return num_cut_faces -end - function neighbor_across_face(f, ex, ey) if f==1 return (ex-1, ey) @@ -62,39 +52,6 @@ function neighbor_across_face(f, ex, ey) end end -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) - num_cartesian_cells = sum(region_flags .== 0) - num_cut_faces = sum(cut_faces_per_cell) - - num_points_per_face = length(rd.rf) ÷ num_faces(rd.element_type) - - face_centroids_x = NamedArrayPartition(cartesian=zeros(num_faces(rd.element_type), num_cartesian_cells), - cut=zeros(num_cut_faces)) - face_centroids_y = similar(face_centroids_x) - - for e in 1:num_cartesian_cells - xf_element = reshape(view(xf.cartesian, :, e), num_points_per_face, num_faces(rd.element_type)) - yf_element = reshape(view(yf.cartesian, :, e), num_points_per_face, num_faces(rd.element_type)) - view(face_centroids_x.cartesian, :, e) .= vec(sum(xf_element, dims=1) / num_points_per_face) - view(face_centroids_y.cartesian, :, e) .= vec(sum(yf_element, dims=1) / num_points_per_face) - end - - for e in 1:num_cut_cells - face_node_ids = (1:(num_points_per_face * cut_faces_per_cell[e])) .+ cut_face_offsets[e] * num_points_per_face - xf_element = reshape(view(xf.cut, face_node_ids), num_points_per_face, cut_faces_per_cell[e]) - yf_element = reshape(view(yf.cut, face_node_ids), num_points_per_face, cut_faces_per_cell[e]) - - face_ids = (1:cut_faces_per_cell[e]) .+ cut_face_offsets[e] - view(face_centroids_x.cut, face_ids) .= vec(sum(xf_element, dims=1) / num_points_per_face) - view(face_centroids_y.cut, face_ids) .= vec(sum(yf_element, dims=1) / num_points_per_face) - end - - return face_centroids_x, face_centroids_y -end - function compute_element_indices(region_flags) cells_per_dimension_x, cells_per_dimension_y = size(region_flags) @@ -133,41 +90,6 @@ function is_inside_domain(ex, ey, regions) end end -# generates at least Np_target sampling points within a cut cell defined by `curve` -# returns both x_sampled, y_sampled (physical points inside the cut cell), as well as -# r_sampled, y_sampled (reference points which correspond to x_sampled, y_sampled). -function generate_sampling_points(objects, elem, rd, Np_target; N_sampled = 4 * rd.N) - - r_sampled, s_sampled = equi_nodes(rd.element_type, 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(rd.element_type, 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 - - -# 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 @@ -198,314 +120,6 @@ function generate_sampling_points(objects, elem, Np_target::Int, N_sampled::Int) 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)` -function map_points_to_cut_cell_faces(points, curve) - num_faces = length(curve.subcurves) - out = ntuple(_ -> similar(points, (length(points) * num_faces)), 5) - return map_points_to_cut_cell_faces!(out, points, curve) -end - -function map_points_to_cut_cell_faces!(out, points, curve) - xf, yf, nxJ, nyJ, Jf = out - stop_points = curve.stop_pts - r1D = points - for f in 1:length(stop_points)-1 - for i in eachindex(r1D) - fid = i + (f-1) * length(r1D) - - s = map_to_interval(r1D[i], stop_points[f], stop_points[f+1]) - - # compute tangent vector at a node, face Jacobian, and normals - x_node, y_node = curve(s) - xf[fid], yf[fid] = x_node, y_node - tangent_vector = PathIntersections.ForwardDiff.derivative(curve, s) - - # the face Jacobian involves scaling between mapped and reference face - # reference face = [-1, 1] - scaling = (stop_points[f+1] - stop_points[f]) / 2 - Jf[fid] = norm(tangent_vector) * scaling - - # we have to flip the sign to get the outward normal. - # note: we compute the scaled normal nxJ for consistency with other meshes. - normal_node = SVector{2}(tangent_vector[2], -tangent_vector[1]) - nxJ[fid], nyJ[fid] = (-normal_node / norm(normal_node)) .* Jf[fid] - end - end - return vec.((xf, yf, nxJ, nyJ, Jf)) -end - - -# Computes face geometric terms from a RefElemData, `quad_rule_face = (r1D, w1D)`, -# the vectors of the 1D vertex nodes `vx` and `vy`, and named tuple -# `cutcell_data is a NamedTuple containing `objects`, `region_flags`, `stop_pts``, `cutcells`. -function compute_geometric_data(rd::RefElemData{2, Quad}, quad_rule_face, - vx, vy, cutcell_data; tol=100 * eps()) - - # domain size and reference face weights - cells_per_dimension_x, cells_per_dimension_y = length(vx) - 1, length(vy) - 1 - LX, LY = (x -> x[2] - x[1]).(extrema.((vx, vy))) - - r1D, w1D = quad_rule_face - - (; objects, region_flags, cutcells, cut_faces_per_cell ) = cutcell_data - - # count number of cells and cut face nodes - num_cartesian_cells = sum(region_flags .== 0) - num_cut_cells = sum(region_flags .> 0) - nodes_per_face = length(r1D) - num_cut_face_nodes = nodes_per_face * sum(cut_faces_per_cell) - - # compute face data - xf, yf, nxJ, nyJ, Jf = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Nfq, num_cartesian_cells), - cut=zeros(num_cut_face_nodes)), 5) - - # the face Jacobian involves scaling between mapped and reference domain - # this is precomputed here since it's needed to compute the normals - face_ids_left_right = 1:(length(rd.rf) ÷ 2) - face_ids_top_bottom = ((length(rd.rf) ÷ 2) + 1):length(rd.rf) - Jf.cartesian[face_ids_top_bottom, :] .= LX / (cells_per_dimension_x * sum(w1D)) - Jf.cartesian[face_ids_left_right, :] .= LY / (cells_per_dimension_y * 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]) - vx_element = SVector(vx[ex], vx[ex + 1], vx[ex], vx[ex + 1]) - vy_element = SVector(vy[ey], vy[ey], vy[ey + 1], vy[ey + 1]) - x_element, y_element = map(x -> rd.V1 * x, (vx_element, vy_element)) - mul!(view(xf.cartesian, :, e), rd.Vf, x_element) - mul!(view(yf.cartesian, :, e), rd.Vf, y_element) - view(nxJ.cartesian, :, e) .= rd.nrJ .* view(Jf.cartesian, :, e) - view(nyJ.cartesian, :, e) .= rd.nsJ .* view(Jf.cartesian, :, e) - e = e + 1 - end - end - - e = 1 - offset = 0 - for (e, curve) in enumerate(cutcells) - # map 1D quadrature points to faces - num_cut_faces = length(curve.subcurves) - fids = (1:length(r1D) * num_cut_faces) .+ offset - out = map(x->view(x, fids), (xf.cut, yf.cut, nxJ.cut, nyJ.cut, Jf.cut)) - map_points_to_cut_cell_faces!(out, r1D, curve) - offset += length(fids) - end - - # compute face points + shifting/scaling coefficients for physical frame cut elements. - physical_frame_elements = PhysicalFrame{2}[] # populate this as we iterate through cut cells - - # store cut-cell scaling/shifting coefficients - (; cut_faces_per_cell, cut_face_offsets ) = cutcell_data - num_points_per_face = length(r1D) - - e = 1 - for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y - if is_cut(region_flags[ex, ey]) - - # here, we evaluate a PhysicalFrame basis by shifting and scaling the - # coordinates on an element back to the reference element [-1, 1]^2. - cut_face_node_ids = (1:num_points_per_face * cut_faces_per_cell[e]) .+ - num_points_per_face * cut_face_offsets[e] - - # store face nodes (extremal) and coordinates of background Cartesian cell - physical_frame_element = - PhysicalFrame(xf.cut[cut_face_node_ids], yf.cut[cut_face_node_ids], - SVector(vx[ex], vx[ex+1]), SVector(vy[ey], vy[ey+1])) - - push!(physical_frame_elements, physical_frame_element) - - e += 1 - end - end - - # interpolation points - x, y = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Np, num_cartesian_cells), - cut=zeros(Np_cut(rd.N), num_cut_cells)), 2) - - # compute interpolation points on cartesian elements - e = 1 - for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y - if is_Cartesian(region_flags[ex, ey]) - vx_element = SVector(vx[ex], vx[ex + 1], vx[ex], vx[ex + 1]) - vy_element = SVector(vy[ey], vy[ey], vy[ey + 1], vy[ey + 1]) - x_element, y_element = map(x -> rd.V1 * x, (vx_element, vy_element)) - view(x.cartesian, :, e) .= x_element - view(y.cartesian, :, e) .= y_element - e = e + 1 - end - end - - # 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, rd, 2 * Np_cut(rd.N)) - V = vandermonde(physical_frame_element, rd.N, x_sampled, y_sampled) - - # use pivoted QR to find good interpolation points - QRfac = qr(V', ColumnNorm()) - ids = QRfac.p[1:Np_cut(rd.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 = 20 * rd.N - iter, maxiters = 0, 100 - while condV > 1e8 && iter < maxiters - x_sampled, y_sampled = - generate_sampling_points(objects, physical_frame_element, rd, 2 * Np_cut(rd.N); - N_sampled = N_sampled) - V = vandermonde(physical_frame_element, rd.N, x_sampled, y_sampled) - - # use pivoted QR to find good interpolation points - QRfac = qr(V', ColumnNorm()) - ids = QRfac.p[1:Np_cut(rd.N)] - condV = cond(V[ids, :]) - N_sampled += 5 * rd.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.cut, :, e) .= x_sampled[ids] - view(y.cut, :, e) .= y_sampled[ids] - end - - # 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)) - - rstxyzJ = SMatrix{2, 2, typeof(rxJ), 4}(rxJ, sxJ, ryJ, syJ) # pack geometric terms together - - return physical_frame_elements, x, y, rstxyzJ, J, xf, yf, nxJ, nyJ, Jf -end - -""" - connect_mesh(rd, face_centroids, region_flags, cutcells; tol = 1e2 * eps()) - -Connects faces of a cut mesh to each other, returns `FToF` such that face -`f` is connected to `FToF[f]`. - -Inputs: -- rd::RefElemData -- face_centroids = (face_centroids_x, face_centroids_y), where `face_centroids_x/y` - are vectors of coordinates of face centroids -- `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::RefElemData, face_centroids, cutcell_data; tol = 1e2 * eps()) - - (; region_flags, cut_faces_per_cell, cut_face_offsets ) = cutcell_data - - cells_per_dimension_x, cells_per_dimension_y = size(region_flags) - num_cartesian_cells = sum(region_flags .== 0) - num_cut_faces = sum(cut_faces_per_cell) - num_total_faces = num_cut_faces + num_faces(rd.element_type) * 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 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y - - 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(rd.element_type)) .+ (e-1) * num_faces(rd.element_type) - 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(rd.element_type)) .+ (e_nbr-1) * num_faces(rd.element_type) - 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 - is_Cartesian(flag) = flag==0 ? true : false is_cut(flag) = flag > 0 @@ -626,13 +240,13 @@ function subtriangulated_cutcell_quadrature(cutcell, rd_tri, return xq, yq, wJq end -struct MomentFitting end +# quadrature type for cut cell meshes struct Subtriangulation end """ MeshData(rd::RefElemData, objects, vx::AbstractVector, vy::AbstractVector, - quadrature_type::MomentFitting; + quadrature_type=Subtriangulation(); quad_rule_face=get_1d_quadrature(rd), precompute_operators=false) @@ -661,252 +275,6 @@ function MeshData(rd::RefElemData, objects, return MeshData(rd, objects, vx, vy, quadrature_type; kwargs...) end -function MeshData(rd::RefElemData, objects, - vx::AbstractVector, vy::AbstractVector, - quadrature_type::MomentFitting; - quad_rule_face=get_1d_quadrature(rd), - precompute_operators=false) - - cells_per_dimension_x = length(vx) - 1 - cells_per_dimension_y = length(vy) - 1 - - # 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) - - # Compute volume, face points, and physical frame element scalings - physical_frame_elements, x, y, rstxyzJ, J, xf, yf, nxJ, nyJ, Jf = - compute_geometric_data(rd, quad_rule_face, vx, vy, cutcell_data) - - # Compute face-to-face connectivity by matching face centroids - face_centroids = compute_face_centroids(rd, xf, yf, cutcell_data) - FToF = connect_mesh(rd, face_centroids, cutcell_data) - - # Compute node-to-node mappings - # !!! Warning: this only works if the same quadrature rule is used for all 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 - - # compute cut-cell surface quadrature - _, w1D = quad_rule_face - wJf = similar(Jf) - wJf.cartesian = reshape(Diagonal(w1D) * reshape(Jf.cartesian, length(w1D), :), size(Jf.cartesian)) - wJf.cut = reshape(Diagonal(w1D) * reshape(Jf.cut, length(w1D), :), size(Jf.cut)) - - # The minimum number of cut cell quadrature points is `Np_cut(2 * rd.N)`. However, - # oversampling slightly seems to improve the conditioning of the quadrature weights. - num_cut_quad_points = Np_cut(2 * rd.N) + 1 - xq, yq, wJq = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Nq, num_cartesian_cells), - cut=zeros(num_cut_quad_points, num_cut_cells)), 3) - - # compute quadrature rules for the Cartesian cells - e = 1 - for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y - if is_Cartesian(region_flags[ex, ey]) - dx = vx[ex+1] - vx[ex] - dy = vy[ey+1] - vy[ey] - J = dx * dy / sum(rd.wq) - @. xq.cartesian[:, e] = dx * 0.5 * (1 + rd.rq) + vx[ex] - @. yq.cartesian[:, e] = dy * 0.5 * (1 + rd.sq) + vy[ey] - @. wJq.cartesian[:, e] = rd.wq * J - e += 1 - end - end - - # polynomial antidifferentiation operators for computing volume integrals - Ix, Iy = StartUpDG.antidiff_operators(2 * rd.N) - - # refine the surface rule used to compute the volume quadrature - if length(quad_rule_face[1]) < 3 * rd.N + 1 - r1D, w1D = gauss_quad(0, 0, 3 * rd.N) - else - r1D, w1D = quad_rule_face - end - # compute quadrature rules for the cut cells. integrate exactly degree 2N basis - for e in eachindex(cutcells) - # compute these quantities using a higher accuracy surface quadrature rule - xf_element, yf_element, nxJ_element, nyJ_element, Jf_element = - map_points_to_cut_cell_faces(r1D, cutcells[e]) - nx_element = nxJ_element ./ Jf_element - ny_element = nyJ_element ./ Jf_element - wJf_element = vec(Diagonal(w1D) * reshape(Jf_element, length(w1D), :)) - - # compute volume integrals using numerical Green's theorem and surface integrals - scaling = physical_frame_elements[e].scaling - Vf = vandermonde(physical_frame_elements[e], 2 * rd.N + 1, xf_element, yf_element) - b = 0.5 * ((Vf * Ix)' * (wJf_element .* nx_element) ./ scaling[1] + - (Vf * Iy)' * (wJf_element .* ny_element) ./ scaling[2]) - - # compute degree 2N basis matrix at sampled points - x_sampled, y_sampled = generate_sampling_points(objects, physical_frame_elements[e], - rd, 2 * Np_cut(rd.N); N_sampled = 8 * rd.N) - Vq = vandermonde(physical_frame_elements[e], 2 * rd.N, x_sampled, y_sampled) - - # if the condition number of the VDM matrix is large, select more points until cond(VDM) < 1e8 - QRfac = qr(Vq', ColumnNorm()) - - # if we have enough points, compute the cond number. Otherwise, assume it's Inf - if length(QRfac.p) >= num_cut_quad_points - ids = view(QRfac.p, 1:num_cut_quad_points) - condV_original = cond(Vq[ids, :]) - else - condV_original = Inf - end - condV = condV_original - N_sampled = 20 * rd.N # start with a lot more points - iter, maxiters = 0, 100 - while condV > 1e8 && iter < maxiters - x_sampled, y_sampled = - generate_sampling_points(objects, physical_frame_elements[e], rd, 2 * Np_cut(rd.N); - N_sampled = N_sampled) - - # if the number of sampled points isn't large enough, add more - if length(x_sampled) < num_cut_quad_points - N_sampled += 5 * rd.N - continue - end - - Vq = vandermonde(physical_frame_elements[e], 2 * rd.N, x_sampled, y_sampled) - - # use pivoted QR to find good interpolation points - QRfac = qr(Vq', ColumnNorm()) - ids = view(QRfac.p, 1:num_cut_quad_points) - condV = cond(Vq[ids, :]) - N_sampled += 5 * rd.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 - - # naive approach to computing quadrature weights; no guarantees of positivity - QRfac = qr(Vq', ColumnNorm()) - ids = view(QRfac.p, 1:num_cut_quad_points) - wq = Vq[ids, :]' \ b - - quadrature_error = norm(Vq[ids, :]' * wq - b) - quadrature_condition_number = sum(abs.(wq)) / sum(wq) - if quadrature_condition_number > 10 || quadrature_error > 1e3 * eps() - @warn "Quadrature error on element $e is $quadrature_error, " * - "quadrature condition number = $quadrature_condition_number. " * - "Condition number of quadrature VDM is $condV." - end - - if e==249 - if isdefined(Main, :Infiltrator) - Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) - end - end - - view(xq.cut, :, e) .= view(x_sampled, ids) - view(yq.cut, :, e) .= view(y_sampled, ids) - view(wJq.cut, :, e) .= wq - end - - # 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 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 - - # 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) - -end - num_cartesian_elements(md::MeshData{DIM, <:CutCellMesh}) where {DIM} = size(md.x.cartesian, 2) num_cut_elements(md::MeshData{DIM, <:CutCellMesh}) where {DIM} = size(md.x.cut, 2) diff --git a/src/cut_cell_moment_fitting.jl b/src/cut_cell_moment_fitting.jl new file mode 100644 index 00000000..23310fe2 --- /dev/null +++ b/src/cut_cell_moment_fitting.jl @@ -0,0 +1,632 @@ +# !!! Warning: MomentFitting cut cell quadrature is experimental +# !!! and may be removed in a future release. +struct MomentFitting end + +function MeshData(rd::RefElemData, objects, + vx::AbstractVector, vy::AbstractVector, + quadrature_type::MomentFitting; + quad_rule_face=get_1d_quadrature(rd), + precompute_operators=false) + + cells_per_dimension_x = length(vx) - 1 + cells_per_dimension_y = length(vy) - 1 + + # 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) + + # Compute volume, face points, and physical frame element scalings + physical_frame_elements, x, y, rstxyzJ, J, xf, yf, nxJ, nyJ, Jf = + compute_geometric_data(rd, quad_rule_face, vx, vy, cutcell_data) + + # Compute face-to-face connectivity by matching face centroids + face_centroids = compute_face_centroids(rd, xf, yf, cutcell_data) + FToF = connect_mesh(rd, face_centroids, cutcell_data) + + # Compute node-to-node mappings + # !!! Warning: this only works if the same quadrature rule is used for all 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 + + # compute cut-cell surface quadrature + _, w1D = quad_rule_face + wJf = similar(Jf) + wJf.cartesian = reshape(Diagonal(w1D) * reshape(Jf.cartesian, length(w1D), :), size(Jf.cartesian)) + wJf.cut = reshape(Diagonal(w1D) * reshape(Jf.cut, length(w1D), :), size(Jf.cut)) + + # The minimum number of cut cell quadrature points is `Np_cut(2 * rd.N)`. However, + # oversampling slightly seems to improve the conditioning of the quadrature weights. + num_cut_quad_points = Np_cut(2 * rd.N) + 1 + xq, yq, wJq = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Nq, num_cartesian_cells), + cut=zeros(num_cut_quad_points, num_cut_cells)), 3) + + # compute quadrature rules for the Cartesian cells + e = 1 + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + if is_Cartesian(region_flags[ex, ey]) + dx = vx[ex+1] - vx[ex] + dy = vy[ey+1] - vy[ey] + J = dx * dy / sum(rd.wq) + @. xq.cartesian[:, e] = dx * 0.5 * (1 + rd.rq) + vx[ex] + @. yq.cartesian[:, e] = dy * 0.5 * (1 + rd.sq) + vy[ey] + @. wJq.cartesian[:, e] = rd.wq * J + e += 1 + end + end + + # polynomial antidifferentiation operators for computing volume integrals + Ix, Iy = StartUpDG.antidiff_operators(2 * rd.N) + + # refine the surface rule used to compute the volume quadrature + if length(quad_rule_face[1]) < 3 * rd.N + 1 + r1D, w1D = gauss_quad(0, 0, 3 * rd.N) + else + r1D, w1D = quad_rule_face + end + # compute quadrature rules for the cut cells. integrate exactly degree 2N basis + for e in eachindex(cutcells) + # compute these quantities using a higher accuracy surface quadrature rule + xf_element, yf_element, nxJ_element, nyJ_element, Jf_element = + map_points_to_cut_cell_faces(r1D, cutcells[e]) + nx_element = nxJ_element ./ Jf_element + ny_element = nyJ_element ./ Jf_element + wJf_element = vec(Diagonal(w1D) * reshape(Jf_element, length(w1D), :)) + + # compute volume integrals using numerical Green's theorem and surface integrals + scaling = physical_frame_elements[e].scaling + Vf = vandermonde(physical_frame_elements[e], 2 * rd.N + 1, xf_element, yf_element) + b = 0.5 * ((Vf * Ix)' * (wJf_element .* nx_element) ./ scaling[1] + + (Vf * Iy)' * (wJf_element .* ny_element) ./ scaling[2]) + + # compute degree 2N basis matrix at sampled points + x_sampled, y_sampled = generate_sampling_points(objects, physical_frame_elements[e], + rd, 2 * Np_cut(rd.N); N_sampled = 8 * rd.N) + Vq = vandermonde(physical_frame_elements[e], 2 * rd.N, x_sampled, y_sampled) + + # if the condition number of the VDM matrix is large, select more points until cond(VDM) < 1e8 + QRfac = qr(Vq', ColumnNorm()) + + # if we have enough points, compute the cond number. Otherwise, assume it's Inf + if length(QRfac.p) >= num_cut_quad_points + ids = view(QRfac.p, 1:num_cut_quad_points) + condV_original = cond(Vq[ids, :]) + else + condV_original = Inf + end + condV = condV_original + N_sampled = 20 * rd.N # start with a lot more points + iter, maxiters = 0, 100 + while condV > 1e8 && iter < maxiters + x_sampled, y_sampled = + generate_sampling_points(objects, physical_frame_elements[e], rd, 2 * Np_cut(rd.N); + N_sampled = N_sampled) + + # if the number of sampled points isn't large enough, add more + if length(x_sampled) < num_cut_quad_points + N_sampled += 5 * rd.N + continue + end + + Vq = vandermonde(physical_frame_elements[e], 2 * rd.N, x_sampled, y_sampled) + + # use pivoted QR to find good interpolation points + QRfac = qr(Vq', ColumnNorm()) + ids = view(QRfac.p, 1:num_cut_quad_points) + condV = cond(Vq[ids, :]) + N_sampled += 5 * rd.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 + + # naive approach to computing quadrature weights; no guarantees of positivity + QRfac = qr(Vq', ColumnNorm()) + ids = view(QRfac.p, 1:num_cut_quad_points) + wq = Vq[ids, :]' \ b + + quadrature_error = norm(Vq[ids, :]' * wq - b) + quadrature_condition_number = sum(abs.(wq)) / sum(wq) + if quadrature_condition_number > 10 || quadrature_error > 1e3 * eps() + @warn "Quadrature error on element $e is $quadrature_error, " * + "quadrature condition number = $quadrature_condition_number. " * + "Condition number of quadrature VDM is $condV." + end + + if e==249 + if isdefined(Main, :Infiltrator) + Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) + end + end + + view(xq.cut, :, e) .= view(x_sampled, ids) + view(yq.cut, :, e) .= view(y_sampled, ids) + view(wJq.cut, :, e) .= wq + end + + # 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 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 + + # 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) + +end + +# Computes face geometric terms from a RefElemData, `quad_rule_face = (r1D, w1D)`, +# the vectors of the 1D vertex nodes `vx` and `vy`, and named tuple +# `cutcell_data is a NamedTuple containing `objects`, `region_flags`, `stop_pts``, `cutcells`. +function compute_geometric_data(rd::RefElemData{2, Quad}, quad_rule_face, + vx, vy, cutcell_data; tol=100 * eps()) + + # domain size and reference face weights + cells_per_dimension_x, cells_per_dimension_y = length(vx) - 1, length(vy) - 1 + LX, LY = (x -> x[2] - x[1]).(extrema.((vx, vy))) + + r1D, w1D = quad_rule_face + + (; objects, region_flags, cutcells, cut_faces_per_cell ) = cutcell_data + + # count number of cells and cut face nodes + num_cartesian_cells = sum(region_flags .== 0) + num_cut_cells = sum(region_flags .> 0) + nodes_per_face = length(r1D) + num_cut_face_nodes = nodes_per_face * sum(cut_faces_per_cell) + + # compute face data + xf, yf, nxJ, nyJ, Jf = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Nfq, num_cartesian_cells), + cut=zeros(num_cut_face_nodes)), 5) + + # the face Jacobian involves scaling between mapped and reference domain + # this is precomputed here since it's needed to compute the normals + face_ids_left_right = 1:(length(rd.rf) ÷ 2) + face_ids_top_bottom = ((length(rd.rf) ÷ 2) + 1):length(rd.rf) + Jf.cartesian[face_ids_top_bottom, :] .= LX / (cells_per_dimension_x * sum(w1D)) + Jf.cartesian[face_ids_left_right, :] .= LY / (cells_per_dimension_y * 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]) + vx_element = SVector(vx[ex], vx[ex + 1], vx[ex], vx[ex + 1]) + vy_element = SVector(vy[ey], vy[ey], vy[ey + 1], vy[ey + 1]) + x_element, y_element = map(x -> rd.V1 * x, (vx_element, vy_element)) + mul!(view(xf.cartesian, :, e), rd.Vf, x_element) + mul!(view(yf.cartesian, :, e), rd.Vf, y_element) + view(nxJ.cartesian, :, e) .= rd.nrJ .* view(Jf.cartesian, :, e) + view(nyJ.cartesian, :, e) .= rd.nsJ .* view(Jf.cartesian, :, e) + e = e + 1 + end + end + + e = 1 + offset = 0 + for (e, curve) in enumerate(cutcells) + # map 1D quadrature points to faces + num_cut_faces = length(curve.subcurves) + fids = (1:length(r1D) * num_cut_faces) .+ offset + out = map(x->view(x, fids), (xf.cut, yf.cut, nxJ.cut, nyJ.cut, Jf.cut)) + map_points_to_cut_cell_faces!(out, r1D, curve) + offset += length(fids) + end + + # compute face points + shifting/scaling coefficients for physical frame cut elements. + physical_frame_elements = PhysicalFrame{2}[] # populate this as we iterate through cut cells + + # store cut-cell scaling/shifting coefficients + (; cut_faces_per_cell, cut_face_offsets ) = cutcell_data + num_points_per_face = length(r1D) + + e = 1 + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + if is_cut(region_flags[ex, ey]) + + # here, we evaluate a PhysicalFrame basis by shifting and scaling the + # coordinates on an element back to the reference element [-1, 1]^2. + cut_face_node_ids = (1:num_points_per_face * cut_faces_per_cell[e]) .+ + num_points_per_face * cut_face_offsets[e] + + # store face nodes (extremal) and coordinates of background Cartesian cell + physical_frame_element = + PhysicalFrame(xf.cut[cut_face_node_ids], yf.cut[cut_face_node_ids], + SVector(vx[ex], vx[ex+1]), SVector(vy[ey], vy[ey+1])) + + push!(physical_frame_elements, physical_frame_element) + + e += 1 + end + end + + # interpolation points + x, y = ntuple(_ -> NamedArrayPartition(cartesian=zeros(rd.Np, num_cartesian_cells), + cut=zeros(Np_cut(rd.N), num_cut_cells)), 2) + + # compute interpolation points on cartesian elements + e = 1 + for ex in 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + if is_Cartesian(region_flags[ex, ey]) + vx_element = SVector(vx[ex], vx[ex + 1], vx[ex], vx[ex + 1]) + vy_element = SVector(vy[ey], vy[ey], vy[ey + 1], vy[ey + 1]) + x_element, y_element = map(x -> rd.V1 * x, (vx_element, vy_element)) + view(x.cartesian, :, e) .= x_element + view(y.cartesian, :, e) .= y_element + e = e + 1 + end + end + + # 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, rd, 2 * Np_cut(rd.N)) + V = vandermonde(physical_frame_element, rd.N, x_sampled, y_sampled) + + # use pivoted QR to find good interpolation points + QRfac = qr(V', ColumnNorm()) + ids = QRfac.p[1:Np_cut(rd.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 = 20 * rd.N + iter, maxiters = 0, 100 + while condV > 1e8 && iter < maxiters + x_sampled, y_sampled = + generate_sampling_points(objects, physical_frame_element, rd, 2 * Np_cut(rd.N); + N_sampled = N_sampled) + V = vandermonde(physical_frame_element, rd.N, x_sampled, y_sampled) + + # use pivoted QR to find good interpolation points + QRfac = qr(V', ColumnNorm()) + ids = QRfac.p[1:Np_cut(rd.N)] + condV = cond(V[ids, :]) + N_sampled += 5 * rd.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.cut, :, e) .= x_sampled[ids] + view(y.cut, :, e) .= y_sampled[ids] + end + + # 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)) + + rstxyzJ = SMatrix{2, 2, typeof(rxJ), 4}(rxJ, sxJ, ryJ, syJ) # pack geometric terms together + + return physical_frame_elements, x, y, rstxyzJ, J, xf, yf, nxJ, nyJ, Jf +end + +""" + connect_mesh(rd, face_centroids, region_flags, cutcells; tol = 1e2 * eps()) + +Connects faces of a cut mesh to each other, returns `FToF` such that face +`f` is connected to `FToF[f]`. + +Inputs: +- rd::RefElemData +- face_centroids = (face_centroids_x, face_centroids_y), where `face_centroids_x/y` + are vectors of coordinates of face centroids +- `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::RefElemData, face_centroids, cutcell_data; tol = 1e2 * eps()) + + (; region_flags, cut_faces_per_cell, cut_face_offsets ) = cutcell_data + + cells_per_dimension_x, cells_per_dimension_y = size(region_flags) + num_cartesian_cells = sum(region_flags .== 0) + num_cut_faces = sum(cut_faces_per_cell) + num_total_faces = num_cut_faces + num_faces(rd.element_type) * 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 1:cells_per_dimension_x, ey in 1:cells_per_dimension_y + + 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(rd.element_type)) .+ (e-1) * num_faces(rd.element_type) + 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(rd.element_type)) .+ (e_nbr-1) * num_faces(rd.element_type) + 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 + +# 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)` +function map_points_to_cut_cell_faces(points, curve) + num_faces = length(curve.subcurves) + out = ntuple(_ -> similar(points, (length(points) * num_faces)), 5) + return map_points_to_cut_cell_faces!(out, points, curve) +end + +function map_points_to_cut_cell_faces!(out, points, curve) + xf, yf, nxJ, nyJ, Jf = out + stop_points = curve.stop_pts + r1D = points + for f in 1:length(stop_points)-1 + for i in eachindex(r1D) + fid = i + (f-1) * length(r1D) + + s = map_to_interval(r1D[i], stop_points[f], stop_points[f+1]) + + # compute tangent vector at a node, face Jacobian, and normals + x_node, y_node = curve(s) + xf[fid], yf[fid] = x_node, y_node + tangent_vector = PathIntersections.ForwardDiff.derivative(curve, s) + + # the face Jacobian involves scaling between mapped and reference face + # reference face = [-1, 1] + scaling = (stop_points[f+1] - stop_points[f]) / 2 + Jf[fid] = norm(tangent_vector) * scaling + + # we have to flip the sign to get the outward normal. + # note: we compute the scaled normal nxJ for consistency with other meshes. + normal_node = SVector{2}(tangent_vector[2], -tangent_vector[1]) + nxJ[fid], nyJ[fid] = (-normal_node / norm(normal_node)) .* Jf[fid] + end + end + return vec.((xf, yf, nxJ, nyJ, Jf)) +end + +# generates at least Np_target sampling points within a cut cell defined by `curve` +# returns both x_sampled, y_sampled (physical points inside the cut cell), as well as +# r_sampled, y_sampled (reference points which correspond to x_sampled, y_sampled). +function generate_sampling_points(objects, elem, rd, Np_target; N_sampled = 4 * rd.N) + + r_sampled, s_sampled = equi_nodes(rd.element_type, 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(rd.element_type, 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 + +function count_cut_faces(cutcells) + num_cut_faces = zeros(Int, length(cutcells)) + for e in eachindex(cutcells) + curve = cutcells[e] + stop_points = curve.stop_pts + num_cut_faces[e] = length(stop_points) - 1 + end + return num_cut_faces +end + +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) + num_cartesian_cells = sum(region_flags .== 0) + num_cut_faces = sum(cut_faces_per_cell) + + num_points_per_face = length(rd.rf) ÷ num_faces(rd.element_type) + + face_centroids_x = NamedArrayPartition(cartesian=zeros(num_faces(rd.element_type), num_cartesian_cells), + cut=zeros(num_cut_faces)) + face_centroids_y = similar(face_centroids_x) + + for e in 1:num_cartesian_cells + xf_element = reshape(view(xf.cartesian, :, e), num_points_per_face, num_faces(rd.element_type)) + yf_element = reshape(view(yf.cartesian, :, e), num_points_per_face, num_faces(rd.element_type)) + view(face_centroids_x.cartesian, :, e) .= vec(sum(xf_element, dims=1) / num_points_per_face) + view(face_centroids_y.cartesian, :, e) .= vec(sum(yf_element, dims=1) / num_points_per_face) + end + + for e in 1:num_cut_cells + face_node_ids = (1:(num_points_per_face * cut_faces_per_cell[e])) .+ cut_face_offsets[e] * num_points_per_face + xf_element = reshape(view(xf.cut, face_node_ids), num_points_per_face, cut_faces_per_cell[e]) + yf_element = reshape(view(yf.cut, face_node_ids), num_points_per_face, cut_faces_per_cell[e]) + + face_ids = (1:cut_faces_per_cell[e]) .+ cut_face_offsets[e] + view(face_centroids_x.cut, face_ids) .= vec(sum(xf_element, dims=1) / num_points_per_face) + view(face_centroids_y.cut, face_ids) .= vec(sum(yf_element, dims=1) / num_points_per_face) + end + + return face_centroids_x, face_centroids_y +end \ No newline at end of file From 7069ce6e7b7d20499d3fff02016b664a85589e5c Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 9 May 2024 22:48:16 -0500 Subject: [PATCH 02/10] fix default target_degree for volume quadrature --- src/cut_cell_meshes.jl | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index f038a9cc..cf97e44b 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -285,18 +285,29 @@ function num_elements(md::MeshData{DIM, <:CutCellMesh}) where {DIM} end """ - construct_cut_volume_quadrature(N, cutcells; target_degree = 2 * N) + construct_cut_volume_quadrature(N, cutcells; target_degree = 2 * N - 1) Constructs volume quadrature using subtriangulations of cut cells and Caratheodory pruning. The resulting quadrature is exact for polynomials of degree `target_degree`. + +We set `target_degree` to `2 * N - 1` by default, which is sufficient to +ensure that ∫du/dx * v is integrated exactly so that integration by parts +holds under the generated cut cell quadrature. """ function construct_cut_volume_quadrature(N, cutcells, physical_frame_elements; - target_degree = 2 * N) + target_degree = 2 * N - 1) + + # We make volume quadrature exact for degree N^2 + N(N-1) + 2(N-1) polynomials on the + # reference element, which ensures that ∫du/dx * v * J is integrated exactly over D̂. + # This is because du/dx * v ∈ P^{2N-1}(D^k) and (du/dx * v) ∘ x(r,s) ∈ P^{(2N-1) * N}). + # + # The minimum exactness of volume quadrature is degree N(N-1) + 2N-2 polynomials, + # which ensures Qh * 1 = 0, where Qh = hybridized SBP operator. - # 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)) + # Integral ∫ v du/dx J over reference element D̂ + # TODO: why can we decrease this and still see exactness? + N_phys_frame_geo = N^2 + N * (N-1) + 2 * (N-1) rd_tri = RefElemData(Tri(), Polynomial(MultidimensionalQuadrature()), N, quad_rule_vol=NodesAndModes.quad_nodes_tri(N_phys_frame_geo)) @@ -323,14 +334,12 @@ function construct_cut_volume_quadrature(N, cutcells, physical_frame_elements; w = vec(wJq) w_pruned, inds = StartUpDG.caratheodory_pruning_qr(Vtarget, w) - if target_degree >= 2*N - # test exactness of the pruned quadrature rule + # test exactness of the pruned quadrature rule if applicable + if target_degree >= 2 * N 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] From 3bcd03430f2caad63cad0e01394f9e12a50baef0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 9 May 2024 22:48:44 -0500 Subject: [PATCH 03/10] formatting and modifying default parameters for cut interpolation nodes --- src/cut_cell_meshes.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index cf97e44b..6acd0bc9 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -351,7 +351,7 @@ end # points (e.g., QR-DEIM). function construct_cut_interpolation_nodes(N, objects, physical_frame_elements) - N_sampled = 4 * N + N_sampled = 4 * N # this is arbitrary num_cut_elements = length(physical_frame_elements) x, y = ntuple(_ -> zeros(Np_cut(N), num_cut_elements), 2) @@ -362,6 +362,7 @@ function construct_cut_interpolation_nodes(N, objects, physical_frame_elements) 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 @@ -372,7 +373,7 @@ function construct_cut_interpolation_nodes(N, objects, physical_frame_elements) # number of sampled points. condV_original = cond(V[ids, :]) condV = condV_original - N_sampled = 10 * N + N_sampled = 8 * N iter, maxiters = 0, 100 while condV > 1e8 && iter < maxiters x_sampled, y_sampled = @@ -384,17 +385,18 @@ function construct_cut_interpolation_nodes(N, objects, physical_frame_elements) QRfac = qr(V', ColumnNorm()) ids = QRfac.p[1:Np_cut(N)] condV = cond(V[ids, :]) - N_sampled += 5 * N + N_sampled += 2 * N if condV < 1e8 - @warn "Conditioning of old VDM for element $e is $condV_original. " * + @info "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." + @warn "Adaptive sampling of cut-cell element $e: " * + "maximum number of iterations $maxiters exceeded." end view(x, :, e) .= x_sampled[ids] From 722db6c57090b5b60e296f7e7dfbd86fd895da9e Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 9 May 2024 22:53:02 -0500 Subject: [PATCH 04/10] setting default boundary quadrature + renaming variables for clarity --- src/cut_cell_meshes.jl | 35 +++++++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 6acd0bc9..29628270 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -539,16 +539,14 @@ 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)) +function construct_cut_surface_quadrature(N, cutcells, + quad_rule_1D = gauss_quad(0, 0, N), + quad_rule_boundary = get_boundary_quadrature(N)) rd_line = RefElemData(Line(), N, quad_rule_vol = quad_rule_1D) + rq_boundary, wq_boundary = quad_rule_boundary - # 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 + interp_to_cut_boundary = vandermonde(Line(), rd_line.N, rq_boundary) / rd_line.VDM xf, yf, nxJ, nyJ, wf = ntuple(_ -> Vector{Float64}[], 5) @@ -584,7 +582,7 @@ function construct_cut_surface_quadrature(N, cutcells, quad_rule_1D = gauss_quad 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 + wf_face = (is_boundary_face) ? wq_boundary : rd_line.wq push!(xf_element, xfq_face) push!(yf_element, yfq_face) @@ -740,6 +738,18 @@ function connect_mesh(face_centroids, region_flags, return FToF end +# If we wish to exactly integrate ∫ u * v * (nx * Jf) and ∫ u * v * (ny * Jf), +# the reference quadrature must be exact for degree 2 * N^2 + (N-1) polynomials +# on the reference line. +# +# Note that we don't perform Caratheodory pruning for boundary faces at the moment. +get_boundary_quadrature(rd::RefElemData) = get_boundary_quadrature(rd.N) +function get_boundary_quadrature(N) + # N_boundary = ceil(Int, (N^2 + (N-1) - 1) / 2) # to ensure weak SBP + N_boundary = ceil(Int, (2 * N^2 + (N - 1) - 1) / 2) # to ensure GSBP property + return gauss_quad(0, 0, N_boundary) +end + """ function MeshData(rd, geometry, vxyz...) @@ -758,7 +768,8 @@ 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, + quad_rule_boundary=get_boundary_quadrature(rd), + cut_quadrature_target_degree = 2 * rd.N-1, precompute_operators=false) cells_per_dimension_x = length(vx) - 1 @@ -793,8 +804,12 @@ function MeshData(rd::RefElemData, objects, construct_cut_volume_quadrature(N, cutcells, physical_frame_elements; target_degree = cut_quadrature_target_degree) + # On the curved cut boundary, 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). quad_rule_boundary defaults to this 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_cut_surface_quadrature(N, cutcells, quad_rule_face, quad_rule_boundary) #################################################### # Construct Cartesian cells # From 1f6383fc53d20ffe95f7ef19fa144802d6da9068 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 9 May 2024 22:53:54 -0500 Subject: [PATCH 05/10] code cleanup --- src/cut_cell_meshes.jl | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 29628270..dc413321 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -923,16 +923,17 @@ function MeshData(rd::RefElemData, objects, # 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) + map(x -> UnitRange(first(x[1]), last(x[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)}[] + volume_interpolation_matrices, face_interpolation_matrices = + ntuple(_ -> Matrix{eltype(x)}[], 3) + + # TODO: remove these? these can be computed from other matrices + mass_matrices, lift_matrices = ntuple(_ -> Matrix{eltype(x)}[], 2) for (e, elem) in enumerate(physical_frame_elements) VDM = vandermonde(elem, rd.N, x.cut[:, e], y.cut[:, e]) @@ -944,20 +945,24 @@ function MeshData(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_node_indices_by_cell[e]], - yf.cut[cut_face_node_indices_by_cell[e]]) / VDM + xf_e = xf.cut[cut_face_node_indices_by_cell[e]] + yf_e = yf.cut[cut_face_node_indices_by_cell[e]] + Vf = vandermonde(elem, rd.N, xf_e, yf_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!(volume_interpolation_matrices, Vq) push!(face_interpolation_matrices, Vf) push!(differentiation_matrices, (Dx_e, Dy_e)) push!(mass_matrices, M) + push!(lift_matrices, M \ (Vf' * diagm(wf))) end - cut_cell_operators = (; differentiation_matrices, face_interpolation_matrices, + cut_cell_operators = (; volume_interpolation_matrices, + face_interpolation_matrices, + differentiation_matrices, mass_matrices, lift_matrices) else From 2e43244ec7d1fbc13c6080aac41d3c3a5e54b54b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 9 May 2024 22:54:13 -0500 Subject: [PATCH 06/10] storing cut_face_node_indices_by_elem_per_face --- src/cut_cell_meshes.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index dc413321..541057b0 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -918,6 +918,7 @@ function MeshData(rd::RefElemData, objects, vxyz=(vx, vy), cells_per_dimension=(cells_per_dimension_x, cells_per_dimension_y), + cut_face_node_indices_by_elem_by_face = cut_face_node_indices, ) # background Cartesian grid info # get flattened indices of cut face nodes. note that From a2bcc552ec5a10c1ae90c93c96ee73f2ccbfcaec Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 9 May 2024 22:54:25 -0500 Subject: [PATCH 07/10] adding hybridized SBP operators on cut cells --- src/cut_cell_meshes.jl | 48 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 541057b0..77d8abc3 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -1005,4 +1005,52 @@ function MeshData(rd::RefElemData, objects, FToF, (x, y), (xf, yf), (xq, yq), wJq, mapM, mapP, mapB, rstxyzJ, J, (nxJ, nyJ), Jf, is_periodic) +end + +function hybridized_SBP_operators(md::MeshData{2, <:CutCellMesh}) + mt = md.mesh_type + (; volume_interpolation_matrices, + face_interpolation_matrices, + differentiation_matrices, + mass_matrices) = mt.cut_cell_operators + (; wJf) = mt.cut_cell_data + wf = wJf ./ md.Jf + (; x, nxJ, nyJ) = md + + hybridized_operators = Tuple{Matrix{eltype(x)}, Matrix{eltype(x)}}[] + interpolation_operators, projection_operators, project_and_interp_operators = + ntuple(_ -> Matrix{eltype(x)}[], 3) + for e in eachindex(differentiation_matrices) + Dx, Dy = differentiation_matrices[e] + Vq = volume_interpolation_matrices[e] + Vf = face_interpolation_matrices[e] + M = mass_matrices[e] + + Pq = M \ (Vq' * Diagonal(md.wJq.cut[:, e])) + E = Vf * Pq + + fids = mt.cut_face_nodes[e] + Bx = Diagonal(wf.cut[fids] .* nxJ.cut[fids]) + By = Diagonal(wf.cut[fids] .* nyJ.cut[fids]) + + Qx = Pq' * M * Dx * Pq + Qy = Pq' * M * Dy * Pq + Qxh = 0.5 * [Qx - Qx' E' * Bx; + -Bx * E Bx ] + Qyh = 0.5 * [Qy - Qy' E' * By; + -By * E By ] + + # interpolation matrices + Vh = [Vq; Vf] + VhP = Vh * Pq + Ph = M \ Vh' + + push!(hybridized_operators, (Qxh, Qyh)) + push!(interpolation_operators, Vh) + push!(projection_operators, Ph) + push!(project_and_interp_operators, VhP) + end + + return hybridized_operators, project_and_interp_operators, + projection_operators, interpolation_operators end \ No newline at end of file From fa6c6b80afe4ec5e5e832322ff3d904b7074611e Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 9 May 2024 22:54:46 -0500 Subject: [PATCH 08/10] adding cut cell hybridized SBP tests --- test/cut_mesh_tests.jl | 52 ++++++++++++++++++++++++++---------------- 1 file changed, 32 insertions(+), 20 deletions(-) diff --git a/test/cut_mesh_tests.jl b/test/cut_mesh_tests.jl index c15de575..aa079fd9 100644 --- a/test/cut_mesh_tests.jl +++ b/test/cut_mesh_tests.jl @@ -109,29 +109,41 @@ end cells_per_dimension = 2 circle = PresetGeometries.Circle(R=0.66, x0=.1, y0=0) - rd = RefElemData(Quad(), N=3) + rd = RefElemData(Quad(), N=2) 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() + mt = md.mesh_type + + (; x, y, xq, yq) = md + u = @. x^2 + x * y + y^2 + uq = @. xq^2 + xq * yq + yq^2 + + (; volume_interpolation_matrices, differentiation_matrices) = md.mesh_type.cut_cell_operators + Qxyh, hybridized_project_interp_matrices, + hybridized_projection_matrices, hybridized_interp_matrices = + hybridized_SBP_operators(md) + + for e in eachindex(Qxyh) + Qxh, Qyh = Qxyh[e] + Dx, Dy = differentiation_matrices[e] + Vq = volume_interpolation_matrices[e] + M = Vq' * diagm(md.wJq.cut[:,e]) * Vq + Pq = M \ (Vq' * diagm(md.wJq.cut[:,e])) + + # test exactness of interpolation, projection + @test norm(Vq * u.cut[:,e] - uq.cut[:,e]) < 100 * eps() + @test norm(Pq * uq.cut[:,e] - u.cut[:,e]) < 100 * eps() + + # test weak SBP property + @test norm(sum(Qxh, dims=2)) + norm(sum(Qyh, dims=2)) < 100 * eps() + + # test accuracy of hybridized SBP operators when rd::RefElemData uses + # full accuracy quadrature. + Ph = hybridized_projection_matrices[e] + Vh = hybridized_interp_matrices[e] + @test norm(Ph * Qxh * Vh - Dx) + + norm(Ph * Qyh * Vh - Dy) < 1e3 * eps() end end \ No newline at end of file From b19c5250c6753cea9efd4742a77e909b0724d890 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 9 May 2024 23:05:00 -0500 Subject: [PATCH 09/10] update test tolerances --- test/cut_mesh_tests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/cut_mesh_tests.jl b/test/cut_mesh_tests.jl index aa079fd9..2a4cdae7 100644 --- a/test/cut_mesh_tests.jl +++ b/test/cut_mesh_tests.jl @@ -96,8 +96,8 @@ using StartUpDG: PathIntersections srd(e) # in-place application of SRD functor # test exactness - @test norm(e .- ecopy) < 1e3 * eps() - @test norm(u .- ucopy) < 1e3 * eps() + @test norm(e .- ecopy) < 1e2 * length(e) * eps() + @test norm(u .- ucopy) < 1e2 * length(u) * eps() end end From 26b985a3d72cfa1c34a62389aff8dd1765eef38a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 9 May 2024 23:17:48 -0500 Subject: [PATCH 10/10] update hybridized SBP test tolerance --- test/cut_mesh_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/cut_mesh_tests.jl b/test/cut_mesh_tests.jl index 2a4cdae7..93c8cf03 100644 --- a/test/cut_mesh_tests.jl +++ b/test/cut_mesh_tests.jl @@ -144,6 +144,6 @@ end Ph = hybridized_projection_matrices[e] Vh = hybridized_interp_matrices[e] @test norm(Ph * Qxh * Vh - Dx) + - norm(Ph * Qyh * Vh - Dy) < 1e3 * eps() + norm(Ph * Qyh * Vh - Dy) < 100 * length(Dx) * eps() end end \ No newline at end of file