InterfaceValues for DG interface integration

benchmark/benchmarks-assembly.jl

@@ -12,7 +12,8 @@ for spatial_dim ∈ 1:3
     for geo_type ∈ FerriteBenchmarkHelper.geo_types_for_spatial_dim(spatial_dim)
         COMMON_LOCAL_ASSEMBLY["spatial-dim",spatial_dim][string(geo_type)] = BenchmarkGroup()
 
-        grid = generate_grid(geo_type, tuple(repeat([1], spatial_dim)...));
+        grid = generate_grid(geo_type, tuple(repeat([2], spatial_dim)...));
+        topology = ExclusiveTopology(grid)
         ip_geo = Ferrite.default_interpolation(geo_type)
         ref_type = FerriteBenchmarkHelper.getrefshape(geo_type)
 
@@ -66,6 +67,20 @@ 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(qr_face, ip, ip_geo; geom_interpol_b = ip_geo);
+                isv2 = InterfaceValues(qr_face, ip, ip_geo; geom_interpol_b = ip_geo);
+                dh = DofHandler(grid)
+                add!(dh, :u, ip)
+                close!(dh)
+
+                LAGRANGE_SUITE["ritz-galerkin"]["interface-{grad}⋅[[val]]"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_interfaces($dh, $isv, shape_gradient_average, shape_value_jump, ⋅)
+                LAGRANGE_SUITE["petrov-galerkin"]["interface-{grad}⋅[[val]]"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_interfaces($dh, $isv, shape_gradient_average, $isv2, shape_value_jump, ⋅)
+
+                LAGRANGE_SUITE["ritz-galerkin"]["interface-interior-penalty"] = @benchmarkable FerriteAssemblyHelper._generalized_ritz_galerkin_assemble_interfaces($dh, $isv, shape_value_jump, shape_value_jump, ⋅)
+                LAGRANGE_SUITE["petrov-galerkin"]["interface-interior-penalty"] = @benchmarkable FerriteAssemblyHelper._generalized_petrov_galerkin_assemble_interfaces($dh, $isv, shape_value_jump, $isv2, shape_value_jump, ⋅)
+
             end
         end
     end

benchmark/helper.jl

@@ -65,6 +65,30 @@ function _generalized_ritz_galerkin_assemble_local_matrix(grid::Ferrite.Abstract
     f
 end
 
+function _generalized_ritz_galerkin_assemble_interfaces(dh::Ferrite.AbstractDofHandler, interfacevalues::InterfaceValues{<: FaceValues{<: Ferrite.InterpolationByDim{dim}}}, f_shape, f_test, op) where {dim}
+    n_basefuncs = getnbasefunctions(interfacevalues)
+
+    K = zeros(ndofs(dh), ndofs(dh))
+
+    for ic in InterfaceIterator(dh)
+        reinit!(interfacevalues, ic)
+        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)
+                f_test == shape_value_jump && (test *= getnormal(interfacevalues, q_point))
+                for j in 1:n_basefuncs
+                    shape = f_shape(interfacevalues, q_point, j)
+                    f_shape == shape_value_jump && (shape *= getnormal(interfacevalues, q_point))
+                    K[interfacedofs(ic)[i], interfacedofs(ic)[j]] += op(test, shape) * dΓ
+                end
+            end
+        end
+    end
+
+    K
+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 +146,32 @@ function _generalized_petrov_galerkin_assemble_local_matrix(grid::Ferrite.Abstra
     f
 end
 
+function _generalized_petrov_galerkin_assemble_interfaces(dh::Ferrite.AbstractDofHandler, interfacevalues_shape::InterfaceValues{<: FaceValues{<: Ferrite.InterpolationByDim{dim}}}, f_shape, interfacevalues_test::InterfaceValues{<: FaceValues{<: Ferrite.InterpolationByDim{dim}}}, f_test, op) where {dim}
+    n_basefuncs_shape = getnbasefunctions(interfacevalues_shape)
+    n_basefuncs_test = getnbasefunctions(interfacevalues_test)
+
+    K = zeros(ndofs(dh), ndofs(dh))
+
+    for ic in InterfaceIterator(dh)
+        reinit!(interfacevalues_shape, ic)
+        reinit!(interfacevalues_test, ic)
+        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)
+                f_test == shape_value_jump && (test *= getnormal(interfacevalues_test, q_point))
+                for j in 1:n_basefuncs_shape
+                    shape = f_shape(interfacevalues_shape, q_point, j)
+                    f_shape == shape_value_jump && (shape *= getnormal(interfacevalues_shape, q_point))
+                    K[interfacedofs(ic)[i], interfacedofs(ic)[j]] += op(test, shape) * dΓ
+                end
+            end
+        end
+    end
+
+    K
+end
+
 function _assemble_mass(dh, cellvalues, sym)
     n_basefuncs = getnbasefunctions(cellvalues)
     Me = zeros(n_basefuncs, n_basefuncs)

docs/src/devdocs/elements.md

@@ -35,5 +35,7 @@ Ferrite.getvertexinstances
 Ferrite.filterfaces
 Ferrite.filteredges
 Ferrite.filtervertices
+Ferrite.element_to_face_transformation
+Ferrite.face_to_element_transformation
+Ferrite.InterfaceTransformation
 ```
-

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_points!
+```

src/FEValues/common_values.jl

@@ -57,7 +57,7 @@ getnquadpoints(fe::FaceValues) = getnquadpoints(fe.qr, fe.current_face[])
     getdetJdV(fe_v::AbstractValues, q_point::Int)
 
 Return the product between the determinant of the Jacobian and the quadrature
-point weight for the given quadrature point: ``\\det(J(\\mathbf{x})) w_q``
+point weight for the given quadrature point: ``\\det(J(\\mathbf{x})) w_q``.
 
 This value is typically used when one integrates a function on a
 finite element cell or face as
@@ -78,6 +78,12 @@ quadrature point `q_point`.
 @propagate_inbounds shape_value(cv::CellValues, q_point::Int, base_func::Int) = cv.N[base_func, q_point]
 @propagate_inbounds shape_value(bv::FaceValues, q_point::Int, base_func::Int) = bv.N[base_func, q_point, bv.current_face[]]
 
+"""
+    geometric_value(fe_v::AbstractValues, q_point::Int, base_function::Int)
+
+Return the value of shape function `base_function` used for geometric interpolation evaluated in
+quadrature point `q_point`.
+"""
 @propagate_inbounds geometric_value(cv::CellValues, q_point::Int, base_func::Int) = cv.M[base_func, q_point]
 @propagate_inbounds geometric_value(bv::FaceValues, q_point::Int, base_func::Int) = bv.M[base_func, q_point, bv.current_face[]]
 
@@ -174,13 +180,14 @@ function function_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector
 end
 
 # TODO: Deprecate this, nobody is using this in practice...
-function function_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector{<:Vec})
+function function_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector{<:Vec}, dof_range = eachindex(u))
     n_base_funcs = getnbasefunctions(fe_v)
