From 73c2809e85d6c3946474b373a1858003b6c54839 Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Wed, 24 Jun 2015 21:03:10 +0300
Subject: [PATCH] fixed a bug with calculating jacobian

---
 src/elasticity_solver.jl       |  2 +-
 test/test_elasticity_solver.jl | 38 +++++++++++++++++++++++++++++++---
 2 files changed, 36 insertions(+), 4 deletions(-)

diff --git a/src/elasticity_solver.jl b/src/elasticity_solver.jl
index da7898f..d8a5b11 100644
--- a/src/elasticity_solver.jl
+++ b/src/elasticity_solver.jl
@@ -81,7 +81,7 @@ function calc_local_matrices!(X, u, R, Kt, N, dNdξ, λ_, μ_, ipoints, iweights
     #@debug("Jt:\n",Jᵀ)
     λ = interpolate(λ_, N, ξ)
     μ = interpolate(μ_, N, ξ)
-    Jᵀ = interpolate(X, dNdξ, ξ)'
+    Jᵀ = interpolate(X, dNdξ, ξ)
         detJ = det(Jᵀ)
         ∇N = inv(Jᵀ)*dNdξ(ξ)'
         ∇u = u*∇N'
diff --git a/test/test_elasticity_solver.jl b/test/test_elasticity_solver.jl
index a7c6c8c..76a1fc7 100644
--- a/test/test_elasticity_solver.jl
+++ b/test/test_elasticity_solver.jl
@@ -5,9 +5,7 @@ using Logging
 
 using JuliaFEM.elasticity_solver: solve_elasticity_increment!
 
-
-facts("test solve elasticity increment") do
-
+function one_elem_fixture()
   X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'
   elmap = [1; 2; 3; 4]
   nodalloads = [0 0; 0 0; 0 -2; 0 0]'
@@ -40,6 +38,15 @@ facts("test solve elasticity increment") do
   ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]
   iweights = [1 1 1 1]
 
+    return (X, u, du, elmap, nodalloads, dirichletbc,
+     la, mu, N, dNdξ, ipoints, iweights)
+end
+
+facts("test solve elasticity increment") do
+
+    (X, u, du, elmap, nodalloads, dirichletbc,
+     la, mu, N, dNdξ, ipoints, iweights) = one_elem_fixture()
+
   for i=1:10
     solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc,
                                 la, mu, N, dNdξ, ipoints, iweights)
@@ -54,6 +61,31 @@ facts("test solve elasticity increment") do
 end
 
 
+facts("test solve elasticity increment rot 30") do
+
+    (X, u, du, elmap, nodalloads, dirichletbc,
+     la, mu, N, dNdξ, ipoints, iweights) = one_elem_fixture()
+
+    phi = 30/180*pi
+    rmat = [cos(phi) -sin(phi); sin(phi) cos(phi)]
+    X = rmat*X
+    nodalloads = rmat*nodalloads
+
+  for i=1:10
+    solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc,
+                                la, mu, N, dNdξ, ipoints, iweights)
+    @debug("increment:\n",du)
+    u += du
+    if norm(du) < 1.0e-9
+      break
+    end
+  end
+    u = rmat'*u
+  @debug("solution\n",u)
+  @fact u[2, 3] => roughly(-2.222244754401764)  # Tested against Elmer solution
+end
+
+
 using JuliaFEM.elasticity_solver: interpolate
 facts("test interpolation of different field variables") do
     N(xi) = [