Assemble several elements at a time

Now assemble! takes a vector of elements as input. This makes it
possible to preallocate memory for common matrices making code super
fast.

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)

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,

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

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))
 

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,

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