Add three-argument reinit! for InterfaceCellValues

[MelanieInky]

Add three-arguments reinit! for InterfaceCellValues. I am currently working on a cohesive zone model using FerriteAssembly and it uses the three argument versions of reinit!. Doing so leads to the error

"
ERROR: LoadError: MethodError: no method matching get_geo_mapping(::InterfaceCellValues{CellValues{Ferrite.FunctionValues{…}, Ferrite.GeometryMapping{…}, QuadratureRule{…}, Vector{…}}})
The function get_geo_mapping exists, but no method is defined for this combination of argument types.

Closest candidates are:
get_geo_mapping(::CellValues)
@ Ferrite ~/.julia/packages/Ferrite/e5O5M/src/FEValues/CellValues.jl:90
get_geo_mapping(::MultiFieldCellValues)
@ Ferrite ~/.julia/packages/Ferrite/e5O5M/src/FEValues/CellValues.jl:265
get_geo_mapping(::FacetValues)
@ Ferrite ~/.julia/packages/Ferrite/e5O5M/src/FEValues/FacetValues.jl:91
"

Implementing a three-argument version solves the issue. Current behavior:

I have added some AI generated tests. Also below an AI generated minimum working example using FerriteAssembly involving a simple cohesive element. I tested with and without the patch. It assembles fine with it and throws a MethodError without.

# MWE: FerriteAssembly cannot assemble a FerriteInterfaceElements domain because
# CellBuffer reinits through the 3-arg `reinit!`, which InterfaceCellValues lacks.
#
# Run top to bottom. It errors, then defines the proposed 3-arg method and succeeds.
 
using Ferrite, FerriteInterfaceElements, FerriteAssembly, Tensors
using SparseArrays: nnz
 
# --- a trivial penalty interface element, so the only thing under test is reinit! ---
struct DummyCohesive{T}
    K::T
end
 
function FerriteAssembly.element_routine!(Ke, re, state, ae, m::DummyCohesive,
                                          cv::InterfaceCellValues, buffer)
    for qp in 1:getnquadpoints(cv)
        dΓ = getdetJdV_average(cv, qp)
        jump_u = function_value_jump(cv, qp, ae)
        for i in 1:getnbasefunctions(cv)
            δN = shape_value_jump(cv, qp, i)
            re[i] += m.K * (δN ⋅ jump_u) * dΓ
            for j in 1:getnbasefunctions(cv)
                Ke[i, j] += m.K * (δN ⋅ shape_value_jump(cv, qp, j)) * dΓ
            end
        end
    end
    return nothing
end
 
function setup()
    grid = generate_grid(Quadrilateral, (2, 1))
    addcellset!(grid, "A", Set([1]))
    addcellset!(grid, "B", Set([2]))
    grid = insert_interfaces(grid, ["A", "B"])
 
    ip_bulk = Lagrange{RefQuadrilateral,1}()^2
    ip_int  = InterfaceCellInterpolation(Lagrange{RefLine,1}())^2
 
    dh = DofHandler(grid)
    sdh_bulk = SubDofHandler(dh, union(getcellset(grid, "A"), getcellset(grid, "B")))
    add!(sdh_bulk, :u, ip_bulk)
    sdh_int = SubDofHandler(dh, getcellset(grid, "interfaces"))
    add!(sdh_int, :u, ip_int)
    close!(dh)
 
    cv_int = InterfaceCellValues(QuadratureRule{RefLine}(2), ip_int)
    db = setup_domainbuffer(DomainSpec(sdh_int, DummyCohesive(1.0e3), cv_int))
    return dh, db
end
 
function assemble(dh, db)
    K = allocate_matrix(dh)
    r = zeros(ndofs(dh))
    a = zeros(ndofs(dh))
    set_time_increment!(db, 0.0)          # Δt defaults to NaN
    work!(start_assemble(K, r), db; a = a)
    return K, r
end
 
dh, db = setup()
 
println("Ferrite               ", pkgversion(Ferrite))
println("FerriteAssembly       ", pkgversion(FerriteAssembly))
println("FerriteInterfaceElements ", pkgversion(FerriteInterfaceElements))
 
assemble(dh,db)

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 92.57%. Comparing base (6e2a63b) to head (416431a).

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #48      +/-   ##
==========================================
+ Coverage   91.66%   92.57%   +0.90%     
==========================================
  Files           6        6              
  Lines         348      350       +2     
==========================================
+ Hits          319      324       +5     
+ Misses         29       26       -3     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[KnutAM]

Nice!

[MelanieInky]

Edit: as of 45a0ce4 my comment is no longer necessary as the issue was fixed upstream.
Adding a comment to say that the line

n_coords_per_side = length(x) ÷ 2

seems to work fine but it leaves me confused as to the distinction when the element is not isoparametric.

We use (from the constructor)

base_indices_here  = collect( get_interface_index(ip, :here,  i) for i in 1:getnbasefunctions(ip.base) )

But the length of x is the number of geometric nodes * 2 (ip_geo), not the shape values (ip). I am not super well versed in non isoparametric elements so I left it as is.

[KnutAM]

Looks good! Just a few small things to fix and test.

[KnutAM]

Suggested change

Ferrite.reinit!(cv::InterfaceCellValues, cc::CellCache) = reinit!(cv, getcells(cc.grid, Ferrite.cellid(cc)), cc.coords)

This can be deleted right since cv::AbstractCellValues and Ferrite.reinit_needs_cell defined above?

[MelanieInky]

You're right it can branch on Ferrite.jl version.

[KnutAM]

Are these tests relevant for this PR? Seems like a different change (adding tests)? Isn't this already tested somewhere else?

