InterfaceValues for DG interface integration

Project.toml

@@ -28,7 +28,7 @@ Metis = "1.3"
 NearestNeighbors = "0.4"
 Preferences = "1"
 Reexport = "1"
-Tensors = "1"
+Tensors = "1.14"
 WriteVTK = "1.13"
 julia = "1.6"
 

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/devdocs/interpolations.md

@@ -37,7 +37,6 @@ Ferrite.edgedof_interior_indices(::Interpolation)
 Ferrite.celldof_interior_indices(::Interpolation)
 Ferrite.getnbasefunctions(::Interpolation)
 Ferrite.reference_coordinates(::Interpolation)
-Ferrite.face_to_element_transformation
 Ferrite.is_discontinuous(::Interpolation)
 Ferrite.adjust_dofs_during_distribution(::Interpolation)
 ```

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[]]
 
@@ -116,7 +122,7 @@ end
 curl_from_gradient(∇v::SecondOrderTensor{3}) = Vec{3}((∇v[3,2] - ∇v[2,3], ∇v[1,3] - ∇v[3,1], ∇v[2,1] - ∇v[1,2]))
 
 """
-    function_value(fe_v::AbstractValues, q_point::Int, u::AbstractVector)
+    function_value(fe_v::AbstractValues, q_point::Int, u::AbstractVector, [dof_range])
 
 Compute the value of the function in a quadrature point. `u` is a vector with values
 for the degrees of freedom. For a scalar valued function, `u` contains scalars.
@@ -146,7 +152,7 @@ shape_value_type(::Union{CellValues{<:Any, N_t}, FaceValues{<:Any, N_t}}) where
 function_value_init(cv::AbstractValues, ::AbstractVector{T}) where {T} = zero(shape_value_type(cv)) * zero(T)
 
 """
-    function_gradient(fe_v::AbstractValues{dim}, q_point::Int, u::AbstractVector)
+    function_gradient(fe_v::AbstractValues{dim}, q_point::Int, u::AbstractVector, [dof_range])
 
 Compute the gradient of the function in a quadrature point. `u` is a vector with values
 for the degrees of freedom. For a scalar valued function, `u` contains scalars.
@@ -196,7 +202,7 @@ function function_gradient_init(cv::AbstractValues, ::AbstractVector{T}) where {
 end
 
 """
-    function_symmetric_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector)
+    function_symmetric_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector, [dof_range])
 
 Compute the symmetric gradient of the function, see [`function_gradient`](@ref).
 Return a `SymmetricTensor`.
@@ -205,19 +211,17 @@ The symmetric gradient of a scalar function is computed as
 ``\\left[ \\mathbf{\\nabla}  \\mathbf{u}(\\mathbf{x_q}) \\right]^\\text{sym} =  \\sum\\limits_{i = 1}^n  \\frac{1}{2} \\left[ \\mathbf{\\nabla} N_i (\\mathbf{x}_q) \\otimes \\mathbf{u}_i + \\mathbf{u}_i  \\otimes  \\mathbf{\\nabla} N_i (\\mathbf{x}_q) \\right]``
 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))
+function function_symmetric_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector, dof_range)
     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})
+function function_symmetric_gradient(fe_v::AbstractValues, q_point::Int, u::AbstractVector)
     grad = function_gradient(fe_v, q_point, u)
     return symmetric(grad)
 end
-
 """
-    function_divergence(fe_v::AbstractValues, q_point::Int, u::AbstractVector)
+    function_divergence(fe_v::AbstractValues, q_point::Int, u::AbstractVector, [dof_range])
 
 Compute the divergence of the vector valued function in a quadrature point.
 
@@ -241,6 +245,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, [dof_range])
+
+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 +270,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 +286,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)