perf_test.jl

# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/PDEAssembler.jl/blob/master/LICENSE

using FEMBase, Base.Threads

struct Poisson <: FieldProblem end

FEMBase.get_unknown_field_name(::Problem{Poisson}) = "u"

function FEMBase.assemble_elements!(problem::Problem{Poisson},
                                    assembly::Assembly,
                                    elements::Vector{Element{E}},
                                    time::Float64) where E

    bi = BasisInfo(E)
    ndofs = length(bi)
    Ke = zeros(ndofs, ndofs)
    fe = zeros(ndofs)

    for element in elements
        fill!(Ke, 0.0)
        fill!(fe, 0.0)
        for ip in get_integration_points(element)
            J, detJ, N, dN = element_info!(bi, element, ip, time)
            k = 1.0
            if haskey(element, "coefficient")
                k = element("coefficient", ip, time)
            end
            Ke += ip.weight * k*dN'*dN * detJ
            if haskey(element, "source")
                f = element("source", ip, time)
                fe += ip.weight * N'*f * detJ
            end
        end
        gdofs = get_gdofs(problem, element)
        add!(assembly.K, gdofs, gdofs, Ke)
        add!(assembly.f, gdofs, fe)
    end

    return nothing

end

function create_element()
    x1 = [2.0, 3.0, 4.0]
    x2 = [6.0, 3.0, 2.0]
    x3 = [2.0, 5.0, 1.0]
    x4 = [4.0, 3.0, 6.0]
    x5 = 0.5*(x1+x2)
    x6 = 0.5*(x2+x3)
    x7 = 0.5*(x3+x1)
    x8 = 0.5*(x1+x4)
    x9 = 0.5*(x2+x4)
    x10 = 0.5*(x3+x4)
    X = Dict(
        1 => x1, 2 => x2, 3 => x3, 4 => x4, 5 => x5,
        6 => x6, 7 => x7, 8 => x8, 9 => x9, 10 => x10)
    u = Dict(i => zeros(3) for i in 1:10)
    element = Element(Tet10, (1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
    update!(element, "geometry", X)
    update!(element, "displacement", u)
    return element
end

function test_assembler_timing()
    problem = Problem(Poisson, "test problem", 1)
    assembly = problem.assembly
    elements = [create_element() for _=1:100000]
    time = 0.0
    assemble_elements!(problem, assembly, elements, time)
    empty!(assembly)
    @time assemble_elements!(problem, assembly, elements, time)
end

test_assembler_timing()