From f48dcf70de0a9dbd2bedcc6bd13e126e51a9c9de Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Sat, 6 Feb 2016 22:00:36 +0200
Subject: [PATCH] continuum elasticity test

---
 test/test_elasticity_continuum.jl | 82 +++++++++++++++++++++++++++++++
 1 file changed, 82 insertions(+)
 create mode 100644 test/test_elasticity_continuum.jl

diff --git a/test/test_elasticity_continuum.jl b/test/test_elasticity_continuum.jl
new file mode 100644
index 0000000..32bd07b
--- /dev/null
+++ b/test/test_elasticity_continuum.jl
@@ -0,0 +1,82 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Test
+
+using JuliaFEM.Core: Quad4, Hex8, Elasticity, Dirichlet, Problem, Node, Solver
+using JuliaFEM.Core: update!
+
+@testset "test simple continuum block with surface traction" begin
+
+    nodes = Dict{Int64, Node}(
+        1 => [0.0, 0.0, 0.0],
+        2 => [1.0, 0.0, 0.0],
+        3 => [1.0, 1.0, 0.0],
+        4 => [0.0, 1.0, 0.0],
+        5 => [0.0, 0.0, 1.0],
+        6 => [1.0, 0.0, 1.0],
+        7 => [1.0, 1.0, 1.0],
+        8 => [0.0, 1.0, 1.0])
+    element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8])
+    element2 = Quad4([5, 6, 7, 8])
+    update!([element1, element2], "geometry", nodes)
+    update!(element1, "youngs modulus", 900.0)
+    update!(element1, "poissons ratio", 0.25)
+    element2["displacement traction force"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]
+
+    problem = Problem(Elasticity, "block", 3)
+    push!(problem, element1, element2)
+
+#=
+    free_dofs = zeros(Bool, 8, 3)
+    x = 1
+    y = 2
+    z = 3
+    free_dofs[2, x] = true
+    free_dofs[3, [x, y]] = true
+    free_dofs[4, y] = true
+    free_dofs[5, z] = true
+    free_dofs[6, [x, z]] = true
+    free_dofs[7, [x, y, z]] = true
+    free_dofs[8, [y, z]] = true
+    free_dofs = find(vec(free_dofs'))
+    info("free dofs: $free_dofs")
+
+    ass = assemble(problem, 0.0)
+    f = full(ass.force_vector)
+    K = full(ass.stiffness_matrix)
+#   info("initial force vector")
+#   dump(reshape(f, 3, 8))
+#   info("initial stiffness matrix")
+#   dump(round(Int, K)[free_dofs, free_dofs])
+    u = zeros(3, 8)
+    u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs]
+    info("result vector")
+    dump(u)
+=#
+
+    dx = Quad4([1, 4, 8, 5])
+    dx["displacement 1"] = 0.0
+    dy = Quad4([1, 5, 6, 2])
+    dy["displacement 2"] = 0.0
+    dz = Quad4([1, 2, 3, 4])
+    dz["displacement 3"] = 0.0
+    bc = Problem(Dirichlet, "symmetries", 3, "displacement")
+    update!([dx, dy, dz], "geometry", nodes)
+    push!(bc, dx, dy, dz)
+
+    solver = Solver()
+    push!(solver, problem, bc)
+    solver()
+
+    X = element1("geometry", [1.0, 1.0, 1.0], 0.0)
+    u = element1("displacement", [1.0, 1.0, 1.0], 0.0)
+    info("displacement at $X = $u")
+
+    # verified using Code Aster.
+    # 2015-12-12-continuum-elasticity/c3d_linear.*
+    # [1/36, 1/36, -1/9]
+    @test isapprox(u, [2.77777777777778E-02, 2.77777777777778E-02, -1.11111111111111E-01])
+end
+