Add Lagrange{RefTetrahedron,3}

[termi-official]

No description provided.

[codecov]

Codecov Report

❌ Patch coverage is 37.64706% with 53 lines in your changes missing coverage. Please review.
✅ Project coverage is 93.47%. Comparing base (65773c2) to head (b01d1b6).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/interpolations.jl	37.64%	53 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1255      +/-   ##
==========================================
- Coverage   94.19%   93.47%   -0.72%     
==========================================
  Files          40       40              
  Lines        6662     6747      +85     
==========================================
+ Hits         6275     6307      +32     
- Misses        387      440      +53     

☔ View full report in Codecov by Sentry.
📢 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]

Seems like this should be a macro or generated function?

[termi-official]

Not 100% sure yet. This is for now just a template to get things started on higher order element generators, as we can simply generate these in a systematic way.

[AbdAlazezAhmed]

Maybe

@generated function do_magic(::Lagrange{RefTetrahedron, order}, ξ::Vec{3, Float64}, i::Int) where {order}

    N = 0
    for p₃ in 0:order
        for p₂ in 0:(order - p₃)
            for p₁ in 0:(order - p₃ - p₂)
                N += 1
            end
        end
    end

    @assert N == 20
    xs = reference_coordinates(Lagrange{RefTetrahedron, order})
    V = zeros(N, N)
    j = 0
    for p₃ in 0:order
        for p₂ in 0:(order - p₃)
            for p₁ in 0:(order - p₃ - p₂)
                j += 1
                for i in 1:N
                    x = xs[i]
                    V[i, j] = x[1]^p₁ * x[2]^p₂ * x[3]^p₃
                end
            end
        end
    end
    C = inv(V)' # These are our coefficients
    ex = :(    ξ_x = ξ[1];
    ξ_y = ξ[2];
    ξ_z = ξ[3];)
    for i in 1:N
        j = 0
        x = xs[i]
        ex_i = :(0.0)
        for p₃ in 0:order
            for p₂ in 0:(order - p₃)
                for p₁ in 0:(order - p₃ - p₂)
                    j += 1
                    ex_i = :($ex_i + $(C[i, j]) * ξ_x^$p₁ * ξ_y^$p₂ * ξ_z^$p₃)
                end
            end
        end
        ex = :($ex; i == $i && return $ex_i)
    end

    return ex
end

julia> @benchmark Ferrite.do_magic($ip, $(Vec(-0.0,-0.0,-0.0)),$(1))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  4.882 ns … 20.719 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     5.008 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.030 ns ±  0.389 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅     █▂                                                   
  ▅█▆▃▂▂▄██▄▂▂▂▂▁▂▂▂▁▂▄▄▃▂▂▂▂▂▂▂▂▁▁▂▂▂▂▁▁▁▂▂▂▂▁▁▂▂▁▁▁▁▁▂▂▁▁▂ ▃
  4.88 ns        Histogram: frequency by time        5.86 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark Ferrite.reference_shape_value($ip, $(Vec(-0.0,-0.0,-0.0)),$(1))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  4.581 ns … 9.742 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     4.619 ns             ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.673 ns ± 0.209 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅█                                                        
  ▂██▂▂▂▂▄▄▂▂▂▂▂▂▂▂▂▂▂▄▂▂▂▂▂▂▂▁▂▂▂▁▂▂▂▂▁▁▁▁▁▁▁▂▁▂▂▂▂▂▁▁▁▁▁▂ ▂
  4.58 ns        Histogram: frequency by time       5.53 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

compiler vroom

[AbdAlazezAhmed]

maybe

Suggested change

	getnbasefunctions(::Lagrange{RefTetrahedron, 3}) = 20
	getnbasefunctions(::Lagrange{RefTetrahedron, order}) where order = (order+1)*(order+2)*(order+3)÷6

would be beneficial in the future

[AbdAlazezAhmed]

Maybe

@generated function do_magic(::Lagrange{RefTetrahedron, order}, ξ::Vec{3, Float64}, i::Int) where {order}

    N = 0
    for p₃ in 0:order
        for p₂ in 0:(order - p₃)
            for p₁ in 0:(order - p₃ - p₂)
                N += 1
            end
        end
    end

    @assert N == 20
    xs = reference_coordinates(Lagrange{RefTetrahedron, order})
    V = zeros(N, N)
    j = 0
    for p₃ in 0:order
        for p₂ in 0:(order - p₃)
            for p₁ in 0:(order - p₃ - p₂)
                j += 1
                for i in 1:N
                    x = xs[i]
                    V[i, j] = x[1]^p₁ * x[2]^p₂ * x[3]^p₃
                end
            end
        end
    end
    C = inv(V)' # These are our coefficients
    ex = :(    ξ_x = ξ[1];
    ξ_y = ξ[2];
    ξ_z = ξ[3];)
    for i in 1:N
        j = 0
        x = xs[i]
        ex_i = :(0.0)
        for p₃ in 0:order
            for p₂ in 0:(order - p₃)
                for p₁ in 0:(order - p₃ - p₂)
                    j += 1
                    ex_i = :($ex_i + $(C[i, j]) * ξ_x^$p₁ * ξ_y^$p₂ * ξ_z^$p₃)
                end
            end
        end
        ex = :($ex; i == $i && return $ex_i)
    end

    return ex
end

julia> @benchmark Ferrite.do_magic($ip, $(Vec(-0.0,-0.0,-0.0)),$(1))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  4.882 ns … 20.719 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     5.008 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.030 ns ±  0.389 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅     █▂                                                   
  ▅█▆▃▂▂▄██▄▂▂▂▂▁▂▂▂▁▂▄▄▃▂▂▂▂▂▂▂▂▁▁▂▂▂▂▁▁▁▂▂▂▂▁▁▂▂▁▁▁▁▁▂▂▁▁▂ ▃
  4.88 ns        Histogram: frequency by time        5.86 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark Ferrite.reference_shape_value($ip, $(Vec(-0.0,-0.0,-0.0)),$(1))
BenchmarkTools.Trial: 10000 samples with 1000 evaluations per sample.
 Range (min … max):  4.581 ns … 9.742 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     4.619 ns             ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.673 ns ± 0.209 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅█                                                        
  ▂██▂▂▂▂▄▄▂▂▂▂▂▂▂▂▂▂▂▄▂▂▂▂▂▂▂▁▂▂▂▁▂▂▂▂▁▁▁▁▁▁▁▂▁▂▂▂▂▂▁▁▁▁▁▂ ▂
  4.58 ns        Histogram: frequency by time       5.53 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

compiler vroom