mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-21 02:18:56 +00:00
lot of tests
This commit is contained in:
Jukka Aho
@@ -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
|
||||
|
||||
+24
-26
@@ -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
|
||||
|
||||
|
||||
+169
-147
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
+208
-22
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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
|
||||
=#
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
Reference in New Issue
Block a user