From cb72cd17fec12e9883ee1a7d2010be1dc202c4b4 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 16 Aug 2015 13:33:02 +0300 Subject: [PATCH] interpolation routines to separate file. --- src/JuliaFEM.jl | 2 ++ src/elasticity_solver.jl | 48 ------------------------------ src/interpolate.jl | 53 ++++++++++++++++++++++++++++++++++ src/types.jl | 28 ++++++++++++++++++ test/test_elasticity_solver.jl | 29 ------------------- test/test_interpolate.jl | 33 +++++++++++++++++++++ 6 files changed, 116 insertions(+), 77 deletions(-) create mode 100644 src/interpolate.jl create mode 100644 src/types.jl create mode 100644 test/test_interpolate.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index a8920eb..a454199 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -6,6 +6,8 @@ module JuliaFEM VERSION < v"0.4-" && using Docile using Lexicon +include("types.jl") # type definitions +include("interpolate.jl") include("elasticity_solver.jl") include("xdmf.jl") include("abaqus_reader.jl") diff --git a/src/elasticity_solver.jl b/src/elasticity_solver.jl index 38580a3..b23ee62 100644 --- a/src/elasticity_solver.jl +++ b/src/elasticity_solver.jl @@ -47,54 +47,6 @@ function dummy(a) end -""" -Interpolate field variable using basis functions f for point ip. -This function tries to be as general as possible and allows interpolating -lot of different fields. - -Parameters ----------- -field :: Array{Number, dim} - Field variable -basis :: Function - Basis functions -ip :: Array{Number, 1} - Point to interpolate -""" -function interpolate{T<:Real}(field::Array{T,1}, basis::Function, ip) - result = dot(field, basis(ip)) - return result -end -function interpolate{T<:Real}(field::Array{T,2}, basis::Function, ip) - m, n = size(field) - bip = basis(ip) - tmp = size(bip) - if length(tmp) == 1 - ndim = 1 - nnodes = tmp[1] - else - ndim, nnodes = size(bip) - end - if ndim == 1 - if n == nnodes - result = field * bip - elseif m == nnodes - result = field' * bip - end - else - if n == nnodes - result = bip' * field - elseif m == nnodes - result = bip' * field' - end - end - if length(result) == 1 - result = result[1] - end - return result -end - - """ Calculate local tangent stiffness matrix and residual force vector R = T - F for elasticity problem. diff --git a/src/interpolate.jl b/src/interpolate.jl new file mode 100644 index 0000000..e0441ac --- /dev/null +++ b/src/interpolate.jl @@ -0,0 +1,53 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +""" +Interpolate field variable using basis functions f for point ip. +This function tries to be as general as possible and allows interpolating +lot of different fields. + +Parameters +---------- +field :: Array{Number, dim} + Field variable +basis :: Function + Basis functions +ip :: Array{Number, 1} + Point to interpolate +""" +function interpolate{T<:Real}(field::Array{T,1}, basis::Function, ip) + result = dot(field, basis(ip)) + return result +end +function interpolate{T<:Real}(field::Array{T,2}, basis::Function, ip) + m, n = size(field) + bip = basis(ip) + tmp = size(bip) + if length(tmp) == 1 + ndim = 1 + nnodes = tmp[1] + else + ndim, nnodes = size(bip) + end + if ndim == 1 + if n == nnodes + result = field * bip + elseif m == nnodes + result = field' * bip + end + else + if n == nnodes + result = bip' * field + elseif m == nnodes + result = bip' * field' + end + end + if length(result) == 1 + result = result[1] + end + return result +end +function interpolate(e::Element, field::ASCIIString, x::Array{Float64,1}; derivative=false) + return interpolate(e.attributes[field], derivative ? e.dbasis : e.basis, x) +end + diff --git a/src/types.jl b/src/types.jl new file mode 100644 index 0000000..6ddedfd --- /dev/null +++ b/src/types.jl @@ -0,0 +1,28 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + + +type Element + id :: Int + element_type :: Int + node_ids :: Array{Int, 1} + basis :: Function + dbasis :: Function + attributes :: Dict{ASCIIString, Any} + ipoints :: Array{Float64, 2} + iweights :: Array{Float64, 1} +end + + +type Assembly + # LHS + I :: Array{Int64, 1} + J :: Array{Int64, 1} + A :: Array{Float64, 1} + # RHS + i :: Array{Int64, 1} + b :: Array{Float64, 1} + # global dofs for each element + gdofs :: Dict{Int64, Array{Int64, 1}} +end + diff --git a/test/test_elasticity_solver.jl b/test/test_elasticity_solver.jl index 167c9e3..57263b3 100644 --- a/test/test_elasticity_solver.jl +++ b/test/test_elasticity_solver.jl @@ -140,35 +140,6 @@ facts("test solve elasticity increment, two elements") do end -using JuliaFEM.elasticity_solver: interpolate -facts("test interpolation of different field variables") do - 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 - -(1+ξ[2])/4.0 (1-ξ[1])/4.0] - F1 = [36.0, 36.0, 36.0, 36.0] - F2 = [36.0 36.0 36.0 36.0] - F3 = F2' - F4 = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]' - F5 = F4' - F6 = [36, 36, 36, 36] - - @fact interpolate(F1, N, [0.0, 0.0]) => 36.0 - @fact interpolate(F2, N, [0.0, 0.0]) => 36.0 - @fact interpolate(F3, N, [0.0, 0.0]) => 36.0 - @fact interpolate(F4, N, [0.0, 0.0]) => [5.0; 0.5] - @fact interpolate(F5, N, [0.0, 0.0]) => [5.0; 0.5] - @fact interpolate(F5, dNdξ, [0.0, 0.0]) => [5.0 0.0; 0.0 0.5] - @fact interpolate(F6, N, [0.0, 0.0]) => 36 -end - - using JuliaFEM.elasticity_solver: assemble! diff --git a/test/test_interpolate.jl b/test/test_interpolate.jl new file mode 100644 index 0000000..37caf1b --- /dev/null +++ b/test/test_interpolate.jl @@ -0,0 +1,33 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM: interpolate + +using FactCheck + +facts("test interpolation of different field variables") do + 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 + -(1+ξ[2])/4.0 (1-ξ[1])/4.0] + F1 = [36.0, 36.0, 36.0, 36.0] + F2 = [36.0 36.0 36.0 36.0] + F3 = F2' + F4 = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]' + F5 = F4' + F6 = [36, 36, 36, 36] + + @fact interpolate(F1, N, [0.0, 0.0]) --> 36.0 + @fact interpolate(F2, N, [0.0, 0.0]) --> 36.0 + @fact interpolate(F3, N, [0.0, 0.0]) --> 36.0 + @fact interpolate(F4, N, [0.0, 0.0]) --> [5.0; 0.5] + @fact interpolate(F5, N, [0.0, 0.0]) --> [5.0; 0.5] + @fact interpolate(F5, dNdξ, [0.0, 0.0]) --> [5.0 0.0; 0.0 0.5] + @fact interpolate(F6, N, [0.0, 0.0]) --> 36 +end