-    length(u) == n_base_funcs || throw_incompatible_dof_length(length(u), n_base_funcs)
+    length(dof_range) == n_base_funcs || throw_incompatible_dof_length(length(dof_range), n_base_funcs)
+    @boundscheck checkbounds(u, dof_range)
     @boundscheck checkquadpoint(fe_v, q_point)
     grad = function_gradient_init(fe_v, u)
-    @inbounds for i in 1:n_base_funcs
-        grad += u[i] ⊗ shape_gradient(fe_v, q_point, i)
+    @inbounds for (i, j) in pairs(dof_range)
+        grad += u[j] ⊗ shape_gradient(fe_v, q_point, i)
     end
     return grad
 end
@@ -206,13 +213,7 @@ The symmetric gradient of a scalar function is computed as
 where ``\\mathbf{u}_i`` are the nodal values of the function.
 """
 function function_symmetric_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector, dof_range = eachindex(u))
-    grad = function_gradient(fe_v, q_point, u, dof_range)
-    return symmetric(grad)
-end
-
-# TODO: Deprecate this, nobody is using this in practice...
-function function_symmetric_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector{<:Vec})
-    grad = function_gradient(fe_v, q_point, u)
+    grad = dof_range == eachindex(u) ? function_gradient(fe_v, q_point, u) : function_gradient(fe_v, q_point, u, dof_range) # Workaround for calling deprecated method
     return symmetric(grad)
 end
 
@@ -241,6 +242,15 @@ function function_divergence(fe_v::AbstractValues, q_point::Int, u::AbstractVect
     return diverg
 end
 
+"""
+    function_curl(fe_v::AbstractValues, q_point::Int, u::AbstractVector)
+
+Compute the curl of the vector valued function in a quadrature point.
+
+The curl of a vector valued functions in the quadrature point ``\\mathbf{x}_q)`` is computed as
+``\\mathbf{\\nabla} \\times \\mathbf{u}(\\mathbf{x_q}) = \\sum\\limits_{i = 1}^n \\mathbf{\\nabla} N_i \\times (\\mathbf{x_q}) \\cdot \\mathbf{u}_i``
+where ``\\mathbf{u}_i`` are the nodal values of the function.
+"""
 function_curl(fe_v::AbstractValues, q_point::Int, u::AbstractVector, dof_range = eachindex(u)) =
     curl_from_gradient(function_gradient(fe_v, q_point, u, dof_range))
 
@@ -257,11 +267,11 @@ coordinates of the cell.
 The coordinate is computed, using the geometric interpolation, as
 ``\\mathbf{x} = \\sum\\limits_{i = 1}^n M_i (\\mathbf{x}) \\mathbf{\\hat{x}}_i``
 """
