analytical tests, hollow sphere and radial displacement + longitudinal vibration of rod (modal analysis)

src/JuliaFEM.jl

@@ -9,7 +9,7 @@ module JuliaFEM
 importall Base
 
 include("fields.jl")
-export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV
+export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
 include("types.jl")  # data types: Point, IntegrationPoint, ...
 export AbstractPoint, Point, IntegrationPoint, IP, Node
 
@@ -132,7 +132,10 @@ include("postprocess_utils.jl")
 export calc_nodal_values!, 
        get_nodal_vector,
        get_nodal_dict,
-       copy_field!
+       copy_field!,
+       calculate_area,
+       calculate_center_of_mass,
+       calculate_second_moment_of_mass
 include("postprocess_xdmf.jl")
 export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save!
 end

src/fields.jl

@@ -410,3 +410,23 @@ function Base.(:*)(grad::Matrix, field::DVTI)
     return sum([kron(grad[:,i], field[i]') for i=1:length(field)])'
 end
 
+function DVTV(data::Pair{Float64, Vector}...)
+    return DVTV([Increment(d[1], d[2]) for d in data])
+end
+
+function start(f::DVTV)
+    return start(f.data)
+end
+
+function next(f::DVTV, state)
+    return next(f.data, state)
+end
+
+function done(f::DVTV, state)
+    return done(f.data, state)
+end
+
+""" Return time vector from time variable field. """
+function keys(field::DVTV)
+    return Float64[increment.time for increment in field]
+end

src/postprocess_utils.jl

@@ -164,3 +164,83 @@ function call(problem::Problem, field_name::AbstractString, X::Vector, time::Flo
     return fillna
 end
 
+""" Calculate area of cross-section. """
+function calculate_area(problem::Problem, X=[0.0, 0.0], time=0.0)
+    A = 0.0
+    for element in get_elements(problem)
+        elsize = size(element)
+        elsize[1] == 2 || error("wrong dimension of problem for area calculation, element size = $elsize")
+        for ip in get_integration_points(element)
+            w = ip.weight*element(ip, time, Val{:detJ})
+            A += w
+        end
+    end
+    return A
+end
+
+""" Calculate volume of body. """
+function calculate_volume(problem::Problem, X=[0.0, 0.0, 0.0], time=0.0)
+    V = 0.0
+    for element in get_elements(problem)
+        elsize = size(element)
+        elsize[1] == 3 || error("wrong dimension of problem for area calculation, element size = $elsize")
+        for ip in get_integration_points(element)
+            w = ip.weight*element(ip, time, Val{:detJ})
+            V += w
+        end
+    end
+    return V
+end
+
+""" Calculate center of mass of body with respect to X.
+https://en.wikipedia.org/wiki/Center_of_mass
+"""
+function calculate_center_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], time=0.0)
+    M = 0.0
+    Xc = zeros(X)
+    for element in get_elements(problem)
+        for ip in get_integration_points(element)
+            w = ip.weight*element(ip, time, Val{:detJ})
+            M += w
+            rho = haskey(element, "density") ? element("density", ip, time) : 1.0
+            Xp = element("geometry", ip, time)
+            Xc += w*rho*(Xp-X)
+        end
+    end
+    return 1.0/M * Xc
+end
+
+""" Calculate second moment of mass with respect to X.
+https://en.wikipedia.org/wiki/Second_moment_of_area
+"""
+function calculate_second_moment_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], time=0.0)
+    n = length(X)
+    I = zeros(n, n)
+    for element in get_elements(problem)
+        for ip in get_integration_points(element)
+            w = ip.weight*element(ip, time, Val{:detJ})
+            rho = haskey(element, "density") ? element("density", ip, time) : 1.0
+            Xp = element("geometry", ip, time) - X
+            I += w*rho*Xp*Xp'
+        end
+    end
+    return I
+end
+
+function getindex(problem::Problem, field_name::AbstractString)
+    info("fetching result $field_name")
+    timeframes = []
+    for frame in first(problem.elements)[field_name].data
+        push!(timeframes, frame.time)
+    end
+    info("time frames: $timeframes")
+    conn = get_connectivity(problem)
+    increments = Increment[]
+    for time in timeframes
+        p = problem(field_name, time)
+        data = [p[id] for id in conn]
+        push!(increments, Increment(time, data))
+    end
+    return DVTV(increments)
+end
+

src/problems.jl

@@ -313,6 +313,10 @@ function push!(problem::Problem, elements_::Vector...)
     end
 end
 
+function get_connectivity(problem::Problem)
+    return union([get_connectivity(element) for element in get_elements(problem)]...)
+end
+
 function get_gdofs(element::Element, dim::Int)
     conn = get_connectivity(element)
     if length(conn) == 0

src/problems_elasticity.jl

