Improve efficiency of cut-cell MeshData (#167)

* diagm -> Diagonal

* use views

* preallocate sorting permutation vector p

* improve efficiency of caratheodory

* remove type instability for PhysicalFrame basis

* bump compat for PathIntersections - should be more type stable

* fix compat for PathIntersections

* add docstring

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.1, 0.2"
 RecipesBase = "1"
 RecursiveArrayTools = "3.4"
 Reexport = "1"

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]
@@ -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])
@@ -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
@@ -939,11 +941,11 @@ 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' * 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 +961,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, 
@@ -1007,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, 

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