Fix + example (#3)

* Fix updating fields when nonzero rows in matrix (after elimination of 
boundary conditions)
* Example, transient response of long bar under self-weight

examples/block.jl

@@ -0,0 +1,96 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/LinearImplicitDynamics.jl/blob/master/LICENSE
+
+using JuliaFEM
+using JuliaFEM.Preprocess
+add_elements! = JuliaFEM.add_elements!
+using LinearImplicitDynamics
+using Logging
+Logging.configure(level=DEBUG)
+
+## Generating structured block mesh
+
+# TODO: this should be in some package!
+function quad_mesh(nel_x=10, nel_y=10; lx=1.0, ly=1.0)
+    nnodes_x = nel_x+1
+    nnodes_y = nel_y+1
+    nnode = nnodes_x*nnodes_y
+    nodemap = reshape(1:nnode, nnodes_x, nnodes_y)
+    nodes_1 = vec(nodemap[1:nnodes_x-1, 1:nnodes_y-1])
+    nodes_2 = vec(nodemap[2:nnodes_x, 1:nnodes_y-1])
+    nodes_3 = vec(nodemap[2:nnodes_x, 2:nnodes_y])
+    nodes_4 = vec(nodemap[1:nnodes_x-1, 2:nnodes_y])
+
+    mesh = Mesh()
+
+    nid = 1
+    for y in linspace(0, ly, nnodes_y)
+        for x in linspace(0, lx, nnodes_x)
+            add_node!(mesh, nid, [x, y])
+            add_node_to_node_set!(mesh, :NALL, nid)
+            nid += 1
+        end
+    end
+
+    elid = 1
+    for c in zip(nodes_1, nodes_2, nodes_3, nodes_4)
+        add_element!(mesh, elid, :Quad4, collect(c))
+        add_element_to_element_set!(mesh, :EALL, elid)
+        elid += 1
+    end
+
+    return mesh
+end
+
+nel_x = 10
+nel_y = 3
+mesh = quad_mesh(nel_x, nel_y; lx=5000.0, ly=90.0)
+for nid=1:nel_x+1:(nel_x+1)*(nel_y+1)
+    add_node_to_node_set!(mesh, :LEFT, nid)
+end
+
+block_elements = create_elements(mesh, "EALL")
+update!(block_elements, "youngs modulus", 210.0e3)
+update!(block_elements, "poissons ratio", 0.3)
+update!(block_elements, "density", 7.80e-7)
+update!(block_elements, "displacement load 2", 7.80e-7*9.81)
+block = Problem(Elasticity, "10x1 block", 2)
+block.properties.formulation = :plane_stress
+add_elements!(block, block_elements)
+
+fixed = Problem(Dirichlet, "fixed", 2, "displacement")
+fixed_nodes = mesh.node_sets[:LEFT]
+fixed_elements = [Element(Poi1, [j]) for j in fixed_nodes]
+update!(fixed_elements, "geometry", mesh.nodes)
+for j=1:2
+    update!(fixed_elements, "displacement $j", 0.0)
+end
+add_elements!(fixed, fixed_elements)
+
+xdmf = Xdmf("block_resulst"; overwrite=true)
+
+# analysis = Analysis(Nonlinear, "modal analysis of block under self-weight")
+# add_results_writer!(analysis, xdmf)
+# add_problems!(analysis, [block, fixed])
+# run!(analysis)
+
+# analysis = Analysis(Modal, "modal analysis of block under self-weight")
+# add_results_writer!(analysis, xdmf)
+# add_problems!(analysis, [block, fixed])
+# run!(analysis)
+
+analysis = Analysis(LinearImplicit, "transient analysis of block under self-weight")
+add_results_writer!(analysis, xdmf)
+add_problems!(analysis, [block, fixed])
+run!(analysis)
+
+close(xdmf.hdf)
+
+max_el_id = maximum(keys(mesh.elements))
+element = first(filter(e -> e.id == max_el_id, block_elements))
+disp = element.fields["displacement"].data
+t = [d.first for d in disp]
+u = [d.second[end][2] for d in disp]
+
+# using Plots
+# plot(t[1:10:end], u[1:10:end])

src/LinearImplicitDynamics.jl

@@ -63,23 +63,30 @@ function get_global_matrices(analysis::Analysis{LinearImplicit})
     K = sparse(K, dim, dim) + sparse(Kg, dim, dim)
     f = sparse(f, dim, 1) + sparse(fg, dim, 1)
 
-    for problem in get_problems(analysis)
-        is_boundary_problem(problem) || continue
-        eliminate_boundary_conditions!(problem, K, M, f)
-    end
-
-    K = 1/2*(K + K')
-    M = 1/2*(M + M')
-    p.K = dropzeros(K)
-    p.M = dropzeros(M)
+    p.K = 1/2*(K + K')
+    p.M = 1/2*(M + M')
     p.f = f
     return p.K, p.M, p.f
 end
 
 function FEMBase.run!(analysis::Analysis{LinearImplicit})
 
+    # Preprocess: create global matrices, eliminate boundary conditions, remove
+    # any nonzero rows and factorize
+
     K, M, f = get_global_matrices(analysis)
     ndofs = size(K, 1)
+    full_vec = zeros(ndofs) # for temporary operations
+    dim = get_unknown_field_dimension(first(get_problems(analysis)))
+    nnodes = Int(ndofs/dim)
+
+    for problem in get_problems(analysis)
+        is_boundary_problem(problem) || continue
+        eliminate_boundary_conditions!(problem, K, M, f)
+    end
+    dropzeros!(K)
+    dropzeros!(M)
+
     if isempty(analysis.properties.u0)
         analysis.properties.u0 = spzeros(ndofs)
     end
@@ -96,6 +103,9 @@ function FEMBase.run!(analysis::Analysis{LinearImplicit})
     u0 = u0[nz]
     v0 = v0[nz]
     M_fact = cholfact(M)
+    ndofs = size(K, 1)
+
+    # Solution
 
     function model(dx, x, p, t)
         info("Solving at time $t")
@@ -106,25 +116,26 @@ function FEMBase.run!(analysis::Analysis{LinearImplicit})
         dx[ndofs+1:end] = M_fact \ full(f - K*u)
     end
 
-    ndofs = size(K, 1)
     x0 = [u0; v0]
     prob = ODEProblem(model, x0, analysis.properties.tspan)
     sol = analysis.properties.sol = solve(prob)
 
-    # FIXME: support for variable number of dofs / node
-    dim = get_unknown_field_dimension(first(get_problems(analysis)))
-    nnodes = Int(ndofs/dim)
-
-    info("Number of time steps = ", length(sol))
+    # postprocess, update displacement and velocity back to elements
 
+    # FIXME: support for variable number of dofs / node
     function to_dict(u, dim, nnodes)
         return Dict(j => [u[dim*(j-1)+k] for k=1:dim] for j=1:nnodes)
     end
 
+    info("Number of nodes = $nnodes, number of time steps = $(length(sol))")
+
     for (time, x) in zip(sol.t, sol.u)
-        u = to_dict(x[1:ndofs], dim, nnodes)
-        v = to_dict(x[ndofs+1:end], dim, nnodes)
+        full_vec[nz] = x[1:ndofs]
+        u = to_dict(full_vec, dim, nnodes)
+        full_vec[nz] = x[ndofs+1:end]
+        v = to_dict(full_vec, dim, nnodes)
         for problem in get_problems(analysis)
+            is_field_problem(problem) || continue
             for elements in get_elements(problem)
                 update!(elements, "displacement", time => u)
                 update!(elements, "velocity", time => v)