From 7cc14c57c8eb4d3e00def61115dbc9b1ea984614 Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Thu, 11 Feb 2016 02:51:27 +0200
Subject: [PATCH] minor fixes to incremental formulation

---
 src/problems.jl                      |  5 +++-
 src/solvers.jl                       | 14 +++++++++---
 test/test_elasticity_surface_load.jl | 34 +++++++++++++++++-----------
 3 files changed, 36 insertions(+), 17 deletions(-)

diff --git a/src/problems.jl b/src/problems.jl
index 0c5931f..f56ec07 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -169,9 +169,11 @@ function update_assembly!(problem, u, la)
     assembly.u_prev = copy(assembly.u)
     assembly.la_prev = copy(assembly.la)
     if get_formulation_type(problem) == :incremental
+        info("incremental formulation, adding increment to solution vector")
         assembly.u += u
         assembly.la += la
     else
+        info("total formulation, replacing solution vector with new values")
         assembly.u = u
         assembly.la = la
     end
@@ -179,7 +181,8 @@ function update_assembly!(problem, u, la)
     # calculate change of norm
     assembly.u_norm_change = norm(assembly.u - assembly.u_prev)
     assembly.la_norm_change = norm(assembly.la - assembly.la_prev)
-    return assembly.u_norm_change, assembly.la_norm_change
+    #return assembly.u_norm_change, assembly.la_norm_change
+    return assembly.u, assembly.la
 end
 
 """ Update solutions to elements.
diff --git a/src/solvers.jl b/src/solvers.jl
index dbc7b1a..f73bdec 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -348,16 +348,24 @@ Notes
 -----
 Default convergence criteria is obtained by checking each sub-problem convergence.
 """
-function has_converged(solver::Solver)
+function has_converged(solver::Solver; check_convergence_for_boundary_problems=false)
     converged = true
     eps = solver.nonlinear_system_convergence_tolerance
     for problem in solver.problems
         has_converged = true
         if is_field_problem(problem)
             has_converged = problem.assembly.u_norm_change/norm(problem.assembly.u) < eps
+            info("Details for problem $(problem.name)")
+            info("Norm: $(norm(problem.assembly.u))")
+            info("Norm change: $(problem.assembly.u_norm_change)")
+            info("Has converged? $(has_converged)")
         end
-        if is_boundary_problem(problem)
+        if is_boundary_problem(problem) && check_convergence_for_boundary_problems
             has_converged = problem.assembly.la_norm_change/norm(problem.assembly.la) < eps
+            info("Details for problem $(problem.name)")
+            info("Norm: $(norm(problem.assembly.la))")
+            info("Norm change: $(problem.assembly.la_norm_change)")
+            info("Has converged? $(has_converged)")
         end
         converged &= has_converged
     end
@@ -394,7 +402,7 @@ function call(solver::Solver)
 
         # 2.3 update solution back to elements
         for problem in solver.problems
-            update_assembly!(problem, u, la)
+            u, la = update_assembly!(problem, u, la)
             update_elements!(problem, u, la)
         end
 
diff --git a/test/test_elasticity_surface_load.jl b/test/test_elasticity_surface_load.jl
index 3f27743..f8afcc7 100644
--- a/test/test_elasticity_surface_load.jl
+++ b/test/test_elasticity_surface_load.jl
@@ -4,8 +4,7 @@
 using JuliaFEM.Test
 using JuliaFEM.Core: Node, update!, Quad4, Seg2, Problem, Elasticity, Solver, Dirichlet
 
-@testset "test 2d linear elasticity with surface load." begin
-
+function get_test_problem()
     nodes = Dict{Int64, Node}(
         1 => [0.0, 0.0],
         2 => [1.0, 0.0],
@@ -39,16 +38,25 @@ using JuliaFEM.Core: Node, update!, Quad4, Seg2, Problem, Elasticity, Solver, Di
     # type, name, dimension, unknown_field_name
     boundary_problem = Problem(Dirichlet, "symmetry boundaries", 2, "displacement")
     push!(boundary_problem, sym13, sym23)
-
-    solver = Solver("solve block problem")
-    push!(solver, elasticity_problem)
-    push!(solver, boundary_problem)
-
-    call(solver)
-
-    u_disp = element1("displacement", [1.0, 1.0], 0.0)
-    info("Displacement = $u_disp")
-
-    @test isapprox(u_disp, expected)
+    return elasticity_problem, boundary_problem
+end
+
+@testset "test 2d linear with surface load." begin
+
+    E = 288.0
+    nu = 1.0/3.0
+    f = -E/10.0
+    expected = f/E*[-nu, 1]
+    elasticity_problem, boundary_problem = get_test_problem()
+    solver = Solver("solve block problem")
+    #solver.is_linear_system = true
+    push!(solver, elasticity_problem)
+    push!(solver, boundary_problem)
+    call(solver)
+    element1 = elasticity_problem.elements[1]
+    u_disp = element1("displacement", [1.0, 1.0], 0.0)
+    info("Displacement = $u_disp")
+    @test isapprox(u_disp, expected)
+
 end