diff --git a/test/test_abaqus_reader.jl b/test/test_abaqus_reader.jl index 49115d9..d4f77d9 100644 --- a/test/test_abaqus_reader.jl +++ b/test/test_abaqus_reader.jl @@ -9,13 +9,13 @@ using JuliaFEM.Test model = open(parse_abaqus, Pkg.dir("JuliaFEM")*"/geometry/3d_beam/palkki.inp") @test length(model["nodes"]) == 298 @test length(model["elements"]) == 120 - @test length(model["elsets"]["BODY1"]) == 120 + @test length(model["elsets"]["Body1"]) == 120 @test length(model["nsets"]["SUPPORT"]) == 9 @test length(model["nsets"]["LOAD"]) == 9 @test length(model["nsets"]["TOP"]) == 83 end -@testset "test that reader throws error when dimension information of elemenet is missing" begin +@testset "test that reader throws error when dimension information of element is missing" begin # *ELEMENT, TYPE=neverseenbefore, ELSET=Body1 data = """ 1, 243, 240, 191, 117, 245, 242, 244, @@ -23,7 +23,7 @@ end """ model = Dict() header = Dict("section"=>"ELEMENT", "options" => Dict("TYPE" => "neverseenbefore", "ELSET"=>"Body1")) - @test_throws parse_element_section(model, header, data) + @test_throws Exception parse_element_section(model, header, data) end @testset "test read element section" begin diff --git a/test/test_assembly.jl b/test/test_assembly.jl index 4d53ae4..e9ca24b 100644 --- a/test/test_assembly.jl +++ b/test/test_assembly.jl @@ -5,39 +5,36 @@ using JuliaFEM using JuliaFEM.Test @testset "test static condensation" begin - nodes = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] - el1 = Quad4([1, 2, 3, 4]) - el1["geometry"] = Vector[nodes[1], nodes[2], nodes[3], nodes[4]] - el1["temperature thermal conductivity"] = 6.0 - el1["temperature load"] = [12.0, 12.0, 12.0, 12.0] - el2 = Seg2([1, 2]) - el2["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]] - el2["temperature flux"] = 6.0 - field_problem = HeatProblem() - push!(field_problem, el1) - push!(field_problem, el1) + K = sparse([ + 4.0 -1.0 -2.0 -1.0 + -1.0 4.0 -1.0 -2.0 + -2.0 -1.0 4.0 -1.0 + -1.0 -2.0 -1.0 4.0]) - el3 = Seg2([3, 4]) - el3["geometry"] = Vector[nodes[3], nodes[4]] - el3["temperature"] = 0.0 - boundary_problem = DirichletProblem("temperature", 1) - push!(boundary_problem, el3) + f = sparse([6.0, 6.0, 3.0, 3.0]) - fass = assemble(field_problem, 0.0) - bass = assemble(boundary_problem, 0.0) + I = [1, 2] + B = [3, 4] -# interior_dofs = [1, 2] - boundary_dofs = [3, 4] - cass = condensate(fass, boundary_dofs) - @test isapprox(full(cass.Kc), [ + Kc, fc = eliminate_interior_dofs(K, f, B, I) + Kc = full(Kc) + fc = full(fc) + dump(Kc) + dump(fc) + + Kc_expected = [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 4.8 -4.8 - 0.0 0.0 -4.8 4.8]) - @test isapprox(full(cass.fc)', [0.0 0.0 12.0 12.0]) - @test cass.interior_dofs == [1, 2] + 0.0 0.0 2.4 -2.4 + 0.0 0.0 -2.4 2.4] + fc_expected = [0.0, 0.0, 9.0, 9.0] + # TODO: needs to check numbers + @test isapprox(Kc, Kc_expected) + @test isapprox(fc, fc_expected) + + #= x = sparse(zeros(4))' la = sparse(zeros(4))' la[3] = la[4] = 24.0 @@ -47,5 +44,6 @@ using JuliaFEM.Test info(x) @test isapprox(x[1], 1.0) @test isapprox(x[2], 1.0) + =# end diff --git a/test/test_basis.jl b/test/test_basis.jl index 0f248f0..28d97a9 100644 --- a/test/test_basis.jl +++ b/test/test_basis.jl @@ -1,15 +1,15 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -module BasisTests - +using JuliaFEM using JuliaFEM.Test -using JuliaFEM.Core: AbstractElement, Element -import JuliaFEM.Core: get_basis, get_dbasis +importall Base +import JuliaFEM: get_basis, get_dbasis -abstract TestElement <: AbstractElement +type TestElement <: AbstractElement +end -function get_basis(::Type{TestElement}, xi::Vector{Float64}) +function get_basis(element::Element{TestElement}, xi, time) 1/4*[ (1-xi[1])*(1-xi[2]) (1+xi[1])*(1-xi[2]) @@ -17,85 +17,104 @@ function get_basis(::Type{TestElement}, xi::Vector{Float64}) (1-xi[1])*(1+xi[2])]' end -function get_dbasis(::Type{TestElement}, xi::Vector{Float64}) +function get_dbasis(element::Element{TestElement}, xi, time) 1/4*[ -(1-xi[2]) (1-xi[2]) (1+xi[2]) -(1+xi[2]) -(1-xi[1]) -(1+xi[1]) (1+xi[1]) (1-xi[1])] end +function length(element::Element{TestElement}) + return 4 +end + +function size(element::Element{TestElement}) + return (2, 4) +end + function get_element() - element = Element{TestElement}([1, 2, 3, 4]) - element["geometry"] = Vector{Float64}[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] - element["temperature"] = Float64[1.0, 2.0, 3.0, 4.0] - element["displacement1"] = Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]] - element["displacement2"] = Vector{Float64}[[0.0, 0.0], [1.0, -1.0], [2.0, 3.0], [0.0, 0.0]] + element = Element(TestElement, [1, 2, 3, 4]) + X = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [1.0, 0.0], + 3 => [1.0, 1.0], + 4 => [0.0, 1.0]) + T = Dict{Int64, Float64}( + 1 => 1.0, + 2 => 2.0, + 3 => 3.0, + 4 => 4.0) + u1 = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [0.0, 0.0], + 3 => [1/4, 0.0], + 4 => [0.0, 0.0]) + u2 = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [1.0, -1.0], + 3 => [2.0, 3.0], + 4 => [0.0, 0.0]) + update!(element, "geometry", X) + update!(element, "temperature", T) + update!(element, "displacement1", u1) + update!(element, "displacement2", u2) return element end -### Test interpolation in spatial domain - -function test_basis_interpolation() +@testset "spatial interpolation in basis" begin element = get_element() - info(element([0.0, 0.0])) - @test element([0.0, 0.0]) == 1/4*[1 1 1 1] - @test element([0.0, 0.0], 1.0) == 1/4*[1 1 1 1] + @test isapprox(element([0.0, 0.0], 0.0), 1/4*[1 1 1 1]) + @test isapprox(element([0.0, 0.0], 1.0), 1/4*[1 1 1 1]) end -function test_basis_gradient_interpolation() +@testset "gradient of shape functions" begin element = get_element() - grad = element([0.0, 0.0], Val{:grad}) - info("grad = \n$grad") - @test grad == 1/2*[-1 1 1 -1; -1 -1 1 1] + grad = element([0.0, 0.0], 0.0, Val{:Grad}) + @test isapprox(grad, 1/2*[-1 1 1 -1; -1 -1 1 1]) end -function test_interpolation_of_scalar_field_in_spatial_domain() +@testset "interpolation of scalar field in spatial domain" begin # in unit square: T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2]) element = get_element() T_known(X) = 1 + X[1] + 3*X[2] - 2*X[1]*X[2] - T_interpolated = element("temperature", [0.0, 0.0]) - @test T_interpolated == T_known([0.5, 0.5]) + T_interpolated = element("temperature", [0.0, 0.0], 0.0) + @test isapprox(T_interpolated, T_known([0.5, 0.5])) end -function test_interpolation_of_gradient_of_scalar_field_in_spatial_domain() +@testset "interpolation of gradient of scalar field in spatial domain" begin # in unit square: grad(T)(X) = [1-2X[2], 3-2*X[1]] element = get_element() - gradT = element("temperature", [0.0, 0.0], Val{:grad}) + gradT = element("temperature", [0.0, 0.0], 0.0, Val{:Grad}) gradT_expected(X) = [1-2*X[2] 3-2*X[1]] - @test gradT == gradT_expected([0.5, 0.5]) + @test isapprox(gradT, gradT_expected([0.5, 0.5])) end -function test_interpolation_of_vector_field() +@testset "test interpolation of vector field" begin # in unit square, u(X,t) = [1/4*t*X[1]*X[2], 0, 0] element = get_element() - u = element("displacement1", [0.0, 0.0]) + u = element("displacement1", [0.0, 0.0], 0.0) # x = X+u u_expected(X) = [1/4*X[1]*X[2], 0] # @test isapprox(x, [9/16, 1/2]) @test isapprox(u, u_expected([0.5, 0.5])) end - -function test_interpolation_of_gradient_of_vector_field() +@testset "interpolation of gradient of vector_field" begin # in unit square, u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)] # => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]] element = get_element() # displacement = Field( # (0.5, Vector[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]), # (1.5, Vector[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]])) - gradu = element("displacement2", [0.0, 0.0], Val{:grad}) + gradu = element("displacement2", [0.0, 0.0], 0.0, Val{:Grad}) gradu_expected(X) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]] @test isapprox(gradu, gradu_expected([0.5, 0.5])) end -### Test interpolation in time domain - -#= - -function test_linear_time_extrapolation_of_field() +@testset "linear time extrapolation of field" begin #T_known(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2]) - T = Field( - (0.0, [0.0, 0.0, 0.0, 0.0]), - (1.0, [1.0, 2.0, 3.0, 4.0])) + T = DVTV() + update!(T, 0.0 => [0.0, 0.0, 0.0, 0.0]) + update!(T, 1.0 => [1.0, 2.0, 3.0, 4.0]) @test T(-1.0) == -1.0*[1.0, 2.0, 3.0, 4.0] @test T( 3.0) == 3.0*[1.0, 2.0, 3.0, 4.0] # when going to \pm infinity, return the last one. @@ -103,150 +122,153 @@ function test_linear_time_extrapolation_of_field() @test T(+Inf) == 1.0*[1.0, 2.0, 3.0, 4.0] end -function test_constant_time_extrapolation_of_field() +@testset "constant time extrapolation of field" begin #T_known(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2]) - T = Field( - (0.0, [0.0, 0.0, 0.0, 0.0]), - (1.0, [1.0, 2.0, 3.0, 4.0])) - @test T(-1.0, :constant) == [0.0, 0.0, 0.0, 0.0] - @test T( 3.0, :constant) == [1.0, 2.0, 3.0, 4.0] + T = DVTV() + update!(T, 0.0 => [0.0, 0.0, 0.0, 0.0]) + update!(T, 1.0 => [1.0, 2.0, 3.0, 4.0]) + @test isapprox(T(-1.0, Val{:constant}), [0.0, 0.0, 0.0, 0.0]) + @test isapprox(T( 3.0, Val{:constant}), [1.0, 2.0, 3.0, 4.0]) end -function test_time_extrapolation_of_field_with_single_timestep() - T = Field([1.0, 2.0, 3.0, 4.0]) - @test T(1.0) == [1.0, 2.0, 3.0, 4.0] +@testset "time extrapolation of field with only one timestep" begin + T = DVTV() + update!(T, 0.0 => [1.0, 2.0, 3.0, 4.0]) + @test isapprox(T(1.0), [1.0, 2.0, 3.0, 4.0]) end -function test_interpolation_in_temporal_basis() - i1 = Increment(0.0) - i2 = Increment(1.0) - i3 = Increment(2.0) - t1 = TimeStep(0.0, Increment[i1]) - t2 = TimeStep(2.0, Increment[i2]) - t3 = TimeStep(4.0, Increment[i3]) - field = Field(TimeStep[t1, t2, t3]) - @test field(-Inf) == [0.0] - @test field( 0.0) == [0.0] - @test field( 1.0) == [0.5] - @test field( 2.0) == [1.0] - @test field( 3.0) == [1.5] - @test field( 4.0) == [2.0] - @test field(+Inf) == [2.0] +@testset "interpolation in temporal direction" begin + field = DCTV() + update!(field, 0.0 => 0.0) + update!(field, 2.0 => 1.0) + update!(field, 4.0 => 2.0) + @test isapprox(field(-Inf), 0.0) + @test isapprox(field( 0.0), 0.0) + @test isapprox(field( 1.0), 0.5) + @test isapprox(field( 2.0), 1.0) + @test isapprox(field( 3.0), 1.5) + @test isapprox(field( 4.0), 2.0) + @test isapprox(field(+Inf), 2.0) end -function test_derivative_interpolation_in_temporal_basis_in_constant_velocity() - i1 = Increment(0.0) - i2 = Increment(1.0) - i3 = Increment(2.0) - t1 = TimeStep(0.0, Increment[i1]) - t2 = TimeStep(2.0, Increment[i2]) - t3 = TimeStep(4.0, Increment[i3]) - field = Field(TimeStep[t1, t2, t3]) - @test field(+Inf, Val{:diff}) == [0.5] - @test field(-Inf, Val{:diff}) == [0.5] - @test field( 0.0, Val{:diff}) == [0.5] - @test field( 0.5, Val{:diff}) == [0.5] - @test field( 1.0, Val{:diff}) == [0.5] - @test field( 1.5, Val{:diff}) == [0.5] - @test field( 2.0, Val{:diff}) == [0.5] +@testset "time derivative interpolation in temporal basis in constant velocity" begin + field = DCTV() + update!(field, 0.0 => 0.0) + update!(field, 2.0 => 1.0) + update!(field, 4.0 => 2.0) + @test isapprox(field(+Inf, Val{:diff}), 0.5) + @test isapprox(field(-Inf, Val{:diff}), 0.5) + @test isapprox(field( 0.0, Val{:diff}), 0.5) + @test isapprox(field( 0.5, Val{:diff}), 0.5) + @test isapprox(field( 1.0, Val{:diff}), 0.5) + @test isapprox(field( 1.5, Val{:diff}), 0.5) + @test isapprox(field( 2.0, Val{:diff}), 0.5) end -function test_derivative_interpolation_in_temporal_basis_in_variable_velocity() - t = linspace(0, 2, 5) - x = 1/2*t.^2 - timesteps = TimeStep[] - for (ti, xi) in zip(t, x) - increment = Increment(xi) - push!(timesteps, TimeStep(ti, increment)) +@testset "time derivative interpolation in temporal basis in variable velocity" begin + pos = DCTV() + for ti in linspace(0, 2, 5) + update!(pos, ti => 1/2*ti^2) end # => ((0.0,0.0),(0.5,0.125),(1.0,0.5),(1.5,1.125),(2.0,2.0)) - pos = Field(timesteps) - - velocity = pos(1.0, Val{:diff})[1] + velocity = pos(1.0, Val{:diff}) v1 = (0.500 - 0.125)/0.5 v2 = (1.125 - 0.500)/0.5 @test isapprox(velocity, mean([v1, v2])) # = 1.00 - - velocity = pos(2.0, Val{:diff})[1] + velocity = pos(2.0, Val{:diff}) @test isapprox(velocity, (2.0-1.125)/0.5) # = 1.75 end -function test_derivative_interpolation_in_temporal_basis_in_variable_velocity_check_type() - t = linspace(0, 2, 5) - x = 1/2*t.^2 - timesteps = TimeStep[] - for (ti, xi) in zip(t, x) - increment = Increment(xi) - push!(timesteps, TimeStep(ti, increment)) - end - # => ((0.0,0.0),(0.5,0.125),(1.0,0.5),(1.5,1.125),(2.0,2.0)) - pos = Field(timesteps) - velocity = pos(1.0, Val{:diff}) - # after interpolation, we are expecting to have same type where we started - @test isa(velocity, Increment) == true -end - function test_time_derivative_gradient_interpolation_of_field() # in unit square, u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)] # => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]] # => d(u_i,j)/dt = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]] - geometry = Field([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]') - displacement = Field( - (0.5, Vector[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]), - (1.5, Vector[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]])) + X = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [1.0, 0.0], + 3 => [1.0, 1.0], + 4 => [0.0, 1.0]) + u1 = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [0.5, -0.5], + 3 => [1.0, 1.5], + 4 => [0.0, 0.0]) + u2 = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [1.5, -1.5], + 3 => [3.0, 4.5], + 4 => [0.0, 0.0]) + element = Element(TestElement, [1, 2, 3, 4]) + update!(element, "geometry", X) + update!(element, "displacement", 0.5 => u1) + update!(element, "displacement", 1.5 => u2) - # wanted - #u = get_basis(element, "displacement") - #L = grad(diff(u)) - #D = 1/2*(L + L') - #@test isapprox(D([0.0, 0.0], 1.0), ...) - - basis, dbasis = get_basis() - N = Basis(basis, dbasis) xi = [0.0, 0.0] time = 1.2 - grad = ElementGradientBasis(N, geometry)(xi, time) - increment = displacement(time, Val{:derivative}) - diffgradu = sum([grad[:,i]*increment[i]' for i=1:length(increment)])' + diffgradu = element("displacement", xi, time, Val{:diff}, Val{:Grad}) diffgradu_expected(X, t) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]] @test diffgradu == diffgradu_expected([0.5, 0.5], 1.2) end -"""basic continuum interpolations""" -function test_basic_interpolations() - - element = Quad4([1, 2, 3, 4]) - - element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] - element["temperature"] = ([0.0, 0.0, 0.0, 0.0], [1.0, 2.0, 3.0, 4.0]) - element["displacement"] = ( - Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]], - Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]]) +@testset "some continuum mechanics interpolations" begin + X = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [1.0, 0.0], + 3 => [1.0, 1.0], + 4 => [0.0, 1.0]) + u1 = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [0.0, 0.0], + 3 => [0.0, 0.0], + 4 => [0.0, 0.0]) + u2 = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [0.0, 0.0], + 3 => [1/4, 0.0], + 4 => [0.0, 0.0]) + element = Element(Quad4, [1, 2, 3, 4]) + update!(element, "geometry", X) + update!(element, "displacement", 0.0 => u1) + update!(element, "displacement", 1.0 => u2) # from my old home works - basis = get_basis(element) - dbasis = grad(basis) - @test isapprox(basis("geometry", [0.0, 0.0], 1.0) + basis("displacement", [0.0, 0.0], 1.0), [9/16, 1/2]) - gradu = dbasis("displacement", [0.0, 0.0], 1.0) + X = element("geometry", [0.0, 0.0], 1.0) + u = element("displacement", [0.0, 0.0], 1.0) + x = X + u + x_expected = [9/16, 1/2] + gradu = element("displacement", [0.0, 0.0], 1.0, Val{:Grad}) epsilon = 1/2*(gradu + gradu') rotation = 1/2*(gradu - gradu') - X = basis("geometry", [0.0, 0.0], 1.0) k = 0.25 - epsilon_wanted = [X[2]*k 1/2*X[1]*k; 1/2*X[1]*k 0] - rotation_wanted = [0 k/2*X[1]; -k/2*X[1] 0] - @test isapprox(epsilon, epsilon_wanted) - @test isapprox(rotation, rotation_wanted) + epsilon_expected = [ + X[2]*k 1/2*X[1]*k + 1/2*X[1]*k 0] + rotation_expected = [ + 0 k/2*X[1] + -k/2*X[1] 0] F = I + gradu - @test isapprox(F, [X[2]*k+1 X[1]*k; 0 1]) + F_expected = [ + X[2]*k+1 X[1]*k + 0 1] C = F'*F - @test isapprox(C, [(X[2]*k+1)^2 (X[2]*k+1)*X[1]*k; (X[2]*k+1)*X[1]*k X[1]^2*k^2+1]) + C_expected = [ + (X[2]*k+1)^2 (X[2]*k+1)*X[1]*k + (X[2]*k+1)*X[1]*k X[1]^2*k^2+1] E = 1/2*(F'*F - I) - @test isapprox(E, [1/2*(X[2]*k + 1)^2-1/2 1/2*(X[2]*k+1)*X[1]*k; 1/2*(X[2]*k + 1)*X[1]*k 1/2*X[1]^2*k^2]) + E_expected = [ + 1/2*(X[2]*k + 1)^2-1/2 1/2*(X[2]*k+1)*X[1]*k + 1/2*(X[2]*k + 1)*X[1]*k 1/2*X[1]^2*k^2] U = 1/sqrt(trace(C) + 2*sqrt(det(C)))*(C + sqrt(det(C))*I) - @test isapprox(U, [1.24235 0.13804; 0.13804 1.02149]) + U_expected = [1.24235 0.13804; 0.13804 1.02149] + + @test isapprox(x, x_expected) + @test isapprox(epsilon, epsilon_expected) + @test isapprox(rotation, rotation_expected) + @test isapprox(F, F_expected) + @test isapprox(C, C_expected) + @test isapprox(E, E_expected) + @test isapprox(U, U_expected) end -=# -end diff --git a/test/test_contact_2d_small_sliding.jl b/test/test_contact_2d_small_sliding.jl index a78b8b4..62406e7 100644 --- a/test/test_contact_2d_small_sliding.jl +++ b/test/test_contact_2d_small_sliding.jl @@ -19,24 +19,24 @@ function get_model(::Type{Val{Symbol("curved 2d contact small sliding")}}) upper = Problem(Elasticity, "upper", 2) upper.properties.formulation = :plane_stress upper.elements = create_elements(mesh, "UPPER") - update!(upper.elements, "youngs modulus", 96.0) - update!(upper.elements, "poissons ratio", 1/3) + update!(upper, "youngs modulus", 96.0) + update!(upper, "poissons ratio", 1/3) lower = Problem(Elasticity, "lower", 2) lower.properties.formulation = :plane_stress lower.elements = create_elements(mesh, "LOWER") - update!(lower.elements, "youngs modulus", 96.0) - update!(lower.elements, "poissons ratio", 1/3) + update!(lower, "youngs modulus", 96.0) + update!(lower, "poissons ratio", 1/3) bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement") bc_upper.elements = create_elements(mesh, "UPPER_TOP") - update!(bc_upper.elements, "displacement 1", 0.0) - update!(bc_upper.elements, "displacement 2", -0.15) + update!(bc_upper, "displacement 1", 0.0) + update!(bc_upper, "displacement 2", -0.15) bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement") bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM") - update!(bc_lower.elements, "displacement 1", 0.0) - update!(bc_lower.elements, "displacement 2", 0.0) + update!(bc_lower, "displacement 1", 0.0) + update!(bc_lower, "displacement 2", 0.0) interface = Problem(Contact, "contact between upper and lower block", 2, "displacement") interface.properties.dimension = 1 @@ -45,6 +45,7 @@ function get_model(::Type{Val{Symbol("curved 2d contact small sliding")}}) interface_master_elements = create_elements(mesh, "UPPER_BOTTOM") update!(interface_slave_elements, "master elements", interface_master_elements) interface.elements = [interface_master_elements; interface_slave_elements] + info("type of list is ", typeof(first(interface_slave_elements)("master elements", 0.0))) solver = Solver(Nonlinear) push!(solver, upper, lower, bc_upper, bc_lower, interface) @@ -77,24 +78,24 @@ function get_model(::Type{Val{Symbol("hertz contact, full 2d model")}}) upper = Problem(Elasticity, "CYLINDER", 2) upper.properties.formulation = :plane_strain upper.elements = create_elements(mesh, "CYLINDER") - update!(upper.elements, "youngs modulus", 70.0e3) - update!(upper.elements, "poissons ratio", 0.3) + update!(upper, "youngs modulus", 70.0e3) + update!(upper, "poissons ratio", 0.3) lower = Problem(Elasticity, "BLOCK", 2) lower.properties.formulation = :plane_strain lower.elements = create_elements(mesh, "BLOCK") - update!(lower.elements, "youngs modulus", 210.0e3) - update!(lower.elements, "poissons ratio", 0.3) + update!(lower, "youngs modulus", 210.0e3) + update!(lower, "poissons ratio", 0.3) # support block to ground bc_fixed = Problem(Dirichlet, "fixed", 2, "displacement") bc_fixed.elements = create_elements(mesh, "FIXED") - update!(bc_fixed.elements, "displacement 2", 0.0) + update!(bc_fixed, "displacement 2", 0.0) # symmetry line bc_sym_23 = Problem(Dirichlet, "symmetry line 23", 2, "displacement") bc_sym_23.elements = create_elements(mesh, "SYM23") - update!(bc_sym_23.elements, "displacement 1", 0.0) + update!(bc_sym_23, "displacement 1", 0.0) nid = find_nearest_nodes(mesh, [0.0, 100.0]) #load = Problem(Dirichlet, "load", 2, "displacement") @@ -102,7 +103,7 @@ function get_model(::Type{Val{Symbol("hertz contact, full 2d model")}}) load.properties.formulation = :plane_strain load.elements = [Element(Poi1, nid)] #update!(load.elements, "displacement 2", -10.0) - update!(load.elements, "displacement traction force 2", -35.0e3) + update!(load, "displacement traction force 2", -35.0e3) contact = Problem(Contact, "contact between block and cylinder", 2, "displacement") contact.properties.rotate_normals = true diff --git a/test/test_elasticity_3d_linear_with_surface_load.jl b/test/test_elasticity_3d_linear_with_surface_load.jl index b3426ec..b811d59 100644 --- a/test/test_elasticity_3d_linear_with_surface_load.jl +++ b/test/test_elasticity_3d_linear_with_surface_load.jl @@ -2,6 +2,7 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md using JuliaFEM +using JuliaFEM.Preprocess using JuliaFEM.Test @testset "test continuum 3d linear elasticity with surface load" begin @@ -38,9 +39,7 @@ using JuliaFEM.Test boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement") push!(boundary_problem, symxy, symxz, symyz) - solver = Solver("solve 3d block") - push!(solver, elasticity_problem) - push!(solver, boundary_problem) + solver = LinearSolver(elasticity_problem, boundary_problem) call(solver) disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) @@ -48,3 +47,39 @@ using JuliaFEM.Test u_expected = 2.0 * [-1/3, -1/3, 1.0] @test isapprox(disp, u_expected) end + +function solve_rod_model_elasticity(eltype) + fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med" + mesh = aster_read_mesh(fn, eltype) + element_sets = join(keys(mesh.element_sets), ", ") + info("element sets: $element_sets") + p1 = Problem(Elasticity, "rod", 3) + p2 = Problem(Elasticity, "trac", 3) + p3 = Problem(Dirichlet, "fixed", 3, "displacement") + p4 = Problem(Dirichlet, "fixed", 3, "displacement") + p5 = Problem(Dirichlet, "fixed", 3, "displacement") + p1.elements = create_elements(mesh, "ROD") + p2.elements = create_elements(mesh, "FACE2") + p3.elements = create_elements(mesh, "FACE1") + p4.elements = create_elements(mesh, "FACE3") + p5.elements = create_elements(mesh, "FACE5") + update!(p1, "youngs modulus", 96.0) + update!(p1, "poissons ratio", 1/3) + update!(p2, "displacement traction force 1", 96.0) + update!(p3, "displacement 1", 0.0) + update!(p4, "displacement 2", 0.0) + update!(p5, "displacement 3", 0.0) + solver = LinearSolver(p1, p2, p3, p4, p5) + call(solver) + u_max = maximum(p1.assembly.u) + info("$eltype, u_max = $u_max") + return u_max +end +@testset "compare 3d rod to CA solution" begin + @test isapprox(solve_rod_model_elasticity("Tet4"), 0.2) + @test isapprox(solve_rod_model_elasticity("Tet10"), 0.2) + @test isapprox(solve_rod_model_elasticity("Hex8"), 0.2) + @test isapprox(solve_rod_model_elasticity("Hex20"), 0.2) + @test isapprox(solve_rod_model_elasticity("Hex27"), 0.2) +end + diff --git a/test/test_elements.jl b/test/test_elements.jl index 8c018c4..57a0187 100644 --- a/test/test_elements.jl +++ b/test/test_elements.jl @@ -108,3 +108,13 @@ end @test isapprox(el("foo1", 1.5), el("foo2", 1.5)) end +@testset "add elements to elements" begin + el1 = Element(Seg2, [1, 2]) + el2 = Element(Seg2, [3, 4]) + update!(el1, "master elements", [el2]) + lst = el1("master elements", 0.0) + info("lst = ", el1["master elements"]) + info("typeof lst = ", typeof(lst)) + @test isa(lst, Vector) +end + diff --git a/test/test_heat.jl b/test/test_heat.jl index 2f2967b..62709f2 100644 --- a/test/test_heat.jl +++ b/test/test_heat.jl @@ -4,8 +4,42 @@ using JuliaFEM using JuliaFEM.Test using JuliaFEM.Preprocess +using JuliaFEM.Postprocess -@testset "test one element heat problem" begin +@testset "Tet10 + convection" begin + # For some reason Tet10 fails, maybe because of convection. + mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med" + mesh = aster_read_mesh(mesh_file, "Tet10") + prob = Problem(Heat, "tet", 1) + face = Problem(Heat, "face 4", 1) + fixed = Problem(Dirichlet, "fixed face 3", 1, "temperature") + prob.elements = create_elements(mesh, "TET") + update!(prob, "temperature thermal conductivity", 50.0) + face.elements = create_elements(mesh, "FACE4") + update!(face, "temperature external temperature", 20.0) + update!(face, "temperature heat transfer coefficient", 60.0) + fixed.elements = create_elements(mesh, "FACE2") + info("# of elements in fixed set: $(length(fixed))") + update!(fixed, "temperature 1", 0.0) + solver = LinearSolver(prob, face, fixed) + call(solver) + T = prob.assembly.u + info("Solution: $T") + T_expected = [ # using code aster + 1.45606533688540E+01 + 5.01315339269860E-17 + 3.02236827927507E-17 + -2.01049663215778E-16 + 1.05228712963739E+01 + 0.00000000000000E+00 + 9.44202309239159E+00 + 1.05228712963739E+01 + 4.44089209850063E-16 + 0.00000000000000E+00] + @test isapprox(T, T_expected; rtol=1.0e-6) +end + +@testset "one element heat problem" begin X = Dict{Int, Vector{Float64}}( 1 => [0.0,0.0], @@ -124,37 +158,189 @@ end end =# -@testset "test 3d heat problem" begin +@testset "compare simple 3d heat problem to code aster solution" begin fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med" - mesh = aster_read_mesh(fn, "SHORT_ROD_RECTANGLE_HEX8") + mesh = aster_read_mesh(fn, "Hex8") + element_sets = join(keys(mesh.element_sets), ", ") + info("element sets: $element_sets") p1 = Problem(Heat, "rod", 1) - push!(p1, create_elements(mesh, "ROD")) - push!(p1, create_elements(mesh, "SIDES")) - push!(p1, create_elements(mesh, "RIGHT")) - update!(p1, "temperature thermal conductivity", 50.0) - update!(p1, "temperature external temperature", 20.0) - update!(p1, "temperature heat transfer coefficient", 10.0) + rod = create_elements(mesh, "ROD") + face2 = create_elements(mesh, "FACE2") + face3 = create_elements(mesh, "FACE3") + face4 = create_elements(mesh, "FACE4") + face5 = create_elements(mesh, "FACE5") + face6 = create_elements(mesh, "FACE6") + update!(rod, "temperature thermal conductivity", 50.0) + update!(face2, "temperature external temperature", 20.0) + update!(face2, "temperature heat transfer coefficient", 60.0) + update!(face3, "temperature external temperature", 30.0) + update!(face3, "temperature heat transfer coefficient", 50.0) + update!(face4, "temperature external temperature", 40.0) + update!(face4, "temperature heat transfer coefficient", 40.0) + update!(face5, "temperature external temperature", 50.0) + update!(face5, "temperature heat transfer coefficient", 30.0) + update!(face6, "temperature external temperature", 60.0) + update!(face6, "temperature heat transfer coefficient", 20.0) + push!(p1, rod, face2, face3, face4, face5, face6) p2 = Problem(Dirichlet, "left support T=100", 1, "temperature") - push!(p2, create_elements(mesh, "LEFT")) + push!(p2, create_elements(mesh, "FACE1")) update!(p2, "temperature 1", 100.0) solver = LinearSolver(p1, p2) call(solver) - T_min = minimum(p1.assembly.u) + # fields extracted from Code Aster .resu file + TEMP = Dict{Int64, Float64}( + 1 => 1.00000000000000E+02, + 2 => 1.00000000000000E+02, + 3 => 1.00000000000000E+02, + 4 => 1.00000000000000E+02, + 5 => 3.01613322896279E+01, + 6 => 3.01263406641066E+01, + 7 => 3.02559777927923E+01, + 8 => 3.02209215997131E+01) + FLUX_ELGA = Dict{Int64, Vector{Float64}}( + 1 => [1.74565160615448E+04, -9.99903237329079E+01, -3.69874201221677E+01], + 2 => [1.74565160615448E+04, -3.73168968436642E+02, -1.38038931136833E+02], + 3 => [1.74428571293096E+04, -9.99903237329079E+01, -3.70268090662933E+01], + 4 => [1.74428571293096E+04, -3.73168968436642E+02, -1.38185932677561E+02], + 5 => [1.74615686370955E+04, -9.99509347888079E+01, -3.69874201221677E+01], + 6 => [1.74615686370955E+04, -3.73021966895897E+02, -1.38038931136833E+02], + 7 => [1.74479150854902E+04, -9.99509347888065E+01, -3.70268090662933E+01], + 8 => [1.74479150854901E+04, -3.73021966895874E+02, -1.38185932677561E+02]) + FLUX_NOEU = Dict{Int64, Vector{Float64}}( + 1 => [1.74596669275930E+04, 7.55555618070503E-11, 3.68594044175552E-12], + 2 => [1.74684148339734E+04, 1.10418341137120E-11, 3.48876483258209E-12], + 3 => [1.74360055518019E+04, 7.91828824731056E-11, 1.95399252334028E-13], + 4 => [1.74447696000717E+04, -3.49587025993969E-12, 3.55271367880050E-13], + 5 => [1.74596669275931E+04, -4.73227515822099E+02, -1.74958127606525E+02], + 6 => [1.74684148339733E+04, -4.72904678032251E+02, -1.74958127606524E+02], + 7 => [1.74360055518019E+04, -4.73227515822118E+02, -1.75280965396335E+02], + 8 => [1.74447696000717E+04, -4.72904678032179E+02, -1.75280965396335E+02]) - # Code Aster solution - T_CA_HEX20 = 4.58158267950429E+01 - T_CA_HEX8 = 3.77215189873436E+01 - info("T_min = $T_min") - info("T_acc = $(T_acc(0.2))") - rtol1 = norm(T_min-T_CA_HEX8)/max(T_min,T_CA_HEX8)*100.0 - rtol2 = norm(T_min-T_acc(0.2))/max(T_min,T_acc(0.2))*100.0 - info("rel. tol to CA solution: $rtol1 %") - info("rel. tol to accurate solution: $rtol2 %") + postprocessor = Postprocessor(p1) + flux = full(call(postprocessor)) + fluxd = Dict{Int64, Vector{Float64}}() + for j=1:8 + fluxd[j] = vec(flux[j,:]) + end + + T = p1("temperature") + + for j in sort(collect(keys(T))) + T1 = T[j][1] + T2 = TEMP[j] + rtol = norm(T1-T2)/max(T1,T2)*100.0 + @printf "node %i temp, JF: %e, CA: %e, rtol: %10.6f %%\n" j T1 T2 rtol + @test rtol < 1.0e-9 + end + + for j=1:8 + q1 = get_integration_points(first(rod))[j]("heat flux", 0.0) + q2 = FLUX_ELGA[j] + rtol = norm(q1-q2)/max(norm(q1),norm(q2))*100.0 + @printf "ip %i flux, JF: (% e,% e,% e), CA: (% e,% e,% e), rtol: %10.6f %%\n" j q1... q2... rtol + # @test rtol < 0.05 + # testing in integration points makes no sense because they are in different order in CA + end + + for j in sort(collect(keys(fluxd))) + q1 = fluxd[j] + q2 = FLUX_NOEU[j] + rtol = norm(q1-q2)/max(norm(q1),norm(q2))*100.0 + @printf "node %i flux, JF: (% e,% e,% e), CA: (% e,% e,% e), rtol: %10.6f %%\n" j q1... q2... rtol + @test rtol < 1.0e-9 + end + +end + +@testset "compare simple 3d heat problem to analytical solution" begin + function calc_3d_heat_model(mesh_name) + fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med" + mesh = aster_read_mesh(fn, mesh_name) + p1 = Problem(Heat, "rod", 1) + p2 = Problem(Dirichlet, "left support T=100", 1, "temperature") + p1.elements = create_elements(mesh, "ROD", "FACE2") + p2.elements = create_elements(mesh, "FACE1") + update!(p1, "temperature thermal conductivity", 100.0) + update!(p1, "temperature external temperature", 0.0) + update!(p1, "temperature heat transfer coefficient", 1000.0) + update!(p2, "temperature 1", 100.0) + solver = LinearSolver(p1, p2) + call(solver) + T_min = minimum(p1.assembly.u) + return T_min + end + for model in ["Tet4", "Tet10", "Hex8", "Hex20", "Hex27"] + Tmin = calc_3d_heat_model(model) + Tacc = 100/3 + rtol = norm(Tmin-Tacc)/max(Tmin,Tacc)*100.0 + @printf "%-10s : Tmin = % g, Tacc = % g, rtol = %g %%\n" model Tmin Tacc rtol + @test isapprox(Tmin, 100/3) + end +end + +@testset "compare simple 3d heat problem to code aster solution" begin + + function calc_3d_heat_model(mesh_name) + fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med" + mesh = aster_read_mesh(fn, mesh_name) + element_sets = join(keys(mesh.element_sets), ", ") + info("element sets: $element_sets") + # x -> FACE1 ... FACE2 + # y -> FACE3 ... FACE4 + # z -> FACE5 ... FACE6 + # rod has longer dimension in x direction, first face comes + # first in corresponding axis direction + p1 = Problem(Heat, "rod", 1) + rod = create_elements(mesh, "ROD") + face2 = create_elements(mesh, "FACE2") + face3 = create_elements(mesh, "FACE3") + face4 = create_elements(mesh, "FACE4") + face5 = create_elements(mesh, "FACE5") + face6 = create_elements(mesh, "FACE6") + update!(rod, "temperature thermal conductivity", 50.0) + update!(face2, "temperature external temperature", 20.0) + update!(face2, "temperature heat transfer coefficient", 60.0) + update!(face3, "temperature external temperature", 30.0) + update!(face3, "temperature heat transfer coefficient", 50.0) + update!(face4, "temperature external temperature", 40.0) + update!(face4, "temperature heat transfer coefficient", 40.0) + update!(face5, "temperature external temperature", 50.0) + update!(face5, "temperature heat transfer coefficient", 30.0) + update!(face6, "temperature external temperature", 60.0) + update!(face6, "temperature heat transfer coefficient", 20.0) + push!(p1, rod, face2, face3, face4, face5, face6) + p2 = Problem(Dirichlet, "left support T=100", 1, "temperature") + p2.elements = create_elements(mesh, "FACE1") + update!(p2, "temperature 1", 100.0) + solver = LinearSolver(p1, p2) + call(solver) + return p1.assembly.u + end + + CA_sol = Dict( + "Tet4" => 3.01872246268290E+01, + "Hex8" => 3.01263406641066E+01, + "Tet10" => 4.38924023356612E+01, + "Hex20" => 4.57539800177123E+01, + "Hex27" => 4.57760386068096E+01) + + models = ["Tet4", "Hex8", "Hex20", "Hex27", "Tet10"] + + for model in models + T = calc_3d_heat_model(model) + T_min = minimum(T) + T_ca = CA_sol[model] + rtol = norm(T_min-T_ca)/max(T_min,T_ca)*100.0 + @printf "%-10s : T_min = % g, T_ca = % g, rtol = %g %%\n" model T_min T_ca rtol + if rtol > 1.0e-9 + info("Solution vector") + dump(T) + end + @test rtol < 1.0e-9 + end - @test isapprox(T_min, T_acc(0.2); rtol=18.0e-2) - @test isapprox(T_min, T_CA_HEX8; rtol=1.0e-9) end diff --git a/test/test_lagrange.jl b/test/test_lagrange.jl new file mode 100644 index 0000000..f50995d --- /dev/null +++ b/test/test_lagrange.jl @@ -0,0 +1,51 @@ +# 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.Test + +ALL_ELEMENTS = [ + Seg2, Seg3, + Tri3, Tri6, Quad4, Quad8, Quad9, + Tet4, Tet10, Hex8, Hex20, Hex27 +] + +@testset "Evaluating basis" begin + for T in ALL_ELEMENTS + el = Element(T) + nnodes = length(el) + for (i, X) in enumerate(get_reference_coordinates(T)) + Ni = vec(el(X)) + expected = zeros(nnodes) + expected[i] = 1.0 + @test isapprox(Ni, expected) + end + end +end + +function get_volume{T<:AbstractElement}(::Type{T}) + X = get_reference_coordinates(T) + element = Element(T) + update!(element, "geometry", X) + V = 0.0 + for ip in get_integration_points(element) + V += ip.weight*element(ip, 0.0, Val{:detJ}) + end + return V +end + +@testset "Calculate reference element length/area/volume" begin + @test isapprox(get_volume(Seg2), 2.0) + @test isapprox(get_volume(Seg3), 2.0) + @test isapprox(get_volume(Tri3), 0.5) + @test isapprox(get_volume(Tri6), 0.5) + @test isapprox(get_volume(Quad4), 2.0^2) + @test isapprox(get_volume(Quad8), 2.0^2) + @test isapprox(get_volume(Quad9), 2.0^2) + @test isapprox(get_volume(Tet4), 1/6) + @test isapprox(get_volume(Tet10), 1/6) + @test isapprox(get_volume(Hex8), 2.0^3) + @test isapprox(get_volume(Hex20), 2.0^3) + @test isapprox(get_volume(Hex27), 2.0^3) +end + diff --git a/test/test_preprocess.jl b/test/test_preprocess.jl new file mode 100644 index 0000000..b94124d --- /dev/null +++ b/test/test_preprocess.jl @@ -0,0 +1,85 @@ +# 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.Test + +@testset "renumber element nodes" begin + mesh = Mesh() + add_element!(mesh, 1, :Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) + mapping = Dict{Symbol, Vector{Int}}( + :Tet10 => [1, 2, 4, 3, 5, 6, 7, 8, 9, 10]) + reorder_element_connectivity!(mesh, mapping) + @test mesh.elements[1] == [1, 2, 4, 3, 5, 6, 7, 8, 9, 10] + invmapping = Dict{Symbol, Vector{Int}}() + invmapping[:Tet10] = invperm(mapping[:Tet10]) + reorder_element_connectivity!(mesh, invmapping) + @test mesh.elements[1] == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] +end + +function get_volume(element::Element, time=0.0) + V = 0.0 + for ip in get_integration_points(element) + V += ip.weight*element(ip, time, Val{:detJ}) + end + return V +end + +function get_volume(elements::Vector{Element}, time=0.0) + return sum([get_volume(element, time) for element in elements]) +end + +#= +@testset "Hex8 element connectivity order" begin + fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med" + mesh = aster_read_mesh(fn, "SHORT_ROD_RECTANGLE_HE8_1ELEM") + # 1. check volume of element + rod = create_elements(mesh, "ROD") + V = get_volume(rod) + V_expected = 0.01^2*0.2 + info("Volume of rod = $V, expected = $V_expected") + @test isapprox(V, V_expected) + # 2. put some field value and calculate flux in gauss points + T = Dict{Int64, Float64}( + 1 => 100.0, 2 => 100.0, 3 => 100.0, 4 => 100.0, + 5 => 200.0, 6 => 300.0, 7 => 400.0, 8 => 500.0) + update!(rod, "temperature", T) + # it has been verified using code aster that flux in integration + # points is + FLUX_ELGA = Dict{Int, Vector{Float64}}( + 1 => [-4.08493649053890E+04, -2.11324865405187E+05, 1.05662432702594E+05], + 2 => [-4.08493649053890E+04, -7.88675134594813E+05, 3.94337567297406E+05], + 3 => [-6.97168783648703E+04, -2.11324865405187E+05, 1.05662432702594E+05], + 4 => [-6.97168783648703E+04, -7.88675134594813E+05, 3.94337567297406E+05], + 5 => [-5.52831216351297E+04, -2.11324865405187E+05, 1.05662432702594E+05], + 6 => [-5.52831216351297E+04, -7.88675134594813E+05, 3.94337567297406E+05], + 7 => [-8.41506350946110E+04, -2.11324865405187E+05, 1.05662432702594E+05], + 8 => [-8.41506350946110E+04, -7.88675134594813E+05, 3.94337567297406E+05]) + # flux is q̄(ξ) = -k∇T + element = first(rod) + k = -50.0 + flux(xi, time) = -k*vec(element("temperature", xi, time, Val{:Grad})) + weights = ones(8) + # code aster integration points (FPG8) + a = -1.0/sqrt(3.0) + points = Vector{Float64}[ + [-a, -a, -a], + [-a, -a, a], + [-a, a, -a], + [-a, a, a], + [ a, -a, -a], + [ a, -a, a], + [ a, a, -a], + [ a, a, a]] + for i=1:8 + q1 = flux(points[i], 0.0) + q2 = FLUX_ELGA[i] + rtol = norm(q1-q2)/max(norm(q1),norm(q2))*100.0 + @printf "ip %i flux, JF: (% e,% e,% e), CA: (% e,% e,% e), rtol: %10.6f %%\n" i q1... q2... rtol + @test rtol < 0.05 + end +end +=# + diff --git a/test/test_preprocess_aster_reader.jl b/test/test_preprocess_aster_reader.jl index dba65e8..c106ed4 100644 --- a/test/test_preprocess_aster_reader.jl +++ b/test/test_preprocess_aster_reader.jl @@ -150,3 +150,33 @@ end @test haskey(mesh2.element_sets, "BLOCK") @test length(mesh2.elements) == 1 end + +function calculate_volume(eltype::Symbol) + fn = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med" + mesh = aster_read_mesh(fn, "$eltype") + elements = create_elements(mesh, eltype) + V = 0.0 + time = 0.0 + for element in elements + for ip in get_integration_points(element) + detJ = element(ip, time, Val{:detJ}) + detJ > 0 || warn("negative determinant for element $eltype !") + V += ip.weight*detJ + end + end + info("volume of $eltype is $V") + return V +end + +@testset "calculate volume for primitives" begin + @test isapprox(calculate_volume(:Tet4), 1/6) + @test isapprox(calculate_volume(:Tet10), 1/6) + @test isapprox(calculate_volume(:Hex8), 2^3) + @test isapprox(calculate_volume(:Hex20), 2^3) + @test isapprox(calculate_volume(:Hex27), 2^3) +# @test isapprox(get_volume("PE6"), V) +# @test isapprox(get_volume("PY5"), V) +# @test isapprox(get_volume("P15"), V) +# @test isapprox(get_volume("P13"), V) +end +