-function spatial_coordinate(fe_v::AbstractValues, q_point::Int, x::AbstractVector{Vec{dim,T}}) where {dim,T}
+function spatial_coordinate(fe_v::AbstractValues, q_point::Int, x::AbstractVector{<:Vec})
     n_base_funcs = getngeobasefunctions(fe_v)
     length(x) == n_base_funcs || throw_incompatible_coord_length(length(x), n_base_funcs)
     @boundscheck checkquadpoint(fe_v, q_point)
-    vec = zero(Vec{dim,T})
+    vec = zero(eltype(x))
     @inbounds for i in 1:n_base_funcs
         vec += geometric_value(fe_v, q_point, i) * x[i]
     end
@@ -273,7 +283,7 @@ function Base.show(io::IO, ::MIME"text/plain", fe_v::AbstractValues)
 end
 
 # copy
-for ValueType in (CellValues, FaceValues)
+for ValueType in (CellValues, FaceValues, InterfaceValues)
     args = [:(copy(cv.$fname)) for fname in fieldnames(ValueType)]
     @eval begin
         function Base.copy(cv::$ValueType)

src/FEValues/face_integrals.jl

@@ -1,11 +1,20 @@
 """
     face_to_element_transformation(point::Vec, ::Type{<:AbstractRefShape}, face::Int)
 
-Transform quadrature point from face's reference (N-1)D coordinates to ND coordinates on the
-cell's face.
+Transform quadrature point from face's reference coordinates to coordinates on the
+cell's face, increasing the number of dimensions by one.
 """
 face_to_element_transformation
 
