diff --git a/src/assembly.jl b/src/assembly.jl index dad60dc..f2ec1c6 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -37,9 +37,7 @@ function assemble!(problem::Problem, time=0.0; auto_initialize=true) if method_exists(assemble_prehook!, Tuple{typeof(problem), Float64}) assemble_prehook!(problem, time) end - for element in get_elements(problem) - assemble!(problem.assembly, problem, element, time) - end + assemble!(get_assembly(problem), problem, get_elements(problem), time) if method_exists(assemble_posthook!, Tuple{typeof(problem), Float64}) assemble_posthook!(problem, time) end @@ -73,7 +71,7 @@ function assemble!(problem::Problem, time::Real, ::Type{Val{:mass_matrix}}; dens end end -function assemble!(assembly::Assembly, problem::Problem, elements::Vector{Element}, time::Real) +function assemble!(assembly::Assembly, problem::Problem, elements::Vector{Element}, time) warn("assemble!() this is default assemble operation, decreased performance can be expected without preallocation of memory!") for element in elements assemble!(assembly, problem, element, time) diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 61c2778..11b5664 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -54,18 +54,23 @@ function get_formulation_type(problem::Problem{Elasticity}) end """ + assemble!(assembly:Assembly, problem::Problem{Elasticity}, elements, time) + Start finite element assembly procedure for Elasticity problem. + +Function groups elements to arrays by their type and assembles one element type +at time. This makes it possible to pre-allocate matrices common to same type +of elements. """ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{Element}, time) - assemble!(assembly, problem, elements, time, Val{problem.properties.formulation}) -end - -""" -This is for backward compatibility, will be removed asap. -""" -function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time) - warn("try to avoid single element assembly function as it's not possible to preallocate causing a slow code") - assemble!(assembly, problem, [element], time, Val{problem.properties.formulation}) + formulation = Val{problem.properties.formulation} + element_types = unique(map(get_element_type, elements)) + for element_type in element_types + elements_subset = filter_by_element_type(element_type, elements) + elements_subset = [element for element in elements_subset] + nelements = length(elements_subset) + assemble!(assembly, problem, elements_subset, time, formulation) + end end include("problems_elasticity_2d.jl") @@ -102,26 +107,6 @@ function get_keys(element) map(x -> all_keys[x], idx) end -""" Continuum elements assembly entry point. - -This splits elements to arrays by their type and assemble one element type -at time. This makes it possible to pre-allocate matrices common to same type -of elements. -""" -function assemble!(assembly::Assembly, problem::Problem{Elasticity}, - all_elements::Vector{Element}, time, ::Type{Val{:continuum}}) - element_types = unique(map(get_element_type, all_elements)) - for element_type in element_types - elements = filter_by_element_type(element_type, all_elements) - # FIXME: there must be better way to do this - # to promote array for certain elemene type - elements = [element for element in elements] - nelements = length(elements) - debug("elasticity 3d: assembling $nelements of type $element_type") - assemble!(assembly, problem, elements, time, Val{:continuum}) - end -end - """ Assemble 3d continuum elements in general solid mechanics problem. """ function assemble!{El<:Elasticity3DVolumeElements}(assembly::Assembly, diff --git a/test/test_common_failures.jl b/test/test_common_failures.jl index eb2b886..774d044 100644 --- a/test/test_common_failures.jl +++ b/test/test_common_failures.jl @@ -7,8 +7,9 @@ using JuliaFEM.Testing @testset "geometry missing" begin el = Element(Quad4, [1, 2, 3, 4]) pr = Problem(Elasticity, "problem", 2) + add_elements!(pr, [el]) # this throws KeyError: geometry not found. # it's descriptive enough to give hint to user # what went wrong - @test_throws KeyError assemble!(pr, el) + @test_throws KeyError assemble!(pr) end diff --git a/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl b/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl index e3f2e7c..d3f6c5c 100644 --- a/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl +++ b/test/test_elasticity_2d_plane_stress_stiffness_matrix.jl @@ -24,8 +24,9 @@ using JuliaFEM.Testing update!(element, "displacement load", DCTI([4.0, 8.0])) problem = Problem(Elasticity, "[0x1] x [0x1] block", 2) - problem.properties.formulation = :plane_stress - assemble!(problem, element) + update!(problem.properties, "formulation" => "plane_stress") + add_elements!(problem, [element]) + assemble!(problem) K = full(problem.assembly.K) f = vec(full(problem.assembly.f)) diff --git a/test/test_elasticity_tet10_stiffness_matrix.jl b/test/test_elasticity_tet10_stiffness_matrix.jl index 34e1f7d..52f97f5 100644 --- a/test/test_elasticity_tet10_stiffness_matrix.jl +++ b/test/test_elasticity_tet10_stiffness_matrix.jl @@ -2,7 +2,6 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md using JuliaFEM -using JuliaFEM.Preprocess using JuliaFEM.Testing @testset "test tet10 stiffness matrix" begin @@ -29,8 +28,9 @@ using JuliaFEM.Testing update!(el, "geometry", X) update!(el, "displacement", u) pr = Problem(Elasticity, "tet10", 3) - ass = Assembly() - assemble!(ass, pr, el, 0.0) + add_elements!(pr, [el]) + assemble!(pr) + ass = pr.assembly Kt = full(ass.K) eigs = real(eigvals(Kt)) eigs_expected = [8809.45, 4936.01, 2880.56, 2491.66, 2004.85, diff --git a/test/test_elasticity_tet4_stiffness_matrix.jl b/test/test_elasticity_tet4_stiffness_matrix.jl index 962ad9d..cd3e89c 100644 --- a/test/test_elasticity_tet4_stiffness_matrix.jl +++ b/test/test_elasticity_tet4_stiffness_matrix.jl @@ -2,7 +2,6 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md using JuliaFEM -using JuliaFEM.Preprocess using JuliaFEM.Testing @testset "test tet4 stiffness matrix" begin @@ -17,8 +16,9 @@ using JuliaFEM.Testing el["geometry"] = Vector{Float64}[x1, x2, x3, x4] u = Vector{Float64}[u1, u2, u3, u4] pr = Problem(Elasticity, "tet4", 3) - as = Assembly() - assemble!(as, pr, el, 0.0) + add_elements!(pr, [el]) + assemble!(pr) + as = pr.assembly Kt = full(as.K) Kt_expected = [ 149.0 108.0 24.0 -1.0 6.0 12.0 -54.0 -48.0 0.0 -94.0 -66.0 -36.0