From 74ed32a96c89ef18f02cccb7df4965504c10dc47 Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Thu, 9 Jun 2016 23:17:57 +0300
Subject: [PATCH] better heat example

---
 src/solvers.jl    |  8 ++++++
 test/test_heat.jl | 65 ++++++++++++++++++++++++++++++-----------------
 2 files changed, 49 insertions(+), 24 deletions(-)

diff --git a/src/solvers.jl b/src/solvers.jl
index 6d553a7..671d08d 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -312,14 +312,22 @@ function call(solver::Solver)
         info("Starting nonlinear iteration #$(solver.iteration)")
 
         # 2.1 update linearized assemblies (if needed)
+        info("Assembling problems ...")
+        tic()
         for problem in solver.problems
             problem.assembly.changed = true  # force reassembly
             assemble!(problem, solver.time)
         end
+        t1 = round(toq(), 2)
+        info("Assembled in $t1 seconds.")
 
         # 2.2 call solver for linearized system (default: direct lu factorization)
+        info("Solve linear system ...")
+        tic()
         u, la = solve_linear_system(solver, Val{solver.linear_system_solver})
         push!(solver.norms, (norm(u), norm(la)))
+        t1 = round(toq(), 2)
+        info("Solved Ax = b in $t1 seconds.")
 
         # 2.3 update solution back to elements
         for problem in solver.problems
diff --git a/test/test_heat.jl b/test/test_heat.jl
index 6d3921b..bbfd036 100644
--- a/test/test_heat.jl
+++ b/test/test_heat.jl
@@ -13,48 +13,65 @@ using JuliaFEM.Test
         3 => [1.0,1.0],
         4 => [0.0,1.0])
 
-    # volume element
-    element = Element(Quad4, [1, 2, 3, 4])
+    # define volume element
+    el1 = Element(Quad4, [1, 2, 3, 4])
 
-    update!(element, "geometry", X)
-    update!(element, "temperature thermal conductivity", 6.0)
-    update!(element, "temperature load", [12.0, 12.0, 12.0, 12.0])
-    update!(element, "density", 36.0)
+    update!(el1, "geometry", X)
+    update!(el1, "temperature thermal conductivity", 6.0)
+    update!(el1, "temperature load", 12.0)
+    update!(el1, "density", 36.0)
 
-    # boundary element 
-    boundary_element = Element(Seg2, [1, 2])
-    update!(boundary_element, "geometry", X)
-    # linear ramp from 0 to 6 in time 0 to 1
-    update!(boundary_element, "temperature flux", 0.0 => 0.0, 1.0 => 6.0)
+    # define boundary element for flux
+    el2 = Element(Seg2, [1, 2])
+    update!(el2, "geometry", X)
+    # linear ramp from 0 -> 6 in time 0 -> 1
+    update!(el2, "temperature flux", 0.0 => 0.0, 1.0 => 6.0)
 
+    # define heat problem and push elements to problem
     problem = Problem(Heat, "one element heat problem", 1)
-    push!(problem, element, boundary_element)
+    push!(problem, el1, el2)
 
-    # Set constant source f=12 with k=6. Accurate solution is
-    # T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x)
+    # define boundary element for dirichlet boundary condition
+    el3 = Element(Seg2, [3, 4])
+    update!(el3, "geometry", X)
+    update!(el3, "temperature 1", 0.0)
+
+    boundary_condition = Problem(Dirichlet, "T=0 on top", 1, "temperature")
+    push!(boundary_condition, el3)
+
+    # manual assembling of problem + solution:
     assemble!(problem, 0.0)
     A = full(problem.assembly.K)
     b = full(problem.assembly.f)
-
     A_expected = [
          4.0 -1.0 -2.0 -1.0
         -1.0  4.0 -1.0 -2.0
         -2.0 -1.0  4.0 -1.0
         -1.0 -2.0 -1.0  4.0]
-
-    @test isapprox(A, A_expected)
-
     free_dofs = [1, 2]
+    @test isapprox(A, A_expected)
     @test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [1.0, 1.0])
 
+    # using Solver
+    solver = Solver("solve heat problem")
+    solver.is_linear_system = true
+    push!(solver, problem, boundary_condition)
+
+    # Set constant source f=12 with k=6. Accurate solution is
+    # T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x)
+    # when boundary flux not active (at t=0)
+    solver.time = 0.0
+    call(solver)
+    # interpolate temperature at middle of element 2 (flux boundary) at time t=0:
+    T = el2("temperature", [0.0], 0.0)
+    @test isapprox(T[1], 1.0)
+
     # Set constant flux g=6 on boundary. Accurate solution is
     # u(x,y) = x which equals T=1 on boundary.
     # at time t=1.0 all loads should be on.
-    empty!(problem.assembly)
-    assemble!(problem, 1.0)
-    A = full(problem.assembly.K)
-    b = full(problem.assembly.f)
-    T = A[free_dofs, free_dofs] \ b[free_dofs]
-    @test isapprox(T, [2.0, 2.0])
+    solver.time = 1.0
+    call(solver)
+    T = el2("temperature", [0.0], 1.0)
+    @test isapprox(T[1], 2.0)
 end