diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl
index 9a543d3..766b6a1 100644
--- a/src/problems_elasticity.jl
+++ b/src/problems_elasticity.jl
@@ -53,17 +53,19 @@ function get_formulation_type(problem::Problem{Elasticity})
     return :incremental
 end
 
-function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time=0.0)
-    props = problem.properties
-    gdofs = get_gdofs(problem, element)
-    formulation = props.formulation
-    if formulation in [:plane_stress, :plane_strain]
-        formulation = :plane
-    end
-    Km, Kg, f = assemble(problem, element, time, Val{formulation})
-    add!(assembly.K, gdofs, gdofs, Km)
-    add!(assembly.Kg, gdofs, gdofs, Kg)
-    add!(assembly.f, gdofs, f)
+"""
+Start finite element assembly procedure for Elasticity problem.
+"""
+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})
 end
 
 include("problems_elasticity_2d.jl")
@@ -100,235 +102,282 @@ function get_keys(element)
     map(x -> all_keys[x], idx)
 end
 
-""" Elasticity equations, 3d nonlinear. """
-function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
+""" 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,
+                                                   problem::Problem{Elasticity},
+                                                   elements::Vector{Element{El}},
+                                                   time, ::Type{Val{:continuum}})
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
-    nnodes = length(element)
-    ndofs = dim*nnodes
-    BL = zeros(6, ndofs)
-    BNL = zeros(9, ndofs)
-    Km = zeros(ndofs, ndofs)
-    Kg = zeros(ndofs, ndofs)
-    f = zeros(ndofs)
 
-    for ip in get_integration_points(element)
-        detJ = element(ip, time, Val{:detJ})
-        w = ip.weight*detJ
-        N = element(ip, time)
-        dN = element(ip, time, Val{:Grad})
+    for element in elements
 
-        # kinematics; calculate deformation gradient and strain
+        nnodes = length(element)
+        ndofs = dim*nnodes
+        BL = zeros(6, ndofs)
+        BNL = zeros(9, ndofs)
+        Km = zeros(ndofs, ndofs)
+        Kg = zeros(ndofs, ndofs)
+        f = zeros(ndofs)
 
-        gradu = zeros(dim, dim)
-        if haskey(element, "displacement")
-            gradu += element("displacement", ip, time, Val{:Grad})
-        end
-        strain = 1/2*(gradu' + gradu)
+        for ip in get_integration_points(element)
+            detJ = element(ip, time, Val{:detJ})
+            w = ip.weight*detJ
+            N = element(ip, time)
+            dN = element(ip, time, Val{:Grad})
 
-        F = eye(dim)
-        if props.finite_strain
-            F += gradu
-            strain += 1/2*gradu'*gradu
-        end
+            # kinematics; calculate deformation gradient and strain
 
-        # material stiffness start
+            gradu = zeros(dim, dim)
+            if haskey(element, "displacement")
+                gradu += element("displacement", ip, time, Val{:Grad})
+            end
+            strain = 1/2*(gradu' + gradu)
 
-        fill!(BL, 0.0)
-        for i=1:nnodes
-            BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i]
-            BL[1, 3*(i-1)+2] = F[2,1]*dN[1,i]
-            BL[1, 3*(i-1)+3] = F[3,1]*dN[1,i]
-            BL[2, 3*(i-1)+1] = F[1,2]*dN[2,i]
-            BL[2, 3*(i-1)+2] = F[2,2]*dN[2,i]
-            BL[2, 3*(i-1)+3] = F[3,2]*dN[2,i]
-            BL[3, 3*(i-1)+1] = F[1,3]*dN[3,i]
-            BL[3, 3*(i-1)+2] = F[2,3]*dN[3,i]
-            BL[3, 3*(i-1)+3] = F[3,3]*dN[3,i]
-            BL[4, 3*(i-1)+1] = F[1,1]*dN[2,i] + F[1,2]*dN[1,i]
-            BL[4, 3*(i-1)+2] = F[2,1]*dN[2,i] + F[2,2]*dN[1,i]
-            BL[4, 3*(i-1)+3] = F[3,1]*dN[2,i] + F[3,2]*dN[1,i]
-            BL[5, 3*(i-1)+1] = F[1,2]*dN[3,i] + F[1,3]*dN[2,i]
-            BL[5, 3*(i-1)+2] = F[2,2]*dN[3,i] + F[2,3]*dN[2,i]
-            BL[5, 3*(i-1)+3] = F[3,2]*dN[3,i] + F[3,3]*dN[2,i]
-            BL[6, 3*(i-1)+1] = F[1,3]*dN[1,i] + F[1,1]*dN[3,i]
-            BL[6, 3*(i-1)+2] = F[2,3]*dN[1,i] + F[2,1]*dN[3,i]
-            BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i]
-        end
+            F = eye(dim)
+            if props.finite_strain
+                F += gradu
+                strain += 1/2*gradu'*gradu
+            end
 
-        strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]]
+            # material stiffness start
 
-        # calculate stress
-        E = element("youngs modulus", ip, time)
-        nu = element("poissons ratio", ip, time)
-        D = E/((1.0+nu)*(1.0-2.0*nu)) * [
-            1.0-nu nu nu 0.0 0.0 0.0
-            nu 1.0-nu nu 0.0 0.0 0.0
-            nu nu 1.0-nu 0.0 0.0 0.0
-            0.0 0.0 0.0 0.5-nu 0.0 0.0
-            0.0 0.0 0.0 0.0 0.5-nu 0.0
-            0.0 0.0 0.0 0.0 0.0 0.5-nu]
+            fill!(BL, 0.0)
+            for i=1:nnodes
+                BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i]
+                BL[1, 3*(i-1)+2] = F[2,1]*dN[1,i]
+                BL[1, 3*(i-1)+3] = F[3,1]*dN[1,i]
+                BL[2, 3*(i-1)+1] = F[1,2]*dN[2,i]
+                BL[2, 3*(i-1)+2] = F[2,2]*dN[2,i]
+                BL[2, 3*(i-1)+3] = F[3,2]*dN[2,i]
+                BL[3, 3*(i-1)+1] = F[1,3]*dN[3,i]
+                BL[3, 3*(i-1)+2] = F[2,3]*dN[3,i]
+                BL[3, 3*(i-1)+3] = F[3,3]*dN[3,i]
+                BL[4, 3*(i-1)+1] = F[1,1]*dN[2,i] + F[1,2]*dN[1,i]
+                BL[4, 3*(i-1)+2] = F[2,1]*dN[2,i] + F[2,2]*dN[1,i]
+                BL[4, 3*(i-1)+3] = F[3,1]*dN[2,i] + F[3,2]*dN[1,i]
+                BL[5, 3*(i-1)+1] = F[1,2]*dN[3,i] + F[1,3]*dN[2,i]
+                BL[5, 3*(i-1)+2] = F[2,2]*dN[3,i] + F[2,3]*dN[2,i]
+                BL[5, 3*(i-1)+3] = F[3,2]*dN[3,i] + F[3,3]*dN[2,i]
+                BL[6, 3*(i-1)+1] = F[1,3]*dN[1,i] + F[1,1]*dN[3,i]
+                BL[6, 3*(i-1)+2] = F[2,3]*dN[1,i] + F[2,1]*dN[3,i]
+                BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i]
+            end
 
-        element_keys = get_keys(element)
+            strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]]
 
-            if "plasticity" in element_keys
-                plastic_def = element("plasticity")[ip.id]
+            # calculate stress
+            E = element("youngs modulus", ip, time)
+            nu = element("poissons ratio", ip, time)
+            D = E/((1.0+nu)*(1.0-2.0*nu)) * [
+                1.0-nu nu nu 0.0 0.0 0.0
+                nu 1.0-nu nu 0.0 0.0 0.0
+                nu nu 1.0-nu 0.0 0.0 0.0
+                0.0 0.0 0.0 0.5-nu 0.0 0.0
+                0.0 0.0 0.0 0.0 0.5-nu 0.0
+                0.0 0.0 0.0 0.0 0.0 0.5-nu]
 
-                calculate_stress! = plastic_def["type"]
-                yield_surface_ = plastic_def["yield_surface"]
-                params = plastic_def["params"]
+            element_keys = get_keys(element)
 
-                initialize_internal_params!(params, ip, Val{:type_3d})
+                if "plasticity" in element_keys
+                    plastic_def = element("plasticity")[ip.id]
 
-                if time == 0.0
-                    error("Given step time = $(time). Please select time > 0.0")
+                    calculate_stress! = plastic_def["type"]
+                    yield_surface_ = plastic_def["yield_surface"]
+                    params = plastic_def["params"]
+
+                    initialize_internal_params!(params, ip, Val{:type_3d})
+
+                    if time == 0.0
+                        error("Given step time = $(time). Please select time > 0.0")
+                    end
+
+                    t_last = ip("prev_time", time)
+                    update!(ip, "prev_time", time => t_last)
+
+                    dt = time - t_last
+
+                    stress_last = ip("stress", t_last)
+                    strain_last = ip("strain", t_last)
+
+                    dstrain_vec = strain_vec - strain_last
+                    stress_vec = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+                    plastic_strain = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+                    Dtan = [0.0 0.0 0.0 0.0 0.0 0.0;
+                            0.0 0.0 0.0 0.0 0.0 0.0;
+                            0.0 0.0 0.0 0.0 0.0 0.0
+                            0.0 0.0 0.0 0.0 0.0 0.0;
+                            0.0 0.0 0.0 0.0 0.0 0.0;
+                            0.0 0.0 0.0 0.0 0.0 0.0]
+                    calculate_stress!(stress_vec, stress_last, dstrain_vec, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d})
+
+
+                # plastic_def = element.dev["plasticity"]
+                # calculate_stress! = plastic_def["stress"]
+                # params = plastic_def["params"]
+                # yield_surface_ = plastic_def["yield_surface"]
+                # (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:type_3d})
+                # dstrain_vec = strain_vec - strain_last
+
+                # calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, Val{:type_3d})
+            else
+
+                stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec)
+                Dtan = D
+            end
+
+            :strain in props.store_fields && update!(ip, "strain", time => strain_vec)
+            :stress in props.store_fields && update!(ip, "stress", time => stress_vec)
+            :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1])
+            :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2])
+            :stress33 in props.store_fields && update!(ip, "stress33", time => stress_vec[3])
+            :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[4])
+            :stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5])
+            :stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6])
+            :plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain)
+
+            Km += w*BL'*Dtan*BL
+            # material stiffness end
+
+            if props.geometric_stiffness
+                # take geometric stiffness into account
+
+                fill!(BNL, 0.0)
+
+                for i=1:size(dN, 2)
+                    BNL[1, 3*(i-1)+1] = dN[1,i]
+                    BNL[2, 3*(i-1)+1] = dN[2,i]
+                    BNL[3, 3*(i-1)+1] = dN[3,i]
+                    BNL[4, 3*(i-1)+2] = dN[1,i]
+                    BNL[5, 3*(i-1)+2] = dN[2,i]
+                    BNL[6, 3*(i-1)+2] = dN[3,i]
+                    BNL[7, 3*(i-1)+3] = dN[1,i]
+                    BNL[8, 3*(i-1)+3] = dN[2,i]
+                    BNL[9, 3*(i-1)+3] = dN[3,i]
                 end
 
-                t_last = ip("prev_time", time)
-                update!(ip, "prev_time", time => t_last)
+                S3 = zeros(3*dim, 3*dim)
+                S3[1,1] = stress_vec[1]
+                S3[2,2] = stress_vec[2]
+                S3[3,3] = stress_vec[3]
+                S3[1,2] = S3[2,1] = stress_vec[4]
+                S3[2,3] = S3[3,2] = stress_vec[5]
+                S3[1,3] = S3[3,1] = stress_vec[6]
+                S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3]
 
-                dt = time - t_last
+                Kg += w*BNL'*S3*BNL
 
-                stress_last = ip("stress", t_last)
-                strain_last = ip("strain", t_last)
+            end
 
-                dstrain_vec = strain_vec - strain_last
-                stress_vec = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-                plastic_strain = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-                Dtan = [0.0 0.0 0.0 0.0 0.0 0.0;
-                        0.0 0.0 0.0 0.0 0.0 0.0;
-                        0.0 0.0 0.0 0.0 0.0 0.0
-                        0.0 0.0 0.0 0.0 0.0 0.0;
-                        0.0 0.0 0.0 0.0 0.0 0.0;
-                        0.0 0.0 0.0 0.0 0.0 0.0]
-                calculate_stress!(stress_vec, stress_last, dstrain_vec, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d})
+            # external load start
 
+            if haskey(element, "displacement load")
+                T = element("displacement load", ip, time)
+                f += w*vec(T*N)
+            end
 
-            # plastic_def = element.dev["plasticity"]
-            # calculate_stress! = plastic_def["stress"]
-            # params = plastic_def["params"]
-            # yield_surface_ = plastic_def["yield_surface"]
-            # (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:type_3d})
-            # dstrain_vec = strain_vec - strain_last
+            for i=1:dim
+                if haskey(element, "displacement load $i")
+                    b = element("displacement load $i", ip, time)
+                    f[i:dim:end] += w*vec(b*N)
+                end
+            end
 
-            # calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, Val{:type_3d})
-        else
+            # external load end
+
+            if get_formulation_type(problem) == :incremental
+                f -= w*BL'*stress_vec
+            end
 
-            stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec)
-            Dtan = D
         end
 
-        :strain in props.store_fields && update!(ip, "strain", time => strain_vec)
-        :stress in props.store_fields && update!(ip, "stress", time => stress_vec)
-        :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1])
-        :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2])
-        :stress33 in props.store_fields && update!(ip, "stress33", time => stress_vec[3])
-        :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[4])
-        :stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5])
-        :stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6])
-        :plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain)
+        gdofs = get_gdofs(problem, element)
 
-        Km += w*BL'*Dtan*BL
-        # material stiffness end
+        # add contributions to K, Kg, f
+        add!(assembly.K, gdofs, gdofs, Km)
 
         if props.geometric_stiffness
-            # take geometric stiffness into account
-
-            fill!(BNL, 0.0)
-
-            for i=1:size(dN, 2)
-                BNL[1, 3*(i-1)+1] = dN[1,i]
-                BNL[2, 3*(i-1)+1] = dN[2,i]
-                BNL[3, 3*(i-1)+1] = dN[3,i]
-                BNL[4, 3*(i-1)+2] = dN[1,i]
-                BNL[5, 3*(i-1)+2] = dN[2,i]
-                BNL[6, 3*(i-1)+2] = dN[3,i]
-                BNL[7, 3*(i-1)+3] = dN[1,i]
-                BNL[8, 3*(i-1)+3] = dN[2,i]
-                BNL[9, 3*(i-1)+3] = dN[3,i]
-            end
-
-            S3 = zeros(3*dim, 3*dim)
-            S3[1,1] = stress_vec[1]
-            S3[2,2] = stress_vec[2]
-            S3[3,3] = stress_vec[3]
-            S3[1,2] = S3[2,1] = stress_vec[4]
-            S3[2,3] = S3[3,2] = stress_vec[5]
-            S3[1,3] = S3[3,1] = stress_vec[6]
-            S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3]
-
-            Kg += w*BNL'*S3*BNL
-
+            add!(assembly.Kg, gdofs, gdofs, Kg)
         end
 
-        # external load start
-
-        if haskey(element, "displacement load")
-            T = element("displacement load", ip, time)
-            f += w*vec(T*N)
-        end
-
-        for i=1:dim
-            if haskey(element, "displacement load $i")
-                b = element("displacement load $i", ip, time)
-                f[i:dim:end] += w*vec(b*N)
-            end
-        end
-
-        # external load end
-
-        if get_formulation_type(problem) == :incremental
-            f -= w*BL'*stress_vec
-        end
+        add!(assembly.f, gdofs, f)
 
     end
 
-    return Km, Kg, f
+    return nothing
 end
 
 """ Elasticity equations, surface traction for continuum formulation. """
-function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
+function assemble!{El<:Elasticity3DSurfaceElements}(assembly::Assembly,
+                                                    problem::Problem{Elasticity},
+                                                    elements::Vector{Element{El}},
+                                                    time, ::Type{Val{:continuum}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
-    nnodes = size(element, 2)
-    Km = zeros(dim*nnodes, dim*nnodes)
-    Kg = zeros(dim*nnodes, dim*nnodes)
-    f = zeros(dim*nnodes)
 
-    has_concentrated_forces = false
-    for ip in get_integration_points(element)
-        detJ = element(ip, time, Val{:detJ})
-        w = ip.weight*detJ
-        N = element(ip, time)
-        if haskey(element, "displacement traction force")
-            T = element("displacement traction force", ip, time)
-            f += w*vec(T*N)
-        end
-        for i in 1:dim
-            if haskey(element, "displacement traction force $i")
-                T = element("displacement traction force $i", ip, time)
-                f[i:dim:end] += w*vec(T*N)
+    for element in elements
+        nnodes = size(element, 2)
+        f = zeros(dim*nnodes)
+
+        has_concentrated_forces = false
+        for ip in get_integration_points(element)
+            detJ = element(ip, time, Val{:detJ})
+            w = ip.weight*detJ
+            N = element(ip, time)
+            if haskey(element, "displacement traction force")
+                T = element("displacement traction force", ip, time)
+                f += w*vec(T*N)
             end
-            if haskey(element, "concentrated force $i")
-                has_concentrated_forces = true
-                T = element("concentrated force $i", ip, time)
-                f[i:dim:end] += w*vec(T*N)
+            for i in 1:dim
+                if haskey(element, "displacement traction force $i")
+                    T = element("displacement traction force $i", ip, time)
+                    f[i:dim:end] += w*vec(T*N)
+                end
+                if haskey(element, "concentrated force $i")
+                    has_concentrated_forces = true
+                    T = element("concentrated force $i", ip, time)
+                    f[i:dim:end] += w*vec(T*N)
+                end
+            end
+            if haskey(element, "surface pressure")
+                J = element(ip, time, Val{:Jacobian})'
+                n = cross(J[:,1], J[:,2])
+                n /= norm(n)
+                # sign convention, positive pressure is towards surface
+                p = -element("surface pressure", ip, time)
+                f += w*p*vec(n*N)
             end
         end
-        if haskey(element, "surface pressure")
-            J = element(ip, time, Val{:Jacobian})'
-            n = cross(J[:,1], J[:,2])
-            n /= norm(n)
-            # sign convention, positive pressure is towards surface
-            p = -element("surface pressure", ip, time)
-            f += w*p*vec(n*N)
+        if has_concentrated_forces
+            update!(element, "concentrated force", time => Any[f])
         end
+
+        gdofs = get_gdofs(problem, element)
+        add!(assembly.f, gdofs, f)
+
     end
-    if has_concentrated_forces
-        update!(element, "concentrated force", time => Any[f])
-    end
-    return Km, Kg, f
 end
 
 """ Return strain tensor. """
diff --git a/src/problems_elasticity_2d.jl b/src/problems_elasticity_2d.jl
index 3eb8559..2be26f9 100644
--- a/src/problems_elasticity_2d.jl
+++ b/src/problems_elasticity_2d.jl
@@ -4,7 +4,9 @@
 const Elasticity2DSurfaceElements = Union{Poi1,Seg2,Seg3}
 const Elasticity2DVolumeElements = Union{Tri3,Tri6,Quad4,Quad8,Quad9}
 
-function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{Element}, time::Real, ::Type{Val{:plane_stress}})
+function assemble!{T}(assembly::Assembly, problem::Problem{Elasticity},
+                      elements::Union{Vector{Element}, Vector{Element{T}}},
+                      time, ::Type{Val{:plane_stress}})
     for element in elements
         gdofs = get_gdofs(problem, element)
         Km, Kg, f = assemble(problem, element, time, Val{:plane})
@@ -14,7 +16,9 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::V
     end
 end
 
-function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{Element}, time::Real, ::Type{Val{:plane_strain}})
+function assemble!{T}(assembly::Assembly, problem::Problem{Elasticity},
+                      elements::Union{Vector{Element}, Vector{Element{T}}},
+                      time, ::Type{Val{:plane_strain}})
     for element in elements
         gdofs = get_gdofs(problem, element)
         Km, Kg, f = assemble(problem, element, time, Val{:plane})
@@ -25,7 +29,9 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::V
 end
 
 """ Plane elasticity equations (plane stress, plane strain). """
-function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}})
+function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
+                                                  element::Element{El}, time,
+                                                  ::Type{Val{:plane}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
@@ -178,7 +184,9 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
     return Km, Kg, f
 end
 
-function assemble{El<:Elasticity2DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:plane}})
+function assemble{El<:Elasticity2DSurfaceElements}(problem::Problem{Elasticity},
+                                                   element::Element{El},
+                                                   time, ::Type{Val{:plane}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)