+"""
+    element_to_face_transformation(point::AbstractVector, cell::AbstractCell{AbstractRefShape}, face::Int)
+
+Transform quadrature point from cell's coordinates to the face's reference coordinates, decreasing the number of dimensions by one.
+This is the inverse of `face_to_element_transformation`.
+This is the inverse of `face_to_element_transformation`.
+"""
+element_to_face_transformation
+
 """
     weighted_normal(J::AbstractTensor, fv::FaceValues, face::Int)
     weighted_normal(J::AbstractTensor, ::Type{<:AbstractRefShape}, face::Int)
@@ -62,12 +71,20 @@ end
 ##################
 
 # Mapping from to 0D node to 1D line vertex.
-function face_to_element_transformation(::Vec{0, T}, ::Type{RefLine}, face::Int) where {T}
+function face_to_element_transformation(::Union{Vec{0, T},Vec{1, T}}, ::Type{RefLine}, face::Int) where {T}
     face == 1 && return Vec{1, T}(( -one(T),))
     face == 2 && return Vec{1, T}((  one(T),))
     throw(ArgumentError("unknown face number"))
 end
 
+# Mapping from 1D line to point.
+function element_to_face_transformation(point::Vec{1, T}, ::Type{RefLine}, face::Int) where T
+    x = point[]
+    face == 1 && return Vec(-x)
+    face == 2 && return Vec( x)
+    throw(ArgumentError("unknown face number"))
+end
+
 function weighted_normal(::Tensor{2,1,T}, ::Type{RefLine}, face::Int) where {T}
     face == 1 && return Vec{1,T}((-one(T),))
     face == 2 && return Vec{1,T}(( one(T),))
@@ -88,6 +105,16 @@ function face_to_element_transformation(point::Vec{1, T}, ::Type{RefQuadrilatera
     throw(ArgumentError("unknown face number"))
 end
 
+# Mapping from 2D face of a quadrilateral to 1D line.
+function element_to_face_transformation(point::Vec{2, T}, ::Type{RefQuadrilateral}, face::Int) where T
+    x, y = point
+    face == 1 && return Vec( x)
+    face == 2 && return Vec( y)
+    face == 3 && return Vec( -x)
+    face == 4 && return Vec( -y)
+    throw(ArgumentError("unknown face number"))
+end
+
 function weighted_normal(J::Tensor{2,2}, ::Type{RefQuadrilateral}, face::Int)
     @inbounds begin
         face == 1 && return Vec{2}(( J[2,1], -J[1,1]))
@@ -111,6 +138,15 @@ function face_to_element_transformation(point::Vec{1, T},  ::Type{RefTriangle},
     throw(ArgumentError("unknown face number"))
 end
 
+# Mapping from 2D face of a triangle to 1D line.
+function element_to_face_transformation(point::Vec{2, T}, ::Type{RefTriangle}, face::Int) where T
+    x, y = point
+    face == 1 && return Vec( one(T) - x * 2)
+    face == 2 && return Vec( one(T) - y * 2 )
+    face == 3 && return Vec( x * 2 - one(T))
+    throw(ArgumentError("unknown face number"))
+end
+
 function weighted_normal(J::Tensor{2,2}, ::Type{RefTriangle}, face::Int)
     @inbounds begin
         face == 1 && return Vec{2}((-(J[2,1] - J[2,2]), J[1,1] - J[1,2]))
@@ -136,6 +172,18 @@ function face_to_element_transformation(point::Vec{2, T}, ::Type{RefHexahedron},
     throw(ArgumentError("unknown face number"))
 end
 
+# Mapping from 3D face of a hexahedron to 2D quadrilateral.
+function element_to_face_transformation(point::Vec{3, T}, ::Type{RefHexahedron}, face::Int) where T
+    x, y, z = point
+    face == 1 && return Vec( y, x)
+    face == 2 && return Vec( x, z)
+    face == 3 && return Vec( y,  z)
+    face == 4 && return Vec( -x,  z)
+    face == 5 && return Vec( z,  y)
+    face == 6 && return Vec( x,  y)
+    throw(ArgumentError("unknown face number"))
+end
+
 function weighted_normal(J::Tensor{2,3}, ::Type{RefHexahedron}, face::Int)
     @inbounds begin
         face == 1 && return J[:,2] × J[:,1]
@@ -162,6 +210,16 @@ function face_to_element_transformation(point::Vec{2, T}, ::Type{RefTetrahedron}
     throw(ArgumentError("unknown face number"))
 end
 
+# Mapping from 3D face of a tetrahedon to 2D triangle.
+function element_to_face_transformation(point::Vec{3, T}, ::Type{RefTetrahedron}, face::Int) where T
+    x, y, z = point
+    face == 1 && return Vec( one(T)-x-y,  y)
+    face == 2 && return Vec( one(T)-z-x,  x)
+    face == 3 && return Vec( x,  y)
+    face == 4 && return Vec( one(T)-y-z,  z)
+    throw(ArgumentError("unknown face number"))
+end
+
 function weighted_normal(J::Tensor{2,3}, ::Type{RefTetrahedron}, face::Int)
     @inbounds begin
         face == 1 && return J[:,2] × J[:,1]
@@ -188,6 +246,17 @@ function face_to_element_transformation(point::Vec{2, T}, ::Type{RefPrism}, face
     throw(ArgumentError("unknown face number"))
 end
 
+# Mapping from 3D face of a wedge to 2D triangle or 2D quadrilateral.
+function element_to_face_transformation(point::Vec{3, T}, ::Type{RefPrism}, face::Int) where T
+    x, y, z = point
+    face == 1 && return Vec( one(T)-x-y,  y)
+    face == 2 && return Vec( 2*x - one(T), 2*z - one(T) )
+    face == 3 && return Vec( 2*(one(T) - y) - one(T), 2*z - one(T) )
+    face == 4 && return Vec( 2*y - one(T), 2*z - one(T) )
+    face == 5 && return Vec( one(T) - x - y,  x)
+    throw(ArgumentError("unknown face number"))
+end
+
 function weighted_normal(J::Tensor{2,3}, ::Type{RefPrism}, face::Int)
     @inbounds begin
         face == 1 && return J[:,2] × J[:,1]
@@ -214,6 +283,17 @@ function face_to_element_transformation(point::Vec{2, T}, ::Type{RefPyramid}, fa
     throw(ArgumentError("unknown face number"))
 end
 
+# Mapping from 3D face of a pyramid to 2D triangle or 2D quadrilateral.
+function element_to_face_transformation(point::Vec{3, T}, ::Type{RefPyramid}, face::Int) where T
+    x, y, z = point
+    face == 1 && return Vec( 2*y - one(T),  2*x - one(T))
+    face == 2 && return Vec( one(T) - z - x, x)
+    face == 3 && return Vec( one(T) - y - z, z)
+    face == 4 && return Vec( x - y, y)
+    face == 5 && return Vec( one(T) - x - z,  z)
+    throw(ArgumentError("unknown face number"))
+end
+
 function weighted_normal(J::Tensor{2,3}, ::Type{RefPyramid}, face::Int)
     @inbounds begin
         face == 1 && return J[:,2] × J[:,1]