@@ -69,7 +69,7 @@ end
 typealias Elasticity2DSurfaceElements Union{Poi1, Seg2, Seg3}
 typealias Elasticity2DVolumeElements Union{Tri3, Tri6, Quad4, Quad8, Quad9}
 typealias Elasticity3DSurfaceElements Union{Poi1, Tri3, Tri6, Quad4, Quad8, Quad9}
-typealias Elasticity3DVolumeElements Union{Tet4, Tet10, Hex8, Hex20, Hex27}
+typealias Elasticity3DVolumeElements Union{Tet4, Wedge6, Hex8, Tet10, Hex20, Hex27}
 
 
 """ Elasticity equations for 2d cases. """
@@ -544,11 +544,12 @@ function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity},
                 f[i:dim:end] += w*vec(T*N)
             end
         end
-        if haskey(element, "displacement traction force n")
+        if haskey(element, "surface pressure")
             J = element(ip, time, Val{:Jacobian})'
             n = cross(J[:,1], J[:,2])
             n /= norm(n)
-            p = element("displacement traction force n", ip, time)
+            # sign convention, positive pressure is towards surface
+            p = -element("surface pressure", ip, time)
             f += w*p*vec(n*N)
         end
     end

src/solvers.jl

@@ -19,6 +19,13 @@ function Solver{S<:AbstractSolver}(::Type{S}, name="solver", properties...)
     return solver
 end
 
+function Solver{S<:AbstractSolver}(::Type{S}, problems::Problem...)
+    variant = S()
+    solver = Solver{S}("$(S)Solver", 0.0, [], [], 0, variant)
+    push!(solver.problems, problems...)
+    return solver
+end
+
 function get_problems(solver::Solver)
     return solver.problems
 end
@@ -334,6 +341,26 @@ function assemble!(solver::Solver; show_info=true)
     show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
 end
 
+function get_unknown_fields(solver::Solver)
+    fields = Dict()
+    for problem in get_field_problems(solver)
+        field_name = get_unknown_field_name(problem)
+        field_dim = get_unknown_field_dimension(problem)
+        fields[field_name] = field_dim
+    end
+    return fields
+end
+
+function get_unknown_field_name(solver::Solver)
+    fields = get_unknown_fields(solver)
+    return join(sort(collect(keys(fields))), ", ")
+end
+
+function get_unknown_field_dimension(solver::Solver)
+    fields = get_unknown_fields(solver)
+    return sum(values(fields))
+end
+
 """ Default initializer for solver. """
 function initialize!(solver::Solver; show_info=true)
     show_info && info("Initializing problems ...")
@@ -342,8 +369,8 @@ function initialize!(solver::Solver; show_info=true)
     t0 = Base.time()
     field_problems = get_field_problems(solver)
     length(field_problems) != 0 || warn("No field problem found from solver, add some..?")
-    field_dim = get_unknown_field_dimension(first(field_problems))
-    field_name = get_unknown_field_name(first(field_problems))
+    field_name = get_unknown_field_name(solver)
+    field_dim = get_unknown_field_dimension(solver)
     info("initialize!(): looks we are solving $field_name, $field_dim dofs/node")
     nodes = Set{Int64}()
     for problem in problems

src/solvers_modal.jl

@@ -56,6 +56,8 @@ function call(solver::Solver{Modal}; show_info=true, debug=false)
     P, h = create_projection(C1, g)
     K_red = P'*K*P
     M_red = P'*M*P