[MelanieInky]

There was a typo that passed through tests before #47 that was fixed in commit 2ec45e7 in this PR.

The commit was a one line change:

- x_there = @view x[cv.base_indices_there]
+ x_there = @view x[cv.base_indices_there[1:n_coords_per_side]]

The tests didn't catch it because there was no appropriate covering when include_R==true with superparametric elements. Now there is so this is partially redundant.
Now the code has been changed upstream so the typo fix is no longer relevant, but the fact that no test caught it is still relevant. I will check again if this is tested since I just appended both this PR and #47 tests.

[MelanieInky]

Trimmed in 416431a, so there should not be any redundant part now

With your comments in mind I am not a fan of this testset I think what it tests should be included in other testsets, or maybe it should be expanded to handle rotations for a variety of geometric order? Something like (Current code passes this)

@testset "midplane_rotation under rigid rotation" begin
    for (shape, cells, dim) in ((RefLine, (Line, QuadraticLine), 2),
                                (RefTriangle, (Triangle, QuadraticTriangle), 3),
                                (RefQuadrilateral, (Quadrilateral, QuadraticQuadrilateral), 3))
        Q = dim == 2 ? rotation_tensor(0.7) :
                       rotation_tensor(Vec{3}((1.0, 2.0, 3.0)) / sqrt(14.0), 0.7)
        for forder in (1, 2), gorder in (1, 2)
            base_fip = Lagrange{shape, forder}()
            base_gip = Lagrange{shape, gorder}()
            fip = InterfaceCellInterpolation(base_fip)
            gip = InterfaceCellInterpolation(base_gip)
            qr  = QuadratureRule{shape}(2)
            cv  = InterfaceCellValues(qr, fip, gip; include_R = true)

            xh = [Vec{dim}(i -> i < dim ? ξ[i] : 0.0) for ξ in Ferrite.reference_coordinates(base_gip)]
            xt = [x + Vec{dim}(i -> i == dim ? 1.0 : 0.0) for x in xh]
            nn = length(xh)
            C  = cells[gorder]
            cell = InterfaceCell(C(Tuple(1:nn)), C(Tuple(nn+1:2nn)))
            x = vcat(xh, xt)[collect(cell.nodes)]

            reinit!(cv, x)
            R0 = [midplane_rotation(cv, qp) for qp in 1:getnquadpoints(cv)]

            reinit!(cv, [Q ⋅ xi for xi in x])
            for qp in 1:getnquadpoints(cv)
                R = midplane_rotation(cv, qp)
                @test tdot(R) ≈ one(R)
                @test det(R) ≈ 1.0
                @test R ≈ Q ⋅ R0[qp]     # equivariance: rotating the cell rotates the frame
            end
        end
    end
end

[KnutAM]

Looks like a good test! Some suggestions to simplify and use the querry functions we have (not tested code so probably some bugs)

@testset "midplane_rotation under rigid rotation" begin
    for base_cell in (Line, QuadraticLine, Triangle, QuadraticTriangle, Quadrilateral, QuadraticQuadrilateral)
        dim = Ferrite.getrefdim(base_cell) + 1
        base_gip = geometric_interpolation(base_cell)
        gip = InterfaceCellInterpolation(base_gip)
        shape = Ferrite.getrefshape(base_cell)
        Q = dim == 2 ? rotation_tensor(0.7) :
                       rotation_tensor(Vec{3}((1.0, 2.0, 3.0)) / sqrt(14.0), 0.7)
        for forder in (1, 2)
            base_fip = Lagrange{shape, forder}()
            fip = InterfaceCellInterpolation(base_fip)
            qr  = QuadratureRule{shape}(2)
            cv  = InterfaceCellValues(qr, fip, gip; include_R = true)

            # Needs a comment to motivate what xh and xt are.
            xh = [Vec{dim}(i -> i < dim ? ξ[i] : 0.0) for ξ in Ferrite.reference_coordinates(base_gip)]
            xt = [x + Vec{dim}(i -> i == dim ? 1.0 : 0.0) for x in xh]

            nn = length(xh)
            cell = InterfaceCell(base_cell(Tuple(1:nn)), base_cell(Tuple(nn+1:2nn)))
            x = vcat(xh, xt)[collect(cell.nodes)]

            reinit!(cv, x)
            R0 = [midplane_rotation(cv, qp) for qp in 1:getnquadpoints(cv)]

            reinit!(cv, [Q ⋅ x_i for x_i in x])
            for qp in 1:getnquadpoints(cv)
                R = midplane_rotation(cv, qp)
                @test tdot(R) ≈ one(R)
                @test det(R) ≈ 1.0
                @test R ≈ Q ⋅ R0[qp]     # equivariance: rotating the cell rotates the frame
            end
        end
    end
end

For nonlinear geometry, we could potentially ensure that also the geometry is nonlinear by doing e.g.

- xt = [x + Vec{dim}(i -> i == dim ? 1.0 : 0.0) for x in xh]
+ xt = [x + (1 + norm(x)^2)/2) * Vec{dim}(i -> i == dim ? 1.0 : 0.0) for x in xh]

[MelanieInky]

Should it be its own PR?

[KnutAM]

Yes, that could be better! Then remove the full test and add in another PR!

[KnutAM]

Should also add a test with CellCache to include the Ferrite.reinit_needs_cell overload.

[MelanieInky]

Done. I am not sure how to handle the case where reinit_needs_cell == true because I don't think we support non identity mapping in this package yet.

[KnutAM]

It would be possible to create a dummy interpolation to test with, but I think it is enough as it is now.