Merge pull request #2 from JuliaFEM/initial_implementation

Initial implementation of dynamics solver using DifferentialEquations.jl

REQUIRE

@@ -1,2 +1,4 @@
 julia 0.6
 FEMBase
+Reexport
+DifferentialEquations

examples/ball.jl

@@ -0,0 +1,139 @@
+# 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
+using JuliaFEM.Postprocess
+using Logging
+Logging.configure(level=INFO)
+add_elements! = JuliaFEM.add_elements!
+using LinearImplicitDynamics
+
+# Model construction starts
+
+datadir = Pkg.dir("LinearImplicitDynamics", "examples", "ball")
+meshfile = joinpath(datadir, "ball.med")
+mesh = aster_read_mesh(meshfile)
+mesh.element_sets[:BALL] = bset = Set{Int64}()
+for (elid, eltype) in mesh.element_types
+    if eltype == :Tet4
+        push!(bset, elid)
+    end
+end
+
+for (elset_name, element_ids) in mesh.element_sets
+    nel = length(element_ids)
+    println("Element set $elset_name contains $nel elements.")
+end
+
+for (nset_name, node_ids) in mesh.node_sets
+    nno = length(node_ids)
+    println("Node set $nset_name contains $nno nodes.")
+end
+
+nnodes = length(mesh.nodes)
+println("Total number of nodes in mesh: $nnodes")
+nelements = length(mesh.elements)
+println("Total number of elements in mesh: $nelements")
+
+ball_elements = create_elements(mesh, "BALL")
+nel_ball_elements = length(ball_elements)
+println("ball contains $nel_ball_elements elements")
+update!(ball_elements, "youngs modulus", 288.0)
+update!(ball_elements, "poissons ratio", 1/3)
+update!(ball_elements, "density", 36.0e-3)
+#update!(ball_elements, "displacement load 1", 0.20)
+update!(ball_elements, "displacement load 2", 0.20)
+
+xmax = -Inf
+nid = 0
+for (j, coord) in mesh.nodes
+    x, y, z = coord
+    if x > xmax
+        xmax = x
+        nid = j
+    end
+end
+println("max_x node = $nid, xval = $xmax")
+
+load_element = Element(Poi1, [nid])
+update!(load_element, "geometry", mesh.nodes)
+
+V_ball = 0.0
+m_ball = 0.0
+time = 0.0
+for element in ball_elements
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        rho = element("density", ip, time)
+        V_ball += ip.weight * detJ
+        m_ball += ip.weight * rho * detJ
+    end
+end
+println("Ball volume: $V_ball")
+println("Ball mass: $m_ball")
+
+ball = Problem(Elasticity, "ball", 3)
+add_elements!(ball, ball_elements)
+update!(load_element, "displacement traction force z", 10e4)
+
+# Model construction end.
+
+analysis = Analysis(LinearImplicit)
+analysis.properties.tspan = (0.0, 10.0)
+add_problems!(analysis, [ball])
+xdmf = Xdmf(joinpath(datadir, "ball_results"); overwrite=true)
+add_results_writer!(analysis, xdmf)
+
+# Start analysis.
+
+run!(analysis)
+
+# Analysis ready
+
+# sol = analysis.properties.sol
+#
+# dim = 3
+# nnodes = length(mesh.nodes)
+# ndofs = dim * nnodes
+#
+# u = hcat([x[1:ndofs] for x in sol.u]...)
+# v = hcat([x[ndofs+1:end] for x in sol.u]...)
+#
+# tspan = analysis.properties.tspan
+# t0, t1 = tspan
+# nsteps = length(sol.u)
+# println("Number of time steps = $nsteps")
+# t = linspace(t0, t1, nsteps)
+#
+# for (time, x) in zip(t, sol.u)
+#     u = x[1:ndofs]
+#     v = x[ndofs+1:end]
+#     ur = reshape(u, dim, nnodes)
+#     vr = reshape(v, dim, nnodes)
+#     ud = Dict(j => ur[:,j] for j=1:nnodes)
+#     vd = Dict(j => vr[:,j] for j=1:nnodes)
+#     update!(ball_elements, "displacement", time => ud)
+#     update!(ball_elements, "velocity", time => vd)
+# end
+
+# nid2 = find_nearest_node(mesh, [0.0, 0.0, 0.0])
+
+# ball("displacement", 0.0)
+# trajectory = [ball("displacement", time)[nid2] for time in t]
+# traj_x = [t[1] for t in trajectory]
+# traj_y = [t[2] for t in trajectory]
+# traj_z = [t[3] for t in trajectory]
+# xmin, xmax = minimum(traj_x), maximum(traj_x)
+# ymin, ymax = minimum(traj_y), maximum(traj_y)
+# zmin, zmax = minimum(traj_z), maximum(traj_z)
+# println("bounding box x = ($xmin, $xmax)")
+# println("bounding box y = ($ymin, $ymax)")
+# println("bounding box z = ($zmin, $zmax)")
+#
+# for time in linspace(0.0, 10.0, 600)
+#     JuliaFEM.write_results!(analysis, time)
+# end
+close(xdmf.hdf) # src
+
+#

src/LinearImplicitDynamics.jl

@@ -4,17 +4,140 @@
 """ Linear, implicit dynamics solver for JuliaFEM. """
 module LinearImplicitDynamics
 
-using FEMBase
+using Reexport
+@reexport using FEMBase
+using DifferentialEquations
 
