From 6f66810ea3f7d8f7587791e5df64d3d1eb6a2021 Mon Sep 17 00:00:00 2001
From: Jukka Aho <ahojukka5@gmail.com>
Date: Tue, 23 Jun 2015 21:36:04 +0300
Subject: [PATCH] added variable la,mu parameters

---
 src/elasticity_solver.jl       | 19 +++++++++++--------
 test/test_elasticity_solver.jl | 24 +++++++++++++++---------
 2 files changed, 26 insertions(+), 17 deletions(-)

diff --git a/src/elasticity_solver.jl b/src/elasticity_solver.jl
index 69a9c9f..406a527 100644
--- a/src/elasticity_solver.jl
+++ b/src/elasticity_solver.jl
@@ -26,8 +26,8 @@ end
 @doc """
 Calculate local tangent stiffness matrix and residual force vector R = T - F
 """ ->
-function calc_local_matrices!(X, u, R, Kt, dNdξ, λ, μ, ipoints, iweights)
-  dim, N = size(X)
+function calc_local_matrices!(X, u, R, Kt, N, dNdξ, λ_, μ_, ipoints, iweights)
+  dim, nnodes = size(X)
   I = eye(dim)
   R[:,:] = 0.0
   Kt[:,:] = 0.0
@@ -35,8 +35,11 @@ function calc_local_matrices!(X, u, R, Kt, dNdξ, λ, μ, ipoints, iweights)
   dF = zeros(dim, dim)
 
   for m = 1:length(iweights)
-        w = iweights[m]
-        ξ = ipoints[m, :]
+    w = iweights[m]
+    ξ = ipoints[m, :]
+    # interpolate material parameters from element node fields
+    λ = (λ_*N(ξ))[1]
+    μ = (μ_*N(ξ))[1]
         Jᵀ = X*dNdξ(ξ)
         detJ = det(Jᵀ)
         ∇N = inv(Jᵀ)*dNdξ(ξ)'
@@ -47,14 +50,14 @@ function calc_local_matrices!(X, u, R, Kt, dNdξ, λ, μ, ipoints, iweights)
         P = F*S  # PK1 stress tensor
         R[:,:] += w*P*∇N*detJ
 
-        for p = 1:N
+        for p = 1:nnodes
             for i = 1:dim
                 dF[:,:] = 0.0
                 dF[i,:] = ∇N[:,p]
                 dE = 1/2*(F'*dF + dF'*F)
                 dS = λ*trace(dE)*I + 2*μ*dE
                 dP = dF*S + F*dS
-                for q = 1:N
+                for q = 1:nnodes
                     for j = 1:dim
                         Kt[dim*(p-1)+i,dim*(q-1)+j] += w*(dP[j,:]*∇N[:,q])[1]*detJ
                     end
@@ -70,9 +73,9 @@ end
 Solve one increment of elasticity problem
 """ ->
 function solve_elasticity_increment!(X, u, du, R, Kt, elmap, nodalloads,
-                                     dirichletbc, λ, μ, dNdξ, ipoints,
+                                     dirichletbc, λ, μ, N, dNdξ, ipoints,
                                      iweights)
-  calc_local_matrices!(X, u, R, Kt, dNdξ, λ, μ, ipoints, iweights)
+  calc_local_matrices!(X, u, R, Kt, N, dNdξ, λ, μ, ipoints, iweights)
   # FIXME: boundary conditions
   free_dofs = find(isnan(dirichletbc))
   R -= nodalloads
diff --git a/test/test_elasticity_solver.jl b/test/test_elasticity_solver.jl
index f892863..3c5564b 100644
--- a/test/test_elasticity_solver.jl
+++ b/test/test_elasticity_solver.jl
@@ -13,18 +13,25 @@ facts("test solve elasticity increment") do
   dirichletbc = [0 0; NaN NaN; NaN NaN; 0 0]'
 
   E = 90
-  ν = 0.25
-  μ = E/(2*(1+ν))
-  λ = E*ν/((1+ν)*(1-2*ν))
-  λ = 2*λ*μ/(λ + 2*μ)
+  nu = 0.25
+  mu = E/(2*(1+nu))
+  la = E*nu/((1+nu)*(1-2*nu))
+  la = 2*la*mu/(la + 2*mu)
 
-  #E = 90.0*ones(2, 4)
-  #nu = 0.25*ones(2, 4)
+  la = la*ones(1, 4)
+  mu = mu*ones(1, 4)
   u = zeros(2, 4)
   du = zeros(2, 4)
   R = zeros(2,4)
   Kt = zeros(8,8)
 
+  N(xi) = [
+      (1-xi[1])*(1-xi[2])/4
+      (1+xi[1])*(1-xi[2])/4
+      (1+xi[1])*(1+xi[2])/4
+      (1-xi[1])*(1+xi[2])/4
+    ]
+
   dNdξ(ξ) = [-(1-ξ[2])/4.0    -(1-ξ[1])/4.0
               (1-ξ[2])/4.0    -(1+ξ[1])/4.0
               (1+ξ[2])/4.0     (1+ξ[1])/4.0
@@ -34,9 +41,8 @@ facts("test solve elasticity increment") do
   iweights = [1 1 1 1]
 
   for i=1:10
-    JuliaFEM.elasticity_solver.solve_elasticity_increment!(X, u, du, R, Kt, elmap, nodalloads,
-                                dirichletbc, λ, μ, dNdξ, ipoints,
-                                iweights)
+    JuliaFEM.elasticity_solver.solve_elasticity_increment!(
+      X, u, du, R, Kt,elmap, nodalloads, dirichletbc, la, mu, N, dNdξ, ipoints, iweights)
     @debug("increment:\n",du)
     u += du
     if norm(du) < 1.0e-9