Question about `IP3()` and `IP6()`

[marinlauber]

@cristophermoen

I've been through the source code of this package and of Ferrite.jl and I've got a few questions. I'll put them here as it might be useful to keep track of those.

I've noticed that the IP3() and IP6() ScalarInterpolation are very similar to the Lagrange{RefTriangle, order} from Ferrite

struct Lagrange{shape, order} <: ScalarInterpolation{shape, order}
    ...
end

The only difference that I could find is that there is a switch in the node numbering, which means that the shape functions are swapped.

function reference_shape_value(ip::Lagrange{RefTriangle, 1}, ξ::Vec{2}, i::Int)
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 1 && return ξ_x
    i == 2 && return ξ_y
    i == 3 && return 1 - ξ_x - ξ_y
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

The same goes for the IP6() and Lagrange{RefTriangle, 2}.

I've tested it on a Cook's membrane case, and I get identical results if I swap them for Ferrite's internal ones

# ip3 = TriShellFiniteElement.IP3()
# ip6 = TriShellFiniteElement.IP6()
ip3 = Lagrange{RefTriangle, 1}()
ip6 = Lagrange{RefTriangle, 2}()

This leads to the two following questions

• Q1. Is it not possible to use the internal ScalarInterpolation from Ferrite?
• Q2. Regardless if the answer to Q1, I am not sure why IP3() works since we have a 3D displacement field which would require IP3()^3 or similar? I guess this crashes due to some internal assumption in the CellValues?

CooksMembrane.jl script?

using TriShellFiniteElement
using Ferrite

function create_cook_grid(nx, ny)
    corners = [Tensors.Vec{2}((0.0, 0.0)),
               Tensors.Vec{2}((48.0, 44.0)),
               Tensors.Vec{2}((48.0, 60.0)),
               Tensors.Vec{2}((0.0, 44.0))]
    return generate_grid(Triangle, (nx, ny), corners)
end

# integrate the traction force
function traction_force_vector!(cell, fv, traction)
    fe = zeros(length(cell.dofs))
    for facet in 1:nfacets(cell)
        if (cellid(cell), facet) ∈ getfacetset(grid, "traction")
            reinit!(fv, cell, facet)
            for q_point in 1:getnquadpoints(fv)
                dΓ = getdetJdV(fv, q_point)
                for i in 1:getnbasefunctions(fv)
                    Nᵢ = shape_value(fv, q_point, i)
                    # @show i, Nᵢ, traction, dΓ
                    fe[i] += (Nᵢ ⋅ traction) * dΓ
                end
            end
        end
    end
    return fe
end

# explicit assembly of the stiffness and force vector
function assemble_shell!(K, F, dh, qr1, qr3, ip3, ip6, E, ν, t, traction)
    assembler = start_assemble(K, F)
    for cell in CellIterator(dh)
        x = getcoordinates(cell)
        ke = TriShellFiniteElement.elastic_stiffness_matrix!(qr1, qr3, ip3, ip6, E, ν, t, x)
        fe = traction_force_vector!(cell, fv, traction)
        assemble!(assembler, celldofs(cell), ke, fe)
    end
    return K, F
end

# Grid, dofhandler, boundary condition
n = 50
grid = create_cook_grid(n, n)

# facesets for boundary conditions
addfacetset!(grid, "clamped", x -> x[1] ≈ 0.0)
addfacetset!(grid, "traction", x -> x[1] ≈ 48.0)

# interpolation order
ip = Lagrange{RefTriangle,1}() #to define fields only
#TODO can be swapped
ip3 = TriShellFiniteElement.IP3()
ip6 = TriShellFiniteElement.IP6()
# ip3 = Lagrange{RefTriangle,1}()
# ip6 = Lagrange{RefTriangle,2}()

qr1 = QuadratureRule{RefTriangle}(1)
qr3 = QuadratureRule{RefTriangle}(2)
facet_qr = FacetQuadratureRule{RefTriangle}(1)

# cell and face value for u
cv = CellValues(qr1, ip3, ip3)
fv = FacetValues(facet_qr, ip^3)

# degrees of freedom for displacements and rotations
dh = DofHandler(grid)
add!(dh, :u, ip^3)
add!(dh, :θ, ip^2)
close!(dh)

# boundary conditions
dbc = ConstraintHandler(dh)
add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "clamped"), x -> zero(x), [1, 2]))
add!(dbc, Dirichlet(:θ, getfacetset(dh.grid, "clamped"), x -> zero(x), [1, 2]))
close!(dbc)

# material properties
E = 0.7 # 70 Pa in N/dm^2
t = 0.5
ν = 0.3333

# Assembly and solve
Ke = allocate_matrix(dh)
f = zeros(ndofs(dh))

# point evaluation
ph = PointEvalHandler(grid, [Tensors.Vec{2}((48.0, 60.0))])

# traction vector in units N/dm/thickness
traction = Tensors.Vec((0.0, 1/16*t, 0.0))

# assemble the system
Ke, f = assemble_shell!(Ke, f, dh, qr1, qr3, ip3, ip6, E, ν, t, traction)

# apply the BCs
apply!(Ke, f, dbc)
# solve
@time u = Ke \ f

# evaluate point value
u_eval = first(evaluate_at_points(ph, dh, u, :u))
@show u_eval
@show ip3, ip6

VTKGridFile("mindlin_shell", dh) do vtk
    write_solution(vtk, dh, u)
end

[cristophermoen]

@marinlauber thanks for thinking through all this

we went with custom interpolations mainly because we wanted to play around with them, hoping to have more control in future element formulations, when we want to tweak things.

Q1. agreed that the Ferrite internal will work for membrane.
Q2. Still trying to figure this out, there are a few things going on. For the element in its local coordinate system, we treat membrane and bending (shear, flexure) separately, and they need different shape functions, thus IP3 for membrane (is a 2D displacement field) and IP6 for bending (this is the trick implemented by others, Tessler, Hughes, and co) and then static condensation to 'absorb' intermediate dof. You can see some implementation details in this paper. We can have a talk this week again if you want.

Also, in general, we are still debugging and verifying Ke and Kg implementations, especially the 3D rotations (from global to local to global) and also the drilling dof. I am working on it about 1 hour a day, so kind of slow.

Moen_et_al_SSRC_2026.pdf

[marinlauber]

@cristophermoen Thanks for the clarification.