-import FEMBase: solve!
+type LinearImplicit <: AbstractAnalysis
+    tspan :: Tuple{Float64,Float64}
+    sol
+    K :: SparseMatrixCSC{Float64}
+    M :: SparseMatrixCSC{Float64}
+    f :: SparseVector{Float64}
+    u0 :: SparseVector{Float64}
+    v0 :: SparseVector{Float64}
+end
+
+function LinearImplicit()
+    tspan = (0.0, 1.0)
+    return LinearImplicit(tspan, nothing, spzeros(0,0), spzeros(0,0),
+                          spzeros(0), spzeros(0), spzeros(0))
+end
+
+"""
+    get_global_matrices(analysis)
+
+Assemble global matrices for problem.
+"""
+function get_global_matrices(analysis::Analysis{LinearImplicit})
+    p = analysis.properties
+    if !isempty(p.K) && !isempty(p.M)
+        return p.K, p.M, p.f
+    end
+
+    time = first(analysis.properties.tspan)
+
+    for problem in get_problems(analysis)
+        assemble!(problem, time)
+        is_field_problem(problem) && assemble_mass_matrix!(problem, time)
+    end
+
+    M = SparseMatrixCOO()
+    K = SparseMatrixCOO()
+    Kg = SparseMatrixCOO()
+    f = SparseMatrixCOO()
+    fg = SparseMatrixCOO()
+
+    for problem in get_problems(analysis)
+        is_field_problem(problem) || continue
+        append!(M, problem.assembly.M)
+        append!(K, problem.assembly.K)
+        append!(Kg, problem.assembly.Kg)
+        append!(f, problem.assembly.f)
+        append!(fg, problem.assembly.fg)
+    end
+
+    dim = size(K, 1)
+
+    M = sparse(M, dim, dim)
+    K = sparse(K, dim, dim) + sparse(Kg, dim, dim)
+    f = sparse(f, dim, 1) + sparse(fg, dim, 1)
 
-type LinearImplicit <: AbstractSolver
+    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.f = f
+    return p.K, p.M, p.f
 end
 
-function solve!(solver::Solver{LinearImplicit}, time)
-    # solve everything
+function FEMBase.run!(analysis::Analysis{LinearImplicit})
+
+    K, M, f = get_global_matrices(analysis)
+    ndofs = size(K, 1)
+    if isempty(analysis.properties.u0)
+        analysis.properties.u0 = spzeros(ndofs)
+    end
+    if isempty(analysis.properties.v0)
+        analysis.properties.v0 = spzeros(ndofs)
+    end
+    u0 = analysis.properties.u0
+    v0 = analysis.properties.v0
+
+    nz = get_nonzero_rows(K)
+    K = K[nz,nz]
+    M = M[nz,nz]
+    f = f[nz]
+    u0 = u0[nz]
+    v0 = v0[nz]
+    M_fact = cholfact(M)
+
+    function model(dx, x, p, t)
+        info("Solving at time $t")
+        ndofs = round(Int, length(dx)/2)
+        u = x[1:ndofs]
+        v = x[ndofs+1:end]
+        dx[1:ndofs] = v
+        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))
+
+    function to_dict(u, dim, nnodes)
+        return Dict(j => [u[dim*(j-1)+k] for k=1:dim] for j=1:nnodes)
+    end
+
+    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)
+        for problem in get_problems(analysis)
+            for elements in get_elements(problem)
+                update!(elements, "displacement", time => u)
+                update!(elements, "velocity", time => v)
+            end
+        end
+    end
+
 end
 
+# function FEMBase.write_results!(analysis::LinearImplicitDynamics, writer::Xdmf)
+#     # Not implemented yet
+# end
+
 export LinearImplicit
 
 end

test/REQUIRE

@@ -1 +1 @@
-TimerOutputs
+JuliaFEM

test/runtests.jl

@@ -1,14 +1,9 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/LinearImplicitDynamics.jl/blob/master/LICENSE
 
-using FEMBase
 using LinearImplicitDynamics
-
-using Base.Test
+using FEMBase.Test
 
 @testset "LinearImplicitDynamics.jl" begin
-    # Test stiffness matrix
-    K = 1.0
-    K_expected = 1.0
-    @test isapprox(K, K_expected)
+    include("test_dropping_element.jl")
 end

test/test_dropping_element.jl

@@ -0,0 +1,28 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/LinearImplicitDynamics.jl/blob/master/LICENSE
+
+using LinearImplicitDynamics
+using JuliaFEM
+using FEMBase.Test
+
+X = Dict(1 => [0.0, 0.0],
+         2 => [2.0, 0.0],
+         3 => [2.0, 2.0],
+         4 => [0.0, 2.0])
+
+element = Element(Quad4, [1, 2, 3, 4])
+update!(element, "geometry", X)
+update!(element, "youngs modulus", 288.0)
+update!(element, "poissons ratio", 1/3)
+update!(element, "density", 36.0)
+update!(element, "displacement load 2", 36.0*9.81)
+problem = Problem(Elasticity, "2x2 element", 2)
+problem.properties.formulation = :plane_stress
+add_elements!(problem, [element])
+analysis = Analysis(LinearImplicit, "dropping element")
+add_problems!(analysis, [problem])
+run!(analysis)
+midpnt_u = element("displacement", (0.0, 0.0), Inf)
+midpnt_v = element("velocity", (0.0, 0.0), Inf)
+@test isapprox(midpnt_u, [0.0, 1/2*9.81])
+@test isapprox(midpnt_v, sqrt.(2.0*9.81*midpnt_u))