+    K_red = 1/2*(K_red + K_red')
+    M_red = 1/2*(M_red + M_red')
     t1 = round(toq(), 2)
     info("Eliminated dirichlet boundaries in $t1 seconds.")
 
@@ -71,7 +73,22 @@ function call(solver::Solver{Modal}; show_info=true, debug=false)
     end
 
     tic()
-    om2, X = eigs(K_red[nz,nz], M_red[nz,nz]; nev=props.nev, which=props.which)
+    om2 = nothing
+    X = nothing
+    try
+        om2, X = eigs(K_red[nz,nz], M_red[nz,nz]; nev=props.nev, which=props.which)
+    catch
+        info("failed to calculate eigenvalues")
+        info("K sym?", issym(K_red[nz,nz]))
+        info("M sym?", issym(M_red[nz,nz]))
+        info("K posdef?", isposdef(K_red[nz,nz]))
+        info("M posdef?", isposdef(M_red[nz,nz]))
+        k1 = maximum(abs(K_red[nz,nz] - K_red[nz,nz]'))
+        m1 = maximum(abs(M_red[nz,nz] - M_red[nz,nz]'))
+        info("K skewness ", k1)
+        info("M skewness ", m1)
+        rethrow()
+    end
     props.eigvals = om2
     props.eigvecs = zeros(ndofs, length(om2))
     v = zeros(ndofs)
@@ -82,6 +99,23 @@ function call(solver::Solver{Modal}; show_info=true, debug=false)
     end
     t1 = round(toq(), 2)
     info("Eigenvalues computed in $t1 seconds. Eigenvalues: $om2")
+
+    for i=1:length(om2)
+        u = props.eigvecs[:,i]
+        field_dim = get_unknown_field_dimension(solver)
+        field_name = get_unknown_field_name(solver)
+        nnodes = round(Int, length(u)/field_dim)
+        solution = reshape(u, field_dim, nnodes)
+        for problem in get_problems(solver)
+            local_sol = Dict{Int64, Vector{Float64}}()
+            for node_id in get_connectivity(problem)
+                local_sol[node_id] = solution[:, node_id]
+            end
+            freq = real(sqrt(om2[i])/(2.0*pi))
+            update!(problem, field_name, freq => local_sol)
+        end
+    end
+    
     return true
 end
 

test/test_elasticity_hollow_sphere_with_surface_pressure.jl

@@ -0,0 +1,60 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Preprocess
+using JuliaFEM.Testing
+
+#=
+test subjects:
+- surface pressure load in curved surface
+- verification of elements wedge6 and wedge15
+
+from Code Aster:
+N38      -8.85861895037377E-01 -3.46944695195361E-18 -3.46944695195361E-18
+coords of N38 = (1.0, 0.0, 0.0)
+
+Analytical solution
+uᵣ(r) = b³p/(2Er²(a³-b³)) * (a³(ν+1) + r³(-4ν+2)), where
+a = inner surface radial distance, b = outer surface ...
+
+if a=0.9, b=1.0, ν=1/3, E = 24580 and p = 7317 equation yields
+-9/10 for radial displacement
+=#
+
+@testset """1/8 hollow sphere with surface load""" begin
+    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh = aster_read_mesh(mesh_file, "HOLLOWSPHERE8_WEDGE6")
+    body = Problem(Elasticity, "hollow sphere 1/8 model", 3)
+    body.elements = create_elements(mesh, "HOLLOWSPHERE8")
+    update!(body, "youngs modulus", 24580.0)
+    update!(body, "poissons ratio", 1/3)
+    bc = Problem(Dirichlet, "symmetry bc", 3, "displacement")
+    el1 = create_elements(mesh, "FACE1")
+    update!(el1, "displacement 3", 0.0)
+    el2 = create_elements(mesh, "FACE2")
+    update!(el2, "displacement 2", 0.0)
+    el3 = create_elements(mesh, "FACE3")
+    update!(el3, "displacement 1", 0.0)
+    bc.elements = [el1; el2; el3]
+    lo = Problem(Elasticity, "pressure load", 3)
+    lo.elements = create_elements(mesh, "OUTER")
+    update!(lo, "surface pressure", 7317.0)
+    solver = LinearSolver(body, bc, lo)
+    solver()
+
+    X = lo("geometry")
+    u = lo("displacement")
+    nids = sort(collect(keys(X)))
+    umag = Float64[norm(u[id]) for id in nids]
+    um = mean(umag)
+    info("mean umag = $um")
+    info("std umag = ", std(umag))
+    rtol = norm(um - 0.9) / max(norm(um), 0.9) * 100.0
+    info("rtol = $rtol")
+    @test rtol < 1.5 # percents
+    u_CA = [-8.85861895037377E-01, -3.46944695195361E-18, -3.46944695195361E-18]
+    rtol = norm(u[38] - u_CA) / max(norm(u[38]), norm(u_CA)) * 100.0
+    info("rel diff to CA = $rtol %")
+    @test isapprox(u[38], u_CA)
+end

test/test_modal_analysis_elasticity.jl

@@ -0,0 +1,142 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Preprocess
+using JuliaFEM.Postprocess
+using JuliaFEM.Testing
+
+@testset "calculate cross-sectional properties" begin
+    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh = aster_read_mesh(mesh_file, "CYLINDER_20_TET4")
+    # calculate cross-sectional properties A and Iₓ
+    fixed1 = Problem(Dirichlet, "left support", 3, "displacement")
+    fixed1.elements = create_elements(mesh, "FACE1")
+    A = calculate_area(fixed1)
+    info("cross-section area: $A")
+    # real area is π
+    @test isapprox(A, pi; rtol=0.1)
+    Xc = calculate_center_of_mass(fixed1)
+    info("center of mass: $Xc")
+    @test isapprox(Xc, [0.0, 0.0, 0.0]; atol=1.0e-12)
+    I = calculate_second_moment_of_mass(fixed1)
+    info("moments:")
+    info(I)
+    I_expected = zeros(3, 3)
+    I_expected[2,2] = I_expected[3,3] = pi/4
+    rtol = norm(I[2,2]-I_expected[2,2]) / max(I[2,2],I_expected[2,2])
+    info("I rtol = $rtol")
+    @test isapprox(I, I_expected; rtol = 0.2)
+end
+
+#=
+test subjects:
+- calculate cross-sectional properties
+- modal analysis with known solution
+
+Fixed-fixed solution is ωᵢ = λᵢ²√(EI/ρA) , where λᵢ = cosh(λᵢℓ)cos(λᵢℓ)
+
+1:  4.730040744862704
+2:  7.853204624095838
+3: 10.995607838001671
+
+[1] De Silva, Clarence W. Vibration: fundamentals and practice. CRC press, 2006, p.355
+=#
+
+@testset "long rod under point load" begin
+    mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
+    mesh = aster_read_mesh(mesh_file, "CYLINDER_20_TET10")
+#   for (id, coords) in mesh.nodes
+#       mesh.nodes[id][1] *= 5.0
+#   end
+    body = Problem(Elasticity, "rod", 3)
+    body.elements = create_elements(mesh, "CYLINDER")
+    E = 50475.44814745859
+    rho = 1.0
+    update!(body.elements, "youngs modulus", E)
+    update!(body.elements, "poissons ratio", 0.3)
+    update!(body.elements, "density", rho)
+    # calculate cross-sectional properties A and Iₓ
+    fixed1 = Problem(Dirichlet, "left support", 3, "displacement")
+    fixed1.elements = create_elements(mesh, "FACE1")
+    update!(fixed1.elements, "displacement 1", 0.0)
+    update!(fixed1.elements, "displacement 2", 0.0)
+    update!(fixed1.elements, "displacement 3", 0.0)
+    fixed2 = Problem(Dirichlet, "right support", 3, "displacement")
+    fixed2.elements = create_elements(mesh, "FACE2")
+    update!(fixed2.elements, "displacement 1", 0.0)
+    update!(fixed2.elements, "displacement 2", 0.0)
+    update!(fixed2.elements, "displacement 3", 0.0)
+    A = calculate_area(fixed1)
+    info("cross-section area: $A")
+    # using SALOME / SMESH, A = 2.82843
+    # real area is π
+    @test isapprox(A, pi; rtol=0.1)
+    Xc = calculate_center_of_mass(fixed1)
+    info("center of mass: $Xc")
+    @test isapprox(Xc, [0.0, 0.0, 0.0]; atol=1.0e-5)
+    I = calculate_second_moment_of_mass(fixed1)
+    info("moments:")
+    info(I)
+    I_expected = zeros(3, 3)
+    r = 1.0
+    I_expected[2,2] = I_expected[3,3] = pi/4*r^2
+    rtol = norm(I[2,2]-I_expected[2,2]) / max(I[2,2],I_expected[2,2])
+    info("I rtol = $rtol")
+    @test isapprox(I, I_expected; rtol = 0.2)
+#=
+    # apply transform Tx + b, in this case move cross-section to
+    # xy-plane from yz-plane, i.e.
+    # x₁ = y₂
+    # y₁ = z₂
+    T = [
+        0.0 1.0 0.0
+        0.0 0.0 1.0]
+    b = [0.0, 0.0]
+    X 1 = first(cross_section)("geometry", [1/3, 1/3], 0.0)
+    apply_affine_transform!(cross_section, T, b)
+    X2 = first(cross_section)("geometry", [1/3, 1/3], 0.0)
+    info("X1 = $X1, X2 = $X2")
+    @test isapprox(T*X1+b, X2)
+=#
+    c = sqrt(E*I[2,2]/(rho*A))
+    info("c = $c")
+
+
+    # analytical solution is
+    l = 20.0
+    r = 1.0
+    la = 4.730040744862704/l
+
+    # semi-analytical (c numerical)
+    freq_sa = (c*la^2)/(2*pi)
+    info("freq_sa = $freq_sa")
+
+    A = pi*r^2
+    I = pi/4*r^4
+    c = sqrt(E*I/(rho*A))
+    info("c analytical = $c")
+    freq_a = (c*la^2)/(2*pi)
+
+    info("freq_a = $freq_a")
+
+    solver = Solver(Modal, body, fixed1, fixed2)
+    solver()
+    freqs = keys(body["displacement"])
+
+    rtol1 = norm(freq_sa - freqs[2])/max(freq_sa, freqs[2])
+    rtol2 = norm(freq_a - freqs[2])/max(freq_a, freqs[2])
+    info("rtol 1 = $rtol1, rtol 2 = $rtol2")
+    @test rtol2 < 1.0e-2
+#=
+    result = XDMF()
+    for (i, freq) in enumerate(freqs)
+        isapprox(freq, 0.0) && continue
+        info("$i freq: $freq")
+        xdmf_new_result!(result, body, freq)
+        xdmf_save_field!(result, body, freq, "displacement"; field_type="Vector")
+    end
+    xdmf_save!(result, "/tmp/rod_nf.xmf")
+=#
+end
+