From c9586ea0a0442028141a4c7cdf6dde646bdbf67b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 10 May 2024 01:16:33 -0500 Subject: [PATCH 1/8] diagm -> Diagonal --- src/cut_cell_meshes.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 77d8abc3..6224fac8 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -337,7 +337,7 @@ function construct_cut_volume_quadrature(N, cutcells, physical_frame_elements; # 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() + @assert norm(V' * Diagonal(w) * V - V' * Diagonal(w_pruned) * V) < 100 * eps() end @. xq_pruned[:, e] = xq[inds] @@ -941,9 +941,9 @@ function MeshData(rd::RefElemData, objects, 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 + M = Vq' * Diagonal(wJq.cut[:, e]) * Vq + Qr = Vq' * Diagonal(wJq.cut[:, e]) * Vrq + Qs = Vq' * Diagonal(wJq.cut[:, e]) * Vsq Dx_e, Dy_e = M \ Qr, M \ Qs xf_e = xf.cut[cut_face_node_indices_by_cell[e]] @@ -959,7 +959,7 @@ function MeshData(rd::RefElemData, objects, push!(face_interpolation_matrices, Vf) push!(differentiation_matrices, (Dx_e, Dy_e)) push!(mass_matrices, M) - push!(lift_matrices, M \ (Vf' * diagm(wf))) + push!(lift_matrices, M \ (Vf' * Diagonal(wf))) end cut_cell_operators = (; volume_interpolation_matrices, face_interpolation_matrices, From c9b94fd4df6a769a29e8daefaab0ba55883e29ea Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 10 May 2024 01:16:49 -0500 Subject: [PATCH 2/8] use views --- src/cut_cell_meshes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 6224fac8..75226fff 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -939,7 +939,7 @@ function MeshData(rd::RefElemData, objects, 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])) + basis(elem, rd.N, view(xq.cut, :, e), view(yq.cut, :, e))) M = Vq' * Diagonal(wJq.cut[:, e]) * Vq Qr = Vq' * Diagonal(wJq.cut[:, e]) * Vrq From 97e81b5862e02f79a460ed3f4189f8dec013431a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 10 May 2024 01:17:02 -0500 Subject: [PATCH 3/8] preallocate sorting permutation vector p --- src/cut_cell_meshes.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 75226fff..a378e8cf 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -883,15 +883,17 @@ function MeshData(rd::RefElemData, objects, cut = collect(eachindex(xf.cut) .+ length(mapM_cartesian))) mapP = copy(mapM) + # allocate a vector for points on non-boundary faces + p = zeros(Int, length(first(quad_rule_face))) 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]]] + idM = view(mapM, face_node_indices_by_face[f]) + idP = view(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 + StartUpDG.match_coordinate_vectors!(p, xyzM, xyzP) + @. mapP[idM[p]] = idP end end mapB = findall(vec(mapM) .== vec(mapP)) # determine boundary nodes From d2d26c92c3e1417b68ce0a18f0425949d81a56bc Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 10 May 2024 01:54:11 -0500 Subject: [PATCH 4/8] improve efficiency of caratheodory --- src/physical_frame_basis.jl | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/physical_frame_basis.jl b/src/physical_frame_basis.jl index 55af85a0..dcff6a46 100644 --- a/src/physical_frame_basis.jl +++ b/src/physical_frame_basis.jl @@ -176,25 +176,23 @@ function caratheodory_pruning_qr(V, w_in) inds = collect(1:M) m = M-N Q, _ = qr(V) - Q = copy(Q) for _ in 1:m - kvec = Q[:,end] + kvec = view(Q, :, size(Q, 2)) # for subtracting the kernel vector idp = findall(@. kvec > 0) - alphap, k0p = findmin(w[inds[idp]] ./ kvec[idp]) + alphap, k0p = findmin(view(w, inds[idp]) ./ view(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]; + idn = findall(@. kvec < 0) + alphan, k0n = findmax(view(w, inds[idn]) ./ view(kvec, idn)) + k0n = idn[k0n] alpha, k0 = abs(alphan) < abs(alphap) ? (alphan, k0n) : (alphap, k0p) - w[inds] = w[inds] - alpha * kvec + @. w[inds] = w[inds] - alpha * kvec deleteat!(inds, k0) Q, _ = qr(V[inds, :]) - Q = copy(Q) end return w, inds end \ No newline at end of file From 3e3ff7def1fe4a4d3fae05cca1e850dd251b7671 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 10 May 2024 02:19:03 -0500 Subject: [PATCH 5/8] remove type instability for PhysicalFrame basis --- src/cut_cell_meshes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index a378e8cf..6a54febd 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -506,7 +506,7 @@ 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 + physical_frame_elements = typeof(PhysicalFrame())[] # 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]) From 8113792fadbd5a16b173f1870e594bc8f95c8e36 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 10 May 2024 14:48:09 -0500 Subject: [PATCH 6/8] bump compat for PathIntersections - should be more type stable --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 37412a69..37b17f3f 100644 --- a/Project.toml +++ b/Project.toml @@ -35,7 +35,7 @@ FillArrays = "0.13, 1" HDF5 = "0.16, 0.17" Kronecker = "0.5" NodesAndModes = "1" -PathIntersections = "0.1" +PathIntersections = "0.2" RecipesBase = "1" RecursiveArrayTools = "3.4" Reexport = "1" From 2a857d9d997f4b08451427b9412ee6d67a884c58 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 10 May 2024 16:48:57 -0500 Subject: [PATCH 7/8] fix compat for PathIntersections --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 37b17f3f..05d9cb6f 100644 --- a/Project.toml +++ b/Project.toml @@ -35,7 +35,7 @@ FillArrays = "0.13, 1" HDF5 = "0.16, 0.17" Kronecker = "0.5" NodesAndModes = "1" -PathIntersections = "0.2" +PathIntersections = "0.1, 0.2" RecipesBase = "1" RecursiveArrayTools = "3.4" Reexport = "1" From 3f717ba7dbf97db72da620d1754d54d2215136cd Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Thu, 16 May 2024 10:50:13 -0500 Subject: [PATCH 8/8] add docstring --- src/cut_cell_meshes.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/cut_cell_meshes.jl b/src/cut_cell_meshes.jl index 6a54febd..c1420a33 100644 --- a/src/cut_cell_meshes.jl +++ b/src/cut_cell_meshes.jl @@ -1009,6 +1009,17 @@ function MeshData(rd::RefElemData, objects, end +""" + hybridized_SBP_operators(md::MeshData{2, <:CutCellMesh}) + +This constructs hybridized SBP operators using the approach taken in Chan (2019), +"Skew-Symmetric Entropy Stable Modal Discontinuous Galerkin Formulations". +[https://doi.org/10.1007/s10915-019-01026-w](https://doi.org/10.1007/s10915-019-01026-w) + +This function returns `hybridized_operators::Vector{Tuple{<:Matrix, <:Matrix}}` and +`project_and_interp_operators, projection_operators, interpolation_operators`, which are +all `Vector{<:Matrix}`, where each entry corresponds to a cut element. +""" function hybridized_SBP_operators(md::MeshData{2, <:CutCellMesh}) mt = md.mesh_type (; volume_interpolation_matrices,