Numerical Evaluation of Order 3,4,5 Lagrange Triangle Inexact

[termi-official]

At least the current implementation of the order 5 Lagrange triangle is inexact, as revealed by improved test coverage of the dirac delta property of the interpolations in #581 .

Edit: Forgot to update after some discussion with Abdulaziz.

[fredrikekre]

Can you post an MWE here?

[termi-official]

using Ferrite
interpolation = Lagrange{2, RefTetrahedron, 5}();
coords = Ferrite.reference_coordinates(interpolation);
Ferrite.value(interpolation, coords[4])

is already slightly inexact. Not sure what is going on.

[termi-official]

Ah, sorry, I see what is happening. Equidistant points for the basis functions. They are inherently broken, because they suffer from Runge's phenomenon.

[AbdAlazezAhmed]

For values being inexact. I think replacing γ = 1. - ξ_x - ξ_y with γ = 1. - (ξ_x + ξ_y) solves the issue.

function value(ip::Lagrange2Tri345, i::Int, ξ::Vec{2})
    order = getorder(ip)
    i = permdof2D[order][i]
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    γ = 1. - (ξ_x + ξ_y)
    i1, i2, i3 = _numlin_basis2D(i, order)
    val = one(γ)
    i1 ≥ 1 && (val *= prod((order * γ - j) / (j + 1) for j = 0:(i1 - 1)))
    i2 ≥ 1 && (val *= prod((order * ξ_x - j) / (j + 1) for j = 0:(i2 - 1)))
    i3 ≥ 1 && (val *= prod((order * ξ_y - j) / (j + 1) for j = 0:(i3 - 1)))
    return val
end

julia> using Ferrite

julia> interpolation = Lagrange{2, RefTetrahedron, 5}();

julia> coords = Ferrite.reference_coordinates(interpolation);

julia> Ferrite.value(interpolation, coords[4])
21-element Vector{Float64}:
  0.0
 -0.0
  0.0
  1.0
  0.0
 -0.0
  0.0
  0.0
  0.0
  0.0
 -0.0
 -0.0
  0.0
 -0.0
  0.0
  0.0
 -0.0
  0.0
 -0.0
  0.0
 -0.0

What do you think?

[AbdAlazezAhmed]

Turns out it's not for all nodes

julia> Ferrite.value(interpolation, coords[15])
21-element Vector{Float64}:
  0.0
  0.0
  1.1102230246251566e-17
  0.0
 -0.0
  0.0
 -0.0
 -0.0
 -0.0
 -0.0
 -0.0
 -7.401486830834378e-17
  2.220446049250313e-16
 -4.440892098500625e-16
  0.9999999999999998
  0.0
 -0.0
  0.0
  0.0
 -0.0
  0.0

From what I've seen the value of γ is the issue here (1. - 0.8 = 0.19999999999999996)

[AbdAlazezAhmed]

Using γ directly as in i1 ≥ 1 && (val *= prod((order *1-(order * ξ_x +order * ξ_y + j)) / (j + 1) for j = 0:(i1 - 1))) Solves the issue. It looks ugly atm, I'll check it again later today. did some testing on all nodes in order 5, it worked.

[termi-official]

Indeed, this could already be the issue. We have automatized some of the tests regarding the interpolations, so you can check if this works by commenting in this line

Ferrite.jl/test/test_interpolations.jl

Line 12 in cc81e6c

#Lagrange{2, RefTetrahedron, 5}(),

in the test battery.