minor fixes to incremental formulation

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.

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
 

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