jump direction consistent with InterfaceValues

src/cellvalues.jl

@@ -101,15 +101,15 @@ end
 get_side_and_baseindex(cv::InterfaceCellValues, i::Integer) = cv.sides_and_baseindices[i]
 
 """
-    get_base_value(get_value::Function, cv::InterfaceCellValues, qp::Int, i::Int, here::Bool)
+    get_base_value(get_value::Function, cv::InterfaceCellValues, qp::Int, i::Int; here::Bool)
 
 Return a value from an `::InterfaceCellValues` by specifing:
 - `get_value`: function specifing which kind of value, e.g. `shape_value`
 - `qp`: index of the quadrature point
 - `i`: index of the base function
 - `here`: side of the interface, where `true` means "here" and `false` means "there".
 """
-function get_base_value(get_value::Function, cv::InterfaceCellValues, qp::Int, i::Int, here::Bool)
+function get_base_value(get_value::Function, cv::InterfaceCellValues, qp::Int, i::Int; here::Bool)
     side, baseindex = cv.sides_and_baseindices[i]
     if side == :here && here
         return get_value(cv.here, qp, baseindex)
@@ -120,27 +120,27 @@ function get_base_value(get_value::Function, cv::InterfaceCellValues, qp::Int, i
 end
 
 """
-    shape_value(cv::InterfaceCellValues, qp::Int, i::Int, here::Bool)
+    shape_value(cv::InterfaceCellValues, qp::Int, i::Int; here::Bool)
 
 Return the value of shape function `i` evaluated in quadrature point `qp`
 on side `here`, where `true` means "here" and `false` means "there".
 """
-function Ferrite.shape_value(cv::InterfaceCellValues, qp::Int, i::Int, here::Bool)
-    val = get_base_value(shape_value, cv, qp, i, here)
+function Ferrite.shape_value(cv::InterfaceCellValues, qp::Int, i::Int; here::Bool)
+    val = get_base_value(shape_value, cv, qp, i; here=here)
     if isnothing(val)
         return zero(shape_value_type(cv))
     end
     return val
 end
 
 """
-    shape_gradient(cv::InterfaceCellValues, qp::Int, i::Int, here::Bool)
+    shape_gradient(cv::InterfaceCellValues, qp::Int, i::Int; here::Bool)
 
 Return the gradient of shape function `i` evaluated in quadrature point `qp`
 on side `here`, where `true` means "here" and `false` means "there".
 """
-function Ferrite.shape_gradient(cv::InterfaceCellValues, qp::Int, i::Int, here::Bool)
-    grad = get_base_value(shape_gradient, cv, qp, i, here)
+function Ferrite.shape_gradient(cv::InterfaceCellValues, qp::Int, i::Int; here::Bool)
+    grad = get_base_value(shape_gradient, cv, qp, i; here=here)
     if isnothing(grad)
         return zero(shape_gradient_type(cv))
     end
@@ -177,7 +177,7 @@ for computing the value jump on an interface.
 """
 function Ferrite.shape_value_jump(cv::InterfaceCellValues, qp::Int, i::Int)
     side, baseindex = cv.sides_and_baseindices[i]
-    return side == :here ? -shape_value(cv.here, qp, baseindex) : shape_value(cv.there, qp, baseindex)
+    return side == :here ? shape_value(cv.here, qp, baseindex) : -shape_value(cv.there, qp, baseindex)
 end
 
 """
@@ -188,43 +188,43 @@ for computing the gradient jump on an interface.
 """
 function Ferrite.shape_gradient_jump(cv::InterfaceCellValues, qp::Int, i::Int)
     side, baseindex = cv.sides_and_baseindices[i]
-    return side == :here ? -shape_gradient(cv.here, qp, baseindex) : shape_gradient(cv.there, qp, baseindex)
+    return side == :here ? shape_gradient(cv.here, qp, baseindex) : -shape_gradient(cv.there, qp, baseindex)
 end
 
 """
-    function_value(cv::InterfaceCellValues, qp::Int, u::AbstractVector, here::Bool)
+    function_value(cv::InterfaceCellValues, qp::Int, u::AbstractVector; here::Bool)
 
 Compute the value of the function in a quadrature point on side `here`,
 where `true` means "here" and `false` means "there".
 `u` is a vector with values for the degrees of freedom.
 """
-function Ferrite.function_value(cv::InterfaceCellValues, qp::Int, u::AbstractVector, here::Bool, dof_range = eachindex(u))
+function Ferrite.function_value(cv::InterfaceCellValues, qp::Int, u::AbstractVector, dof_range=eachindex(u); here::Bool)
     nbf = getnbasefunctions(cv)
     length(dof_range) == nbf || throw_incompatible_dof_length(length(dof_range), nbf)
     @boundscheck checkbounds(u, dof_range)
     @boundscheck checkquadpoint(cv, qp)
     val = function_value_init(cv, u)
     @inbounds for (i, j) in pairs(dof_range)
-        val += shape_value(cv, qp, i, here) * u[j]
+        val += shape_value(cv, qp, i; here=here) * u[j]
     end
     return val
 end
 
 """
-    function_gradient(cv::InterfaceCellValues, qp::Int, u::AbstractVector, here::Bool)
+    function_gradient(cv::InterfaceCellValues, qp::Int, u::AbstractVector; here::Bool)
 
 Compute the gradient of the function in a quadrature point on side `here`,
 where `true` means "here" and `false` means "there".
 `u` is a vector with values for the degrees of freedom.
 """
-function Ferrite.function_gradient(cv::InterfaceCellValues, qp::Int, u::AbstractVector, here::Bool, dof_range = eachindex(u))
+function Ferrite.function_gradient(cv::InterfaceCellValues, qp::Int, u::AbstractVector, dof_range=eachindex(u); here::Bool)
     nbf = getnbasefunctions(cv)
     length(dof_range) == nbf || throw_incompatible_dof_length(length(dof_range), nbf)
     @boundscheck checkbounds(u, dof_range)
     @boundscheck checkquadpoint(cv, qp)
     grad = function_gradient_init(cv, u)
     @inbounds for (i, j) in pairs(dof_range)
-        grad += shape_gradient(cv, qp, i, here) * u[j]
+        grad += shape_gradient(cv, qp, i; here=here) * u[j]
     end
     return grad
 end
@@ -235,7 +235,7 @@ end
 Compute the average value of the function in a quadrature point.
 """
 function Ferrite.function_value_average(cv::InterfaceCellValues, qp::Int, u::AbstractVector)
-    return (function_value(cv, qp, u, true) + function_value(cv, qp, u, false))/2
+    return (function_value(cv, qp, u; here=true) + function_value(cv, qp, u; here=false))/2
 end
 
 """
@@ -244,7 +244,7 @@ end
 Compute the average gradient of the function in a quadrature point.
 """
 function Ferrite.function_gradient_average(cv::InterfaceCellValues, qp::Int, u::AbstractVector)
-    return (function_gradient(cv, qp, u, true) + function_gradient(cv, qp, u, false)) / 2
+    return (function_gradient(cv, qp, u; here=true) + function_gradient(cv, qp, u; here=false)) / 2
 end
 
 """
@@ -265,7 +265,7 @@ end
 Compute the jump of the function value in a quadrature point.
 """
 function Ferrite.function_value_jump(cv::InterfaceCellValues, qp::Int, u::AbstractVector)
-    return function_value(cv, qp, u, false) - function_value(cv, qp, u, true)
+    return function_value(cv, qp, u; here=true) - function_value(cv, qp, u; here=false)
 end
 
 """
@@ -274,7 +274,7 @@ end
 Compute the jump of the function gradient in a quadrature point.
 """
 function Ferrite.function_gradient_jump(cv::InterfaceCellValues, qp::Int, u::AbstractVector)
-    return function_gradient(cv, qp, u, false) - function_gradient(cv, qp, u, true)
+    return function_gradient(cv, qp, u; here=true) - function_gradient(cv, qp, u; here=false)
 end
 
 function Base.show(io::IO, d::MIME"text/plain", cv::InterfaceCellValues)

test/test_cellvalues.jl

@@ -26,14 +26,45 @@
     here, there = rand(2)
     u = vcat(ones(nbf÷2).*here, ones(nbf÷2).*there)
     for qp in 1:getnquadpoints(cv)
-        @test function_value(cv, qp, u, true) ≈ here
-        @test function_value(cv, qp, u, false) ≈ there
-        @test all(abs.(function_gradient(cv, qp, u, true)) .≤ 1e-14)
-        @test all(abs.(function_gradient(cv, qp, u, false)) .≤ 1e-14)
+        @test function_value(cv, qp, u; here=true) ≈ here
+        @test function_value(cv, qp, u; here=false) ≈ there
+        @test all(abs.(function_gradient(cv, qp, u; here=true)) .≤ 1e-14)
+        @test all(abs.(function_gradient(cv, qp, u; here=false)) .≤ 1e-14)
         @test function_value_average(cv, qp, u) ≈ (here + there)/2
-        @test function_value_jump(cv, qp, u) ≈ there - here
+        @test function_value_jump(cv, qp, u) ≈ here - there
         @test all(abs.(function_gradient_average(cv, qp, u)) .≤ 1e-14)
         @test all(abs.(function_gradient_jump(cv, qp, u)) .≤ 1e-14)
         @test getdetJdV_average(cv, qp) == (getdetJdV(cv.here, qp) + getdetJdV(cv.there, qp)) / 2
     end
+end
+
+@testset "Consistency with InterfaceValues" begin
+    qr = QuadratureRule{RefTriangle}(1)
+    ip = InterfaceCellInterpolation(Lagrange{RefTriangle, 1}())
+    cv = InterfaceCellValues(qr, ip)
+
+    ip = DiscontinuousLagrange{RefTetrahedron, 1}()
+    qr = FacetQuadratureRule{RefTetrahedron}(1)
+    iv = InterfaceValues(qr, ip)
+
+    x₁ = [ Vec((0.0, 0.0, 0.0)), Vec((rand(), 0.0, 0.0)), Vec((0.0, rand(), 0.0)), Vec((0.0, 0.0, rand())) ]
+    x₂ = [ Vec((1.0, 1.0, 1.0)), x₁[2], x₁[3], x₁[4] ]
+    xᵢ = repeat([x₁[2], x₁[3], x₁[4]], 2)
+    u₁, u₂ = rand(4), rand(4)
+    uᵢ = [ u₁[2], u₁[3], u₁[4], u₂[2], u₂[3], u₂[4]]
+
+    c₁ = Tetrahedron((1,2,3,4))
+    c₂ = Tetrahedron((5,2,3,4))
+
+    reinit!(cv, xᵢ)
+    reinit!(iv, c₁, x₁, 3, c₂, x₂, 3)
+
+    for qp in 1:getnquadpoints(cv)
+        #@test getdetJdV(cv, qp) == getdetJdV(iv, qp) <- TODO
+        @test function_value(cv, qp, uᵢ; here=true)  == function_value(iv, qp, vcat(u₁, u₂); here=true)
+        @test function_value(cv, qp, uᵢ; here=false) == function_value(iv, qp, vcat(u₁, u₂); here=false)
+        @test function_value_jump(cv, qp, uᵢ) == function_value_jump(iv, qp, vcat(u₁, u₂))
+        #@test function_gradient_jump(cv, qp, uᵢ) == function_gradient_jump(iv, qp, vcat(u₁, u₂)) <- TODO
+        #@test getnormal(cv, qp) == getnormal(iv, qp) <- TODO
+    end
 end