InterfaceValues for DG interface integration

benchmark/benchmarks-assembly.jl

@@ -13,6 +13,7 @@ for spatial_dim ∈ 1:3
         COMMON_LOCAL_ASSEMBLY["spatial-dim",spatial_dim][string(geo_type)] = BenchmarkGroup()
 
         grid = generate_grid(geo_type, tuple(repeat([1], spatial_dim)...));
+        topology = ExclusiveTopology(grid)
         ip_geo = Ferrite.default_interpolation(geo_type)
         ref_type = FerriteBenchmarkHelper.getrefshape(geo_type)
 
@@ -66,6 +67,15 @@ for spatial_dim ∈ 1:3
 
                 LAGRANGE_SUITE["ritz-galerkin"]["face-flux"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_local_matrix($grid, $fsv, shape_gradient, shape_value, *)
                 LAGRANGE_SUITE["petrov-galerkin"]["face-flux"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $fsv, shape_gradient, $fsv2, shape_value, *)
+                ip = DiscontinuousLagrange{ref_type, order}()
+                isv = InterfaceValues(grid, qr_face, ip, ip_geo; geom_interpol_b = ip_geo);
+                isv2 = InterfaceValues(grid, qr_face, ip, ip_geo; geom_interpol_b = ip_geo);
+                LAGRANGE_SUITE["ritz-galerkin"]["interface-{grad}⋅[[val]]"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_local_matrix($grid, $topology, $isv, shape_gradient_average, shape_value_jump, ⋅)
+                LAGRANGE_SUITE["petrov-galerkin"]["interface-{grad}⋅[[val]]"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $topology, $isv, shape_gradient_average, $isv2, shape_value_jump, ⋅)
+
+                LAGRANGE_SUITE["ritz-galerkin"]["interface-interior-penalty"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_local_matrix($grid, $topology, $isv, shape_value_jump, shape_value_jump, ⋅)
+                LAGRANGE_SUITE["petrov-galerkin"]["interface-interior-penalty"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_local_matrix($grid, $topology, $isv, shape_value_jump, $isv2, shape_value_jump, ⋅)
+
             end
         end
     end

benchmark/helper.jl

@@ -65,6 +65,33 @@ function _generalized_ritz_galerkin_assemble_local_matrix(grid::Ferrite.Abstract
     f
 end
 
+function _generalized_ritz_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, topology::ExclusiveTopology, interfacevalues::InterfaceValues{<: Ferrite.InterpolationByDim{dim}}, f_shape, f_test, op) where {dim}
+    n_basefuncs = getnbasefunctions(interfacevalues)
+
+    Ke = zeros(n_basefuncs, n_basefuncs)
+
+    cell_a_coords = get_cell_coordinates(grid, 1)
+    for face in 1:nfaces(getcells(grid)[1])
+        face_a = FaceIndex(1, face)
+        for face_b in topology.face_neighbor[1, face]
+            cell_b_coords = get_cell_coordinates(grid, face_b[1])
+            reinit!(interfacevalues, face_a, face_b, cell_a_coords, cell_b_coords, grid)
+            for q_point in 1:getnquadpoints(interfacevalues)
+                dΓ = getdetJdV(interfacevalues, q_point)
+                for i in 1:n_basefuncs
+                    test = f_test(interfacevalues, q_point, i)
+                    for j in 1:n_basefuncs
+                        shape = f_shape(interfacevalues, q_point, j)
+                        Ke[i, j] += op(test, shape) * dΓ
+                    end
+                end
+            end
+        end
+    end
+
+    Ke
+end
+
 # Minimal Petrov-Galerkin type local assembly loop. We assume that both function spaces share the same integration rule. Test is applied from the left.
 function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, cellvalues_shape::CellValues{<: Ferrite.InterpolationByDim{dim}}, f_shape, cellvalues_test::CellValues{<: Ferrite.InterpolationByDim{dim}}, f_test, op) where {dim}
     n_basefuncs_shape = getnbasefunctions(cellvalues_shape)
@@ -122,6 +149,35 @@ function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.Abstra
     f
 end
 
+function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.AbstractGrid, topology::ExclusiveTopology,  interfacevalues_shape::InterfaceValues{<: Ferrite.InterpolationByDim{dim}}, f_shape, interfacevalues_test::InterfaceValues{<: Ferrite.InterpolationByDim{dim}}, f_test, op) where {dim}
+    n_basefuncs_shape = getnbasefunctions(interfacevalues_shape)
+    n_basefuncs_test = getnbasefunctions(interfacevalues_test)
+
+    Ke = zeros(n_basefuncs_test, n_basefuncs_shape)
+
+    cell_a_coords = get_cell_coordinates(grid, 1)
+    for face in 1:nfaces(getcells(grid)[1])
+        face_a = FaceIndex(1, face)
+        for face_b in topology.face_neighbor[1, face]
+            cell_b_coords = get_cell_coordinates(grid, face_b[1])
+            reinit!(interfacevalues_shape, face_a, face_b, cell_a_coords, cell_b_coords, grid)
+            reinit!(interfacevalues_test, face_a, face_b, cell_a_coords, cell_b_coords, grid)
+            for q_point in 1:getnquadpoints(interfacevalues_shape)
+                dΓ = getdetJdV(interfacevalues_test, q_point)
+                for i in 1:n_basefuncs_test
+                    test = f_test(interfacevalues_test, q_point, i)
+                    for j in 1:n_basefuncs_shape
+                        shape = f_shape(interfacevalues_shape, q_point, j)
+                        Ke[i, j] += op(test, shape) * dΓ
+                    end
+                end
+            end
+        end
+    end
+
+    Ke
+end
+
 function _assemble_mass(dh, cellvalues, sym)
     n_basefuncs = getnbasefunctions(cellvalues)
     Me = zeros(n_basefuncs, n_basefuncs)

docs/src/devdocs/elements.md

@@ -34,5 +34,7 @@ Ferrite.getvertexinstances
 Ferrite.filterfaces
 Ferrite.filteredges
 Ferrite.filtervertices
+Ferrite.transfer_point_cell_to_face
+Ferrite.transfer_point_face_to_cell
+Ferrite.InterfaceOrientationInfo
 ```
-

docs/src/reference/fevalues.md

@@ -35,3 +35,21 @@ FaceValues
 getcurrentface
 getnquadpoints(::FaceValues)
 ```
+
+## [InterfaceValues](@id reference-interfacevalues)
+
+All of the methods for [`FaceValues`](@ref) apply for `InterfaceValues` as well.
+In addition, there are some methods that are unique for `InterfaceValues`:
+
+```@docs
+InterfaceValues
+shape_value_average
+shape_value_jump
+shape_gradient_average
+shape_gradient_jump
+function_value_average
+function_value_jump
+function_gradient_average
+function_gradient_jump
+transform_interface_point
+```