From 8ec056ec6212aef9dd7d75f35ea408140ae10ced Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Tue, 23 Jun 2015 20:41:25 +0300 Subject: [PATCH] added elasticity solver --- src/JuliaFEM.jl | 4 ++ src/elasticity_solver.jl | 83 ++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + test/test_elasticity_solver.jl | 69 ++++++++++++++++------------ 4 files changed, 128 insertions(+), 29 deletions(-) create mode 100644 src/elasticity_solver.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index dbc00d8..e098cab 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -3,6 +3,10 @@ module JuliaFEM VERSION < v"0.4-" && using Docile +# import solvers +include("elasticity_solver.jl") +#using elasticity_solver + export Model, new_model, new_field, get_field, add_nodes, get_nodes, add_elements, get_elements diff --git a/src/elasticity_solver.jl b/src/elasticity_solver.jl new file mode 100644 index 0000000..69a9c9f --- /dev/null +++ b/src/elasticity_solver.jl @@ -0,0 +1,83 @@ +# This file is a part of JuliaFEM. License is MIT: https://github.com/ovainola/JuliaFEM/blob/master/README.md +module elasticity_solver + +VERSION < v"0.4-" && using Docile + +export solve_elasticity_interface! + +@doc """ +This is generic interface that reads data from data model, solves elasticity +problem and updates model. + +Parameters +---------- +model : to be defined +""" -> +function solve_elasticity_interface!() + return 0 +end + + + +# Below this line is internal functions related to solver. They can be used +# directly if needed or using general interface combining data model and +# solver. + +@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) + I = eye(dim) + R[:,:] = 0.0 + Kt[:,:] = 0.0 + + dF = zeros(dim, dim) + + for m = 1:length(iweights) + w = iweights[m] + ξ = ipoints[m, :] + Jᵀ = X*dNdξ(ξ) + detJ = det(Jᵀ) + ∇N = inv(Jᵀ)*dNdξ(ξ)' + ∇u = u*∇N' + F = I + ∇u # Deformation gradient + E = 1/2*(∇u' + ∇u + ∇u'*∇u) # Green-Lagrange strain tensor + S = λ*trace(E)*I + 2*μ*E # PK2 stress tensor + P = F*S # PK1 stress tensor + R[:,:] += w*P*∇N*detJ + + for p = 1:N + 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 j = 1:dim + Kt[dim*(p-1)+i,dim*(q-1)+j] += w*(dP[j,:]*∇N[:,q])[1]*detJ + end + end + end + end + + end +end + + +@doc """ +Solve one increment of elasticity problem +""" -> +function solve_elasticity_increment!(X, u, du, R, Kt, elmap, nodalloads, + dirichletbc, λ, μ, dNdξ, ipoints, + iweights) + calc_local_matrices!(X, u, R, Kt, dNdξ, λ, μ, ipoints, iweights) + # FIXME: boundary conditions + free_dofs = find(isnan(dirichletbc)) + R -= nodalloads + du[free_dofs] = Kt[free_dofs, free_dofs] \ -reshape(R, 8)[free_dofs] +end + + +end diff --git a/test/runtests.jl b/test/runtests.jl index c063cb7..dddfe88 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -59,3 +59,4 @@ end # include("solver_tests/test_elasticity_solver.jl") # include("test_model.jl") +include("test_elasticity_solver.jl") diff --git a/test/test_elasticity_solver.jl b/test/test_elasticity_solver.jl index 3c6c927..f892863 100644 --- a/test/test_elasticity_solver.jl +++ b/test/test_elasticity_solver.jl @@ -1,37 +1,48 @@ -using JuliaFEM -using Base.Test +using JuliaFEM.elasticity_solver -function test_one_element_solution() - model = new_model() +using FactCheck +using Logging +@Logging.configure(level=DEBUG) - # 1. create mesh structure - # add nodes to model - nodes = Dict(1 => [0.0, 0.0], 2 => [10.0, 0.0], 3 => [10.0, 1.0], 4 => [0.0, 1.0]) - add_nodes(model, nodes) - # 2. add elements to model - const QUAD4 = 0x4 - element = Dict("element_type" => QUAD4, "node_ids" => [1, 2, 3, 4]) - elements = Dict() - elements[1] = element - add_elements(model, elements) +facts("test solve elasticity increment") do + X = [0 0; 10 0; 10 1; 0 1]' + elmap = [1; 2; 3; 4] + nodalloads = [0 0; 0 0; 0 -2; 0 0]' + @debug("nodal loads:\n", nodalloads) + dirichletbc = [0 0; NaN NaN; NaN NaN; 0 0]' - # 3. add boundary conditions - # add dirichlet boundary condition u=0 - boundaries = Dict("SUPPORT" => [1, 4]) - add_dirichlet_boundary_condition(model, boundaries) + E = 90 + ν = 0.25 + μ = E/(2*(1+ν)) + λ = E*ν/((1+ν)*(1-2*ν)) + λ = 2*λ*μ/(λ + 2*μ) - # add nodal load to node 3 - nodal_loads = Dict(3 => [-2, 0]) - add_nodal_loads(mode, nodal_loads) + #E = 90.0*ones(2, 4) + #nu = 0.25*ones(2, 4) + u = zeros(2, 4) + du = zeros(2, 4) + R = zeros(2,4) + Kt = zeros(8,8) - # 4. model is defined, solve it - JuliaFEM.solvers.solve_elasticity!(model; parallel=false) + 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 + -(1+ξ[2])/4.0 (1-ξ[1])/4.0] - # 5. extract results - #write_xdmf(m, "results", ["displacement"]) - disp = get_field(model, "displacement") - @assert_eq disp[2, 4] == -2.22224475 + ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1] + 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) + @debug("increment:\n",du) + u += du + if norm(du) < 1.0e-9 + break + end + end + @debug("solution\n",u) + @fact u[2, 3] => roughly(-2.222244754401764) # Tested against Elmer solution end - -test_one_element_solution()