Fix tests

* Fix deprecation warnings from tests
* Refactor tests so that ´@testset` is usually called in master file
  `runtests.jl`, not inside test file. Later on we can convert tests
  to examples.
* Syntax of tests now follow more closely syntax used currently in
  JuliaFEM. We have had earlier studies with different kind of syntaxes,
  now we have kind of explicit way to do things.
This commit is contained in:
Jukka Aho
2018-09-06 12:08:16 +03:00
parent 52be9e546b
commit ca7e2904cf
53 changed files with 2285 additions and 2581 deletions
+64 -52
View File
@@ -1,8 +1,7 @@
# 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.Testing
using JuliaFEM, SparseArrays, Test
#=
In [36]: C = Matrix([[0], [30], [15]]) # node coordinates
@@ -11,13 +10,13 @@ In [38]: N = P.T*A.inv()
In [39]: Me = integrate(N.T*N, (x, 0, 30))
In [40]: De = diag(*integrate(N, (x, 0, 30)))
In [41]: Me
Out[41]:
Out[41]:
Matrix([
[ 4, -1, 2],
[-1, 4, 2],
[ 2, 2, 16]])
In [42]: De
Out[42]:
Out[42]:
Matrix([
[5, 0, 0],
[0, 5, 0],
@@ -25,39 +24,55 @@ Matrix([
=#
@testset "dirichlet problem in 1 dimension" begin
element = Element(Seg2, [1, 2])
element["geometry"] = Vector{Float64}[[0.0, 0.0], [6.0, 0.0]]
element["temperature 1"] = 0.0
p1 = Problem(Dirichlet, "test problem 1", 1, "temperature")
p1.properties.dual_basis = false
p2 = Problem(Dirichlet, "test problem 2", 1, "temperature")
p2.properties.dual_basis = true
assemble!(p1, element)
assemble!(p2, element)
C1 = full(p1.assembly.C1)
C2 = full(p1.assembly.C2)
time = 0.0
element = Element(Seg2, (1, 2))
X = Dict(1 => [0.0, 0.0], 2 => [6.0, 0.0])
update!(element, "geometry", X)
update!(element, "temperature 1", 0.0)
problem1 = Problem(Dirichlet, "test problem 1", 1, "temperature")
problem1.properties.variational = true
problem1.properties.dual_basis = false
add_element!(problem1, element)
assemble!(problem1, time)
C1 = problem1.assembly.C1
C2 = problem1.assembly.C2
@test isapprox(C1, C2)
@test isapprox(C1, [2.0 1.0; 1.0 2.0])
C1 = full(p2.assembly.C1)
C2 = full(p2.assembly.C2)
problem2 = Problem(Dirichlet, "test problem 2", 1, "temperature")
problem2.properties.variational = true
problem2.properties.dual_basis = true
add_element!(problem2, element)
assemble!(problem2, time)
C1 = problem2.assembly.C1
C2 = problem2.assembly.C2
@test isapprox(C1, C2)
@test isapprox(C1, [3.0 0.0; 0.0 3.0])
element = Element(Seg3, [1, 2, 3])
element["geometry"] = Vector{Float64}[[0.0, 0.0], [30.0, 0.0], [15.0, 0.0]]
element["temperature 1"] = 0.0
p1 = Problem(Dirichlet, "quadratic 1", 1, "temperature")
p1.properties.dual_basis = false
p2 = Problem(Dirichlet, "quadratic 1", 1, "temperature")
p2.properties.dual_basis = true
assemble!(p1, element)
assemble!(p2, element)
C1 = full(p1.assembly.C1)
C2 = full(p1.assembly.C2)
element = Element(Seg3, (1, 2, 3))
X = Dict(1 => [0.0, 0.0], 2 => [30.0, 0.0], 3 => [15.0, 0.0])
update!(element, "geometry", X)
update!(element, "temperature 1", 0.0)
problem3 = Problem(Dirichlet, "quadratic 1", 1, "temperature")
problem3.properties.variational = true
problem3.properties.dual_basis = false
add_element!(problem3, element)
assemble!(problem3, time)
C1 = problem3.assembly.C1
C2 = problem3.assembly.C2
@test isapprox(C1, C2)
@test isapprox(C1, [4.0 -1.0 2.0; -1.0 4.0 2.0; 2.0 2.0 16.0])
C1 = full(p2.assembly.C1)
C2 = full(p2.assembly.C2)
problem4 = Problem(Dirichlet, "quadratic 2", 1, "temperature")
problem4.properties.variational = true
problem4.properties.dual_basis = true
add_element!(problem4, element)
assemble!(problem4, time)
C1 = problem4.assembly.C1
C2 = problem4.assembly.C2
@test isapprox(C1, C2)
@test isapprox(C1, [5.0 0.0 0.0; 0.0 5.0 0.0; 0.0 0.0 20.0])
end
@@ -115,32 +130,29 @@ end
=#
@testset "test analytical boundary condition" begin
X = Dict{Int64, Vector{Float64}}(
1 => [0.0, 0.0],
2 => [1.0, 0.0])
element = Element(Seg2, [1, 2])
X = Dict(1 => [0.0, 0.0],
2 => [1.0, 0.0])
element = Element(Seg2, (1, 2))
update!(element, "geometry", X)
update!(element, "displacement 1", 0.0)
f(xi, time) = begin
info("function called at xi = $xi, time = $time")
X = element("geometry", xi, time)
info("geometry at xi, X = $X")
val = X[1]*time
info("result for field at xi = $val")
function f(element, ip, time)
x, y = element("geometry", ip, time)
val = x*time
@debug("analytical function called", ip, time, x, y, val)
return val
end
update!(element, "displacement 1", 0.0)
update!(element, "displacement 2", f)
p = Problem(Dirichlet, "test boundary", 2, "displacement")
push!(p, element)
assemble!(p, 0.0)
g1 = full(p.assembly.g, 4, 1)
@test isapprox(g1, [0.0, 0.0, 0.0, 0.0])
empty!(p.assembly)
assemble!(p, 1.0)
g2 = full(p.assembly.g, 4, 1)
C2 = full(p.assembly.C2, 4, 4)
problem = Problem(Dirichlet, "test boundary", 2, "displacement")
add_element!(problem, element)
time = 0.0
assemble!(problem, time)
@test isapprox(problem.assembly.g, [0.0, 0.0, 0.0, 0.0])
empty!(problem.assembly)
time = 1.0
assemble!(problem, time)
g2 = Vector(problem.assembly.g, 4)
C2 = Matrix(problem.assembly.C2, 4, 4)
u = C2 \ g2
info("u = $u")
@debug("displacement vector", u)
@test isapprox(u, [0.0, 0.0, 0.0, 1.0])
end
+15 -22
View File
@@ -1,26 +1,19 @@
# 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.Testing
using JuliaFEM, Test
@testset "1d strain" begin
X = Dict{Int64, Vector{Float64}}(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 1.0, 1.0])
u = Dict{Int64, Vector{Float64}}(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 1.0, 1.0])
element = Element(Seg2, [1, 2])
update!(element, "geometry", X)
detJ = element([0.0], 0.0, Val{:detJ})
info("detJ = $detJ")
@test isapprox(detJ, sqrt(3)/2)
J = element([0.0], 0.0, Val{:Jacobian})
info("J = $J")
@test isapprox(J, [0.5 0.5 0.5])
update!(element, "displacement", u)
# FIXME
# gradu = element("displacement", [0.0], 0.0, Val{:Grad})
# info("1d bar: ∇u = $gradu")
end
# 1d strain
X = Dict(1 => [0.0, 0.0, 0.0], 2 => [1.0, 1.0, 1.0])
u = Dict(1 => [0.0, 0.0, 0.0], 2 => [1.0, 1.0, 1.0])
element = Element(Seg2, (1, 2))
update!(element, "geometry", X)
update!(element, "displacement", u)
xi, time = (0.0,), 0.0
detJ = element(xi, time, Val{:detJ})
J = element(xi, time, Val{:Jacobian})
# gradu = element("displacement", xi, time, Val{:Grad})
@debug("1d seg2 info", xi ,time, detJ, J)
@test isapprox(detJ, sqrt(3)/2)
@test isapprox(J, [0.5 0.5 0.5])
@@ -1,52 +1,51 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using Base.Test
using JuliaFEM
using JuliaFEM, Test
@testset "2d linear elasticity + volume load + surface load" begin
# Example of 2d linear elasticity + volume load + surface load
X = Dict(1 => [0.0, 0.0],
2 => [1.0, 0.0],
3 => [1.0, 1.0],
4 => [0.0, 1.0])
# Dictionary containing node coordinates
X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [1.0, 1.0], 4 => [0.0, 1.0])
props = ("formulation" => "plane_stress",
"finite_strain" => "false",
"geometric_stiffness" => "false")
# Create new problem of type `Elasticity`
block = Problem(Elasticity, "block", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = false
block.properties.geometric_stiffness = false
# field problem
block = Problem(Elasticity, "BLOCK", 2)
block.elements = [Element(Quad4, [1, 2, 3, 4])]
update!(block.properties, props...)
update!(block.elements, "geometry", X)
update!(block.elements, "youngs modulus", 288.0)
update!(block.elements, "poissons ratio", 1/3)
update!(block.elements, "displacement load 2", 576.0)
# Add volume element
element = Element(Quad4, (1, 2, 3, 4))
update!(element, "geometry", X)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
update!(element, "displacement load 2", 576.0)
add_element!(block, element)
# traction
traction = Problem(Elasticity, "TRACTION", 2)
traction.elements = [Element(Seg2, [3, 4])]
update!(traction.properties, props...)
update!(traction.elements, "geometry", X)
update!(traction.elements, "displacement traction force 2", 288.0)
# Add boundary element for tractoin force
traction_element = Element(Seg2, (3, 4))
update!(traction_element, "geometry", X)
update!(traction_element, "displacement traction force 2", 288.0)
add_element!(block, traction_element)
# boundary conditions
bc = Problem(Dirichlet, "symmetry boundary conditions", 2, "displacement")
bc.elements = [Element(Seg2, [1, 2]), Element(Seg2, [4, 1])]
update!(bc.elements, "geometry", X)
update!(bc.elements[1], "displacement 2", 0.0)
update!(bc.elements[2], "displacement 1", 0.0)
# Define boundary conditions
bc = Problem(Dirichlet, "symmetry boundary conditions", 2, "displacement")
bc_elements = [Element(Seg2, (1, 2)), Element(Seg2, (4, 1))]
update!(bc_elements, "geometry", X)
update!(bc_elements[1], "displacement 2", 0.0)
update!(bc_elements[2], "displacement 1", 0.0)
add_elements!(bc, bc_elements)
solver = Solver(Linear, "solve 2d linear elasticity problem")
add_problems!(solver, [block, traction, bc])
solve!(solver, 0.0)
# Create analysis, add problems to it and run analysis
analysis = Analysis(Linear)
add_problems!(analysis, block, bc)
run!(analysis)
f = 288.0
g = 576.0
E = 288.0
nu = 1/3
u3 = block("displacement", 0.0)[3]
u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1]
@test isapprox(u3, u3_expected)
end
# Analytical solution is known:
f = 288.0
g = 576.0
E = 288.0
nu = 1/3
u3 = block("displacement", 0.0)[3]
u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1]
@test isapprox(u3, u3_expected)
@@ -1,72 +1,43 @@
# 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
using JuliaFEM, Test
@testset "2d nonlinear elasticity: test nonhomogeneous boundary conditions and stress calculation" begin
# Example of 2d nonlinear elasticity with non-homogeneous boundary contitions
# field problem
block = Problem(Elasticity, "BLOCK", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = true
block.properties.geometric_stiffness = true
# Geometry of nodes
X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [1.0, 1.0], 4 => [0.0, 1.0])
nodes = Dict{Int, Vector{Float64}}(
1 => [0.0, 0.0],
2 => [1.0, 0.0],
3 => [1.0, 1.0],
4 => [0.0, 1.0])
block = Problem(Elasticity, "block", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = true
block.properties.geometric_stiffness = true
element = Element(Quad4, [1, 2, 3, 4])
update!(element, "geometry", nodes)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
push!(block, element)
element = Element(Quad4, (1, 2, 3, 4))
update!(element, "geometry", X)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
add_element!(block, element)
# boundary conditions
bc = Problem(Dirichlet, "bc", 2, "displacement")
bel1 = Element(Seg2, [1, 2])
bel2 = Element(Seg2, [3, 4])
bel3 = Element(Seg2, [4, 1])
update!([bel1, bel2, bel3], "geometry", nodes)
update!(bel1, "displacement 2", 0.0)
update!(bel2, "displacement 2", 0.5)
update!(bel3, "displacement 1", 0.0)
push!(bc, bel1, bel2, bel3)
# boundary conditions
bc = Problem(Dirichlet, "bc", 2, "displacement")
bel1 = Element(Seg2, (1, 2))
bel2 = Element(Seg2, (3, 4))
bel3 = Element(Seg2, (4, 1))
update!((bel1, bel2, bel3), "geometry", X)
update!(bel1, "displacement 2", 0.0)
update!(bel2, "displacement 2", 0.5)
update!(bel3, "displacement 1", 0.0)
add_elements!(bc, bel1, bel2, bel3)
solver = NonlinearSolver("solve block problem")
push!(solver, block, bc)
solver()
analysis = Analysis(Nonlinear)
add_problems!(analysis, block, bc)
run!(analysis)
# from code aster
eps_expected = [-2.08333312468287E-01, 6.25000000000000E-01, 0.0]
sig_expected = [ 4.50685020821470E-06, 4.62857140373777E+02, 0.0]
u3_expected = [-2.36237356855269E-01, 5.00000000000000E-01]
# results are verified using Code Aster
eps_expected = [-2.08333312468287E-01, 6.25000000000000E-01, 0.0]
sig_expected = [ 4.50685020821470E-06, 4.62857140373777E+02, 0.0]
u3_expected = [-2.36237356855269E-01, 5.00000000000000E-01]
u3 = reshape(block.assembly.u, 2, 4)[:, 3]
info("u3 = $u3")
@test isapprox(u3, u3_expected, atol=1.0e-5)
#= TODO: to postprocess
info("strain")
for ip in get_integration_points(element)
eps = ip("strain")
@printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] eps[1] eps[2] eps[3]
@test isapprox(eps, eps_expected, atol=1.0e-5)
end
info("cauchy stress")
for ip in get_integration_points(element)
sig = ip("cauchy stress")
@printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3]
@test isapprox(sig, sig_expected)
end
info("pk2 stress")
for ip in get_integration_points(element)
sig = ip("pk2 stress")
@printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3]
@test isapprox(sig, sig_expected)
end
=#
end
u3 = block("displacement", 0.0)[3]
@test isapprox(u3, u3_expected, atol=1.0e-5)
@@ -1,65 +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.Preprocess
using JuliaFEM.Testing
using JuliaFEM, Test
@testset "test 2d nonlinear elasticity with surface load" begin
meshfile = "/geometry/2d_block/BLOCK_1elem.med"
mesh = aster_read_mesh(dirname(@__DIR__)*meshfile)
# Example, 2d nonlinear elasticity with surface load (mesh is read from file)
# field problem
block = Problem(Elasticity, "BLOCK", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = true
block.properties.geometric_stiffness = true
meshfile = joinpath("test_elasticity_2d_nonlinear_with_surface_load", "BLOCK_1elem.med")
mesh = aster_read_mesh(meshfile)
block.elements = create_elements(mesh, "BLOCK")
update!(block.elements, "youngs modulus", 288.0)
update!(block.elements, "poissons ratio", 1/3)
update!(block.elements, "displacement load 2", 576.0)
# define field problem
block = Problem(Elasticity, "BLOCK", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = true
block.properties.geometric_stiffness = true
traction = create_elements(mesh, "TOP")
update!(traction, "displacement traction force 2", 288.0)
push!(block, traction...)
# Add volume elements
block_elements = create_elements(mesh, "BLOCK")
update!(block_elements, "youngs modulus", 288.0)
update!(block_elements, "poissons ratio", 1/3)
update!(block_elements, "displacement load 2", 576.0)
add_elements!(block, block_elements)
# boundary conditions
bc_sym = Problem(Dirichlet, "symmetry bc", 2, "displacement")
bc_elements_left = create_elements(mesh, "LEFT")
bc_elements_bottom = create_elements(mesh, "BOTTOM")
update!(bc_elements_left, "displacement 1", 0.0)
update!(bc_elements_bottom, "displacement 2", 0.0)
push!(bc_sym, bc_elements_left..., bc_elements_bottom...)
# Add surface elements
traction_elements = create_elements(mesh, "TOP")
update!(traction_elements, "displacement traction force 2", 288.0)
add_elements!(block, traction_elements)
solver = Solver(Nonlinear, "solve block problem")
add_problems!(solver, [block, bc_sym])
solve!(solver, 0.0)
# Create boundary problem, add boundary conditions:
bc_sym = Problem(Dirichlet, "symmetry bc", 2, "displacement")
bc_elements_left = create_elements(mesh, "LEFT")
bc_elements_bottom = create_elements(mesh, "BOTTOM")
update!(bc_elements_left, "displacement 1", 0.0)
update!(bc_elements_bottom, "displacement 2", 0.0)
add_elements!(bc_sym, bc_elements_left)
add_elements!(bc_sym, bc_elements_bottom)
# from code aster
u3_expected = [-4.92316106779943E-01, 7.96321884292103E-01]
eps_zz = -3.71128811855451E-01
eps_expected = [-3.71128532282463E-01, 1.11338615599337E+00, 0.0]
sig_expected = [ 3.36174888827909E-05, 2.23478729403118E+03, 0.0]
# Create analysis, add problems to analysis and run analysis:
analysis = Analysis(Nonlinear)
add_problems!(analysis, block, bc_sym)
run!(analysis)
u3 = reshape(block.assembly.u, 2, 4)[:, 3]
info("u3 = $u3")
@test isapprox(u3, u3_expected, atol=1.0e-5)
#= TODO: Test postprocessing in separate test
info("strain")
for ip in get_integration_points(block.elements[1])
eps = ip("strain")
@printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] eps[1] eps[2] eps[3]
@test isapprox(eps, eps_expected)
end
info("stress")
for ip in get_integration_points(block.elements[1])
sig = ip("stress")
@printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3]
#@test isapprox(sig, sig_expected)
end
=#
end
# Results from Code Aster:
u3_expected = [-4.92316106779943E-01, 7.96321884292103E-01]
eps_zz = -3.71128811855451E-01
eps_expected = [-3.71128532282463E-01, 1.11338615599337E+00, 0.0]
sig_expected = [ 3.36174888827909E-05, 2.23478729403118E+03, 0.0]
u3 = block("displacement", 0.0)[3]
@debug("displacement", u3)
@test isapprox(u3, u3_expected, atol=1.0e-5)
@@ -1,49 +1,38 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
# http://ahojukka5.github.io/posts/finite-element-solution-for-one-element-problem/
using JuliaFEM, Test
using JuliaFEM
using JuliaFEM: add_elements!
using Base.Test
# Local stiffness matrix of plane stress element
@testset "test 2d linear elasticity local matrices" begin
element = Element(Quad4, [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])
u = Dict{Int64, Vector{Float64}}(
1 => [0.0, 0.0],
2 => [0.0, 0.0],
3 => [0.0, 0.0],
4 => [0.0, 0.0])
update!(element, "geometry", X)
update!(element, "displacement", u)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
update!(element, "displacement load", [4.0, 8.0])
element = Element(Quad4, (1, 2, 3, 4))
X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [1.0, 1.0], 4 => [0.0, 1.0])
u = Dict(1 => [0.0, 0.0], 2 => [0.0, 0.0], 3 => [0.0, 0.0], 4 => [0.0, 0.0])
update!(element, "geometry", X)
update!(element, "displacement", u)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
update!(element, "displacement load", [4.0, 8.0])
problem = Problem(Elasticity, "[0x1] x [0x1] block", 2)
update!(problem.properties, "formulation" => "plane_stress")
add_elements!(problem, [element])
assemble!(problem)
K = full(problem.assembly.K)
f = vec(full(problem.assembly.f))
problem = Problem(Elasticity, "[0x1] x [0x1] block", 2)
problem.properties.formulation = :plane_stress
add_element!(problem, element)
assemble!(problem, 0.0)
K_expected = [
144 54 -90 0 -72 -54 18 0
54 144 0 18 -54 -72 0 -90
-90 0 144 -54 18 0 -72 54
0 18 -54 144 0 -90 54 -72
-72 -54 18 0 144 54 -90 0
-54 -72 0 -90 54 144 0 18
18 0 -72 54 -90 0 144 -54
0 -90 54 -72 0 18 -54 144]
K = Matrix(problem.assembly.K)
f = Vector(problem.assembly.f)
f_expected = [1, 2, 1, 2, 1, 2, 1, 2]
K_expected = [
144 54 -90 0 -72 -54 18 0
54 144 0 18 -54 -72 0 -90
-90 0 144 -54 18 0 -72 54
0 18 -54 144 0 -90 54 -72
-72 -54 18 0 144 54 -90 0
-54 -72 0 -90 54 144 0 18
18 0 -72 54 -90 0 144 -54
0 -90 54 -72 0 18 -54 144]
@test isapprox(K, K_expected)
@test isapprox(f, f_expected)
end
f_expected = [1, 2, 1, 2, 1, 2, 1, 2]
@test isapprox(K, K_expected)
@test isapprox(f, f_expected)
+42 -52
View File
@@ -1,59 +1,49 @@
# 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
using JuliaFEM, Test
@testset "test 2d nonlinear residual" 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])
u = Dict{Int64, Vector{Float64}}(
1 => [0.1, 0.2],
2 => [0.3, 0.4],
3 => [0.5, 0.6],
4 => [0.7, 0.8])
T = Dict{Int64, Vector{Float64}}(
3 => [0.0, 288.0],
4 => [0.0, 288.0])
element = Element(Quad4, [1, 2, 3, 4])
update!(element, "geometry", X)
update!(element, "displacement", u)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
traction = Element(Seg2, [3, 4])
update!(traction, "geometry", X)
update!(traction, "displacement", u)
update!(traction, "displacement traction force", T)
# Test stiffness matrix of geometrically nonlinear problem
# field problem
block = Problem(Elasticity, "BLOCK", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = true
block.properties.geometric_stiffness = true
push!(block, element)
#push!(block, traction)
assemble!(block, 0.0)
Km = full(block.assembly.K)
Kg = full(block.assembly.Kg)
K = Km + Kg
f = full(block.assembly.f)
X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0], 3 => [1.0, 1.0], 4 => [0.0, 1.0])
u = Dict(1 => [0.1, 0.2], 2 => [0.3, 0.4], 3 => [0.5, 0.6], 4 => [0.7, 0.8])
T = Dict(3 => [0.0, 288.0], 4 => [0.0, 288.0])
K_expected = [
401.76 200.88 -123.12 -5.76 -191.52 -117.36 -87.12 -77.76
200.88 473.76 -5.76 28.08 -117.36 -205.92 -77.76 -295.92
-123.12 -5.76 197.28 -2.16 -12.24 -25.92 -61.92 33.84
-5.76 28.08 -2.16 298.08 -25.92 -163.44 33.84 -162.72
-191.52 -117.36 -12.24 -25.92 240.48 120.24 -36.72 23.04
-117.36 -205.92 -25.92 -163.44 120.24 312.48 23.04 56.88
-87.12 -77.76 -61.92 33.84 -36.72 23.04 185.76 20.88
-77.76 -295.92 33.84 -162.72 23.04 56.88 20.88 401.76]
element = Element(Quad4, (1, 2, 3, 4))
update!(element, "geometry", X)
update!(element, "displacement", u)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
traction = Element(Seg2, (3, 4))
update!(traction, "geometry", X)
update!(traction, "displacement", u)
update!(traction, "displacement traction force", T)
# field problem
block = Problem(Elasticity, "block", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = true
block.properties.geometric_stiffness = true
add_element!(block, element)
#add_elements!(block, traction)
assemble!(block, 0.0)
Km = Matrix(block.assembly.K)
Kg = Matrix(block.assembly.Kg)
K = Km + Kg
f = Vector(block.assembly.f)
K_expected = [
401.76 200.88 -123.12 -5.76 -191.52 -117.36 -87.12 -77.76
200.88 473.76 -5.76 28.08 -117.36 -205.92 -77.76 -295.92
-123.12 -5.76 197.28 -2.16 -12.24 -25.92 -61.92 33.84
-5.76 28.08 -2.16 298.08 -25.92 -163.44 33.84 -162.72
-191.52 -117.36 -12.24 -25.92 240.48 120.24 -36.72 23.04
-117.36 -205.92 -25.92 -163.44 120.24 312.48 23.04 56.88
-87.12 -77.76 -61.92 33.84 -36.72 23.04 185.76 20.88
-77.76 -295.92 33.84 -162.72 23.04 56.88 20.88 401.76]
# f_expected = [142.272, 214.272, -13.824, 58.176, -75.456, 19.584, -52.992, -4.032]
f_expected = [142.272, 214.272, -13.824, 58.176, -75.456, -124.416, -52.992, -148.032]
@test isapprox(K, K_expected)
@test isapprox(f, f_expected)
end
f_expected = [142.272, 214.272, -13.824, 58.176, -75.456, -124.416, -52.992, -148.032]
@test isapprox(K, K_expected)
@test isapprox(f, f_expected)
@@ -1,85 +1,90 @@
# 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
using JuliaFEM, Test
@testset "test continuum 3d linear elasticity with surface load" begin
nodes = Dict{Int64, Node}(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [1.0, 1.0, 0.0],
4 => [0.0, 1.0, 0.0],
5 => [0.0, 0.0, 1.0],
6 => [1.0, 0.0, 1.0],
7 => [1.0, 1.0, 1.0],
8 => [0.0, 1.0, 1.0])
# test continuum 3d linear elasticity with surface load
element1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8])
element2 = Element(Quad4, [5, 6, 7, 8])
update!([element1, element2], "geometry", nodes)
update!([element1], "youngs modulus", 288.0)
update!([element1], "poissons ratio", 1/3)
update!([element2], "displacement traction force 3", 288.0)
update!([element1], "displacement load 3", 576.0)
X = Dict(1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [1.0, 1.0, 0.0],
4 => [0.0, 1.0, 0.0],
5 => [0.0, 0.0, 1.0],
6 => [1.0, 0.0, 1.0],
7 => [1.0, 1.0, 1.0],
8 => [0.0, 1.0, 1.0])
elasticity_problem = Problem(Elasticity, "solve continuum block", 3)
elasticity_problem.properties.finite_strain = false
push!(elasticity_problem, element1)
push!(elasticity_problem, element2)
element1 = Element(Hex8, (1, 2, 3, 4, 5, 6, 7, 8))
element2 = Element(Quad4, (5, 6, 7, 8))
update!((element1, element2), "geometry", X)
update!(element1, "youngs modulus", 288.0)
update!(element1, "poissons ratio", 1/3)
update!(element2, "displacement traction force 3", 288.0)
update!(element1, "displacement load 3", 576.0)
symxy = Element(Quad4, [1, 2, 3, 4])
symxz = Element(Quad4, [1, 2, 6, 5])
symyz = Element(Quad4, [1, 4, 8, 5])
update!([symxy, symxz, symyz], "geometry", nodes)
symyz["displacement 1"] = 0.0
symxz["displacement 2"] = 0.0
symxy["displacement 3"] = 0.0
boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement")
push!(boundary_problem, symxy, symxz, symyz)
problem = Problem(Elasticity, "solve continuum block", 3)
problem.properties.finite_strain = false
add_elements!(problem, element1, element2)
solver = LinearSolver(elasticity_problem, boundary_problem)
solver()
symxy = Element(Quad4, (1, 2, 3, 4))
symxz = Element(Quad4, (1, 2, 6, 5))
symyz = Element(Quad4, (1, 4, 8, 5))
update!([symxy, symxz, symyz], "geometry", X)
update!(symyz, "displacement 1", 0.0)
update!(symxz, "displacement 2", 0.0)
update!(symxy, "displacement 3", 0.0)
disp = element1("displacement", [1.0, 1.0, 1.0], 0.0)
info("displacement at tip: $disp")
u_expected = 2.0 * [-1/3, -1/3, 1.0]
@test isapprox(disp, u_expected)
end
bc = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement")
add_elements!(bc, symxy, symxz, symyz)
analysis = Analysis(Linear)
add_problems!(analysis, problem, bc)
run!(analysis)
u = element1("displacement", [1.0, 1.0, 1.0], 0.0)
u_expected = 2.0 * [-1/3, -1/3, 1.0]
@debug("displacement at tip", u, u_expected)
@test isapprox(u, u_expected)
# run similar analysis for different meshes
function solve_rod_model_elasticity(eltype)
fn = @__DIR__() * "/testdata/rod_short.med"
mesh = aster_read_mesh(fn, eltype)
element_sets = join(keys(mesh.element_sets), ", ")
info("element sets: $element_sets")
@debug("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)
solver()
u_max = maximum(p1.assembly.u)
info("$eltype, u_max = $u_max")
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_elements, "youngs modulus", 96.0)
update!(p1_elements, "poissons ratio", 1/3)
update!(p2_elements, "displacement traction force 1", 96.0)
update!(p3_elements, "displacement 1", 0.0)
update!(p4_elements, "displacement 2", 0.0)
update!(p5_elements, "displacement 3", 0.0)
add_elements!(p1, p1_elements)
add_elements!(p2, p2_elements)
add_elements!(p3, p3_elements)
add_elements!(p4, p4_elements)
add_elements!(p5, p5_elements)
analysis = Analysis(Linear)
add_problems!(analysis, p1, p2, p3, p4, p5)
run!(analysis)
u = p1("displacement", 0.0)
u_max = maximum(maximum(ui for ui in values(u)))
@debug("analysis of block", eltype, 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
@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)
@@ -1,52 +1,50 @@
# 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.Testing
using JuliaFEM, Test
@testset "test continuum nonlinear elasticity with surface load" begin
# test continuum nonlinear elasticity with surface load
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [1.0, 1.0, 0.0],
4 => [0.0, 1.0, 0.0],
5 => [0.0, 0.0, 1.0],
6 => [1.0, 0.0, 1.0],
7 => [1.0, 1.0, 1.0],
8 => [0.0, 1.0, 1.0])
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [1.0, 1.0, 0.0],
4 => [0.0, 1.0, 0.0],
5 => [0.0, 0.0, 1.0],
6 => [1.0, 0.0, 1.0],
7 => [1.0, 1.0, 1.0],
8 => [0.0, 1.0, 1.0])
element1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8])
element2 = Element(Quad4, [5, 6, 7, 8])
update!([element1, element2], "geometry", X)
update!([element1], "youngs modulus", 900.0)
update!([element1], "poissons ratio", 0.25)
update!([element2], "displacement traction force", [0.0, 0.0, -100.0])
element1 = Element(Hex8, (1, 2, 3, 4, 5, 6, 7, 8))
element2 = Element(Quad4, (5, 6, 7, 8))
update!((element1, element2), "geometry", X)
update!(element1, "youngs modulus", 900.0)
update!(element1, "poissons ratio", 0.25)
update!(element2, "displacement traction force", [0.0, 0.0, -100.0])
elasticity_problem = Problem(Elasticity, "solve continuum block", 3)
elasticity_problem.properties.finite_strain = true
push!(elasticity_problem, element1)
push!(elasticity_problem, element2)
problem = Problem(Elasticity, "solve continuum block", 3)
problem.properties.finite_strain = true
add_elements!(problem, element1, element2)
symxy = Element(Quad4, [1, 2, 3, 4])
symxz = Element(Quad4, [1, 2, 6, 5])
symyz = Element(Quad4, [1, 4, 8, 5])
update!([symxy, symxz, symyz], "geometry", X)
symxy["displacement 3"] = 0.0
symxz["displacement 2"] = 0.0
symyz["displacement 1"] = 0.0
boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement")
push!(boundary_problem, symxy, symxz, symyz)
symxy = Element(Quad4, (1, 2, 3, 4))
symxz = Element(Quad4, (1, 2, 6, 5))
symyz = Element(Quad4, (1, 4, 8, 5))
update!((symxy, symxz, symyz), "geometry", X)
update!(symxy, "displacement 3", 0.0)
update!(symxz, "displacement 2", 0.0)
update!(symyz, "displacement 1", 0.0)
solver = NonlinearSolver("solve 3d block")
push!(solver, elasticity_problem)
push!(solver, boundary_problem)
solver()
bc = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement")
add_elements!(bc, symxy, symxz, symyz)
disp = element1("displacement", [1.0, 1.0, 1.0], 0.0)
info("displacement at tip: $disp")
# verified using Code Aster.
# 2015-12-12-continuum-elasticity/vim c3d_grot_gdep_traction_force.comm
@test isapprox(disp, [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]; rtol=1.0e-4)
end
analysis = Analysis(Nonlinear)
add_problems!(analysis, problem, bc)
run!(analysis)
u = element1("displacement", [1.0, 1.0, 1.0], 0.0)
u_expected = [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]
@debug("displacement at tip", u, u_expected)
# verified using Code Aster:
# 2015-12-12-continuum-elasticity/vim c3d_grot_gdep_traction_force.comm
@test isapprox(u, u_expected; rtol=1.0e-4)
+27 -44
View File
@@ -1,28 +1,25 @@
# 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
using JuliaFEM, Test
function calc_model(mesh_name; with_volume_load=false, debug_print=false)
function calc_model(mesh_name; with_volume_load=false)
meshfile = @__DIR__()*"/testdata/3d_block.med"
mesh = aster_read_mesh(meshfile, mesh_name)
block = Problem(Elasticity, "BLOCK", 3)
block.properties.finite_strain = false
block.properties.geometric_stiffness = false
block.elements = create_elements(mesh, "BLOCK")
update!(block, "youngs modulus", 288.0)
update!(block, "poissons ratio", 1/3)
with_volume_load && update!(block, "displacement load 3", 576.0)
traction = Problem(Elasticity, "traction force", 3)
traction.properties.finite_strain = false
traction.properties.geometric_stiffness = false
traction.elements = create_elements(mesh, "LOAD")
update!(traction, "displacement traction force 3", 288.0)
block_elements = create_elements(mesh, "BLOCK")
update!(block_elements, "youngs modulus", 288.0)
update!(block_elements, "poissons ratio", 1/3)
with_volume_load && update!(block_elements, "displacement load 3", 576.0)
add_elements!(block, block_elements)
traction_elements = create_elements(mesh, "LOAD")
update!(traction_elements, "displacement traction force 3", 288.0)
add_elements!(block, traction_elements)
bc = Problem(Dirichlet, "symmetry boundary condition", 3, "displacement")
symyz = create_elements(mesh, "SYMYZ")
@@ -31,40 +28,26 @@ function calc_model(mesh_name; with_volume_load=false, debug_print=false)
update!(symyz, "displacement 1", 0.0)
update!(symxz, "displacement 2", 0.0)
update!(symxy, "displacement 3", 0.0)
push!(bc, symyz, symxz, symxy)
add_elements!(bc, symyz, symxz, symxy)
solver = LinearSolver("Solver block problem")
push!(solver, block, traction, bc)
solver()
analysis = Analysis(Linear)
add_problems!(analysis, block, bc)
run!(analysis)
max_u = maximum(block.assembly.u)
nu = round(Int, length(block.assembly.u)/3)
u = reshape(block.assembly.u, 3, nu)
if debug_print
f = reshape(full(block.assembly.f), 3, nu)
dump(round(u', 5))
dump(round(f', 5))
info("max |u| = $max_u")
end
return block, u
u = block("displacement", 0.0)
max_u = maximum(maximum(ui for ui in values(u)))
return max_u
end
@testset "test 3d block HEX8" begin
block, u = calc_model("BLOCK_HEX8"; with_volume_load=true)
@test isapprox(maximum(u), 2.0)
end
# 3d block HEX8
max_u = calc_model("BLOCK_HEX8"; with_volume_load=true)
@test isapprox(max_u, 2.0)
@testset "test 3d block TET4" begin
# block, u = calc_model("BLOCK_TET4", :TE4, :TR3; with_volume_load=true)
# @test isapprox(maximum(u), 2.1329516539440205)
block, u = calc_model("BLOCK_TET4"; with_volume_load=false)
@test isapprox(maximum(u), 1.0)
end
@testset "test 3d block TET10" begin
block, u = calc_model("BLOCK_TET10"; with_volume_load=false)
# @test isapprox(maximum(u), 2.13656216413056)
@test isapprox(maximum(u), 1.0)
end
# 3d block TET4
max_u = calc_model("BLOCK_TET4"; with_volume_load=false)
@test isapprox(max_u, 1.0)
# 3d block TET10
max_u = calc_model("BLOCK_TET10"; with_volume_load=false)
@test isapprox(max_u, 1.0)
+2 -6
View File
@@ -1,10 +1,7 @@
# 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
using JuliaFEM, Test
#= TODO: Fix test.
@testset "test forwarddiff version + volume load." begin
@@ -43,7 +40,7 @@ using JuliaFEM.Testing
push!(solver, body, bc)
solver()
disp = element("displacement", [1.0, 1.0], 0.0)
info("displacement at tip: $disp")
@info("displacement at tip: $disp")
# verified using Code Aster, verification/2015-10-22-plane-stress/cplan_grot_gdep_volume_force.resu
@test isapprox(disp[2], -8.77303119819776)
end
@@ -85,4 +82,3 @@ end
@test isapprox(f1, f2)
end
=#
@@ -1,9 +1,7 @@
# 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
using JuliaFEM, Test, SparseArrays, LinearAlgebra, Statistics
function get_stress_tensor(element, ip, time)
haskey(element, "displacement") || return nothing
@@ -13,7 +11,7 @@ function get_stress_tensor(element, ip, time)
nu = element("poissons ratio", ip, time)
mu = E/(2.0*(1.0+nu))
la = E*nu/((1.0+nu)*(1.0-2.0*nu))
S = la*trace(eps)*I + 2.0*mu*eps
S = la*tr(eps)*I + 2.0*mu*eps
return S
end
@@ -49,15 +47,15 @@ function lsq_fit(elements, field, time)
volume += w
end
end
info("Mass matrix for least-squares fit is assembled. Total volume to fit: $volume")
@info("Mass matrix for least-squares fit is assembled. Total volume to fit: $volume")
A = sparse(A)
b = sparse(b)
A = 1/2*(A + A')
SparseArrays.droptol!(A, 1.0e-6)
SparseArrays.dropzeros!(A)
nz = get_nonzero_rows(A)
F = ldltfact(A[nz,nz])
F = ldlt(A[nz,nz])
x = zeros(size(b)...)
x[nz, :] = F \ b[nz, :]
@@ -90,10 +88,13 @@ http://mms2.ensmp.fr/emms_paris/plasticite3D/exercices/eSpherePress.pdf
function test_wedge_sphere(model, u_CA, S_CA)
mesh_file = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(mesh_file, model)
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)
body_elements = create_elements(mesh, "HOLLOWSPHERE8")
update!(body_elements, "youngs modulus", 24580.0)
update!(body_elements, "poissons ratio", 1/3)
add_elements!(body, body_elements)
bc = Problem(Dirichlet, "symmetry bc", 3, "displacement")
el1 = create_elements(mesh, "FACE1")
update!(el1, "displacement 3", 0.0)
@@ -101,40 +102,47 @@ function test_wedge_sphere(model, u_CA, S_CA)
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()
add_elements!(bc, el1, el2, el3)
X = lo("geometry", 0.0)
u = lo("displacement", 0.0)
load = Problem(Elasticity, "pressure load", 3)
load_elements = create_elements(mesh, "OUTER")
update!(load_elements, "surface pressure", 7317.0)
add_elements!(load, load_elements)
analysis = Analysis(Linear)
add_problems!(analysis, body, load, bc)
run!(analysis)
X = load("geometry", 0.0)
u = load("displacement", 0.0)
nids = sort(collect(keys(X)))
umag = Float64[norm(u[id]) for id in nids]
um = mean(umag)
us = std(umag)
rtol = norm(um - 0.9) / max(norm(um), 0.9) * 100.0
info("mean umag = $um, std umag = $us, rtol = $rtol")
@debug("Displacement field statistics", mean=um, std=us, rtol=rtol)
@test rtol < 1.5 # percents
info("Verifying displacement against Code Aster solution.. ")
@debug("Verifying displacement against Code Aster solution.. ")
pass = true
for nid in keys(u_CA)
rtol = norm(u[nid] - u_CA[nid]) / max(norm(u[nid]), norm(u_CA[nid])) * 100.0
info("Node id $nid, rel diff to CA = $rtol %")
@test isapprox(u[nid], u_CA[nid])
@debug("Node id $nid, rel diff to CA = $rtol %")
pass &= isapprox(u[nid], u_CA[nid])
end
@test pass
S = lsq_fit(body.elements, get_stress, 0.0)
info("Verifying stress against Code Aster solution.. ")
@debug("Verifying stress against Code Aster solution.. ")
for nid in keys(S_CA)
rtol = norm(S[nid] - S_CA[nid]) / max(norm(S[nid]), norm(S_CA[nid])) * 100.0
info("Node id $nid, S=$(S[nid]), S_CA=$(S_CA[nid]), rel diff to CA = $rtol %")
#@test isapprox(S[nid], S_CA[nid])
@debug("Node id $nid, S=$(S[nid]), S_CA=$(S_CA[nid]), rel diff to CA = $rtol %")
# @test isapprox(S[nid], S_CA[nid])
# http://code-aster.org/doc/default/en/man_r/r3/r3.06.03.pdf
@test rtol < 10.0 # percents
pass &= (rtol < 10.0) # percents
end
@test pass
# Calculate principal stresses in nodes
Sp = Dict()
@@ -146,12 +154,12 @@ function test_wedge_sphere(model, u_CA, S_CA)
s[6] s[5] s[3]]
Sp[nid] = sort(eigvals(stress_tensor))
end
node_ids = sort(collect(keys(Sp)))
for (i, nid) in enumerate(node_ids)
println("$nid -> $(Sp[nid])")
@debug("$nid -> $(Sp[nid])")
if i > 9
println("...")
@debug("...")
break
end
end
@@ -170,20 +178,19 @@ function test_wedge_sphere(model, u_CA, S_CA)
for (i, ip) in enumerate(get_integration_points(element))
S = get_stress(element, ip, time)
rtol = norm(S - S_CA_gp[i]) / max(norm(S), norm(S_CA_gp[i]))
info("$i $S, rtol=$rtol")
@test isapprox(S, S_CA_gp[i])
@debug("$i $S, rtol=$rtol")
pass &= isapprox(S, S_CA_gp[i])
end
end
@test pass
end
@testset """1/8 hollow sphere with surface load""" begin
u_CA = Dict()
u_CA[38] = [-8.85861895037377E-01, -3.46944695195361E-18, -3.46944695195361E-18]
u_CA[49] = [-4.50934684240566E-01, -4.46908405333021E-01, -5.92329377847994E-01]
S_CA = Dict()
S_CA[38] = [-1.94504819940510E+03, -3.47014655438479E+04, -3.40876135857114E+04, 1.70379805227866E+03, 2.16707698388864E+03, 4.32009983342141E-12]
S_CA[49] = [-2.36067010168718E+04, -2.33820775692753E+04, -1.80606705820104E+04, 8.85783922092890E+03, 1.12122431934777E+04, 1.14035796957343E+04]
test_wedge_sphere("HOLLOWSPHERE8_WEDGE6", u_CA, S_CA)
#test_wedge_sphere("HOLLOWSPHERE8_WEDGE15")
end
# 1/8 hollow sphere with surface load
u_CA = Dict()
u_CA[38] = [-8.85861895037377E-01, -3.46944695195361E-18, -3.46944695195361E-18]
u_CA[49] = [-4.50934684240566E-01, -4.46908405333021E-01, -5.92329377847994E-01]
S_CA = Dict()
S_CA[38] = [-1.94504819940510E+03, -3.47014655438479E+04, -3.40876135857114E+04, 1.70379805227866E+03, 2.16707698388864E+03, 4.32009983342141E-12]
S_CA[49] = [-2.36067010168718E+04, -2.33820775692753E+04, -1.80606705820104E+04, 8.85783922092890E+03, 1.12122431934777E+04, 1.14035796957343E+04]
test_wedge_sphere("HOLLOWSPHERE8_WEDGE6", u_CA, S_CA)
#test_wedge_sphere("HOLLOWSPHERE8_WEDGE15")
+35 -37
View File
@@ -1,46 +1,44 @@
# 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
using JuliaFEM, Test
@testset "test Pyr5 elasticity with point load" begin
fn = dirname(@__DIR__) * "/geometry/3d_pyr/Pyr5.med"
mesh = aster_read_mesh(fn)
element_sets = join(keys(mesh.element_sets), ", ")
info("element sets: $element_sets")
element1 = create_elements(mesh,"Pyr5")
baseQuad = create_elements(mesh,"baseQuad")
tipPoint = Element(Poi1, collect(mesh.node_sets[:tipPoint]))
update!([element1,baseQuad,tipPoint], "geometry", mesh.nodes)
update!([element1], "youngs modulus", 288.0)
update!([element1], "poissons ratio", 1/3)
# test Pyr5 elasticity with point load
update!([tipPoint], "displacement traction force 1", 5.0)
update!([tipPoint], "displacement traction force 2", -7.0)
update!([tipPoint], "displacement traction force 3", 3.0)
fn = dirname(@__DIR__) * "/geometry/3d_pyr/Pyr5.med"
mesh = aster_read_mesh(fn)
element_sets = join(keys(mesh.element_sets), ", ")
@debug("element sets: $element_sets")
elasticity_problem = Problem(Elasticity, "solve continuum block", 3)
elasticity_problem.properties.finite_strain = false
push!(elasticity_problem, element1)
push!(elasticity_problem, tipPoint)
element1 = create_elements(mesh,"Pyr5")
baseQuad = create_elements(mesh,"baseQuad")
update!(element1, "youngs modulus", 288.0)
update!(element1, "poissons ratio", 1/3)
baseQuad[1]["displacement 1"] = 0.0
baseQuad[1]["displacement 2"] = 0.0
baseQuad[1]["displacement 3"] = 0.0
boundary_problem = Problem(Dirichlet, "Boundary conditions", 3, "displacement")
push!(boundary_problem, baseQuad)
tip_node_id = first(mesh.node_sets[:tipPoint])
tipPoint = Element(Poi1, [tip_node_id])
update!(tipPoint, "geometry", mesh.nodes)
update!(tipPoint, "displacement traction force 1", 5.0)
update!(tipPoint, "displacement traction force 2", -7.0)
update!(tipPoint, "displacement traction force 3", 3.0)
solver = LinearSolver(elasticity_problem, boundary_problem)
solver()
problem = Problem(Elasticity, "solve continuum block", 3)
problem.properties.finite_strain = false
add_elements!(problem, element1, tipPoint)
disp = element1[1]("displacement", [0.0, 0.0, 1.0], 0.0)
info("########################################################")
info("displacement at tip: $disp")
# Code_Aster Result in verification/2017-05-27-pyramids/Pyr5_displacement.txt
u_expected = [6.9444444444427100E-02,-9.7222222222197952E-02,1.0416666666679683E-02]
@test isapprox(disp, u_expected)
end
update!(first(baseQuad), "displacement 1", 0.0)
update!(first(baseQuad), "displacement 2", 0.0)
update!(first(baseQuad), "displacement 3", 0.0)
bc = Problem(Dirichlet, "Boundary conditions", 3, "displacement")
add_elements!(bc, baseQuad)
analysis = Analysis(Linear)
add_problems!(analysis, problem, bc)
run!(analysis)
xi, time = (0.0, 0.0, 1.0), 0.0
u = first(element1)("displacement", xi, time)
# Code_Aster Result in verification/2017-05-27-pyramids/Pyr5_displacement.txt
u_expected = [6.9444444444427100E-02,-9.7222222222197952E-02,1.0416666666679683E-02]
@debug("displacement at tip", u, u_expected)
@test isapprox(u, u_expected)
+27 -36
View File
@@ -1,45 +1,36 @@
# 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
using JuliaFEM, Test
@testset "test 2d linear elasticity with surface + volume load" begin
meshfile = "/geometry/2d_block/BLOCK_1elem.med"
mesh = aster_read_mesh(dirname(@__DIR__)*meshfile)
meshfile = "/geometry/2d_block/BLOCK_1elem.med"
mesh = aster_read_mesh(dirname(@__DIR__)*meshfile)
# field problem
block = Problem(Elasticity, "BLOCK", 2)
block.properties.formulation = :plane_strain
block.properties.finite_strain = false
block.properties.geometric_stiffness = false
block.elements = create_elements(mesh, "BLOCK")
update!(block.elements, "youngs modulus", 288.0)
update!(block.elements, "poissons ratio", 1/3)
# field problem
block = Problem(Elasticity, "BLOCK", 2)
block.properties.formulation = :plane_strain
block.properties.finite_strain = false
block.properties.geometric_stiffness = false
# traction
traction = Problem(Elasticity, "TRACTION", 2)
traction.properties.formulation = :plane_strain
traction.properties.finite_strain = false
traction.properties.geometric_stiffness = false
traction.elements = create_elements(mesh, "TOP")
update!(traction, "displacement traction force 2", 288.0*9/8)
block_elements = create_elements(mesh, "BLOCK")
update!(block_elements, "youngs modulus", 288.0)
update!(block_elements, "poissons ratio", 1/3)
traction_elements = create_elements(mesh, "TOP")
update!(traction_elements, "displacement traction force 2", 288.0*9/8)
add_elements!(block, block_elements, traction_elements)
# boundary conditions
bc_sym_23 = Problem(Dirichlet, "symmetry bc 23", 2, "displacement")
bc_sym_23.elements = create_elements(mesh, "LEFT")
update!(bc_sym_23, "displacement 1", 0.0)
bc_sym_13 = Problem(Dirichlet, "symmetry bc 13", 2, "displacement")
bc_sym_13.elements = create_elements(mesh, "BOTTOM")
update!(bc_sym_13, "displacement 2", 0.0)
# boundary conditions
bc = Problem(Dirichlet, "symmetry bc 23", 2, "displacement")
bc_sym_23_elements = create_elements(mesh, "LEFT")
bc_sym_13_elements = create_elements(mesh, "BOTTOM")
update!(bc_sym_23_elements, "displacement 1", 0.0)
update!(bc_sym_13_elements, "displacement 2", 0.0)
add_elements!(bc, bc_sym_23_elements, bc_sym_13_elements)
solver = LinearSolver(block, traction, bc_sym_23, bc_sym_13)
solver()
info("u = ", block.assembly.u)
info("λ = ", block.assembly.la)
end
analysis = Analysis(Linear, block, bc)
run!(analysis)
u = block("displacement", 0.0)
u3_expected = [-0.5, 1.0]
@debug("displacement", u, u3_expected)
@test isapprox(u[3], u3_expected)
+35 -36
View File
@@ -1,47 +1,46 @@
# 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
using JuliaFEM, Test
@testset "test Pyr5 elasticity with point load" begin
nodes = Dict{Int64, Node}(
1 => [-1.0,-1.0,-1.0],
2 => [ 1.0,-1.0,-1.0],
3 => [ 1.0, 1.0,-1.0],
4 => [-1.0, 1.0,-1.0],
5 => [ 0.0, 0.0, 1.0])
# Pyr5 elasticity with point load
element1 = Element(Pyr5, [1, 2, 3, 4, 5])
baseQuad = Element(Quad4, [1, 2, 3, 4])
tipPoint = Element(Poi1, [5,])
nodes = Dict(
1 => [-1.0,-1.0,-1.0],
2 => [ 1.0,-1.0,-1.0],
3 => [ 1.0, 1.0,-1.0],
4 => [-1.0, 1.0,-1.0],
5 => [ 0.0, 0.0, 1.0])
update!([element1,baseQuad,tipPoint], "geometry", nodes)
update!([element1], "youngs modulus", 288.0)
update!([element1], "poissons ratio", 1/3)
element1 = Element(Pyr5, (1, 2, 3, 4, 5))
baseQuad = Element(Quad4, (1, 2, 3, 4))
tipPoint = Element(Poi1, (5,))
update!([tipPoint], "displacement traction force 1", 5.0)
update!([tipPoint], "displacement traction force 2", -7.0)
update!([tipPoint], "displacement traction force 3", 3.0)
update!((element1, baseQuad, tipPoint), "geometry", nodes)
update!(element1, "youngs modulus", 288.0)
update!(element1, "poissons ratio", 1/3)
elasticity_problem = Problem(Elasticity, "solve continuum block", 3)
elasticity_problem.properties.finite_strain = false
push!(elasticity_problem, element1, baseQuad, tipPoint)
update!(tipPoint, "displacement traction force 1", 5.0)
update!(tipPoint, "displacement traction force 2", -7.0)
update!(tipPoint, "displacement traction force 3", 3.0)
baseQuad["displacement 1"] = 0.0
baseQuad["displacement 2"] = 0.0
baseQuad["displacement 3"] = 0.0
boundary_problem = Problem(Dirichlet, "Boundary conditions", 3, "displacement")
push!(boundary_problem, baseQuad)
elasticity_problem = Problem(Elasticity, "solve continuum block", 3)
elasticity_problem.properties.finite_strain = false
add_elements!(elasticity_problem, element1, baseQuad, tipPoint)
solver = LinearSolver(elasticity_problem, boundary_problem)
solver()
update!(baseQuad, "displacement 1", 0.0)
update!(baseQuad, "displacement 2", 0.0)
update!(baseQuad, "displacement 3", 0.0)
boundary_problem = Problem(Dirichlet, "Boundary conditions", 3, "displacement")
add_elements!(boundary_problem, baseQuad)
disp = element1("displacement", [0.0, 0.0, 1.0], 0.0)
info("########################################################")
info("displacement at tip: $disp")
# Code_Aster Result in verification/2017-05-27-pyramids/Pyr5_displacement.txt
u_expected = [6.9444444444427100E-02,-9.7222222222197952E-02,1.0416666666679683E-02]
@test isapprox(disp, u_expected)
end
analysis = Analysis(Linear)
add_problems!(analysis, elasticity_problem, boundary_problem)
run!(analysis)
xi, time = (0.0, 0.0, 1.0), 0.0
u = element1("displacement", xi, time)
u_expected = [6.9444444444427100E-02,-9.7222222222197952E-02,1.0416666666679683E-02]
@debug("displacement at tip", xi, time, u, u_expected)
# Code_Aster Result in verification/2017-05-27-pyramids/Pyr5_displacement.txt
@test isapprox(u, u_expected)
+7 -9
View File
@@ -1,9 +1,7 @@
# 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: assemble_mass_matrix!, add_elements!
using Base.Test
using JuliaFEM, Test
X = Dict(
1 => [2.0, 3.0, 4.0],
@@ -17,16 +15,16 @@ X[8] = 1/2*(X[1] + X[4])
X[9] = 1/2*(X[2] + X[4])
X[10] = 1/2*(X[3] + X[4])
element = Element(Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
element = Element(Tet10, (1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
update!(element, "youngs modulus", 480.0)
update!(element, "poissons ratio", 1/3)
update!(element, "geometry", X)
update!(element, "density", 105.0)
body = Problem(Heat, "TET", 1)
add_elements!(body, [element])
assemble_mass_matrix!(body, 0.0)
M = full(body.assembly.M)
problem = Problem(Heat, "tet10", 1)
add_element!(problem, element)
time = 0.0
assemble_mass_matrix!(problem, time)
M_expected = [
6 1 1 1 -4 -6 -4 -4 -6 -6
@@ -40,4 +38,4 @@ M_expected = [
-6 -4 -6 -4 16 16 8 16 32 16
-6 -6 -4 -4 8 16 16 16 16 32]
@test isapprox(M, M_expected)
@test isapprox(problem.assembly.M, M_expected)
+33 -40
View File
@@ -1,44 +1,37 @@
# 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: add_elements!
using Base.Test
using JuliaFEM, Test, LinearAlgebra
@testset "test tet10 stiffness matrix" begin
el = Element(Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
el["youngs modulus"] = 480.0
el["poissons ratio"] = 1/3
x1 = [2.0, 3.0, 4.0]
x2 = [6.0, 3.0, 2.0]
x3 = [2.0, 5.0, 1.0]
x4 = [4.0, 3.0, 6.0]
x5 = 0.5*(x1+x2)
x6 = 0.5*(x2+x3)
x7 = 0.5*(x3+x1)
x8 = 0.5*(x1+x4)
x9 = 0.5*(x2+x4)
x10 = 0.5*(x3+x4)
X = Dict{Int64, Vector{Float64}}(
1 => x1, 2 => x2, 3 => x3, 4 => x4, 5 => x5,
6 => x6, 7 => x7, 8 => x8, 9 => x9, 10 => x10)
u = Dict{Int64, Vector{Float64}}()
for i=1:10
u[i] = [0.0, 0.0, 0.0]
end
update!(el, "geometry", X)
update!(el, "displacement", u)
pr = Problem(Elasticity, "tet10", 3)
add_elements!(pr, [el])
assemble!(pr)
ass = pr.assembly
Kt = full(ass.K)
eigs = real(eigvals(Kt))
eigs_expected = [8809.45, 4936.01, 2880.56, 2491.66, 2004.85,
1632.49, 1264.32, 1212.42, 817.905,
745.755, 651.034, 517.441, 255.1, 210.955,
195.832, 104.008, 72.7562, 64.4376, 53.8515,
23.8417, 16.6354, 9.54682, 6.93361, 2.22099,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
@test isapprox(eigs, eigs_expected; atol=1.0e-2)
end
x1 = [2.0, 3.0, 4.0]
x2 = [6.0, 3.0, 2.0]
x3 = [2.0, 5.0, 1.0]
x4 = [4.0, 3.0, 6.0]
x5 = 0.5*(x1+x2)
x6 = 0.5*(x2+x3)
x7 = 0.5*(x3+x1)
x8 = 0.5*(x1+x4)
x9 = 0.5*(x2+x4)
x10 = 0.5*(x3+x4)
X = Dict(
1 => x1, 2 => x2, 3 => x3, 4 => x4, 5 => x5,
6 => x6, 7 => x7, 8 => x8, 9 => x9, 10 => x10)
u = Dict(i => zeros(3) for i in 1:10)
element = Element(Tet10, (1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
update!(element, "youngs modulus", 480.0)
update!(element, "poissons ratio", 1/3)
update!(element, "geometry", X)
update!(element, "displacement", u)
problem = Problem(Elasticity, "tet10", 3)
add_element!(problem, element)
time = 0.0
assemble!(problem, time)
eigs = real(eigvals(Matrix(problem.assembly.K)))
eigs_expected = [8809.45, 4936.01, 2880.56, 2491.66, 2004.85,
1632.49, 1264.32, 1212.42, 817.905,
745.755, 651.034, 517.441, 255.1, 210.955,
195.832, 104.008, 72.7562, 64.4376, 53.8515,
23.8417, 16.6354, 9.54682, 6.93361, 2.22099,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
@test isapprox(eigs, eigs_expected; atol=1.0e-2)
+29 -41
View File
@@ -1,45 +1,33 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM: add_elements!
using Base.Test
using JuliaFEM, Test
@testset "test tet4 stiffness matrix" begin
el = Element(Tet4, [1, 2, 3, 4])
el["youngs modulus"] = 96.0
el["poissons ratio"] = 1/3
x1 = [2.0, 3.0, 4.0]
x2 = [6.0, 3.0, 2.0]
x3 = [2.0, 5.0, 1.0]
x4 = [4.0, 3.0, 6.0]
u1 = u2 = u3 = u4 = zeros(3)
el["geometry"] = Vector{Float64}[x1, x2, x3, x4]
u = Vector{Float64}[u1, u2, u3, u4]
pr = Problem(Elasticity, "tet4", 3)
add_elements!(pr, [el])
assemble!(pr)
as = pr.assembly
Kt = full(as.K)
Kt_expected = [
149.0 108.0 24.0 -1.0 6.0 12.0 -54.0 -48.0 0.0 -94.0 -66.0 -36.0
108.0 344.0 54.0 -24.0 104.0 42.0 -24.0 -216.0 -12.0 -60.0 -232.0 -84.0
24.0 54.0 113.0 0.0 30.0 35.0 0.0 -24.0 -54.0 -24.0 -60.0 -94.0
-1.0 -24.0 0.0 29.0 -18.0 -12.0 -18.0 24.0 0.0 -10.0 18.0 12.0
6.0 104.0 30.0 -18.0 44.0 18.0 12.0 -72.0 -12.0 0.0 -76.0 -36.0
12.0 42.0 35.0 -12.0 18.0 29.0 0.0 -24.0 -18.0 0.0 -36.0 -46.0
-54.0 -24.0 0.0 -18.0 12.0 0.0 36.0 0.0 0.0 36.0 12.0 0.0
-48.0 -216.0 -24.0 24.0 -72.0 -24.0 0.0 144.0 0.0 24.0 144.0 48.0
0.0 -12.0 -54.0 0.0 -12.0 -18.0 0.0 0.0 36.0 0.0 24.0 36.0
-94.0 -60.0 -24.0 -10.0 0.0 0.0 36.0 24.0 0.0 68.0 36.0 24.0
-66.0 -232.0 -60.0 18.0 -76.0 -36.0 12.0 144.0 24.0 36.0 164.0 72.0
-36.0 -84.0 -94.0 12.0 -36.0 -46.0 0.0 48.0 36.0 24.0 72.0 104.0]
if !isapprox(Kt, Kt_expected)
info("Test failed")
info("Kt_expected")
dump(Kt_expected)
info("Kt")
dump(Kt)
end
@test isapprox(Kt, Kt_expected)
end
# test tet4 stiffness matrix
X = Dict(1 => [2.0, 3.0, 4.0],
2 => [6.0, 3.0, 2.0],
3 => [2.0, 5.0, 1.0],
4 => [4.0, 3.0, 6.0])
element = Element(Tet4, (1, 2, 3, 4))
update!(element, "youngs modulus", 96.0)
update!(element, "poissons ratio", 1/3)
update!(element, "geometry", X)
problem = Problem(Elasticity, "tet4", 3)
add_element!(problem, element)
time = 0.0
assemble!(problem, time)
K_expected = [
149 108 24 -1 6 12 -54 -48 0 -94 -66 -36
108 344 54 -24 104 42 -24 -216 -12 -60 -232 -84
24 54 113 0 30 35 0 -24 -54 -24 -60 -94
-1 -24 0 29 -18 -12 -18 24 0 -10 18 12
6 104 30 -18 44 18 12 -72 -12 0 -76 -36
12 42 35 -12 18 29 0 -24 -18 0 -36 -46
-54 -24 0 -18 12 0 36 0 0 36 12 0
-48 -216 -24 24 -72 -24 0 144 0 24 144 48
0 -12 -54 0 -12 -18 0 0 36 0 24 36
-94 -60 -24 -10 0 0 36 24 0 68 36 24
-66 -232 -60 18 -76 -36 12 144 24 36 164 72
-36 -84 -94 12 -36 -46 0 48 36 24 72 104]
@test isapprox(problem.assembly.K, K_expected)
+38 -61
View File
@@ -1,68 +1,45 @@
# 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
using JuliaFEM.Preprocess
using JuliaFEM.Testing
using SparseArrays, Test
@testset "test tet4 + volume load" begin
X1 = [2.0, 3.0, 4.0]
X2 = [6.0, 3.0, 2.0]
X3 = [2.0, 5.0, 1.0]
X4 = [4.0, 3.0, 6.0]
e1 = Element(Tet4, [1, 2, 3, 4])
e1["geometry"] = Node[X1, X2, X3, X4]
e1["youngs modulus"] = 96.0
e1["poissons ratio"] = 1/3
e1["displacement load 3"] = 784.0/110.0
e2 = Element(Tri3, [1, 2, 3])
e2["geometry"] = Node[X2, X1, X3]
e2["displacement 1"] = 0.0
e2["displacement 2"] = 0.0
e2["displacement 3"] = 0.0
p1 = Problem(Elasticity, "tetra", 3)
p2 = Problem(Dirichlet, "bc", 3, "displacement")
push!(p1, e1)
push!(p2, e2)
s = Solver(Linear)
push!(s, p1, p2)
s()
u_4 = p1.assembly.u[10:end]
u_expected = [-3.0/220.0, -9.0/220.0, 1.0/10.0]
@test isapprox(u_4, u_expected)
end
# test tet4 + volume load
@testset "test tet4 + surface load" begin
nodes = Dict{Int, Vector{Float64}}(
1 => [2.0, 3.0, 4.0],
2 => [6.0, 3.0, 2.0],
3 => [2.0, 5.0, 1.0],
4 => [4.0, 3.0, 6.0])
e1 = Element(Tet4, [1, 2, 3, 4])
update!(e1, "geometry", nodes)
e1["youngs modulus"] = 96.0
e1["poissons ratio"] = 1/3
e2 = Element(Tri3, [1, 2, 3])
update!(e2, "geometry", nodes)
e2["displacement 1"] = 0.0
e2["displacement 2"] = 0.0
e2["displacement 3"] = 0.0
e3 = Element(Tri3, [4, 3, 2])
update!(e3, "geometry", nodes)
update!(e3, "surface pressure", -96.0)
p1 = Problem(Elasticity, "tetra", 3)
p2 = Problem(Dirichlet, "bc", 3, "displacement")
push!(p1, e1, e3)
push!(p2, e2)
s = Solver(Linear)
push!(s, p1, p2)
s()
u_4 = p1.assembly.u[10:end]
u_expected = [-17.0/14.0, -27.0/14.0, 1.0]
info("u_4 = $(u_4)")
info("u_expected = $(u_expected)")
@test isapprox(u_4, u_expected)
end
X = Dict(1 => [2.0, 3.0, 4.0], 2 => [6.0, 3.0, 2.0],
3 => [2.0, 5.0, 1.0], 4 => [4.0, 3.0, 6.0])
element1 = Element(Tet4, (1, 2, 3, 4))
element2 = Element(Tri3, (1, 2, 3))
update!((element1, element2), "geometry", X)
update!(element1, "youngs modulus", 96.0)
update!(element1, "poissons ratio", 1/3)
update!(element1, "displacement load 3", 784.0/110.0)
update!(element2, "displacement 1", 0.0)
update!(element2, "displacement 2", 0.0)
update!(element2, "displacement 3", 0.0)
problem1 = Problem(Elasticity, "tetra", 3)
problem2 = Problem(Dirichlet, "bc", 3, "displacement")
add_element!(problem1, element1)
add_element!(problem2, element2)
analysis = Analysis(Linear)
add_problems!(analysis, problem1, problem2)
run!(analysis)
u4 = problem1("displacement", 0.0)[4]
u4_expected = [-3.0/220.0, -9.0/220.0, 1.0/10.0]
@test isapprox(u4, u4_expected)
# tet4 + surface load
element3 = Element(Tri3, (4, 3, 2))
update!(element3, "geometry", X)
update!(element3, "surface pressure", -96.0)
update!(element1, "displacement load 3", 0.0)
add_element!(problem1, element3)
run!(analysis)
u4 = problem1("displacement", 0.0)[4]
u4_expected = [-17.0/14.0, -27.0/14.0, 1.0]
@debug("displacement at node 4", u4, u4_expected)
#= TODO: Fix test. Make linear perturbation solver.
@testset "test tet4 + buckling" begin
@@ -92,8 +69,8 @@ end
dump(full(Kg))
la = sort(eigs(Km, -Kg)[1])
la_expected = [1.0, 4.0]
info("la = $la")
info("la_expected = $(la_expected)")
@info("la = $la")
@info("la_expected = $(la_expected)")
@test isapprox(la, la_expected)
end
=#
@@ -1,9 +1,7 @@
# 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
using JuliaFEM, Test
#=
@testset "2d nonlinear elasticity: test nonhomogeneous boundary conditions and stress calculation" begin
@@ -55,7 +53,7 @@ using JuliaFEM.Testing
u3_expected = [-2.36237356855269E-01, 5.00000000000000E-01]
u3 = reshape(block.assembly.u, 2, 4)[:, 3]
info("u3 = $u3")
@info("u3 = $u3")
@test isapprox(u3, u3_expected, atol=1.0e-5)
end
@@ -1,9 +1,7 @@
# 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
using JuliaFEM, Test
#=
@@ -58,7 +56,7 @@ using JuliaFEM.Testing
solver()
disp = element("displacement", [1.0, 1.0, 1.0], 1.0)
info("displacement at tip: $disp")
@info("displacement at tip: $disp")
u_expected = 2.0 * [-1/3, -1/3, 1.0]
@test isapprox(disp, u_expected)
end
@@ -69,7 +67,7 @@ end
# fn = @__DIR__() * "/testdata/rod_short.med"
# mesh = aster_read_mesh(fn, eltype)
# element_sets = join(keys(mesh.element_sets), ", ")
# info("element sets: $element_sets")
# @info("element sets: $element_sets")
# p1 = Problem(Elasticity, "rod", 3)
# p2 = Problem(Elasticity, "trac", 3)
# p3 = Problem(Dirichlet, "fixed", 3, "displacement")
@@ -89,7 +87,7 @@ end
# solver = LinearSolver(p1, p2, p3, p4, p5)
# solver()
# u_max = maximum(p1.assembly.u)
# info("$eltype, u_max = $u_max")
# @info("$eltype, u_max = $u_max")
# return u_max
# end
# @testset "compare 3d rod to CA solution" begin
+40 -232
View File
@@ -1,112 +1,28 @@
# 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.Testing
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM, LinearAlgebra, Test
@testset "Tet10 + convection" begin
# For some reason Tet10 fails, maybe because of convection.
mesh_file = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(mesh_file, "TETRA_TET10_1")
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, "thermal conductivity", 50.0)
face.elements = create_elements(mesh, "FACE4")
update!(face, "external temperature", 20.0)
update!(face, "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)
solver()
Temp = prob.assembly.u
info("Solution: $Temp")
Temp_expected = [ # using code aster
1.45606533688540E+01
0.0
0.0
0.0
1.05228712963739E+01
0.0
9.44202309239159E+00
1.05228712963739E+01
0.0
0.0]
info("Expected: $Temp_expected")
rtol = norm(Temp-Temp_expected)/max(norm(Temp), norm(Temp_expected))
info("rtol = $rtol")
@test isapprox(Temp, Temp_expected; rtol=1.0e-6)
end
# compare simple 3d heat problem to code aster solution
@testset "2d heat problem (one element)" begin
X = Dict{Int, Vector{Float64}}(
1 => [0.0,0.0],
2 => [1.0,0.0],
3 => [1.0,1.0],
4 => [0.0,1.0])
# define volume element
el1 = Element(Quad4, [1, 2, 3, 4])
update!(el1, "geometry", X)
update!(el1, "thermal conductivity", 6.0)
update!(el1, "heat source", 12.0)
# define boundary element for flux
el2 = Element(Seg2, [1, 2])
update!(el2, "geometry", X)
# linear ramp from 0 -> 6 in time 0 -> 1
update!(el2, "heat flux", 0.0 => 0.0)
update!(el2, "heat flux", 1.0 => 6.0)
# define heat problem and push elements to problem
problem = Problem(PlaneHeat, "one element heat problem", 1)
push!(problem, el1, el2)
# Set constant source f=12 with k=6. Accurate solution is
# T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x)
# when boundary flux not active (at t=0)
assemble!(problem, 0.0)
A = full(problem.assembly.K)
b = full(problem.assembly.f)
A_expected = [
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]
free_dofs = [1, 2]
@test isapprox(A, A_expected)
@test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [1.0, 1.0])
# Set constant flux g=6 on boundary. Accurate solution is
# u(x,y) = x which equals T=1 on boundary.
# at time t=1.0 all loads should be on.
empty!(problem.assembly)
assemble!(problem, 1.0)
A = full(problem.assembly.K)
b = full(problem.assembly.f)
@test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [2.0, 2.0])
end
@testset "compare simple 3d heat problem to code aster solution" begin
function calc_3d_heat_model(mesh_name)
fn = @__DIR__() * "/testdata/rod_short.med"
mesh = aster_read_mesh(fn, "Hex8")
mesh = aster_read_mesh(fn, mesh_name)
element_sets = join(keys(mesh.element_sets), ", ")
info("element sets: $element_sets")
p1 = Problem(Heat, "rod", 1)
rod = create_elements(mesh, "ROD")
@debug("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
rod = Problem(Heat, "rod", 1)
rod_elements = 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, "thermal conductivity", 50.0)
update!(rod_elements, "thermal conductivity", 50.0)
update!(face2, "external temperature", 20.0)
update!(face2, "heat transfer coefficient", 60.0)
update!(face3, "external temperature", 30.0)
@@ -117,141 +33,33 @@ end
update!(face5, "heat transfer coefficient", 30.0)
update!(face6, "external temperature", 60.0)
update!(face6, "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, "FACE1"))
update!(p2, "temperature 1", 100.0)
solver = LinearSolver(p1, p2)
solver()
# 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])
Temp = p1("temperature", 0.0)
for j in sort(collect(keys(Temp)))
T1 = Temp[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
push!(rod, rod_elements, face2, face3, face4, face5, face6)
bc = Problem(Dirichlet, "left support T=100", 1, "temperature")
bc_elements = create_elements(mesh, "FACE1")
update!(bc_elements, "temperature 1", 100.0)
add_elements!(bc, bc_elements)
analysis = Analysis(Linear)
add_problems!(analysis, rod, bc)
run!(analysis)
time = 0.0
temperature = rod("temperature", time)
minimum_temperature = minimum(values(temperature))
return minimum_temperature
end
@testset "compare simple 3d heat problem to analytical solution" begin
function calc_3d_heat_model(mesh_name)
fn = @__DIR__() * "/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, "thermal conductivity", 100.0)
update!(p1, "external temperature", 0.0)
update!(p1, "heat transfer coefficient", 1000.0)
update!(p2, "temperature 1", 100.0)
solver = LinearSolver(p1, p2)
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 = @__DIR__() * "/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, "thermal conductivity", 50.0)
update!(face2, "external temperature", 20.0)
update!(face2, "heat transfer coefficient", 60.0)
update!(face3, "external temperature", 30.0)
update!(face3, "heat transfer coefficient", 50.0)
update!(face4, "external temperature", 40.0)
update!(face4, "heat transfer coefficient", 40.0)
update!(face5, "external temperature", 50.0)
update!(face5, "heat transfer coefficient", 30.0)
update!(face6, "external temperature", 60.0)
update!(face6, "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)
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
Temp = calc_3d_heat_model(model)
T_min = minimum(Temp)
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
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_min = calc_3d_heat_model(model)
T_ca = CA_sol[model]
rtol = norm(T_min-T_ca)/max(T_min,T_ca)*100.0
@debug("Results for model $model", model, T_min, T_ca, rtol)
@test rtol < 1.0e-9
end
+21 -20
View File
@@ -1,27 +1,29 @@
# 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
using JuliaFEM, LinearAlgebra, Statistics, Test
mesh = aster_read_mesh(@__DIR__()*"/testdata/primitives.med", "CYLINDER_20_TET4")
problem = Problem(Heat, "rod of length 20", 1)
problem.elements = create_elements(mesh, "CYLINDER")
update!(problem, "thermal conductivity", 200.0)
outer = Problem(Heat, "outer surface", 1)
outer.elements = create_elements(mesh, "FACE2", "FACE3")
update!(outer, "external temperature", 20.0)
update!(outer, "heat transfer coefficient", 1.0)
problem_elements = create_elements(mesh, "CYLINDER")
update!(problem_elements, "thermal conductivity", 200.0)
outer_elements = create_elements(mesh, "FACE2", "FACE3")
update!(outer_elements, "external temperature", 20.0)
update!(outer_elements, "heat transfer coefficient", 1.0)
#midline = Problem(Heat, "midline of rod", 1)
#midline.elements = create_elements(mesh, "INNER_LINE")
boundary_elements = create_elements(mesh, "FACE1")
update!(boundary_elements, "temperature 1", 100.0)
problem = Problem(Heat, "rod of length 20", 1)
add_elements!(problem, problem_elements)
outer = Problem(Heat, "outer surface", 1)
add_elements!(outer, outer_elements)
boundary = Problem(Dirichlet, "homogeneous dirichlet boundary", 1, "temperature")
boundary.elements = create_elements(mesh, "FACE1")
update!(boundary, "temperature 1", 100.0)
#solver = LinearSolver(problem, outer, boundary, midline)
solver = LinearSolver(problem, outer, boundary)
solver()
add_elements!(boundary, boundary_elements)
analysis = Analysis(Linear)
add_problems!(analysis, problem, outer, boundary)
run!(analysis)
# Analytical solution
L = 20
@@ -36,12 +38,11 @@ T0 = 100.0
C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T0-Tu, 0]
T_diff = []
for x in linspace(0, 20)
T_FEM = problem("temperature", [x, 0.0, 0.0])[1]
for x in range(0, stop=20)
T_FEM = problem("temperature", [x, 0.0, 0.0], 0.0)[1]
T_ACC = dot(C, [exp(β*x), exp(-β*x)]) + Tu
push!(T_diff, norm(T_FEM - T_ACC))
info("x = $x, T_FEM = $T_FEM, T_ACC = $T_ACC")
end
info("mean diff = ", mean(T_diff))
@debug("mean diff on heat problem", mean(T_diff))
@test mean(T_diff) < 1.2 # mean diff = 1.14
+56
View File
@@ -0,0 +1,56 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM, LinearAlgebra, Test
# 2d heat problem (one element)
X = Dict(
1 => [0.0,0.0],
2 => [1.0,0.0],
3 => [1.0,1.0],
4 => [0.0,1.0])
# define volume element
element1 = Element(Quad4, (1, 2, 3, 4))
update!(element1, "geometry", X)
update!(element1, "thermal conductivity", 6.0)
update!(element1, "heat source", 12.0)
# define boundary element for flux
element2 = Element(Seg2, (1, 2))
update!(element2, "geometry", X)
# linear ramp from 0 -> 6 in time 0 -> 1
update!(element2, "heat flux", 0.0 => 0.0)
update!(element2, "heat flux", 1.0 => 6.0)
# define heat problem and add elements to problem
problem = Problem(PlaneHeat, "one element heat problem", 1)
add_elements!(problem, element1, element2)
# Set constant source f=12 with k=6. Accurate solution is
# T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x)
# when boundary flux not active (at t=0)
time = 0.0
assemble!(problem, time)
A = Matrix(problem.assembly.K)
b = Vector(problem.assembly.f)
A_expected = [
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]
free_dofs = [1, 2]
@test isapprox(A, A_expected)
@test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [1.0, 1.0])
# Set constant flux g=6 on boundary. Accurate solution is
# u(x,y) = x which equals T=1 on boundary.
# at time t=1.0 all loads should be on.
empty!(problem.assembly)
time = 1.0
assemble!(problem, time)
A = Matrix(problem.assembly.K)
b = Vector(problem.assembly.f)
@test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [2.0, 2.0])
+39 -44
View File
@@ -1,52 +1,47 @@
# 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
using JuliaFEM, Test
@testset "2d poisson problem with known analytical solution" begin
# from FENiCS tutorial, u(x,y) = 1 + x² + 2y² on [0x1]×[0,1]
# and u₀(x,y) = 1 + x² + 2y², f(x,y) = -6
# 2d poisson problem with known analytical solution
# from FENiCS tutorial, u(x,y) = 1 + x² + 2y² on [0x1]×[0,1]
# and u₀(x,y) = 1 + x² + 2y², f(x,y) = -6
mesh_file = @__DIR__()*"/testdata/primitives.med"
mesh = aster_read_mesh(mesh_file, "UNITSQUARE_6X4")
mesh_file = @__DIR__()*"/testdata/primitives.med"
mesh = aster_read_mesh(mesh_file, "UNITSQUARE_6X4")
field = Problem(PlaneHeat, "unit square, 6x4 triangular mesh", 1)
field.elements = create_elements(mesh, "UNITSQUARE")
update!(field, "thermal conductivity", 1.0)
update!(field, "heat source", -6.0)
square = Problem(PlaneHeat, "unit square, 6x4 triangular mesh", 1)
square_elements = create_elements(mesh, "UNITSQUARE")
update!(square_elements, "thermal conductivity", 1.0)
update!(square_elements, "heat source", -6.0)
add_elements!(square, square_elements)
bc = Problem(Dirichlet, "u₀(x,y) = 1 + x² + 2y²", 1, "temperature")
#bc.properties.order = 2
#bc.properties.dual_basis = true
bc.properties.variational = false
bc.elements = create_elements(mesh, "FACE1", "FACE2", "FACE3", "FACE4")
function u0(element, ip, time)
x, y = element("geometry", ip, time)
return 1 + x^2 + 2*y^2
end
update!(bc, "temperature 1", u0)
solver = LinearSolver(field, bc)
solver()
T_fem = Float64[]
T_acc = Float64[]
for (nid, X) in field("geometry", 0.0)
push!(T_fem, field("temperature", X)[1])
push!(T_acc, 1.0 + X[1]^2 + 2*X[2]^2)
end
@test maximum(abs.(T_fem-T_acc)) < 1.0e-12
# gradient of field is
gradT(X) = [2*X[1] 4*X[2]]
X = [0.5, 0.5]
gradT1 = gradT(X)
gradT2 = field("temperature", X, 0.0, Val{:Grad})
info("gradT1 = $gradT1, gradT2 = $gradT2")
# [1.1666666666666625 1.5000000000000018] quite big difference ..?
@test isapprox(gradT1, gradT2; rtol=25.0e-2)
bc = Problem(Dirichlet, "u₀(x,y) = 1 + x² + 2y²", 1, "temperature")
bc_elements = create_elements(mesh, "FACE1", "FACE2", "FACE3", "FACE4")
function u0(element, ip, time)
x, y = element("geometry", ip, time)
return 1 + x^2 + 2*y^2
end
update!(bc_elements, "temperature 1", u0)
add_elements!(bc, bc_elements)
analysis = Analysis(Linear)
add_problems!(analysis, square, bc)
run!(analysis)
T_fem = Float64[]
T_acc = Float64[]
for (nid, Xi) in square("geometry", 0.0)
push!(T_fem, square("temperature", Xi, 0.0)[1])
push!(T_acc, 1.0 + Xi[1]^2 + 2*Xi[2]^2)
end
@test maximum(abs.(T_fem-T_acc)) < 1.0e-12
# gradient of field is
X = [0.5, 0.5]
gradT1 = [2*X[1] 4*X[2]]
gradT2 = square("temperature", X, 0.0, Val{:Grad})
@debug("gradT1 = $gradT1, gradT2 = $gradT2")
# [1.1666666666666625 1.5000000000000018] quite big difference ..?
@test isapprox(gradT1, gradT2; rtol=25.0e-2)
+81
View File
@@ -0,0 +1,81 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM, LinearAlgebra, Test
# compare simple 3d heat problem to code aster solution
fn = @__DIR__() * "/testdata/rod_short.med"
mesh = aster_read_mesh(fn, "Hex8")
element_sets = join(keys(mesh.element_sets), ", ")
@debug("element sets: $element_sets")
problem = Problem(Heat, "rod", 1)
volume = 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!(volume, "thermal conductivity", 50.0)
update!(face2, "external temperature", 20.0)
update!(face2, "heat transfer coefficient", 60.0)
update!(face3, "external temperature", 30.0)
update!(face3, "heat transfer coefficient", 50.0)
update!(face4, "external temperature", 40.0)
update!(face4, "heat transfer coefficient", 40.0)
update!(face5, "external temperature", 50.0)
update!(face5, "heat transfer coefficient", 30.0)
update!(face6, "external temperature", 60.0)
update!(face6, "heat transfer coefficient", 20.0)
add_elements!(problem, volume, face2, face3, face4, face5, face6)
bc = Problem(Dirichlet, "left support T=100", 1, "temperature")
bc_elements = create_elements(mesh, "FACE1")
update!(bc_elements, "temperature 1", 100.0)
add_elements!(bc, bc_elements)
analysis = Analysis(Linear)
add_problems!(analysis, problem, bc)
run!(analysis)
# fields extracted from Code Aster .resu file
temp_ca = Dict(
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(
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(
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])
time = 0.0
temp_jf = problem("temperature", time)
@debug("Temperature comparison between JuliaFEM and Code Aster for 3D model",
temp_ca, temp_jf)
T1 = [temp_ca[i] for i in 1:8]
T2 = [temp_jf[i] for i in 1:8]
@test isapprox(T1, T2; rtol=1.0e-9)
+36
View File
@@ -0,0 +1,36 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM, LinearAlgebra, Test
# compare simple 3d heat problem to analytical solution
function calc_3d_heat_model(mesh_name)
fn = @__DIR__() * "/testdata/rod_short.med"
mesh = aster_read_mesh(fn, mesh_name)
problem = Problem(Heat, "rod", 1)
bc = Problem(Dirichlet, "left support T=100", 1, "temperature")
problem_elements = create_elements(mesh, "ROD", "FACE2")
bc_elements = create_elements(mesh, "FACE1")
update!(problem_elements, "thermal conductivity", 100.0)
update!(problem_elements, "external temperature", 0.0)
update!(problem_elements, "heat transfer coefficient", 1000.0)
update!(bc_elements, "temperature 1", 100.0)
add_elements!(problem, problem_elements)
add_elements!(bc, bc_elements)
analysis = Analysis(Linear)
add_problems!(analysis, problem, bc)
run!(analysis)
time = 0.0
T = problem("temperature", time)
T_min = minimum(map(first, values(T)))
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
@debug("Temperature for model $model", model, Tmin, Tacc, rtol)
@test isapprox(Tmin, 100/3)
end
+37 -42
View File
@@ -1,51 +1,46 @@
# 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
using JuliaFEM, Test
using AsterReader: RMEDFile, aster_read_nodes, aster_read_data
#=
Two rings, RING1 = inner, RING2 = outer, RINGS combined mesh. Set T=1.0 for
inner ring and T=2.0 for outer ring, measure temperature from middle of ring.
Results are calculated using Code Aster for comparison.
=#
@testset "test 3d heat, two rings, and compare to CA solution" begin
meshfile = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(meshfile, "RINGS_UNION")
## Two rings
# RING1 = inner, RING2 = outer, RINGS combined mesh. Set T=1.0 for
# inner ring and T=2.0 for outer ring, measure temperature from middle of ring.
# Results are calculated using Code Aster for comparison.
rings = Problem(Heat, "RINGS", 1)
# rings.elements = create_elements(mesh; element_type=:Tet4)
rings.elements = create_elements(mesh, "RING1", "RING2")
update!(rings.elements, "thermal conductivity", 1.0)
bc_inner = Problem(Dirichlet, "INNER SURFACE", 1, "temperature")
bc_inner.elements = create_elements(mesh, "RING1_INNER")
bc_outer = Problem(Dirichlet, "OUTER SURFACE", 1, "temperature")
bc_outer.elements = create_elements(mesh, "RING2_OUTER")
update!(bc_inner, "temperature 1", 1.0)
update!(bc_outer, "temperature 1", 2.0)
info("# of elements in RING1_INNER = ", length(bc_inner.elements))
info("# of elements in RING2_OUTER = ", length(bc_outer.elements))
solver = LinearSolver(rings, bc_inner, bc_outer)
solver()
# test 3d heat, two rings, and compare to CA solution
meshfile = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(meshfile, "RINGS_UNION")
temp_jf = rings("temperature", 0.0)
rings = Problem(Heat, "RINGS", 1)
rings_elements_1 = create_elements(mesh, "RING1")
rings_elements_2 = create_elements(mesh, "RING2")
update!(rings_elements_1, "thermal conductivity", 1.0)
update!(rings_elements_2, "thermal conductivity", 1.0)
add_elements!(rings, rings_elements_1)
add_elements!(rings, rings_elements_2)
fn = @__DIR__() * "/testdata/rings.rmed"
results = RMEDFile(fn)
nodes = aster_read_nodes(results)
temp_ca = aster_read_data(results, "TEMP")
bc_inner = Problem(Dirichlet, "INNER SURFACE", 1, "temperature")
bc_inner_elements = create_elements(mesh, "RING1_INNER")
update!(bc_inner_elements, "temperature 1", 1.0)
add_elements!(bc_inner, bc_inner_elements)
passed = true
for j in sort(collect(keys(temp_jf)))
X = nodes[j]
T1 = temp_jf[j]
T2 = temp_ca[j]
rtol = norm(T1-T2) / max(T1,T2)
@printf "% 5i : %8.5f %8.5f %8.5f | %8.5f %8.5f | %8.5e\n" j X... T1 T2 rtol
passed &= rtol < 1.0e-12
end
@test passed
end
bc_outer = Problem(Dirichlet, "OUTER SURFACE", 1, "temperature")
bc_outer_elements = create_elements(mesh, "RING2_OUTER")
update!(bc_outer_elements, "temperature 1", 2.0)
add_elements!(bc_outer, bc_outer_elements)
analysis = Analysis(Linear)
add_problems!(analysis, rings, bc_inner, bc_outer)
run!(analysis)
temp_jf = rings("temperature", 0.0)
fn = @__DIR__() * "/testdata/rings.rmed"
results = RMEDFile(fn)
nodes = aster_read_nodes(results)
temp_ca = aster_read_data(results, "TEMP")
rtol = [norm(temp_jf[j]-temp_ca[j]) / max(temp_jf[j],temp_ca[j]) for j in keys(temp_jf)]
@test maximum(rtol) < 1.0e-12
+46 -66
View File
@@ -1,81 +1,61 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using Base.Test
using JuliaFEM, Test
@testset "two increments, linear solver" begin
X = Dict{Int, Vector{Float64}}(
1 => [0.0,0.0],
2 => [1.0,0.0],
3 => [1.0,1.0],
4 => [0.0,1.0])
element = Element(Quad4, [1, 2, 3, 4])
update!(element, "geometry", X)
update!(element, "thermal conductivity", 6.0)
update!(element, "heat source", 0.0 => 12.0)
update!(element, "heat source", 1.0 => 24.0)
problem = Problem(PlaneHeat, "one element heat problem", 1)
push!(problem, element)
boundary_element = Element(Seg2, [1, 2])
update!(boundary_element, "geometry", X)
update!(boundary_element, "temperature 1", 0.0)
bc = Problem(Dirichlet, "fixed", 1, "temperature")
push!(bc, boundary_element)
solver = Solver(Linear, problem, bc)
X = Dict(1 => [0.0,0.0], 2 => [1.0,0.0], 3 => [1.0,1.0], 4 => [0.0,1.0])
element = Element(Quad4, (1, 2, 3, 4))
update!(element, "geometry", X)
update!(element, "thermal conductivity", 6.0)
update!(element, "heat source", 0.0 => 12.0)
update!(element, "heat source", 1.0 => 24.0)
problem = Problem(PlaneHeat, "one element heat problem", 1)
add_element!(problem, element)
boundary_element = Element(Seg2, [1, 2])
update!(boundary_element, "geometry", X)
update!(boundary_element, "temperature 1", 0.0)
bc = Problem(Dirichlet, "fixed", 1, "temperature")
add_element!(bc, boundary_element)
analysis = Analysis(Linear)
add_problems!(analysis, problem, bc)
empty!(problem.assembly)
solve!(solver, 0.0)
@test isapprox(solver("temperature", 0.0)[3], 1.0)
# two increments, linear solver
empty!(problem.assembly)
solver()
@test isapprox(solver("temperature", 0.0)[3], 1.0)
run!(analysis)
@test isapprox(analysis("temperature", 0.0)[3], 1.0)
empty!(problem.assembly)
solve!(solver, 1.0)
@test isapprox(solver("temperature", 1.0)[3], 2.0)
empty!(problem.assembly)
run!(analysis)
@test isapprox(analysis("temperature", 0.0)[3], 1.0)
empty!(problem.assembly)
solver()
@test isapprox(solver("temperature", 1.0)[3], 2.0)
empty!(problem.assembly)
analysis.properties.time = 1.0
run!(analysis)
@test isapprox(analysis("temperature", 1.0)[3], 2.0)
end
empty!(problem.assembly)
run!(analysis)
@test isapprox(analysis("temperature", 1.0)[3], 2.0)
@testset "two increments, nonlinear solver" begin
X = Dict{Int, Vector{Float64}}(
1 => [0.0,0.0],
2 => [1.0,0.0],
3 => [1.0,1.0],
4 => [0.0,1.0])
element = Element(Quad4, [1, 2, 3, 4])
update!(element, "geometry", X)
update!(element, "thermal conductivity", 6.0)
update!(element, "heat source", 0.0 => 12.0)
update!(element, "heat source", 1.0 => 24.0)
problem = Problem(PlaneHeat, "one element heat problem", 1)
push!(problem, element)
boundary_element = Element(Seg2, [1, 2])
update!(boundary_element, "geometry", X)
update!(boundary_element, "temperature 1", 0.0)
bc = Problem(Dirichlet, "fixed", 1, "temperature")
push!(bc, boundary_element)
solver = Solver(Nonlinear, problem, bc)
# two increments, nonlinear solver
empty!(problem.assembly)
solve!(solver, 0.0)
@test isapprox(solver("temperature", 0.0)[3], 1.0)
delete!(element.fields, "temperature")
analysis = Analysis(Nonlinear)
add_problems!(analysis, problem, bc)
empty!(problem.assembly)
solver()
@test isapprox(solver("temperature", 0.0)[3], 1.0)
empty!(problem.assembly)
run!(analysis)
@test isapprox(analysis("temperature", 0.0)[3], 1.0)
empty!(problem.assembly)
solve!(solver, 1.0)
@test isapprox(solver("temperature", 1.0)[3], 2.0)
empty!(problem.assembly)
run!(analysis)
@test isapprox(analysis("temperature", 0.0)[3], 1.0)
empty!(problem.assembly)
solve!(solver, 1.0)
@test isapprox(solver("temperature", 1.0)[3], 2.0)
empty!(problem.assembly)
analysis.properties.time = 1.0
run!(analysis)
@test isapprox(analysis("temperature", 1.0)[3], 2.0)
end
empty!(problem.assembly)
run!(analysis)
@test isapprox(analysis("temperature", 1.0)[3], 2.0)
+46
View File
@@ -0,0 +1,46 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM, LinearAlgebra, Test
mesh_file = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(mesh_file, "TETRA_TET10_1")
problem = Problem(Heat, "tet", 1)
problem_elements = create_elements(mesh, "TET")
update!(problem_elements, "thermal conductivity", 50.0)
add_elements!(problem, problem_elements)
face = Problem(Heat, "face 4", 1)
face_elements = create_elements(mesh, "FACE4")
update!(face_elements, "external temperature", 20.0)
update!(face_elements, "heat transfer coefficient", 60.0)
add_elements!(face, face_elements)
fixed = Problem(Dirichlet, "fixed face 3", 1, "temperature")
fixed_elements = create_elements(mesh, "FACE2")
update!(fixed_elements, "temperature 1", 0.0)
@debug("number of elements in fixed set", length(fixed_elements))
add_elements!(fixed, fixed_elements)
analysis = Analysis(Linear)
add_problems!(analysis, problem, face, fixed)
run!(analysis)
Temp = problem("temperature", 0.0)
Temp_expected = Dict(
1 => 1.45606533688540E+01,
2 => 0.0,
3 => 0.0,
4 => 0.0,
5 => 1.05228712963739E+01,
6 => 0.0,
7 => 9.44202309239159E+00,
8 => 1.05228712963739E+01,
9 => 0.0,
10 => 0.0)
@debug("Temperature solution using JuliaFEM", Temp)
@debug("Temperature solution using Code Aster", Temp_expected)
Temp_diff = collect(Temp[i]-Temp_expected[i] for i in 1:10)
@debug("Difference between solutions", Temp_diff)
@test isapprox(Temp_diff, zeros(10); atol=1.0e-9)
+64 -75
View File
@@ -1,85 +1,74 @@
# 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.Testing
using JuliaFEM, Test
function get_model()
X = Dict{Int, Vector{Float64}}(
1 => [2.0, 3.0, 4.0],
2 => [6.0, 3.0, 2.0],
3 => [2.0, 5.0, 1.0],
4 => [4.0, 3.0, 6.0])
u = Dict{Int, Vector{Float64}}(
1 => [0.0, 0.0, 0.0],
2 => [0.0, 0.0, 0.0],
3 => [0.0, 0.0, 0.0],
4 => [0.25, 0.25, 0.25])
e1 = Element(Tet4, [1, 2, 3, 4])
e2 = Element(Tri3, [1, 2, 3])
update!([e1, e2], "geometry", X)
update!([e1, e2], "displacement", 0.0 => u)
update!(e1, "youngs modulus", 96.0)
update!(e1, "poissons ratio", 1.0/3.0)
update!(e1, "density", 420.0)
update!(e2, "displacement 1", 0.0)
update!(e2, "displacement 2", 0.0)
update!(e2, "displacement 3", 0.0)
p1 = Problem(Elasticity, "test problem", 3)
p1.properties.finite_strain = false
p1.properties.geometric_stiffness = false
p2 = Problem(Dirichlet, "boundary condition", 3, "displacement")
push!(p1, e1)
push!(p2, e2)
solver = Solver(Modal)
solver.properties.which = :LM
push!(solver, p1, p2)
return solver
end
X = Dict(
1 => [2.0, 3.0, 4.0],
2 => [6.0, 3.0, 2.0],
3 => [2.0, 5.0, 1.0],
4 => [4.0, 3.0, 6.0])
u = Dict(
1 => [0.0, 0.0, 0.0],
2 => [0.0, 0.0, 0.0],
3 => [0.0, 0.0, 0.0],
4 => [0.25, 0.25, 0.25])
element1 = Element(Tet4, (1, 2, 3, 4))
element2 = Element(Tri3, (1, 2, 3))
update!((element1, element2), "geometry", X)
update!((element1, element2), "displacement", 0.0 => u)
update!(element1, "youngs modulus", 96.0)
update!(element1, "poissons ratio", 1.0/3.0)
update!(element1, "density", 420.0)
update!(element2, "displacement 1", 0.0)
update!(element2, "displacement 2", 0.0)
update!(element2, "displacement 3", 0.0)
problem1 = Problem(Elasticity, "test problem", 3)
problem1.properties.finite_strain = false
problem1.properties.geometric_stiffness = false
problem2 = Problem(Dirichlet, "boundary condition", 3, "displacement")
add_elements!(problem1, element1)
add_elements!(problem2, element2)
analysis = Analysis(Modal)
analysis.properties.which = :LM
add_problems!(analysis, problem1, problem2)
@testset "test eigenvalues for single tet4 element" begin
solver = get_model()
solver()
@test isapprox(solver.properties.eigvals, [4/3, 1/3])
end
# test eigenvalues for single tet4 element
run!(analysis)
@test isapprox(analysis.properties.eigvals, [4/3, 1/3])
@testset "test eigenvalues for single tet4 element, with geometric stiffness" begin
solver = get_model()
problem = first(solver.problems)
# problem.properties.finite_strain = true
problem.properties.geometric_stiffness = true
solver.properties.geometric_stiffness = true
solver()
@test isapprox(solver.properties.eigvals, [5/3, 2/3])
end
# test eigenvalues for single tet4 element, with geometric stiffness
problem1.properties.geometric_stiffness = true
analysis.properties.geometric_stiffness = true
run!(analysis)
@test isapprox(analysis.properties.eigvals, [5/3, 2/3])
@testset "test poisson problem modal analysis without tie" begin
X = Dict{Int64, Vector{Float64}}(
1 => [0.0, 0.0],
2 => [1.0, 0.0],
3 => [1.0, 3.0],
4 => [0.0, 3.0],
5 => [0.0, 3.0],
6 => [1.0, 3.0],
7 => [1.0, 9.0],
8 => [0.0, 9.0])
el1 = Element(Quad4, [1, 2, 3, 4])
el2 = Element(Quad4, [4, 3, 7, 8])
el3 = Element(Seg2, [1, 2])
el4 = Element(Seg2, [7, 8])
update!([el1, el2, el3, el4], "geometry", X)
update!([el1, el2], "density", 6.0)
update!([el1, el2], "thermal conductivity", 36.0)
update!([el3, el4], "temperature 1", 0.0)
p1 = Problem(PlaneHeat, "combined body", 1)
p2 = Problem(Dirichlet, "fixed ends", 1, "temperature")
push!(p1, el1, el2)
push!(p2, el3, el4)
solver = Solver(Modal)
push!(solver, p1, p2)
solver()
@test isapprox(solver.properties.eigvals[1], 1.0)
end
# test poisson problem modal analysis without tie
X = Dict(
1 => [0.0, 0.0],
2 => [1.0, 0.0],
3 => [1.0, 3.0],
4 => [0.0, 3.0],
5 => [0.0, 3.0],
6 => [1.0, 3.0],
7 => [1.0, 9.0],
8 => [0.0, 9.0])
element1 = Element(Quad4, (1, 2, 3, 4))
element2 = Element(Quad4, (4, 3, 7, 8))
element3 = Element(Seg2, (1, 2))
element4 = Element(Seg2, (7, 8))
update!((element1, element2, element3, element4), "geometry", X)
update!((element1, element2), "density", 6.0)
update!((element1, element2), "thermal conductivity", 36.0)
update!((element3, element4), "temperature 1", 0.0)
problem1 = Problem(PlaneHeat, "combined body", 1)
problem2 = Problem(Dirichlet, "fixed ends", 1, "temperature")
add_elements!(problem1, element1, element2)
add_elements!(problem2, element3, element4)
analysis = Analysis(Modal)
add_problems!(analysis, problem1, problem2)
run!(analysis)
@test isapprox(first(analysis.properties.eigvals), 1.0)
#=
@testset "test poisson modal problem with mesh tie" begin
+41 -146
View File
@@ -1,10 +1,7 @@
# 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
using JuliaFEM, Test, LinearAlgebra
#= this has nothing to do here
@testset "calculate cross-sectional properties" begin
@@ -14,19 +11,19 @@ using JuliaFEM.Testing
fixed1 = Problem(Dirichlet, "left support", 3, "displacement")
fixed1.elements = create_elements(mesh, "FACE1")
A = calculate_area(fixed1)
info("cross-section area: $A")
@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")
@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)
@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")
@info("I rtol = $rtol")
@test isapprox(I, I_expected; rtol = 0.2)
end
=#
@@ -64,144 +61,42 @@ numéro fréquence (HZ) norme d'erreur
[1] De Silva, Clarence W. Vibration: fundamentals and practice. CRC press, 2006, p.355
=#
@testset "long rod natural frequencies" begin
mesh_file = @__DIR__() * "/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
E = 50475.5
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")
# long rod natural frequencies
mesh_file = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(mesh_file, "CYLINDER_20_TET10")
body = Problem(Elasticity, "rod", 3)
body_elements = create_elements(mesh, "CYLINDER")
#E = 50475.44814745859
E = 50475.5
rho = 1.0
update!(body_elements, "youngs modulus", E)
update!(body_elements, "poissons ratio", 0.3)
update!(body_elements, "density", rho)
add_elements!(body, body_elements)
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)
add_elements!(fixed1, fixed1_elements)
# analytical solution is
l = 20.0
r = 1.0
la = 4.730040744862704/l
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)
add_elements!(fixed2, fixed2_elements)
# 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.properties.nev = 5
solver()
freqs_jf = sqrt.(solver.properties.eigvals)/(2.0*pi)
# with Tet4 elements
#freqs_ca = [1.19789E+00, 1.20179E+00, 3.07391E+00, 3.08813E+00, 4.87370E+00]
# with Tet10 elements
freqs_ca = [9.65942E-01, 9.66160E-01, 2.52127E+00, 2.52187E+00, 3.48584E+00]
# looks that juliafem results are more close to 1.0, maybe different integration order
rtol1 = norm(freq_sa - freqs_jf[1])/max(freq_sa, freqs_jf[1])
rtol2 = norm(freq_a - freqs_jf[1])/max(freq_a, freqs_jf[1])
info("rtol 1 = $rtol1, rtol 2 = $rtol2")
passed = true
for (f1, f2) in zip(freqs_jf, freqs_ca)
rtol = norm(f1-f2) / max(f1,f2)
@printf "JF: %8.5e | CA: %8.5e | rtol: %8.5e\n" f1 f2 rtol
passed &= (rtol < 3.0e-2)
end
@test rtol2 < 3.5e-2
@test passed
#=
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
@testset "eigenvalues of cube (tet4)" begin
meshfile = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(meshfile, "CUBE_TET4")
cube = Problem(mesh, Elasticity, "CUBE", 3)
update!(cube.elements, "youngs modulus", 10000.0)
update!(cube.elements, "poissons ratio", 0.3)
update!(cube.elements, "density", 10.0)
sym23 = create_elements(mesh, "FACE231")
update!(sym23, "displacement 1", 0.0)
sym13 = create_elements(mesh, "FACE131")
update!(sym13, "displacement 2", 0.0)
sym12 = create_elements(mesh, "FACE121")
update!(sym12, "displacement 3", 0.0)
bcs = Problem(Dirichlet, "bcs", 3, "displacement")
bcs.elements = [sym23; sym13; sym12]
solver = Solver(Modal)
solver.properties.nev = 5
push!(solver, cube, bcs)
solver()
freqs_jf = sqrt.(solver.properties.eigvals)/(2.0*pi)
freqs_ca = [3.73724E+00, 3.73724E+00, 4.93519E+00, 6.59406E+00, 7.65105E+00]
for (f1, f2) in zip(freqs_jf, freqs_ca)
rtol = norm(f1-f2) / max(f1,f2)
@printf "JF: %8.5e | CA: %8.5e | rtol: %8.5e\n" f1 f2 rtol
@test rtol < 1.0e-5
end
end
analysis = Analysis(Modal)
add_problems!(analysis, body, fixed1, fixed2)
analysis.properties.nev = 5
run!(analysis)
freqs_jf = sqrt.(analysis.properties.eigvals)/(2.0*pi)
# with Tet4 elements
#freqs_ca = [1.19789E+00, 1.20179E+00, 3.07391E+00, 3.08813E+00, 4.87370E+00]
# with Tet10 elements
freqs_ca = [9.65942E-01, 9.66160E-01, 2.52127E+00, 2.52187E+00, 3.48584E+00]
rtol = [norm(f1-f2) / max(f1,f2) for (f1, f2) in zip(freqs_jf, freqs_ca)]
@test maximum(rtol) < 3.0e-2
+53 -62
View File
@@ -1,68 +1,59 @@
# 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
using JuliaFEM, Test
@testset "eigenvalues of CYLINDER1" begin
meshfile = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(meshfile, "CYLINDER_1_TET4")
cylinder = Problem(mesh, Elasticity, "CYLINDER", 3)
update!(cylinder.elements, "youngs modulus", 10000.0)
update!(cylinder.elements, "poissons ratio", 0.3)
update!(cylinder.elements, "density", 10.0)
bc1 = create_elements(mesh, "FACE_YZ1")
update!(bc1, "displacement 1", 0.0)
update!(bc1, "displacement 2", 0.0)
update!(bc1, "displacement 3", 0.0)
bcs = Problem(Dirichlet, "bcs", 3, "displacement")
bcs.elements = bc1
solver = Solver(Modal)
solver.properties.nev = 3
push!(solver, cylinder, bcs)
solver()
freqs_jf = sqrt.(solver.properties.eigvals)/(2.0*pi)
freqs_ca = [4.84532E+00, 4.90698E+00, 8.33813E+00]
passed = []
for (f1, f2) in zip(freqs_jf, freqs_ca)
rtol = norm(f1-f2) / max(f1,f2)
@printf "JF: %8.5e | CA: %8.5e | rtol: %8.5e\n" f1 f2 rtol
push!(passed, rtol < 1.0e-5)
end
@test reduce(&, passed)
end
# eigenvalues of CYLINDER1
@testset "eigenvalues of CYLINDER20" begin
meshfile = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(meshfile, "CYLINDER_20_TET4")
cylinder = Problem(mesh, Elasticity, "CYLINDER", 3)
#update!(cylinder.elements, "youngs modulus", 10.0e6)
update!(cylinder.elements, "youngs modulus", 50475.5)
update!(cylinder.elements, "poissons ratio", 0.3)
#update!(cylinder.elements, "density", 10.0)
update!(cylinder.elements, "density", 1.0)
bc1 = create_elements(mesh, "FACE1", "FACE2")
update!(bc1, "displacement 1", 0.0)
update!(bc1, "displacement 2", 0.0)
update!(bc1, "displacement 3", 0.0)
bcs = Problem(Dirichlet, "bcs", 3, "displacement")
bcs.elements = bc1
solver = Solver(Modal)
solver.properties.nev = 3
push!(solver, cylinder, bcs)
solver()
freqs_jf = sqrt.(solver.properties.eigvals)/(2.0*pi)
#freqs_ca = [8.82848E-01, 8.85353E-01, 5.30286E+00] # only face1 fixed
#freqs_ca = [5.33185E+00, 5.34920E+00, 1.36820E+01] # face1 and face2 fixed
freqs_ca = [1.19789E+00, 1.20179E+00, 3.07391E+00]
passed = []
for (f1, f2) in zip(freqs_jf, freqs_ca)
rtol = norm(f1-f2) / max(f1,f2)
@printf "JF: %8.5e | CA: %8.5e | rtol: %8.5e\n" f1 f2 rtol
push!(passed, rtol < 1.0e-5)
end
@test reduce(&, passed)
end
meshfile = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(meshfile, "CYLINDER_1_TET4")
cylinder = Problem(Elasticity, "CYLINDER", 3)
cylinder_elements = create_elements(mesh, "CYLINDER")
update!(cylinder_elements, "youngs modulus", 10000.0)
update!(cylinder_elements, "poissons ratio", 0.3)
update!(cylinder_elements, "density", 10.0)
add_elements!(cylinder, cylinder_elements)
bc = Problem(Dirichlet, "bc", 3, "displacement")
bc_elements = create_elements(mesh, "FACE_YZ1")
update!(bc_elements, "displacement 1", 0.0)
update!(bc_elements, "displacement 2", 0.0)
update!(bc_elements, "displacement 3", 0.0)
add_elements!(bc, bc_elements)
analysis = Analysis(Modal)
analysis.properties.nev = 3
add_problems!(analysis, cylinder, bc)
run!(analysis)
freqs_jf = sqrt.(analysis.properties.eigvals)/(2.0*pi)
freqs_ca = [4.84532E+00, 4.90698E+00, 8.33813E+00]
@test isapprox(freqs_jf, freqs_ca; rtol=1.0e-5)
# eigenvalues of CYLINDER20
meshfile = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(meshfile, "CYLINDER_20_TET4")
cylinder = Problem(Elasticity, "CYLINDER", 3)
cylinder_elements = create_elements(mesh, "CYLINDER")
#update!(cylinder_elements, "youngs modulus", 10.0e6)
update!(cylinder_elements, "youngs modulus", 50475.5)
update!(cylinder_elements, "poissons ratio", 0.3)
#update!(cylinder_elements, "density", 10.0)
update!(cylinder_elements, "density", 1.0)
add_elements!(cylinder, cylinder_elements)
bc = Problem(Dirichlet, "bc", 3, "displacement")
bc_elements = create_elements(mesh, "FACE1", "FACE2")
update!(bc_elements, "displacement 1", 0.0)
update!(bc_elements, "displacement 2", 0.0)
update!(bc_elements, "displacement 3", 0.0)
add_elements!(bc, bc_elements)
analysis = Analysis(Modal)
analysis.properties.nev = 3
add_problems!(analysis, cylinder, bc)
run!(analysis)
freqs_jf = sqrt.(analysis.properties.eigvals)/(2.0*pi)
freqs_ca = [1.19789E+00, 1.20179E+00, 3.07391E+00]
@test isapprox(freqs_jf, freqs_ca; rtol=1.0e-5)
+27 -31
View File
@@ -1,37 +1,33 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using Base.Test
using JuliaFEM
using JuliaFEM, Test
@testset "zero eigenmode model" begin
X = Dict(
1 => [2.0, 3.0, 4.0],
2 => [6.0, 3.0, 2.0],
3 => [2.0, 5.0, 1.0],
4 => [4.0, 3.0, 6.0])
X[5] = 1/2*(X[1] + X[2])
X[6] = 1/2*(X[2] + X[3])
X[7] = 1/2*(X[3] + X[1])
X[8] = 1/2*(X[1] + X[4])
X[9] = 1/2*(X[2] + X[4])
X[10] = 1/2*(X[3] + X[4])
# zero eigenmode model
X = Dict(
1 => [2.0, 3.0, 4.0],
2 => [6.0, 3.0, 2.0],
3 => [2.0, 5.0, 1.0],
4 => [4.0, 3.0, 6.0])
X[5] = 1/2*(X[1] + X[2])
X[6] = 1/2*(X[2] + X[3])
X[7] = 1/2*(X[3] + X[1])
X[8] = 1/2*(X[1] + X[4])
X[9] = 1/2*(X[2] + X[4])
X[10] = 1/2*(X[3] + X[4])
element = Element(Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
update!(element, "youngs modulus", 480.0)
update!(element, "poissons ratio", 1/3)
update!(element, "geometry", X)
update!(element, "density", 105.0)
body = Problem(Elasticity, "TET", 3)
add_elements!(body, [element])
solver = Solver(Modal, body)
solver.properties.nev = 30
assemble!(body, 0.0)
assemble!(body, 0.0, Val{:mass_matrix})
solver()
w2 = solver.properties.eigvals
w2_expected = [5.66054e-16, 1.08925e-15, 1.95035e-15, -2.18372e-15, -8.01069e-15, 9.26902e-15, 0.401407, 0.88248, 1.3185, 2.55833, 3.3538, 5.02371, 7.71933, 8.43733, 11.0521, 19.3262, 26.0944, 28.7033, 51.5814, 59.1422, 83.3597, 91.1507, 121.689, 126.753, 157.938, 160.344, 205.634, 324.608, 468.331]
@test isapprox(w2, w2; atol=1.0e-4)
end
element = Element(Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
update!(element, "youngs modulus", 480.0)
update!(element, "poissons ratio", 1/3)
update!(element, "geometry", X)
update!(element, "density", 105.0)
body = Problem(Elasticity, "TET", 3)
add_elements!(body, element)
analysis = Analysis(Modal)
add_problems!(analysis, body)
analysis.properties.nev = 30
run!(analysis)
w2 = analysis.properties.eigvals
w2_expected = [5.66054e-16, 1.08925e-15, 1.95035e-15, -2.18372e-15, -8.01069e-15, 9.26902e-15, 0.401407, 0.88248, 1.3185, 2.55833, 3.3538, 5.02371, 7.71933, 8.43733, 11.0521, 19.3262, 26.0944, 28.7033, 51.5814, 59.1422, 83.3597, 91.1507, 121.689, 126.753, 157.938, 160.344, 205.634, 324.608, 468.331]
@test isapprox(w2, w2; atol=1.0e-4)
+28 -28
View File
@@ -2,7 +2,7 @@
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Testing
using SparseArrays, Test
function test_auxiliary_plane_transforms()
nodes = Vector{Float64}[
@@ -18,19 +18,19 @@ function test_auxiliary_plane_transforms()
e1["normal-tangential coordinates"] = Matrix{Float64}[R, R, R]
time::Real = 0.0
x0, Q = create_auxiliary_plane(e1, time)
info("x0 = $x0")
info("Q = $Q")
@info("x0 = $x0")
@info("Q = $Q")
@test isapprox(x0, [1.0/3.0, 1.0/3.0, 0.0])
@test isapprox(Q, R)
p1 = Float64[1.0/3.0+0.1, 1.0/3.0+0.1, 1.0]
p2 = project_point_to_auxiliary_plane(p1, x0, Q)
info("point in auxiliary plane p2 = $p2")
@info("point in auxiliary plane p2 = $p2")
@test isapprox(p2, [0.1, 0.1])
theta = project_point_from_plane_to_surface(p2, x0, Q, e1, time)
info("theta = $theta")
@info("theta = $theta")
@test isapprox(theta[1], 0.0)
X = e1("geometry", theta[2:3], time)
info("projected point = $X")
@info("projected point = $X")
@test isapprox(X, Float64[1.0/3.0+0.1, 1.0/3.0+0.1, 0.0])
end
@@ -126,7 +126,7 @@ function test_calculate_polygon_centerpoint()
0.0 1.0 2.0 2.0 1.25 0.0
0.5 0.0 0.0 1.0 1.75 1.33333]
C = calculate_polygon_centerpoint(P)
info("Polygon centerpoint: $C")
@info("Polygon centerpoint: $C")
@test isapprox(C, [1.0397440690338993, 0.8047003412233396])
end
@@ -154,10 +154,10 @@ function test_assemble_3d_problem_tri3()
push!(prob, sel)
stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)
info("stiffness matrix for this problem:\n$stiffness_matrix")
@info("stiffness matrix for this problem:\n$stiffness_matrix")
M = D = 1/24*[2 1 1; 1 2 1; 1 1 2]
B = [D -M] # slave dofs are first in this.
info("expected matrix for this problem:\n$B")
@info("expected matrix for this problem:\n$B")
@test isapprox(stiffness_matrix, B)
# rotate and translate surface and check that we are still having same results
@@ -183,15 +183,15 @@ function test_assemble_3d_problem_tri3()
end
calculate_normal_tangential_coordinates!(sel, 0.0)
stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)
info("sel midpnt: ", sel("geometry", [1/3, 1/3], 0.0))
info("nt basis: ", sel("normal-tangential coordinates", [1/3, 1/3], 0.0))
@info("sel midpnt: ", sel("geometry", [1/3, 1/3], 0.0))
@info("nt basis: ", sel("normal-tangential coordinates", [1/3, 1/3], 0.0))
@test isapprox(stiffness_matrix, B)
end
function test_assemble_3d_problem_quad4()
info("assemble 3d problem in quad4-quad4")
@info("assemble 3d problem in quad4-quad4")
nodes = Vector{Float64}[
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
@@ -215,10 +215,10 @@ function test_assemble_3d_problem_quad4()
D = [16 8 4 8; 8 16 8 4; 4 8 16 8; 8 4 8 16]
M = [25 5 1 5; 20 10 2 4; 16 8 4 8; 20 4 2 10]
B = [D -M] # slave dofs are first in this.
info("expected matrix for this problem:")
@info("expected matrix for this problem:")
dump(round(B, 3))
info("stiffness matrix for this problem:")
@info("stiffness matrix for this problem:")
dump(round(stiffness_matrix, 3))
@test isapprox(stiffness_matrix, B)
@@ -226,7 +226,7 @@ end
function test_assemble_3d_problem_quad4_2()
info("assemble 3d problem in quad4-quad4")
@info("assemble 3d problem in quad4-quad4")
nodes = Vector{Float64}[
[0.0, 0.0, 0.0],
[1/4, 0.0, 0.0],
@@ -261,10 +261,10 @@ function test_assemble_3d_problem_quad4_2()
3456 1152 1152 3456
]
B = [D -M] # slave dofs are first in this.
info("expected matrix for this problem:")
@info("expected matrix for this problem:")
dump(round(B, 3))
info("stiffness matrix for this problem:")
@info("stiffness matrix for this problem:")
dump(round(stiffness_matrix, 3))
@test isapprox(stiffness_matrix, B)
@@ -272,7 +272,7 @@ end
function test_assemble_3d_problem_quad4_3()
info("assemble 3d problem in quad4-quad4")
@info("assemble 3d problem in quad4-quad4")
a = 1/4
b = 1/3
nodes = Vector{Float64}[
@@ -310,10 +310,10 @@ function test_assemble_3d_problem_quad4_3()
]
B = [D -M] # slave dofs are first in this.
info("expected matrix for this problem:")
@info("expected matrix for this problem:")
dump(round(B, 3))
info("stiffness matrix for this problem:")
@info("stiffness matrix for this problem:")
dump(round(stiffness_matrix, 3))
@test isapprox(stiffness_matrix, B)
@@ -375,7 +375,7 @@ function test_3d_problem()
solver()
X = el2("geometry", [1.0, 1.0, 1.0], 0.0)
u = el2("displacement", [1.0, 1.0, 1.0], 0.0)
info("displacement at $X = $u")
@info("displacement at $X = $u")
@test isapprox(u, 1/36*[1, 1, -4])
end
@@ -422,7 +422,7 @@ end
sel1["master elements"] = [mel1, mel2, mel4, mel5]
sel2["master elements"] = [mel2, mel3, mel5, mel6]
stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)*2592*6
info("interface matrix:")
@info("interface matrix:")
dump(round(stiffness_matrix, 3))
B = [
864 432 0 432 216 0 -420 -375 -15 0 -504 -450 -18 0 -84 -75 -3 0
@@ -432,7 +432,7 @@ end
216 864 216 432 1728 432 -24 -138 -138 -24 -144 -828 -828 -144 -120 -690 -690 -120
0 216 432 0 432 864 0 -3 -75 -84 0 -18 -450 -504 0 -15 -375 -420
]
info("expected interface matrix:")
@info("expected interface matrix:")
dump(round(B, 3))
@test isapprox(stiffness_matrix, B)
end
@@ -504,9 +504,9 @@ end
B = full(B)
D = B[1:9,1:9]
M = B[1:9,10:end]
info("interface matrix D:")
@info("interface matrix D:")
dump(round(D, 3))
info("interface matrix M:")
@info("interface matrix M:")
dump(round(M, 3))
D_expected = [
1296 648 0 648 324 0 0 0 0
@@ -531,9 +531,9 @@ end
0 0 0 0 0 -1 -25 -28 0 -25 -625 -700 0 -28 -700 -784]
info("D - D_expected")
@info("D - D_expected")
dump(D - D_expected)
info("M - M_expected")
@info("M - M_expected")
dump(M - M_expected)
@test isapprox(D, D_expected)
@@ -547,7 +547,7 @@ end
216 864 216 432 1728 432 -24 -138 -138 -24 -144 -828 -828 -144 -120 -690 -690 -120
0 216 432 0 432 864 0 -3 -75 -84 0 -18 -450 -504 0 -15 -375 -420
]
info("expected interface matrix:")
@info("expected interface matrix:")
dump(round(B, 3))
@test isapprox(stiffness_matrix, B)
=#
+2 -2
View File
@@ -2,7 +2,7 @@
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Testing
using Test
function get_test_model()
@@ -88,7 +88,7 @@ end
solver()
X = e3("geometry", [1.0, 1.0], 0.0)
u = e3("displacement", [1.0, 1.0], 0.0)
info("displacement at $X: $u")
@info("displacement at $X: $u")
u_expected = [-1/3, 1.0]
@test isapprox(u, u_expected)
end
+2 -2
View File
@@ -2,7 +2,7 @@
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Testing
using Test
function get_test_2d_model()
X = Dict{Int64, Vector{Float64}}(
@@ -107,6 +107,6 @@ end
u = e2("displacement", [1.0, 1.0], 0.0)
u_expected = [-1/3, 1.0]
info("displacement at tip: $u")
@info("displacement at tip: $u")
@test isapprox(u, u_expected)
end
+21 -20
View File
@@ -1,9 +1,7 @@
# 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
using JuliaFEM, Test, LinearAlgebra, Statistics
mesh = Mesh()
add_node!(mesh, 1, [0.0, 0.0])
@@ -30,27 +28,31 @@ add_element_to_element_set!(mesh, :UPPER_BOTTOM, 6)
upper = Problem(Elasticity, "UPPER", 2)
upper.properties.formulation = :plane_stress
upper.elements = create_elements(mesh, "UPPER")
update!(upper.elements, "youngs modulus", 288.0)
update!(upper.elements, "poissons ratio", 1/3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)
lower = Problem(Elasticity, "LOWER", 2)
lower.properties.formulation = :plane_stress
lower.elements = create_elements(mesh, "LOWER")
update!(lower.elements, "youngs modulus", 288.0)
update!(lower.elements, "poissons ratio", 1/3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 2, "displacement")
bc_upper.elements = create_elements(mesh, "UPPER_TOP")
bc_upper_elements = create_elements(mesh, "UPPER_TOP")
#update!(bc_upper.elements, "displacement 1", -17/90)
#update!(bc_upper.elements, "displacement 1", -17/90)
update!(bc_upper.elements, "displacement 1", -0.2)
update!(bc_upper.elements, "displacement 2", -0.2)
update!(bc_upper_elements, "displacement 1", -0.2)
update!(bc_upper_elements, "displacement 2", -0.2)
add_elements!(bc_upper, bc_upper_elements)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 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)
bc_lower_elements = create_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 1", 0.0)
update!(bc_lower_elements, "displacement 2", 0.0)
add_elements!(bc_lower, bc_lower_elements)
contact = Problem(Contact2D, "LOWER_TO_UPPER", 2, "displacement")
contact_slave_elements = create_elements(mesh, "LOWER_TOP")
@@ -58,17 +60,16 @@ contact_master_elements = create_elements(mesh, "UPPER_BOTTOM")
add_slave_elements!(contact, contact_slave_elements)
add_master_elements!(contact, contact_master_elements)
solver = Solver(Nonlinear)
push!(solver, upper, lower, bc_upper, bc_lower, contact)
solver()
analysis = Analysis(Nonlinear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, contact)
run!(analysis)
master = first(contact_master_elements)
slave = first(contact_slave_elements)
um = master("displacement", (0.0,), 0.0)
us = slave("displacement", (0.0,), 0.0)
la = slave("lambda", (0.0,), 0.0)
info("um = $um, us = $us, la = $la")
@debug("um = $um, us = $us, la = $la")
@test isapprox(um, [-0.20, -0.15])
@test isapprox(us, [0.0, -0.05])
@test isapprox(la, [0.0, 30.375])
+123 -116
View File
@@ -1,14 +1,10 @@
# 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 Base.Test
using JuliaFEM, Test
function get_test_model()
X = Dict{Int64, Vector{Float64}}(
X = Dict(
1 => [0.0, 0.0],
2 => [1.0, 0.0],
3 => [1.0, 0.5],
@@ -17,143 +13,154 @@ function get_test_model()
6 => [1.0, 0.6],
7 => [1.0, 1.1],
8 => [0.0, 1.1])
el1 = Element(Quad4, [1, 2, 3, 4])
el2 = Element(Quad4, [5, 6, 7, 8])
el3 = Element(Seg2, [1, 2])
el4 = Element(Seg2, [7, 8])
el5 = Element(Seg2, [4, 3])
el6 = Element(Seg2, [5, 6])
update!([el1, el2, el3, el4, el5, el6], "geometry", X)
update!([el1, el2], "youngs modulus", 96.0)
update!([el1, el2], "poissons ratio", 1/3)
update!([el3], "displacement 1", 0.0)
update!([el3], "displacement 2", 0.0)
update!([el4], "displacement 1", 0.0)
update!([el4], "displacement 2", 0.0)
el1 = Element(Quad4, (1, 2, 3, 4))
el2 = Element(Quad4, (5, 6, 7, 8))
el3 = Element(Seg2, (1, 2))
el4 = Element(Seg2, (7, 8))
el5 = Element(Seg2, (4, 3))
el6 = Element(Seg2, (5, 6))
update!((el1, el2, el3, el4, el5, el6), "geometry", X)
update!((el1, el2), "youngs modulus", 96.0)
update!((el1, el2), "poissons ratio", 1/3)
update!(el3, "displacement 1", 0.0)
update!(el3, "displacement 2", 0.0)
update!(el4, "displacement 1", 0.0)
update!(el4, "displacement 2", 0.0)
update!(el6, "master elements", [el5])
p1 = Problem(Elasticity, "body1", 2)
p2 = Problem(Elasticity, "body2", 2)
p3 = Problem(Dirichlet, "fixed", 2, "displacement")
p4 = Problem(Mortar2D, "interface", 2, "displacement")
push!(p1, el1)
push!(p2, el2)
push!(p3, el3, el4)
push!(p4, el5, el6)
add_elements!(p1, el1)
add_elements!(p2, el2)
add_elements!(p3, el3, el4)
add_elements!(p4, el5, el6)
return p1, p2, p3, p4
end
@testset "test adjust setting in 2d tie contact" begin
p1, p2, p3, p4 = get_test_model()
p1.properties.formulation = :plane_stress
p2.properties.formulation = :plane_stress
#p4.properties.adjust = true
#p4.properties.rotate_normals = false
solver = Solver(Linear)
push!(solver, p1, p2, p3, p4)
solver()
el5 = p4.elements[1]
u = el5("displacement", [0.0], 0.0)
info("u = $u")
@test_broken isapprox(u, [0.0, 0.05])
end
## test adjust setting in 2d tie contact
p1, p2, p3, p4 = get_test_model()
p1.properties.formulation = :plane_stress
p2.properties.formulation = :plane_stress
#p4.properties.adjust = true
#p4.properties.rotate_normals = false
analysis = Analysis(Linear)
add_problems!(analysis, p1, p2, p3, p4)
run!(analysis)
el5 = first(get_elements(p4))
xi, time = (0.0, ), 0.0
u = el5("displacement", xi, time)
@debug("u = $u")
@test_broken isapprox(u, [0.0, 0.05])
@testset "test that interface transfers constant field without error" begin
meshfile = @__DIR__() * "/testdata/block_2d.med"
mesh = aster_read_mesh(meshfile)
upper = Problem(PlaneHeat, "upper", 1)
upper.elements = create_elements(mesh, "UPPER")
update!(upper.elements, "thermal conductivity", 1.0)
## test that interface transfers constant field without error
meshfile = @__DIR__() * "/testdata/block_2d.med"
mesh = aster_read_mesh(meshfile)
lower = Problem(PlaneHeat, "lower", 1)
lower.elements = create_elements(mesh, "LOWER")
update!(lower.elements, "thermal conductivity", 1.0)
upper = Problem(PlaneHeat, "upper", 1)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "thermal conductivity", 1.0)
add_elements!(upper, upper_elements)
bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
bc_upper.elements = create_elements(mesh, "UPPER_TOP")
update!(bc_upper.elements, "temperature 1", 0.0)
lower = Problem(PlaneHeat, "lower", 1)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "thermal conductivity", 1.0)
add_elements!(lower, lower_elements)
bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature")
bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower.elements, "temperature 1", 1.0)
bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
bc_upper_elements = create_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "temperature 1", 0.0)
add_elements!(bc_upper, bc_upper_elements)
interface = Problem(Mortar2D, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
add_master_elements!(interface, interface_master_elements)
add_slave_elements!(interface, interface_slave_elements)
bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature")
bc_lower_elements = create_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "temperature 1", 1.0)
add_elements!(bc_lower, bc_lower_elements)
solver = Solver(Linear)
push!(solver, upper, lower, bc_upper, bc_lower, interface)
solver()
interface = Problem(Mortar2D, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
add_master_elements!(interface, interface_master_elements)
add_slave_elements!(interface, interface_slave_elements)
#interface_norm = norm(interface.assembly)
# for bi-orthogonal:
#interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.44870723441585775, 0.44870723441585775, 0.0, 0.0, 0.0]
#interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.39361633468943247, 0.39361633468943247, 0.0, 0.0, 0.0]
#info("Interface norm: $interface_norm")
#info("Interface norm expected: $interface_norm_expected")
#@test isapprox(interface_norm, interface_norm_expected)
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface)
run!(analysis)
T_upper = first(bc_upper.elements)("temperature", [0.0], 0.0)
T_lower = first(bc_lower.elements)("temperature", [0.0], 0.0)
T_middle = first(interface.elements)("temperature", [0.0], 0.0)
info("T upper: $T_upper, T lower: $T_lower, T interface: $T_middle")
#interface_norm = norm(interface.assembly)
# for bi-orthogonal:
#interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.44870723441585775, 0.44870723441585775, 0.0, 0.0, 0.0]
#interface_norm_expected = [0.0, 0.0, 0.0, 0.0, 0.0, 0.39361633468943247, 0.39361633468943247, 0.0, 0.0, 0.0]
#@info("Interface norm: $interface_norm")
#@info("Interface norm expected: $interface_norm_expected")
#@test isapprox(interface_norm, interface_norm_expected)
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
T = [t[1] for t in temperature]
minT = minimum(T)
maxT = maximum(T)
info("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)
end
xi, time = (0.0, ), 0.0
T_upper = first(bc_upper_elements)("temperature", xi, time)
T_lower = first(bc_lower_elements)("temperature", xi, time)
T_middle = first(interface_slave_elements)("temperature", xi, time)
@debug("T upper: $T_upper, T lower: $T_lower, T interface: $T_middle")
@testset "test mesh tie with splitted block and plane stress elasticity" begin
meshfile = @__DIR__() * "/testdata/block_2d.med"
mesh = aster_read_mesh(meshfile)
node_ids, temperature = get_nodal_vector(interface_slave_elements, "temperature", 0.0)
T = [t[1] for t in temperature]
minT = minimum(T)
maxT = maximum(T)
@debug("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)
upper = Problem(Elasticity, "upper", 2)
upper.properties.formulation = :plane_stress
upper.elements = create_elements(mesh, "UPPER")
update!(upper.elements, "youngs modulus", 100.0)
update!(upper.elements, "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", 100.0)
update!(lower.elements, "poissons ratio", 1/3)
## test mesh tie with splitted block and plane stress elasticity
bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement")
bc_upper.elements = create_elements(mesh, "UPPER_TOP")
meshfile = @__DIR__() * "/testdata/block_2d.med"
mesh = aster_read_mesh(meshfile)
upper = Problem(Elasticity, "upper", 2)
upper.properties.formulation = :plane_stress
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 100.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)
lower = Problem(Elasticity, "lower", 2)
lower.properties.formulation = :plane_stress
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 100.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)
bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement")
bc_upper_elements = create_elements(mesh, "UPPER_TOP")
# update!(bc_upper.elements, "displacement 1", 0.1)
update!(bc_upper.elements, "displacement 2", -0.1)
update!(bc_upper_elements, "displacement 2", -0.1)
add_elements!(bc_upper, bc_upper_elements)
bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement")
bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
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_elements, "displacement 2", 0.0)
add_elements!(bc_lower, bc_lower_elements)
bc_corner = Problem(Dirichlet, "fix model from lower left corner to prevent singularity", 2, "displacement")
node_ids = find_nearest_nodes(mesh, [0.0, 0.0])
bc_corner.elements = [Element(Poi1, node_ids)]
update!(bc_corner.elements, "geometry", mesh.nodes)
update!(bc_corner.elements, "displacement 1", 0.0)
bc_corner = Problem(Dirichlet, "fix model from lower left corner to prevent singularity", 2, "displacement")
node_ids = find_nearest_nodes(mesh, [0.0, 0.0])
bc_corner_elements = [Element(Poi1, node_ids)]
update!(bc_corner_elements, "geometry", mesh.nodes)
update!(bc_corner_elements, "displacement 1", 0.0)
add_elements!(bc_corner, bc_corner_elements)
interface = Problem(Mortar2D, "interface between upper and lower block", 2, "displacement")
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
add_master_elements!(interface, interface_master_elements)
add_slave_elements!(interface, interface_slave_elements)
interface = Problem(Mortar2D, "interface between upper and lower block", 2, "displacement")
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
add_master_elements!(interface, interface_master_elements)
add_slave_elements!(interface, interface_slave_elements)
solver = Solver(Linear)
push!(solver, upper, lower, bc_upper, bc_lower, interface, bc_corner)
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface, bc_corner)
run!(analysis)
solver()
slave_elements = get_slave_elements(interface)
node_ids, la = get_nodal_vector(slave_elements, "lambda", 0.0)
for lai in la
@test isapprox(lai, [0.0, 10.0])
end
slave_elements = get_slave_elements(interface)
node_ids, la = get_nodal_vector(slave_elements, "lambda", 0.0)
for lai in la
@test isapprox(lai, [0.0, 10.0])
end
+2 -2
View File
@@ -1,7 +1,7 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM.Testing
using Test
#= TODO: Fix test.
@testset "calculate mortar matrices and weighted gap vector for 2d model" begin
@@ -23,7 +23,7 @@ using JuliaFEM.Testing
M = C1[:, [5, 6, 7, 8]]*24
X = calculate_nodal_vector("geometry", 2, [sel, mel], 0.0)
@debug begin
info("nodal vector")
@info("nodal vector")
dump(round(full(X), 3))
end
g = -C1*X
+40 -39
View File
@@ -1,10 +1,7 @@
# 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
using JuliaFEM, Test
#=
test subjects:
@@ -120,41 +117,45 @@ IMPR_RESU(
FIN()
"""
@testset "splitted rod with tie contact" begin
# CYLINDER_20_1_FACE1 -- CYLINDER_20_1_FACE2 -- CYLINDER_20_2_FACE_1 -- CYLINDER_20_2_FACE_2
mesh_file = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(mesh_file, "CYLINDER_20_SPLITTED")
body1 = Problem(mesh, Elasticity, "CYLINDER_20_1", 3)
body2 = Problem(mesh, Elasticity, "CYLINDER_20_2", 3)
for body in [body1, body2]
update!(body.elements, "youngs modulus", 54475.45)
update!(body.elements, "poissons ratio", 0.3)
update!(body.elements, "density", 1.0)
end
bc1 = Problem(mesh, Dirichlet, "CYLINDER_20_1_FACE1", 3, "displacement")
bc2 = Problem(mesh, Dirichlet, "CYLINDER_20_2_FACE2", 3, "displacement")
for bc in [bc1, bc2]
update!(bc.elements, "displacement 1", 0.0)
update!(bc.elements, "displacement 2", 0.0)
update!(bc.elements, "displacement 3", 0.0)
end
interface = Problem(Mortar, "interface between bodies", 3, "displacement")
slave = create_elements(mesh, "CYLINDER_20_1_FACE2")
master = create_elements(mesh, "CYLINDER_20_2_FACE1")
update!(slave, "master elements", master)
interface.elements = [slave; master]
solver = Solver(Modal, body1, body2, bc1, bc2, interface)
solver.properties.nev = 5
solver.properties.which = :SM
solver()
freqs_jf = sqrt.(solver.properties.eigvals)/(2*pi)
freqs_ca = [1.12946E+00, 1.13141E+00, 2.93779E+00, 2.94143E+00, 4.51684E+00]
for (i, freq) in enumerate(freqs_jf)
@printf "mode %i | freq JuliaFEM %8.3f | freq Code Aster %8.3f\n" i freqs_jf[i] freqs_ca[i]
end
@test isapprox(freqs_ca, freqs_jf; rtol=0.04)
# CYLINDER_20_1_FACE1 -- CYLINDER_20_1_FACE2 -- CYLINDER_20_2_FACE_1 -- CYLINDER_20_2_FACE_2
mesh_file = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(mesh_file, "CYLINDER_20_SPLITTED")
body1 = Problem(Elasticity, "CYLINDER_20_1", 3)
body2 = Problem(Elasticity, "CYLINDER_20_2", 3)
body1_elements = create_elements(mesh, "CYLINDER_20_1")
body2_elements = create_elements(mesh, "CYLINDER_20_2")
for element_set in [body1_elements, body2_elements]
update!(element_set, "youngs modulus", 54475.45)
update!(element_set, "poissons ratio", 0.3)
update!(element_set, "density", 1.0)
end
add_elements!(body1, body1_elements)
add_elements!(body2, body2_elements)
bc1 = Problem(Dirichlet, "CYLINDER_20_1_FACE1", 3, "displacement")
bc2 = Problem(Dirichlet, "CYLINDER_20_2_FACE2", 3, "displacement")
bc1_elements = create_elements(mesh, "CYLINDER_20_1_FACE1")
bc2_elements = create_elements(mesh, "CYLINDER_20_2_FACE2")
for element_set in [bc1_elements, bc2_elements]
update!(element_set, "displacement 1", 0.0)
update!(element_set, "displacement 2", 0.0)
update!(element_set, "displacement 3", 0.0)
end
add_elements!(bc1, bc1_elements)
add_elements!(bc2, bc2_elements)
interface = Problem(Mortar, "interface between bodies", 3, "displacement")
slave = create_elements(mesh, "CYLINDER_20_1_FACE2")
master = create_elements(mesh, "CYLINDER_20_2_FACE1")
update!(slave, "master elements", master)
interface.elements = [slave; master]
analysis = Analysis(Modal)
add_problems!(analysis, body1, body2, bc1, bc2, interface)
analysis.properties.nev = 5
analysis.properties.which = :SM
run!(analysis)
freqs_jf = sqrt.(analysis.properties.eigvals)/(2*pi)
freqs_ca = [1.12946E+00, 1.13141E+00, 2.93779E+00, 2.94143E+00, 4.51684E+00]
@test isapprox(freqs_ca, freqs_jf; rtol=0.04)
+39 -42
View File
@@ -1,11 +1,7 @@
# 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
using JuliaFEM, Test
using AsterReader: RMEDFile, aster_read_nodes, aster_read_data
#=
@@ -15,48 +11,49 @@ Put constant temperature 1.0 for inner surface of inner ring and 2.0 for outer
surface of outer ring. We should expect constant temperature in contact surface.
This is conforming mesh so result should match to the conforming situation.
=#
@testset "test that curved interface transfers constant field without error, two rings problem" begin
meshfile = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(meshfile, "RINGS")
ring1 = Problem(Heat, "RING1", 1)
ring1.elements = create_elements(mesh, "RING1")
update!(ring1.elements, "thermal conductivity", 1.0)
## test that curved interface transfers constant field without error, two rings problem
ring2 = Problem(Heat, "RING2", 1)
ring2.elements = create_elements(mesh, "RING2")
update!(ring2.elements, "thermal conductivity", 1.0)
meshfile = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(meshfile, "RINGS")
bc_inner = Problem(Dirichlet, "INNER SURFACE", 1, "temperature")
bc_inner.elements = create_elements(mesh, "RING1_INNER")
update!(bc_inner, "temperature 1", 1.0)
ring1 = Problem(Heat, "RING1", 1)
ring1_elements = create_elements(mesh, "RING1")
update!(ring1_elements, "thermal conductivity", 1.0)
add_elements!(ring1, ring1_elements)
bc_outer = Problem(Dirichlet, "OUTER SURFACE", 1, "temperature")
bc_outer.elements = create_elements(mesh, "RING2_OUTER")
update!(bc_outer, "temperature 1", 2.0)
ring2 = Problem(Heat, "RING2", 1)
ring2_elements = create_elements(mesh, "RING2")
update!(ring2_elements, "thermal conductivity", 1.0)
add_elements!(ring2, ring2_elements)
interface = Problem(Mortar, "interface between rings", 1, "temperature")
interface_slave = create_elements(mesh, "RING1_OUTER")
interface_master = create_elements(mesh, "RING2_INNER")
interface.elements = [interface_slave; interface_master]
update!(interface_slave, "master elements", interface_master)
bc_inner = Problem(Dirichlet, "INNER SURFACE", 1, "temperature")
bc_inner_elements = create_elements(mesh, "RING1_INNER")
update!(bc_inner_elements, "temperature 1", 1.0)
add_elements!(bc_inner, bc_inner_elements)
solver = LinearSolver(ring1, ring2, bc_inner, bc_outer, interface)
solver()
bc_outer = Problem(Dirichlet, "OUTER SURFACE", 1, "temperature")
bc_outer_elements = create_elements(mesh, "RING2_OUTER")
update!(bc_outer_elements, "temperature 1", 2.0)
add_elements!(bc_outer, bc_outer_elements)
fn = @__DIR__() * "/testdata/rings.rmed"
results = RMEDFile(fn)
nodes = aster_read_nodes(results)
temp_ca = aster_read_data(results, "TEMP")
interface = Problem(Mortar, "interface between rings", 1, "temperature")
interface_slave = create_elements(mesh, "RING1_OUTER")
interface_master = create_elements(mesh, "RING2_INNER")
update!(interface_slave, "master elements", interface_master)
add_elements!(interface, interface_slave, interface_master)
passed = true
for j in sort(collect(keys(nodes)))
X = nodes[j]
T1 = solver("temperature", X, 0.0)
T2 = temp_ca[j]
rtol = norm(T1-T2) / max(T1,T2)
@printf "% 5i : %8.5f %8.5f %8.5f | %8.5f %8.5f | %8.5f\n" j X... T1 T2 rtol
passed = passed && (rtol < 1.0e-12)
end
@test passed
end
analysis = Analysis(Linear)
add_problems!(analysis, ring1, ring2, bc_inner, bc_outer, interface)
run!(analysis)
fn = @__DIR__() * "/testdata/rings.rmed"
results = RMEDFile(fn)
nodes = aster_read_nodes(results)
temp_ca = aster_read_data(results, "TEMP")
sorted_node_ids = sort(collect(keys(nodes)))
node_coords = [nodes[j] for j in sorted_node_ids]
node_temps_jf = [analysis("temperature", node_coords[j], 0.0) for j in sorted_node_ids]
node_temps_ca = [temp_ca[j] for j in sorted_node_ids]
@test isapprox(node_temps_jf, node_temps_ca; rtol=1.0e-12)
+23 -36
View File
@@ -1,42 +1,29 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM, Test
using JuliaFEM: get_polygon_clip, calculate_polygon_area
using JuliaFEM.Testing
@testset "polygon clipping" begin
S = Vector[
[0.375, 0.0, 0.5],
[0.6, 0.0, 0.5],
[0.5, 0.25, 0.5]]
M = Vector[
[0.50, 0.0, 0.5],
[0.25, 0.0, 0.5],
[0.375, 0.25, 0.5]]
n0 = [0.0, 0.0, 1.0]
P = get_polygon_clip(S, M, n0)
@test length(P) == 3
@test isapprox(calculate_polygon_area(P), 1/128)
S = Vector[
[0.25, 0.0, 0.5],
[0.75, 0.0, 0.5],
[0.50, 0.25, 0.5]]
M = Vector[
[0.50, 0.0, 0.5],
[0.25, 0.0, 0.5],
[0.375, 0.25, 0.5]]
n0 = [0.0, 0.0, 1.0]
P = get_polygon_clip(S, M, n0)
@test length(P) == 3
@test isapprox(calculate_polygon_area(P), 1/48)
## polygon clipping
# visually inspected
Xs = Vector[[0.0, 0.0, 0.5], [1.0, 0.0, 0.5], [0.0, 1.0, 0.5]]
Xm = Vector[[-0.25, 0.50, 0.5], [0.50, -0.25, 0.5], [0.75,0.75, 0.5]]
P_ = Vector[[0.65,0.35,0.0], [0.5625,0.0,0.0], [0.25,0.0,0.0],
[0.0,0.25,0.0], [0.0,0.5625,0.0], [0.35,0.65,0.0]]
P = get_polygon_clip(Xs, Xm, [0.0, 0.0, 1.0])
@test length(P) == length(P_)
end
S = [[0.375, 0.000, 0.500], [0.600, 0.000, 0.500], [0.500, 0.250, 0.500]]
M = [[0.500, 0.000, 0.500], [0.250, 0.000, 0.500], [0.375, 0.250, 0.500]]
n0 = [0.0, 0.0, 1.0]
P = get_polygon_clip(S, M, n0)
@test length(P) == 3
@test isapprox(calculate_polygon_area(P), 1/128)
S = [[0.250, 0.000, 0.500], [0.750, 0.000, 0.500], [0.500, 0.250, 0.500]]
M = [[0.500, 0.000, 0.500], [0.250, 0.000, 0.500], [0.375, 0.250, 0.500]]
n0 = [0.0, 0.0, 1.0]
P = get_polygon_clip(S, M, n0)
@test length(P) == 3
@test isapprox(calculate_polygon_area(P), 1/48)
# visually inspected
Xs = [[ 0.00, 0.00, 0.50], [1.00, 0.00, 0.50], [0.00, 1.00, 0.50]]
Xm = [[-0.25, 0.50, 0.50], [0.50, -0.25, 0.50], [0.75, 0.75, 0.50]]
P_ = [[0.65, 0.35, 0.00], [0.5625, 0.0, 0.0], [0.25, 0.00, 0.0],
[0.00, 0.25, 0.00], [0.0, 0.5625, 0.0], [0.35, 0.65, 0.0]]
P = get_polygon_clip(Xs, Xm, [0.0, 0.0, 1.0])
@test length(P) == length(P_)
+18 -50
View File
@@ -1,55 +1,23 @@
# 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.Testing
using JuliaFEM, Test
@testset "get nodal values" begin
el1 = Element(Seg2, [1, 2])
el2 = Element(Seg2, [2, 3])
X = Dict{Int64, Vector{Float64}}(
1 => [0.0],
2 => [1.0],
3 => [2.0])
T = Dict{Int64, Vector{Float64}}(
1 => [0.0],
2 => [1.0],
3 => [0.0])
P = Problem(Heat, "foo", 1)
push!(P, el1, el2)
update!(P, "geometry", X)
update!(P, "temperature", T)
@test isnan(P("temperature", [-0.1]))
@test isapprox(P("temperature", [0.0]), [0.0])
@test isapprox(P("temperature", [0.5]), [0.5])
@test isapprox(P("temperature", [1.0]), [1.0])
@test isapprox(P("temperature", [1.5]), [0.5])
@test isapprox(P("temperature", [2.0]), [0.0])
@test isnan(P("temperature", [ 2.1]))
end
@testset "interpolate from set of elements" begin
el1 = Element(Seg2, [1, 2])
el2 = Element(Seg2, [2, 3])
X = Dict{Int64, Vector{Float64}}(
1 => [0.0],
2 => [1.0],
3 => [2.0])
T = Dict{Int64, Vector{Float64}}(
1 => [0.0],
2 => [1.0],
3 => [0.0])
P = Problem(Heat, "foo", 1)
push!(P, el1, el2)
update!(P, "geometry", X)
update!(P, "temperature", T)
@test isnan(P("temperature", [-0.1]))
@test isapprox(P("temperature", [0.0]), [0.0])
@test isapprox(P("temperature", [0.5]), [0.5])
@test isapprox(P("temperature", [1.0]), [1.0])
@test isapprox(P("temperature", [1.5]), [0.5])
@test isapprox(P("temperature", [2.0]), [0.0])
@test isnan(P("temperature", [ 2.1]))
end
## interpolate from a set of elements
element1 = Element(Seg2, (1, 2))
element2 = Element(Seg2, (2, 3))
X = Dict(1 => [0.0], 2 => [1.0], 3 => [2.0])
T = Dict(1 => [0.0], 2 => [1.0], 3 => [0.0])
problem = Problem(Heat, "foo", 1)
add_elements!(problem, element1, element2)
problem_elements = get_elements(problem)
update!(problem_elements, "geometry", X)
update!(problem_elements, "temperature", T)
@test isnan(problem("temperature", [-0.1], 0.0))
@test isapprox(problem("temperature", [0.0], 0.0), [0.0])
@test isapprox(problem("temperature", [0.5], 0.0), [0.5])
@test isapprox(problem("temperature", [1.0], 0.0), [1.0])
@test isapprox(problem("temperature", [1.5], 0.0), [0.5])
@test isapprox(problem("temperature", [2.0], 0.0), [0.0])
@test isnan(problem("temperature", [ 2.1], 0.0))
+3 -3
View File
@@ -2,7 +2,7 @@
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Testing
using Test
abstract type HeatProblem <: AbstractProblem
@@ -63,7 +63,7 @@ function test_potential_energy_method()
solve!(problem, [1, 2], 0.0)
temp = element("temperature", [0.0, -1.0], 0.0)
err = temp - 2/3
info("error: $err")
@info("error: $err")
@test isapprox(err, 0.0)
end
@@ -93,7 +93,7 @@ function test_potential_energy_method_2()
temp = element1("temperature", [0.0, -1.0], 0.0)
err = temp - 0.5
info("error: $err")
@info("error: $err")
@test isapprox(err, 0.0, atol=1.0e-6)
# @test isapprox(temp, 2.93509690572300E+00) # tested using Code Aster
+98 -90
View File
@@ -1,35 +1,36 @@
# 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
using JuliaFEM, Statistics, Test
tet4_meshfile = "test_problems_contact_3d/tet4.inp"
tet10_meshfile = "test_problems_contact_3d/tet10.inp"
function get_model(meshfile)
mesh = abaqus_read_mesh(meshfile)
upper = Problem(Elasticity, "UPPER", 3)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "youngs modulus", 3*288.0)
update!(upper, "poissons ratio", 1/3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 3*288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)
lower = Problem(Elasticity, "LOWER", 3)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "youngs modulus", 288.0)
update!(lower, "poissons ratio", 1/3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "displacement 3", -0.4)
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", -0.4)
add_elements!(bc_upper, bc_upper_elements)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "displacement 3", 0.0)
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)
# point-wise boundary conditions to prevent free body move
nid1 = find_nearest_nodes(mesh, [0.0, 0.0, 0.0])[1]
nid2 = find_nearest_nodes(mesh, [1.0, 0.0, 0.0])[1]
@@ -40,97 +41,104 @@ function get_model(meshfile)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
# nodes in X2=0 plane
bc_sym13.elements = [Element(Poi1, [j]) for j in [nid1, nid2, nid4, nid5]]
update!(bc_sym13, "geometry", mesh.nodes)
update!(bc_sym13, "displacement 2", 0.0)
bc_sym13_elements = [Element(Poi1, [j]) for j in [nid1, nid2, nid4, nid5]]
update!(bc_sym13_elements, "geometry", mesh.nodes)
update!(bc_sym13_elements, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
# nodes in X1=0 plane
bc_sym23.elements = [Element(Poi1, [j]) for j in [nid1, nid3, nid4, nid6]]
update!(bc_sym23, "geometry", mesh.nodes)
update!(bc_sym23, "displacement 1", 0.0)
bc_sym23_elements = [Element(Poi1, [j]) for j in [nid1, nid3, nid4, nid6]]
update!(bc_sym23_elements, "geometry", mesh.nodes)
update!(bc_sym23_elements, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements)
interface = Problem(Contact, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
add_elements!(interface, interface_slave_elements, interface_master_elements)
interface.properties.contact_state_in_first_iteration = :AUTO
#append!(bc_sym13.assembly.removed_dofs, [1316, 1319, 1388, 1358])
#append!(bc_sym23.assembly.removed_dofs, [1492, 1627, 1387, 1597])
#append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])
solver = NonlinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
analysis = Analysis(Nonlinear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
# solver.properties.max_iterations = 5
return solver
return analysis
end
@testset "small sliding contact patch test, tet4 + standard basis" begin
solver = get_model(tet4_meshfile)
add_results_writer!(solver, Xdmf("contact_sl_lin_disp_results"; overwrite=true))
interface = solver["LOWER_TO_UPPER"]
interface.properties.dual_basis = false
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
# test postprocess of fields
postprocess!(interface, 0.0, Val{Symbol("contact pressure")})
node_ids, contact_pressure = get_nodal_vector(interface.elements, "contact pressure", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
maxpres = maximum(contact_pressure)
stdpres = std(contact_pressure)
info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
info("max(contact_pressure) = $maxpres, std(contact_pressure) = $stdpres")
@test isapprox(stdabsu3, 0.0; atol=1.0e-12)
@test isapprox(maxpres, 172.8; atol=1.0e-6)
end
## small sliding contact patch test, tet4 + standard basis
analysis = get_model(tet4_meshfile)
xdmf = Xdmf("contact_sl_lin_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
interface = get_problem(analysis, "LOWER_TO_UPPER")
interface.properties.dual_basis = false
run!(analysis)
u = interface("displacement", 0.0)
X = interface("geometry", 0.0)
# test postprocess of fields
postprocess!(interface, 0.0, Val{Symbol("contact pressure")})
contact_pressure = interface("contact pressure", 0.0)
@debug("Contact pressure in interface", contact_pressure)
u3 = [u[3] for u in values(u)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
maxpres = maximum(values(contact_pressure))
stdpres = std(values(contact_pressure))
@debug("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@debug("max(contact_pressure) = $maxpres, std(contact_pressure) = $stdpres")
@test isapprox(stdabsu3, 0.0; atol=1.0e-12)
@test isapprox(maxpres, 172.8; atol=1.0e-6)
@testset "small sliding contact patch test, tet4 + dual basis" begin
solver = get_model(tet4_meshfile)
add_results_writer!(solver, Xdmf("contact_dl_lin_disp_results"; overwrite=true))
interface = solver["LOWER_TO_UPPER"]
interface.properties.dual_basis = true
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-12)
end
@testset "small sliding contact patch test, tet10 + standard basis" begin
solver = get_model(tet10_meshfile)
add_results_writer!(solver, Xdmf("contact_sl_quad_disp_results"; overwrite=true))
interface = solver["LOWER_TO_UPPER"]
interface.properties.dual_basis = false
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)
end
## small sliding contact patch test, tet4 + dual basis
analysis = get_model(tet4_meshfile)
xdmf = Xdmf("contact_dl_lin_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
interface = get_problem(analysis, "LOWER_TO_UPPER")
interface.properties.dual_basis = true
run!(analysis)
u = interface("displacement", 0.0)
X = interface("geometry", 0.0)
u3 = [u[3] for u in values(u)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-12)
@testset "small sliding contact patch test, tet10 + dual basis, alpha=0.2" begin
solver = get_model(tet10_meshfile)
add_results_writer!(solver, Xdmf("contact_dl_quad_disp_results"; overwrite=true))
interface = solver["LOWER_TO_UPPER"]
interface.properties.dual_basis = true
interface.properties.alpha = 0.2
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)
end
## small sliding contact patch test, tet10 + standard basis
analysis = get_model(tet10_meshfile)
xdmf = Xdmf("contact_sl_quad_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
interface = get_problem(analysis, "LOWER_TO_UPPER")
interface.properties.dual_basis = false
run!(analysis)
u = interface("displacement", 0.0)
X = interface("geometry", 0.0)
u3 = [u[3] for u in values(u)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)
## small sliding contact patch test, tet10 + dual basis, alpha=0.2
analysis = get_model(tet10_meshfile)
xdmf = Xdmf("contact_dl_quad_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
interface = get_problem(analysis, "LOWER_TO_UPPER")
interface.properties.dual_basis = true
interface.properties.alpha = 0.2
run!(analysis)
u = interface("displacement", 0.0)
X = interface("geometry", 0.0)
u3 = [u[3] for u in values(u)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)
+20 -27
View File
@@ -1,13 +1,11 @@
# 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
using JuliaFEM, Test
@testset "test postprocessing of strain and stress" begin
# test postprocessing of strain and stress
X = Dict(
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [1.0, 1.0, 0.0],
@@ -16,8 +14,8 @@ using JuliaFEM.Testing
6 => [1.0, 0.0, 1.0],
7 => [1.0, 1.0, 1.0],
8 => [0.0, 1.0, 1.0])
u = Dict(
u = Dict(
1 => [0.0, 0.0, 0.0],
2 => [-1/3, 0.0, 0.0],
3 => [-1/3, -1/3, 0.0],
@@ -27,27 +25,22 @@ using JuliaFEM.Testing
7 => [-1/3, -1/3, 1.0],
8 => [0.0, -1/3, 1.0])
element = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8])
update!(element, "geometry", 0.0 => X)
update!(element, "displacement", 0.0 => u)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
element = Element(Hex8, (1, 2, 3, 4, 5, 6, 7, 8))
update!(element, "geometry", 0.0 => X)
update!(element, "displacement", 0.0 => u)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
body = Problem(Elasticity, "[0,1]³ elastic unit block", 3)
body.elements= [element]
postprocess!(body, 0.0, Val{:strain})
postprocess!(body, 0.0, Val{:stress})
body = Problem(Elasticity, "[0,1]³ elastic unit block", 3)
add_elements!(body, element)
postprocess!(body, 0.0, Val{:strain})
postprocess!(body, 0.0, Val{:stress})
geom = element("geometry", 0.0)
strain = element("strain", 0.0)
stress = element("stress", 0.0)
for i in 1:length(geom)
a = geom[i]
b = strain[i]
c = stress[i]
@test isapprox(b, [-1/3, -1/3, 1.0, 0.0, 0.0, 0.0])
@test isapprox(c, [0.0, 0.0, 288.0, 0.0, 0.0, 0.0])
end
geom = element("geometry", 0.0)
strain = element("strain", 0.0)
stress = element("stress", 0.0)
for i in 1:length(geom)
@test isapprox(strain[i], [-1/3, -1/3, 1.0, 0.0, 0.0, 0.0])
@test isapprox(stress[i], [0.0, 0.0, 288.0, 0.0, 0.0, 0.0])
end
+412 -367
View File
@@ -1,460 +1,505 @@
# 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
using JuliaFEM, Test, Statistics
### temperature patch tests, sl tet4, dl tet4, sl tet10, dl tet 10
tet4_meshfile = "test_problems_mortar_3d/tet4.inp"
tet10_meshfile = "test_problems_mortar_3d/tet10.inp"
tet4_meshfile = joinpath("test_problems_mortar_3d", "tet4.inp")
tet10_meshfile = joinpath("test_problems_mortar_3d", "tet10.inp")
@testset "patch test temperature + abaqus inp + tet4" begin
mesh = abaqus_read_mesh(tet4_meshfile)
# patch test temperature + abaqus inp + tet4
mesh = abaqus_read_mesh(tet4_meshfile)
upper = Problem(Heat, "UPPER", 1)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "thermal conductivity", 1.0)
upper = Problem(Heat, "UPPER", 1)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "thermal conductivity", 1.0)
add_elements!(upper, upper_elements)
lower = Problem(Heat, "LOWER", 1)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "thermal conductivity", 1.0)
lower = Problem(Heat, "LOWER", 1)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "thermal conductivity", 1.0)
add_elements!(lower, lower_elements)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "temperature 1", 0.0)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "temperature 1", 0.0)
add_elements!(bc_upper, bc_upper_elements)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "temperature 1", 1.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "temperature 1", 1.0)
add_elements!(bc_lower, bc_lower_elements)
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_master_elements; interface_slave_elements]
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)
JuliaFEM.diagnose_interface(interface, 0.0)
solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
add_results_writer!(solver, Xdmf("sl_lin_temp_results"; overwrite=true))
JuliaFEM.diagnose_interface(interface, 0.0)
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface)
xdmf = Xdmf("sl_lin_temp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)
solver()
T = values(interface("temperature", 0.0))
minT = minimum(T)
maxT = maximum(T)
@info("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
T = [t[1] for t in temperature]
minT = minimum(T)
maxT = maximum(T)
info("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)
#=
initialize!(solver)
assemble!(solver)
M, K, Kg, f, fg = get_field_assembly(solver)
Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
K = K + Kg + Kb
f = f + fg + fb
K = 1/2*(K + K')
M = 1/2*(M + M')
=#
#=
initialize!(solver)
assemble!(solver)
M, K, Kg, f, fg = get_field_assembly(solver)
Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
K = K + Kg + Kb
f = f + fg + fb
K = 1/2*(K + K')
M = 1/2*(M + M')
=#
end
# patch test temperature + abaqus inp + tet4 + dual basis + adjust
mesh = abaqus_read_mesh(tet4_meshfile)
@testset "patch test temperature + abaqus inp + tet4 + dual basis + adjust" begin
mesh = abaqus_read_mesh(tet4_meshfile)
upper = Problem(Heat, "UPPER", 1)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "thermal conductivity", 1.0)
add_elements!(upper, upper_elements)
upper = Problem(Heat, "UPPER", 1)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "thermal conductivity", 1.0)
lower = Problem(Heat, "LOWER", 1)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "thermal conductivity", 1.0)
add_elements!(lower, lower_elements)
lower = Problem(Heat, "LOWER", 1)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "thermal conductivity", 1.0)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "temperature 1", 0.0)
add_elements!(bc_upper, bc_upper_elements)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "temperature 1", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "temperature 1", 1.0)
add_elements!(bc_lower, bc_lower_elements)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "temperature 1", 1.0)
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_master_elements, interface_slave_elements)
interface.properties.dual_basis = true
#interface.properties.adjust = true
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_master_elements; interface_slave_elements]
interface.properties.dual_basis = true
#interface.properties.adjust = true
JuliaFEM.diagnose_interface(interface, 0.0)
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface)
xdmf = Xdmf("dl_lin_temp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)
JuliaFEM.diagnose_interface(interface, 0.0)
solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
add_results_writer!(solver, Xdmf("dl_lin_temp_results"; overwrite=true))
T = values(interface("temperature", 0.0))
minT = minimum(T)
maxT = maximum(T)
@debug("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)
solver()
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
T = [t[1] for t in temperature]
minT = minimum(T)
maxT = maximum(T)
info("minT = $minT, maxT = $maxT")
@test isapprox(minT, 0.5)
@test isapprox(maxT, 0.5)
# patch test temperature + abaqus inp + tet10, quadratic surface elements
mesh = abaqus_read_mesh(tet10_meshfile)
end
upper = Problem(Heat, "UPPER", 1)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "thermal conductivity", 1.0)
add_elements!(upper, upper_elements)
@testset "patch test temperature + abaqus inp + tet10, quadratic surface elements" begin
mesh = abaqus_read_mesh(tet10_meshfile)
lower = Problem(Heat, "LOWER", 1)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "thermal conductivity", 1.0)
add_elements!(lower, lower_elements)
upper = Problem(Heat, "UPPER", 1)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "thermal conductivity", 1.0)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "temperature 1", 0.0)
add_elements!(bc_upper, bc_upper_elements)
lower = Problem(Heat, "LOWER", 1)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "thermal conductivity", 1.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "temperature 1", 1.0)
add_elements!(bc_lower, bc_lower_elements)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "temperature 1", 0.0)
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_master_elements, interface_slave_elements)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "temperature 1", 1.0)
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_master_elements; interface_slave_elements]
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.alpha = 0.0
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.alpha = 0.0
# JuliaFEM.diagnose_interface(interface, 0.0)
solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
add_results_writer!(solver, Xdmf("sl_quad_temp_results"; overwrite=true))
solver()
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface)
xdmf = Xdmf("sl_quad_temp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
T = [t[1] for t in temperature]
T = values(interface("temperature", 0.0))
minT = minimum(T)
maxT = maximum(T)
stdT = std(T)
@info("minT = $minT, maxT = $maxT, stdT = $stdT")
@test maxT - minT < 1.0e-6
@test isapprox(stdT, 0.0; atol=1.0e-6)
minT = minimum(T)
maxT = maximum(T)
stdT = std(T)
info("minT = $minT, maxT = $maxT, stdT = $stdT")
@test maxT - minT < 1.0e-6
@test isapprox(stdT, 0.0; atol=1.0e-6)
end
# patch test temperature + abaqus inp + tet10 + quadratic surface elements + dual basis + alpha=0.2
mesh = abaqus_read_mesh(tet10_meshfile)
@testset "patch test temperature + abaqus inp + tet10 + quadratic surface elements + dual basis + alpha=0.2" begin
mesh = abaqus_read_mesh(tet10_meshfile)
upper = Problem(Heat, "UPPER", 1)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "thermal conductivity", 1.0)
add_elements!(upper, upper_elements)
upper = Problem(Heat, "UPPER", 1)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "thermal conductivity", 1.0)
lower = Problem(Heat, "LOWER", 1)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "thermal conductivity", 1.0)
add_elements!(lower, lower_elements)
lower = Problem(Heat, "LOWER", 1)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "thermal conductivity", 1.0)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "temperature 1", 0.0)
add_elements!(bc_upper, bc_upper_elements)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "temperature 1", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "temperature 1", 1.0)
add_elements!(bc_lower, bc_lower_elements)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 1, "temperature")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "temperature 1", 1.0)
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_master_elements; interface_slave_elements]
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.dual_basis = true
interface.properties.alpha = 0.2
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.dual_basis = true
interface.properties.alpha = 0.2
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, interface)
xdmf = Xdmf("dl_quad_temp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)
solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
add_results_writer!(solver, Xdmf("dl_quad_temp_results"; overwrite=true))
solver()
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
#node_ids, temperature = get_nodal_vector(interface_slave_elements, "temperature", 0.0)
#=
for (j, (nid, T)) in enumerate(zip(node_ids, temperature))
info("$j: $nid -> $(T[1])")
j == 10 && break
end
=#
T = [t[1] for t in temperature]
minT = minimum(T)
maxT = maximum(T)
stdT = std(T)
info("minT = $minT, maxT = $maxT, stdT = $stdT")
@test maxT - minT < 1.0e-10
@test isapprox(stdT, 0.0; atol=1.0e-10)
end
T = values(interface("temperature", 0.0))
minT = minimum(T)
maxT = maximum(T)
stdT = std(T)
@info("minT = $minT, maxT = $maxT, stdT = $stdT")
@test maxT - minT < 1.0e-10
@test isapprox(stdT, 0.0; atol=1.0e-10)
### displacement patch tests, sl tet4, dl tet4, sl tet10, dl tet 10
@testset "patch test displacement + abaqus inp + tet4 + adjust" begin
# patch test displacement + abaqus inp + tet4 + adjust
mesh = abaqus_read_mesh(tet4_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
upper = Problem(Elasticity, "UPPER", 3)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "youngs modulus", 288.0)
update!(upper, "poissons ratio", 1/3)
lower = Problem(Elasticity, "LOWER", 3)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "youngs modulus", 288.0)
update!(lower, "poissons ratio", 1/3)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "displacement 3", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "displacement 3", 0.0)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")]
update!(bc_sym13, "displacement 2", 0.0)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")]
update!(bc_sym23, "displacement 1", 0.0)
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = false
solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
add_results_writer!(solver, Xdmf("sl_lin_disp_results"; overwrite=true))
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)
mesh = abaqus_read_mesh(tet4_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
@testset "patch test displacement + abaqus inp + tet4 + adjust + dual basis" begin
upper = Problem(Elasticity, "UPPER", 3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)
mesh = abaqus_read_mesh(tet4_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
lower = Problem(Elasticity, "LOWER", 3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)
upper = Problem(Elasticity, "UPPER", 3)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "youngs modulus", 288.0)
update!(upper, "poissons ratio", 1/3)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", 0.0)
add_elements!(bc_upper, bc_upper_elements)
lower = Problem(Elasticity, "LOWER", 3)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "youngs modulus", 288.0)
update!(lower, "poissons ratio", 1/3)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "displacement 3", 0.0)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13_elements_1 = create_surface_elements(mesh, "LOWER_SYM13")
bc_sym13_elements_2 = create_surface_elements(mesh, "UPPER_SYM13")
update!(bc_sym13_elements_1, "displacement 2", 0.0)
update!(bc_sym13_elements_2, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements_1, bc_sym13_elements_2)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "displacement 3", 0.0)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23_elements_1 = create_surface_elements(mesh, "LOWER_SYM23")
bc_sym23_elements_2 = create_surface_elements(mesh, "UPPER_SYM23")
update!(bc_sym23_elements_1, "displacement 1", 0.0)
update!(bc_sym23_elements_2, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements_1, bc_sym23_elements_2)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")]
update!(bc_sym13, "displacement 2", 0.0)
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")]
update!(bc_sym23, "displacement 1", 0.0)
append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = false
append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
xdmf = Xdmf("sl_lin_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = true
displacement = interface("displacement", 0.0)
u3 = [u[3] for u in values(displacement)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)
solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
add_results_writer!(solver, Xdmf("dl_lin_disp_results"; overwrite=true))
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)
# patch test displacement + abaqus inp + tet4 + adjust + dual basis
mesh = abaqus_read_mesh(tet4_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
@testset "patch test displacement + abaqus inp + tet10 + adjust" begin
upper = Problem(Elasticity, "UPPER", 3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)
mesh = abaqus_read_mesh(tet10_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
lower = Problem(Elasticity, "LOWER", 3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)
upper = Problem(Elasticity, "UPPER", 3)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "youngs modulus", 288.0)
update!(upper, "poissons ratio", 1/3)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", 0.0)
add_elements!(bc_upper, bc_upper_elements)
lower = Problem(Elasticity, "LOWER", 3)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "youngs modulus", 288.0)
update!(lower, "poissons ratio", 1/3)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "displacement 3", 0.0)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13_elements_1 = create_surface_elements(mesh, "LOWER_SYM13")
bc_sym13_elements_2 = create_surface_elements(mesh, "UPPER_SYM13")
update!(bc_sym13_elements_1, "displacement 2", 0.0)
update!(bc_sym13_elements_2, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements_1, bc_sym13_elements_2)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "displacement 3", 0.0)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23_elements_1 = create_surface_elements(mesh, "LOWER_SYM23")
bc_sym23_elements_2 = create_surface_elements(mesh, "UPPER_SYM23")
update!(bc_sym23_elements_1, "displacement 1", 0.0)
update!(bc_sym23_elements_2, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements_1, bc_sym23_elements_2)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")]
update!(bc_sym13, "displacement 2", 0.0)
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")]
update!(bc_sym23, "displacement 1", 0.0)
append!(interface.assembly.removed_dofs, [1316, 1319, 1358, 1387, 1388, 1492, 1597, 1627])
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = true
removed_dofs = [1316, 1319, 1325, 1358, 1361, 1387, 1388, 1391, 1492, 1597, 1600, 1627, 1630, 1657]
append!(interface.assembly.removed_dofs, removed_dofs)
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
xdmf = Xdmf("dl_lin_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = false
displacement = interface("displacement", 0.0)
u3 = [u[3] for u in values(displacement)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-10)
solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
add_results_writer!(solver, Xdmf("sl_quad_disp_results"; overwrite=true))
solver()
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-6)
# patch test displacement + abaqus inp + tet10 + adjust
mesh = abaqus_read_mesh(tet10_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
upper = Problem(Elasticity, "UPPER", 3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)
@testset "patch test displacement + abaqus inp + tet10 + adjust + dual basis + alpha=0.2" begin
lower = Problem(Elasticity, "LOWER", 3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)
mesh = abaqus_read_mesh(tet10_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.Preprocess.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", 0.0)
add_elements!(bc_upper, bc_upper_elements)
upper = Problem(Elasticity, "UPPER", 3)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "youngs modulus", 288.0)
update!(upper, "poissons ratio", 1/3)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)
lower = Problem(Elasticity, "LOWER", 3)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "youngs modulus", 288.0)
update!(lower, "poissons ratio", 1/3)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13_elements_1 = create_surface_elements(mesh, "LOWER_SYM13")
bc_sym13_elements_2 = create_surface_elements(mesh, "UPPER_SYM13")
update!(bc_sym13_elements_1, "displacement 2", 0.0)
update!(bc_sym13_elements_2, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements_1, bc_sym13_elements_2)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "displacement 3", 0.0)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23_elements_1 = create_surface_elements(mesh, "LOWER_SYM23")
bc_sym23_elements_2 = create_surface_elements(mesh, "UPPER_SYM23")
update!(bc_sym23_elements_1, "displacement 1", 0.0)
update!(bc_sym23_elements_2, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements_1, bc_sym23_elements_2)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "displacement 3", 0.0)
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13.elements = [create_surface_elements(mesh, "LOWER_SYM13"); create_surface_elements(mesh, "UPPER_SYM13")]
update!(bc_sym13, "displacement 2", 0.0)
removed_dofs = [1316, 1319, 1325, 1358, 1361, 1387, 1388, 1391, 1492, 1597, 1600, 1627, 1630, 1657]
append!(interface.assembly.removed_dofs, removed_dofs)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23.elements = [create_surface_elements(mesh, "LOWER_SYM23"); create_surface_elements(mesh, "UPPER_SYM23")]
update!(bc_sym23, "displacement 1", 0.0)
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = false
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
xdmf = Xdmf("sl_quad_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)
removed_dofs = [1316, 1319, 1325, 1358, 1361, 1387, 1388, 1391, 1492, 1597, 1600, 1627, 1630, 1657]
append!(interface.assembly.removed_dofs, removed_dofs)
displacement = interface("displacement", 0.0)
u3 = [u[3] for u in values(displacement)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-6)
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = true
interface.properties.alpha = 0.2
solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
add_results_writer!(solver, Xdmf("dl_quad_disp_results"; overwrite=true))
solver()
# patch test displacement + abaqus inp + tet10 + adjust + dual basis + alpha=0.2
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-6)
mesh = abaqus_read_mesh(tet10_meshfile)
# modify mesh a bit, find all nodes in elements in element set UPPER and put 0.2 to X3 to test adjust
JuliaFEM.create_node_set_from_element_set!(mesh, :UPPER)
for nid in mesh.node_sets[:UPPER]
mesh.nodes[nid][3] += 0.2
end
upper = Problem(Elasticity, "UPPER", 3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)
lower = Problem(Elasticity, "LOWER", 3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", 0.0)
add_elements!(bc_upper, bc_upper_elements)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
bc_sym13_elements_1 = create_surface_elements(mesh, "LOWER_SYM13")
bc_sym13_elements_2 = create_surface_elements(mesh, "UPPER_SYM13")
update!(bc_sym13_elements_1, "displacement 2", 0.0)
update!(bc_sym13_elements_2, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements_1, bc_sym13_elements_2)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
bc_sym23_elements_1 = create_surface_elements(mesh, "LOWER_SYM23")
bc_sym23_elements_2 = create_surface_elements(mesh, "UPPER_SYM23")
update!(bc_sym23_elements_1, "displacement 1", 0.0)
update!(bc_sym23_elements_2, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements_1, bc_sym23_elements_2)
interface = Problem(Mortar, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
add_elements!(interface, interface_slave_elements, interface_master_elements)
removed_dofs = [1316, 1319, 1325, 1358, 1361, 1387, 1388, 1391, 1492, 1597, 1600, 1627, 1630, 1657]
append!(interface.assembly.removed_dofs, removed_dofs)
interface.properties.linear_surface_elements = false
interface.properties.split_quadratic_slave_elements = false
interface.properties.split_quadratic_master_elements = false
interface.properties.adjust = true
interface.properties.dual_basis = true
interface.properties.alpha = 0.2
analysis = Analysis(Linear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
xdmf = Xdmf("dl_quad_disp_results"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)
displacement = interface("displacement", 0.0)
u3 = [u[3] for u in values(displacement)]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-6)
+124 -166
View File
@@ -1,28 +1,7 @@
# 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 "forget to add elements to problem" begin
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [0.0, 1.0, 0.0],
4 => [-0.25, 0.50, 0.00],
5 => [0.50, -0.25, 0.00],
6 => [0.75, 0.75, 0.00])
s = Element(Tri3, [1, 2, 3])
m = Element(Tri3, [4, 5, 6])
update!([s, m], "geometry", X)
update!(s, "master elements", [m])
p = Problem(Mortar, "two elements", 1, "temperature")
initialize!(p)
assemble!(p)
@test true
end
using JuliaFEM, SparseArrays, LinearAlgebra, Test
""" Calculate mortar projection matrix P = D^-1*M from mortar assembly. """
function calculate_mortar_projection_matrix(problem::Problem{Mortar}, ndim::Int)
@@ -30,17 +9,18 @@ function calculate_mortar_projection_matrix(problem::Problem{Mortar}, ndim::Int)
C1 = sparse(problem.assembly.C1, ndim, ndim)
C2 = sparse(problem.assembly.C2, ndim, ndim)
@assert nnz(sparse(problem.assembly.K)) == 0
@assert nnz(sparse(problem.assembly.D)) == 0
@assert nnz(sparse(problem.assembly.Kg)) == 0
@assert nnz(sparse(problem.assembly.fg)) == 0
@assert nnz(sparse(problem.assembly.f)) == 0
@assert nnz(sparse(problem.assembly.g)) == 0
@debug("problem.assembly", problem.assembly.K, problem.assembly.D,
problem.assembly.Kg, problem.assembly.fg, problem.assembly.f,
problem.assembly.g)
@assert isempty(problem.assembly.K)
@assert isempty(problem.assembly.D)
@assert isempty(problem.assembly.Kg)
@assert isempty(problem.assembly.fg)
@assert isempty(problem.assembly.f)
@assert C1 == C2
#@assert problem.properties.dual_basis == true
@assert problem.properties.adjust == false
S = get_nonzero_rows(C2)
M = setdiff(get_nonzero_columns(C2), S)
@@ -48,148 +28,126 @@ function calculate_mortar_projection_matrix(problem::Problem{Mortar}, ndim::Int)
D_ = C2[S,S]
M_ = -C2[S,M]
#=
P = nothing
if !isdiag(D_)
warn("D is not diagonal, is dual basis used? This might take a long time.")
P = ldltfact(1/2*(D_ + D_')) \ M_
else
P = D_ \ M_
end
=#
P = lufact(D_) \ full(M_)
P = lu(D_) \ Matrix(M_)
return S, M, P
end
@testset "two linear element clipping, calculation of projection matrix P for standard and dual basis" begin
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [0.0, 1.0, 0.0],
4 => [-0.25, 0.50, 0.00],
5 => [0.50, -0.25, 0.00],
6 => [0.75, 0.75, 0.00])
s = Element(Tri3, [1, 2, 3])
m = Element(Tri3, [4, 5, 6])
update!([s, m], "geometry", X)
update!(s, "master elements", [m])
p = Problem(Mortar, "two elements", 1, "temperature")
p.properties.dual_basis = false
p.elements = [s; m]
initialize!(p)
assemble!(p)
C1 = sparse(p.assembly.C1)
C2 = sparse(p.assembly.C2)
D = sparse(p.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(p, 6)
@test S == [1, 2, 3]
@test M == [4, 5, 6]
# visually inspected to be ok result
P_expected = 1/15*[9 9 -3; -7 13 9; 13 -7 9]
@test isapprox(P, P_expected)
um = [7.5, 15.0, 22.5]
@test isapprox(P*um, [9.0, 23.0, 13.0])
# two linear element clipping, calculation of projection matrix P for
# standard and dual basis
empty!(p.assembly)
p.properties.dual_basis = true
assemble!(p)
C1 = sparse(p.assembly.C1)
C2 = sparse(p.assembly.C2)
D = sparse(p.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(p, 6)
@test S == [1, 2, 3]
@test M == [4, 5, 6]
@test isapprox(P, P_expected)
um = [7.5, 15.0, 22.5]
@test isapprox(P*um, [9.0, 23.0, 13.0])
end
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [0.0, 1.0, 0.0],
4 => [-0.25, 0.50, 0.00],
5 => [0.50, -0.25, 0.00],
6 => [0.75, 0.75, 0.00])
slave = Element(Tri3, (1, 2, 3))
master = Element(Tri3, (4, 5, 6))
update!((slave, master), "geometry", X)
update!(slave, "master elements", [master])
problem = Problem(Mortar, "two elements", 1, "temperature")
problem.properties.dual_basis = false
initialize!(problem, 0.0)
assemble!(problem, 0.0)
# forget to add elements to problem
add_elements!(problem, slave, master)
initialize!(problem, 0.0)
assemble!(problem, 0.0)
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
D = sparse(problem.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(problem, 6)
@test S == [1, 2, 3]
@test M == [4, 5, 6]
# visually inspected to be ok result
P_expected = 1/15*[9 9 -3; -7 13 9; 13 -7 9]
@test isapprox(P, P_expected)
um = [7.5, 15.0, 22.5]
@test isapprox(P*um, [9.0, 23.0, 13.0])
empty!(problem.assembly)
problem.properties.dual_basis = true
assemble!(problem, 0.0)
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
D = sparse(problem.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(problem, 6)
@test S == [1, 2, 3]
@test M == [4, 5, 6]
@test isapprox(P, P_expected)
um = [7.5, 15.0, 22.5]
@test isapprox(P*um, [9.0, 23.0, 13.0])
@testset "two quadratic element clipping, calculation of projection matrix P for standard basis" begin
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [0.0, 1.0, 0.0],
7 => [-0.25, 0.50, 0.00],
8 => [0.50, -0.25, 0.00],
9 => [0.75, 0.75, 0.00])
# middle nodes
X[4] = 1/2*(X[1] + X[2])
X[5] = 1/2*(X[2] + X[3])
X[6] = 1/2*(X[3] + X[1])
X[10] = 1/2*(X[7] + X[8])
X[11] = 1/2*(X[8] + X[9])
X[12] = 1/2*(X[9] + X[7])
s = Element(Tri6, [1, 2, 3, 4, 5, 6])
m = Element(Tri6, [7, 8, 9, 10, 11, 12])
update!([s, m], "geometry", X)
update!(s, "master elements", [m])
p = Problem(Mortar, "two elements", 1, "temperature")
p.properties.dual_basis = false
p.properties.alpha = 0.2
p.elements = [s; m]
initialize!(p)
assemble!(p)
C1 = sparse(p.assembly.C1)
C2 = sparse(p.assembly.C2)
D = sparse(p.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(p, 12)
@test S == [1, 2, 3, 4, 5, 6]
@test M == [7, 8, 9, 10, 11, 12]
println(full(P))
# visually inspected to be ok result
P_expected = 1/675*[81 81 189 972 -324 -324; 609 429 81 -1092 1404 -756; 429 609 81 -1092 -756 1404; -39 231 -81 132 396 36; -81 -81 81 108 324 324; 231 -39 -81 132 36 396]
@test isapprox(P, P_expected)
um = 15/2*[1, 2, 3]
um = [um[1], um[2], um[3], 0.5*(um[1]+um[2]), 0.5*(um[2]+um[3]), 0.5*(um[3]+um[1])]
us = P*um
@test isapprox(us, [9.0, 23.0, 13.0, 16.0, 18.0, 11.0])
end
# two quadratic element clipping, calculation of projection matrix P
# for standard basis
@testset "two quadratic element clipping, calculation of projection matrix P for dual lagrange basis" begin
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [0.0, 1.0, 0.0],
7 => [-0.25, 0.50, 0.00],
8 => [0.50, -0.25, 0.00],
9 => [0.75, 0.75, 0.00])
# middle nodes
X[4] = 1/2*(X[1] + X[2])
X[5] = 1/2*(X[2] + X[3])
X[6] = 1/2*(X[3] + X[1])
X[10] = 1/2*(X[7] + X[8])
X[11] = 1/2*(X[8] + X[9])
X[12] = 1/2*(X[9] + X[7])
s = Element(Tri6, [1, 2, 3, 4, 5, 6])
m = Element(Tri6, [7, 8, 9, 10, 11, 12])
update!([s, m], "geometry", X)
update!(s, "master elements", [m])
p = Problem(Mortar, "two elements", 1, "temperature")
p.properties.dual_basis = true
p.properties.alpha = 0.2
p.elements = [s; m]
initialize!(p)
assemble!(p)
C1 = sparse(p.assembly.C1)
C2 = sparse(p.assembly.C2)
D = sparse(p.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(p, 12)
P_expected = 1/675*[81 81 189 972 -324 -324; 609 429 81 -1092 1404 -756; 429 609 81 -1092 -756 1404; -39 231 -81 132 396 36; -81 -81 81 108 324 324; 231 -39 -81 132 36 396]
@test isapprox(P, P_expected)
um = 15/2*[1, 2, 3]
um = [um[1], um[2], um[3], 0.5*(um[1]+um[2]), 0.5*(um[2]+um[3]), 0.5*(um[3]+um[1])]
us = P*um
@test isapprox(us, [9.0, 23.0, 13.0, 16.0, 18.0, 11.0])
end
X = Dict(
1 => [0.0, 0.0, 0.0],
2 => [1.0, 0.0, 0.0],
3 => [0.0, 1.0, 0.0],
7 => [-0.25, 0.50, 0.00],
8 => [0.50, -0.25, 0.00],
9 => [0.75, 0.75, 0.00])
# middle nodes
X[4] = 1/2*(X[1] + X[2])
X[5] = 1/2*(X[2] + X[3])
X[6] = 1/2*(X[3] + X[1])
X[10] = 1/2*(X[7] + X[8])
X[11] = 1/2*(X[8] + X[9])
X[12] = 1/2*(X[9] + X[7])
slave = Element(Tri6, (1, 2, 3, 4, 5, 6))
master = Element(Tri6, (7, 8, 9, 10, 11, 12))
update!((slave, master), "geometry", X)
update!(slave, "master elements", [master])
problem = Problem(Mortar, "two elements", 1, "temperature")
problem.properties.dual_basis = false
problem.properties.alpha = 0.2
add_elements!(problem, slave, master)
initialize!(problem, 0.0)
assemble!(problem, 0.0)
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
D = sparse(problem.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(problem, 12)
@test S == [1, 2, 3, 4, 5, 6]
@test M == [7, 8, 9, 10, 11, 12]
@debug("Projection matrix P", Matrix(P))
# visually inspected to be ok result
P_expected = 1/675*[81 81 189 972 -324 -324; 609 429 81 -1092 1404 -756; 429 609 81 -1092 -756 1404; -39 231 -81 132 396 36; -81 -81 81 108 324 324; 231 -39 -81 132 36 396]
@test isapprox(P, P_expected)
um = 15/2*[1, 2, 3]
um = [um[1], um[2], um[3], 0.5*(um[1]+um[2]), 0.5*(um[2]+um[3]), 0.5*(um[3]+um[1])]
us = P*um
@test isapprox(us, [9.0, 23.0, 13.0, 16.0, 18.0, 11.0])
# two quadratic element clipping, calculation of projection matrix P for
# dual lagrange basis
problem.properties.dual_basis = true
problem.properties.alpha = 0.2
empty!(problem.assembly)
initialize!(problem, 0.0)
assemble!(problem, 0.0)
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
D = sparse(problem.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(problem, 12)
P_expected = 1/675*[81 81 189 972 -324 -324; 609 429 81 -1092 1404 -756; 429 609 81 -1092 -756 1404; -39 231 -81 132 396 36; -81 -81 81 108 324 324; 231 -39 -81 132 36 396]
@test isapprox(P, P_expected)
um = 15/2*[1, 2, 3]
um = [um[1], um[2], um[3], 0.5*(um[1]+um[2]), 0.5*(um[2]+um[3]), 0.5*(um[3]+um[1])]
us = P*um
@test isapprox(us, [9.0, 23.0, 13.0, 16.0, 18.0, 11.0])
+80 -77
View File
@@ -1,94 +1,97 @@
# 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
using JuliaFEM, Test, Statistics
datadir = first(splitext(basename(@__FILE__)))
@testset "test postprocessing of secondary fields" begin
mesh = abaqus_read_mesh(joinpath(datadir, "tet4.inp"))
upper = Problem(Elasticity, "UPPER", 3)
upper.elements = create_elements(mesh, "UPPER")
update!(upper, "youngs modulus", 3*288.0)
update!(upper, "poissons ratio", 1/3)
push!(upper.postprocess_fields, "strain", "stress")
info("upper postprocess: $(upper.postprocess_fields)")
mesh = abaqus_read_mesh(joinpath(datadir, "tet4.inp"))
lower = Problem(Elasticity, "LOWER", 3)
lower.elements = create_elements(mesh, "LOWER")
update!(lower, "youngs modulus", 288.0)
update!(lower, "poissons ratio", 1/3)
push!(lower.postprocess_fields, "strain", "stress")
info("lower postprocess: $(lower.postprocess_fields)")
upper = Problem(Elasticity, "UPPER", 3)
upper_elements = create_elements(mesh, "UPPER")
update!(upper_elements, "youngs modulus", 3*288.0)
update!(upper_elements, "poissons ratio", 1/3)
add_elements!(upper, upper_elements)
push!(upper.postprocess_fields, "strain", "stress")
@info("upper postprocess: $(upper.postprocess_fields)")
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper, "displacement 3", -0.4)
lower = Problem(Elasticity, "LOWER", 3)
lower_elements = create_elements(mesh, "LOWER")
update!(lower_elements, "youngs modulus", 288.0)
update!(lower_elements, "poissons ratio", 1/3)
add_elements!(lower, lower_elements)
push!(lower.postprocess_fields, "strain", "stress")
@info("lower postprocess: $(lower.postprocess_fields)")
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower.elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower, "displacement 3", 0.0)
# point-wise boundary conditions to prevent free body move
nid1 = find_nearest_nodes(mesh, [0.0, 0.0, 0.0])[1]
nid2 = find_nearest_nodes(mesh, [1.0, 0.0, 0.0])[1]
nid3 = find_nearest_nodes(mesh, [0.0, 1.0, 0.0])[1]
nid4 = find_nearest_nodes(mesh, [0.0, 0.0, 1.0])[1]
nid5 = find_nearest_nodes(mesh, [1.0, 0.0, 1.0])[1]
nid6 = find_nearest_nodes(mesh, [0.0, 1.0, 1.0])[1]
bc_upper = Problem(Dirichlet, "UPPER_TOP", 3, "displacement")
bc_upper_elements = create_surface_elements(mesh, "UPPER_TOP")
update!(bc_upper_elements, "displacement 3", -0.4)
add_elements!(bc_upper, bc_upper_elements)
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
# nodes in X2=0 plane
bc_sym13.elements = [Element(Poi1, [j]) for j in [nid1, nid2, nid4, nid5]]
update!(bc_sym13, "geometry", mesh.nodes)
update!(bc_sym13, "displacement 2", 0.0)
bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 3, "displacement")
bc_lower_elements = create_surface_elements(mesh, "LOWER_BOTTOM")
update!(bc_lower_elements, "displacement 3", 0.0)
add_elements!(bc_lower, bc_lower_elements)
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
# nodes in X1=0 plane
bc_sym23.elements = [Element(Poi1, [j]) for j in [nid1, nid3, nid4, nid6]]
update!(bc_sym23, "geometry", mesh.nodes)
update!(bc_sym23, "displacement 1", 0.0)
# point-wise boundary conditions to prevent free body move
nid1 = find_nearest_nodes(mesh, [0.0, 0.0, 0.0])[1]
nid2 = find_nearest_nodes(mesh, [1.0, 0.0, 0.0])[1]
nid3 = find_nearest_nodes(mesh, [0.0, 1.0, 0.0])[1]
nid4 = find_nearest_nodes(mesh, [0.0, 0.0, 1.0])[1]
nid5 = find_nearest_nodes(mesh, [1.0, 0.0, 1.0])[1]
nid6 = find_nearest_nodes(mesh, [0.0, 1.0, 1.0])[1]
for bc in [bc_upper, bc_lower, bc_sym13, bc_sym23]
push!(bc.postprocess_fields, "reaction force")
end
bc_sym13 = Problem(Dirichlet, "SYM13", 3, "displacement")
# nodes in X2=0 plane
bc_sym13_elements = [Element(Poi1, [j]) for j in [nid1, nid2, nid4, nid5]]
update!(bc_sym13_elements, "geometry", mesh.nodes)
update!(bc_sym13_elements, "displacement 2", 0.0)
add_elements!(bc_sym13, bc_sym13_elements)
interface = Problem(Contact, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
interface.properties.dual_basis = true
interface.properties.contact_state_in_first_iteration = :AUTO
push!(interface.postprocess_fields, "contact pressure")
bc_sym23 = Problem(Dirichlet, "SYM23", 3, "displacement")
# nodes in X1=0 plane
bc_sym23_elements = [Element(Poi1, [j]) for j in [nid1, nid3, nid4, nid6]]
update!(bc_sym23_elements, "geometry", mesh.nodes)
update!(bc_sym23_elements, "displacement 1", 0.0)
add_elements!(bc_sym23, bc_sym23_elements)
solver = NonlinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface)
push!(bc_upper.postprocess_fields, "reaction force")
push!(bc_lower.postprocess_fields, "reaction force")
push!(bc_sym13.postprocess_fields, "reaction force")
push!(bc_sym23.postprocess_fields, "reaction force")
xdmf = Xdmf("contact_two_blocks_postprocess"; overwrite=true)
add_results_writer!(solver, xdmf)
solver()
interface = Problem(Contact, "LOWER_TO_UPPER", 3, "displacement")
interface_slave_elements = create_surface_elements(mesh, "LOWER_TO_UPPER")
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
update!(interface_slave_elements, "master elements", interface_master_elements)
interface.elements = [interface_slave_elements; interface_master_elements]
interface.properties.dual_basis = true
interface.properties.contact_state_in_first_iteration = :AUTO
push!(interface.postprocess_fields, "contact pressure")
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
# node_ids, pressure = get_nodal_vector(interface.elements, "contact pressure", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
info("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-12)
upper_X = [0.5, 0.5, 0.75]
lower_X = [0.5, 0.5, 0.25]
strain_upper = upper("strain", upper_X)
stress_upper = upper("stress", upper_X)
strain_lower = lower("strain", lower_X)
stress_lower = lower("stress", lower_X)
info("strain at $upper_X = $strain_upper")
info("stress at $upper_X = $stress_upper")
info("strain at $lower_X = $strain_lower")
info("stress at $lower_X = $stress_lower")
end
analysis = Analysis(Nonlinear)
add_problems!(analysis, upper, lower, bc_upper, bc_lower,
bc_sym13, bc_sym23, interface)
xdmf = Xdmf("contact_two_blocks_postprocess"; overwrite=true)
add_results_writer!(analysis, xdmf)
run!(analysis)
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
# node_ids, pressure = get_nodal_vector(interface.elements, "contact pressure", 0.0)
u3 = [u[3] for u in displacement]
maxabsu3 = maximum(abs.(u3))
stdabsu3 = std(abs.(u3))
@debug("max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
@test isapprox(stdabsu3, 0.0; atol=1.0e-12)
upper_X = [0.5, 0.5, 0.75]
lower_X = [0.5, 0.5, 0.25]
strain_upper = upper("strain", upper_X, 0.0)
stress_upper = upper("stress", upper_X, 0.0)
strain_lower = lower("strain", lower_X, 0.0)
stress_lower = lower("stress", lower_X, 0.0)
@debug("strain at $upper_X = $strain_upper")
@debug("stress at $upper_X = $stress_upper")
@debug("strain at $lower_X = $strain_lower")
@debug("stress at $lower_X = $stress_lower")
+2 -4
View File
@@ -1,8 +1,7 @@
# 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.Testing
using JuliaFEM, Test
abstract type PlaneStressElasticityProblem <: AbstractProblem end
@@ -61,9 +60,8 @@ function test_residual_form()
free_dofs = [3, 4, 5, 6]
solve!(problem, free_dofs, 0.0) # launch a newton solver for single element
disp = element("displacement", [1.0, 1.0], 0.0)
info("displacement at tip: $disp")
@info("displacement at tip: $disp")
# verified using Code Aster.
@test isapprox(disp[2], -8.77303119819776E+00)
end
+4 -4
View File
@@ -3,9 +3,9 @@
#using PyPlot
#using JuliaFEM
#using JuliaFEM.Testing
#using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress, State
#using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress
using Test
#=
function test_von_mises_3D_basic()
@@ -59,7 +59,7 @@ function test_von_mises_3D_basic()
a[3, 2] = b[4]
end
info("Starting calculation")
@info("Starting calculation")
tic()
params = Dict("yield_stress" => stress_y)
stress_new = zeros(Float64, 6)
@@ -117,7 +117,7 @@ function test_von_mises_3D_basic()
end
info("Calculation finished")
@info("Calculation finished")
# plot3D(ee, ss)
@@ -224,7 +224,7 @@ function test_von_mises_planestress_basic()
eig_stress = zeros(Float64, (3, 3))
eig_vals = zeros(Float64, (steps, 3))
info("Starting calculation")
@info("Starting calculation")
tic()
stress_new = zeros(Float64, 3)