Added interpolation Lagrange{2,RefCube,3}

[JanM12]

introduced interpolation Lagrange{2,RefCube,3} to src/interpolations.jl

[koehlerson]

need to add the new interpolation to the test suite, see: https://github.com/Ferrite-FEM/Ferrite.jl/blob/master/test/test_interpolations.jl#L3-L36

[fredrikekre]

Looks like the tolerance here

Ferrite.jl/test/test_interpolations.jl

Line 106 in ade9292

@test N_dof[k] ≈ 0.0 atol=4eps(Float64) #broken=typeof(interpolation)==Lagrange{2, RefTetrahedron, 5}&&dof==4&&k==18

needs to be adjusted (or possibly the long expressions here can be reordered to give a smaller error like what was done in #633, but maybe not possible).

[termi-official]

Give me some more time to think about this adjustment, because I am not really sure about what could break downstream due to such spurious evaluations and what the related thresholds to prevent numerical oscillations are. Also I do not know how these affect the overall assembly precision. Generally speaking we should still be able to reorder the expressions to eliminate errors in the tested delta property (even though it is not trivial).

[fredrikekre]

Perhaps we can BigFloats and truncate to Float64 after evaluation.

[termi-official]

But this will be a big performance hit for the value call then.

[fredrikekre]

Maybe, but it is only called in the constructor.

[AbdAlazezAhmed]

How about using Rational ?

[termi-official]

Thanks for this suggestions!

Maybe, but it is only called in the constructor.

For order 3 (and maybe order 4) this should be fine, but when moving to higher order one quickyl runs into the issue that assembly degrades solver performance quite severly. Here usually matrix-free techniques are deployed which require reevaluation of the basis functions. I think BigFloat could hurt performance in this scenario quite a bit. To not delay this PW we can go either way and discuss this in a separate issue/discussion, investigating this in more detail (and maybe gathering some literature on this).

How about using Rational ?

This would certainly fix the tests, but IIRC then general quadrature points are unfortunately Irrational. However, I am happy to get corrected if I am wrong on this.

[AbdAlazezAhmed]

This would certainly fix the tests, but IIRC then general quadrature points are unfortunately Irrational. However, I am happy to get corrected if I am wrong on this.

You're right. AccurateArithmetic.jl can be helpful?

[termi-official]

Good point. Kahan summation could solve the problem already.

[termi-official]

Okay maybe we can bring this over the finish line as follows. We can increase the number of points in the partition of unity test here (points outside the domain should be rejected!)

Ferrite.jl/test/test_interpolations.jl

Lines 50 to 55 in 1fadc1f

	n_basefuncs = getnbasefunctions(interpolation)
	x = rand(Tensor{1, ref_dim})
	f = (x) -> Ferrite.value(interpolation, Tensor{1, ref_dim}(x))
	@test vec(ForwardDiff.jacobian(f, Array(x))') ≈
	reinterpret(Float64, Ferrite.derivative(interpolation, x))
	@test sum(Ferrite.value(interpolation, x)) ≈ 1.0

and increase the tolerance here

Ferrite.jl/test/test_interpolations.jl

Line 107 in 1fadc1f

@test N_dof[k] ≈ 0.0 atol=4eps(Float64) #broken=typeof(interpolation)==Lagrange{2, RefTetrahedron, 5}&&dof==4&&k==18

to 1e-14 . This should be reasonably small for now.

Alternatively we can try to overdub the value calls in the tests with Kahan summations, but this might be a bit more work. It think the linked package fails with AD (just skipped over the implementation).