mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-30 13:29:22 +00:00
linear elasticity. looks we are having problems with results in 3d
This commit is contained in:
Jukka Aho
@@ -2,6 +2,34 @@
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
|
||||
function aster_parse_nodes(section::ASCIIString; strip_characters=true)
|
||||
nodes = Dict{Any, Vector{Float64}}()
|
||||
has_started = false
|
||||
for line in split(section, '\n')
|
||||
m = matchall(r"[\w.-]+", line)
|
||||
if (length(m) != 1) && (!has_started)
|
||||
continue
|
||||
end
|
||||
if length(m) == 1
|
||||
if (m[1] == "COOR_2D") || (m[1] == "COOR_3D")
|
||||
has_started = true
|
||||
continue
|
||||
end
|
||||
if m[1] == "FINSF"
|
||||
break
|
||||
end
|
||||
end
|
||||
if length(m) == 4
|
||||
nid = m[1]
|
||||
if strip_characters
|
||||
nid = matchall(r"\d", nid)
|
||||
nid = parse(Int, nid[1])
|
||||
end
|
||||
nodes[nid] = float(m[2:end])
|
||||
end
|
||||
end
|
||||
return nodes
|
||||
end
|
||||
|
||||
function parse(mesh::ASCIIString, ::Type{Val{:CODE_ASTER_MAIL}})
|
||||
model = Dict{ASCIIString, Any}()
|
||||
|
||||
@@ -80,6 +80,7 @@ include("equations.jl")
|
||||
include("dirichlet.jl")
|
||||
include("heat.jl")
|
||||
include("elasticity.jl")
|
||||
include("linear_elasticity.jl")
|
||||
|
||||
### ASSEMBLY + SOLVE ###
|
||||
include("assembly.jl")
|
||||
|
||||
+1
-1
@@ -79,10 +79,10 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element:
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
u = element("displacement", time, variation)
|
||||
grad = element(ip, time, Val{:grad})
|
||||
# gradu = element("displacement", ip, time, Val{:grad}, variation)
|
||||
gradu = grad*u
|
||||
|
||||
F = I + gradu # deformation gradient
|
||||
# info("gradu = \n$(ForwardDiff.get_value(gradu))")
|
||||
|
||||
young = element("youngs modulus", ip, time)
|
||||
poisson = element("poissons ratio", ip, time)
|
||||
|
||||
+27
-10
@@ -6,20 +6,22 @@
|
||||
|
||||
### 1d elements
|
||||
|
||||
function get_integration_points(::Type{Seg2}, ::Type{Val{1}})
|
||||
typealias LineElement Union{Type{Seg2}, Type{Seg3}}
|
||||
|
||||
function get_integration_points(::LineElement, ::Type{Val{1}})
|
||||
[
|
||||
IntegrationPoint([0.0], 2.0)
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::Type{Seg2}, ::Type{Val{2}})
|
||||
function get_integration_points(::LineElement, ::Type{Val{2}})
|
||||
[
|
||||
IntegrationPoint([-sqrt(1/3)], 1)
|
||||
IntegrationPoint([+sqrt(1/3)], 1)
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::Type{Seg3}, ::Type{Val{3}})
|
||||
function get_integration_points(::LineElement, ::Type{Val{3}})
|
||||
[
|
||||
IntegrationPoint([0.0], 8/9),
|
||||
IntegrationPoint([-sqrt(3/5)], 5/9),
|
||||
@@ -27,7 +29,7 @@ function get_integration_points(::Type{Seg3}, ::Type{Val{3}})
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::Type{Seg3}, ::Type{Val{4}})
|
||||
function get_integration_points(::LineElement, ::Type{Val{4}})
|
||||
[
|
||||
IntegrationPoint([+sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
|
||||
IntegrationPoint([-sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
|
||||
@@ -36,7 +38,7 @@ function get_integration_points(::Type{Seg3}, ::Type{Val{4}})
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::Union{Type{Seg2}, Type{Seg3}}, ::Type{Val{5}})
|
||||
function get_integration_points(::LineElement, ::Type{Val{5}})
|
||||
[
|
||||
IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900),
|
||||
IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),
|
||||
@@ -58,16 +60,16 @@ end
|
||||
|
||||
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
|
||||
|
||||
typealias TriangularElements Union{Type{Tri3}, Type{Tri6}}
|
||||
typealias TriangularElement Union{Type{Tri3}, Type{Tri6}}
|
||||
|
||||
function get_integration_points(::TriangularElements, ::Type{Val{1}})
|
||||
function get_integration_points(::TriangularElement, ::Type{Val{1}})
|
||||
# http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
|
||||
[
|
||||
IntegrationPoint([1.0/3.0, 1.0/3.0], 0.5)
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::TriangularElements, ::Type{Val{2}})
|
||||
function get_integration_points(::TriangularElement, ::Type{Val{2}})
|
||||
# http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
|
||||
[
|
||||
IntegrationPoint([2.0/3.0, 1.0/6.0], 1.0/6.0),
|
||||
@@ -76,7 +78,7 @@ function get_integration_points(::TriangularElements, ::Type{Val{2}})
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::TriangularElements, ::Type{Val{5}})
|
||||
function get_integration_points(::TriangularElement, ::Type{Val{5}})
|
||||
# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
|
||||
# FIXME: something wrong here with weights ..?
|
||||
[
|
||||
@@ -95,7 +97,7 @@ function get_integration_points(::Type{Tri3})
|
||||
end
|
||||
|
||||
|
||||
function get_integration_points(::Type{Quad4})
|
||||
function get_integration_points(::Type{Quad4}, ::Type{Val{2}})
|
||||
[
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
|
||||
IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
|
||||
@@ -104,8 +106,23 @@ function get_integration_points(::Type{Quad4})
|
||||
]
|
||||
end
|
||||
|
||||
function get_integration_points(::Type{Quad4})
|
||||
return get_integration_points(Quad4, Val{2})
|
||||
end
|
||||
|
||||
### 3d elements
|
||||
|
||||
|
||||
function get_integration_points(::Type{Hex8}, ::Type{Val{2}})
|
||||
p = 1.0/sqrt(3.0)*[-1.0, 1.0]
|
||||
w = 1.0
|
||||
return vec([IntegrationPoint([p[i], p[j], p[k]], w) for i=1:2, j=1:2, k=1:2])
|
||||
end
|
||||
|
||||
function get_integration_points(::Type{Hex8})
|
||||
return get_integration_points(Hex8, Val{2})
|
||||
end
|
||||
|
||||
function get_integration_points(::Type{Tet4})
|
||||
# http://libmesh.github.io/doxygen/quadrature__gauss__3D_8C_source.html
|
||||
[
|
||||
|
||||
@@ -113,6 +113,20 @@ end
|
||||
|
||||
# 3d Lagrange elements
|
||||
|
||||
@create_lagrange_element(Hex8, "8 node hexahedra",
|
||||
[-1.0 1.0 1.0 -1.0 -1.0 1.0 1.0 -1.0
|
||||
-1.0 -1.0 1.0 1.0 -1.0 -1.0 1.0 1.0
|
||||
-1.0 -1.0 -1.0 -1.0 1.0 1.0 1.0 1.0],
|
||||
(xi) -> [
|
||||
1.0,
|
||||
xi[1],
|
||||
xi[2],
|
||||
xi[1]*xi[2],
|
||||
xi[3],
|
||||
xi[1]*xi[3],
|
||||
xi[2]*xi[3],
|
||||
xi[1]*xi[2]*xi[3]])
|
||||
|
||||
@create_lagrange_element(Tet4, "4 node tetrahedron",
|
||||
[0.0 1.0 0.0 0.0
|
||||
0.0 0.0 1.0 0.0
|
||||
|
||||
@@ -0,0 +1,105 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# Linear elasticity
|
||||
|
||||
abstract LinearElasticityProblem <: ElasticityProblem
|
||||
|
||||
function LinearElasticityProblem(dim::Int=3, elements=[])
|
||||
return Problem{LinearElasticityProblem}(dim, elements)
|
||||
end
|
||||
|
||||
""" Elasticity equations, general 3D case. """
|
||||
function assemble!{E<:CG, P<:LinearElasticityProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real)
|
||||
|
||||
gdofs = get_gdofs(element, problem.dim)
|
||||
ndim, nnodes = size(E)
|
||||
B = zeros(6, 3*nnodes)
|
||||
for ip in get_integration_points(element)
|
||||
w = ip.weight*det(element, ip, time)
|
||||
N = element(ip, time)
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
v = element("poissons ratio", ip, time)
|
||||
E_ = element("youngs modulus", ip, time)
|
||||
a = 1 - v
|
||||
b = 1 - 2*v
|
||||
c = 1 + v
|
||||
C = E_/(b*c) .* [
|
||||
a v v 0 0 0
|
||||
v a v 0 0 0
|
||||
v v a 0 0 0
|
||||
0 0 0 b 0 0
|
||||
0 0 0 0 b 0
|
||||
0 0 0 0 0 b]
|
||||
dN = element(ip, time, Val{:grad})
|
||||
fill!(B, 0.0)
|
||||
for i=1:size(dN, 2)
|
||||
B[1, 3*(i-1)+1] = dN[1,i]
|
||||
B[2, 3*(i-1)+2] = dN[2,i]
|
||||
B[3, 3*(i-1)+3] = dN[3,i]
|
||||
B[4, 3*(i-1)+1] = dN[2,i]
|
||||
B[4, 3*(i-1)+2] = dN[1,i]
|
||||
B[5, 3*(i-1)+2] = dN[3,i]
|
||||
B[5, 3*(i-1)+3] = dN[2,i]
|
||||
B[6, 3*(i-1)+1] = dN[3,i]
|
||||
B[6, 3*(i-1)+3] = dN[1,i]
|
||||
end
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B)
|
||||
end
|
||||
if haskey(element, "displacement load")
|
||||
b = element("displacement load", ip, time)
|
||||
add!(assembly.force_vector, gdofs, w*N'*b)
|
||||
end
|
||||
if haskey(element, "displacement traction force")
|
||||
T = element("displacement traction force", ip, time)
|
||||
L = w*T*N
|
||||
# dump(L)
|
||||
add!(assembly.force_vector, gdofs, vec(L))
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
abstract PlaneStressLinearElasticityProblem <: LinearElasticityProblem
|
||||
|
||||
function PlaneStressLinearElasticityProblem(dim::Int=2, elements=[])
|
||||
return Problem{PlaneStressLinearElasticityProblem}(dim, elements)
|
||||
end
|
||||
|
||||
""" Elasticity equations, plane stress. """
|
||||
function assemble!{E<:CG, P<:PlaneStressLinearElasticityProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real)
|
||||
|
||||
gdofs = get_gdofs(element, problem.dim)
|
||||
ndim, nnodes = size(E)
|
||||
B = zeros(3, 2*nnodes)
|
||||
for ip in get_integration_points(element)
|
||||
w = ip.weight*det(element, ip, time)
|
||||
N = element(ip, time)
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
nu = element("poissons ratio", ip, time)
|
||||
E_ = element("youngs modulus", ip, time)
|
||||
C = E_/(1.0 - nu^2) .* [
|
||||
1.0 nu 0.0
|
||||
nu 1.0 0.0
|
||||
0.0 0.0 (1.0-nu)/2.0]
|
||||
dN = element(ip, time, Val{:grad})
|
||||
fill!(B, 0.0)
|
||||
for i=1:size(dN, 2)
|
||||
B[1, 2*(i-1)+1] = dN[1,i]
|
||||
B[2, 2*(i-1)+2] = dN[2,i]
|
||||
B[3, 2*(i-1)+1] = dN[2,i]
|
||||
B[3, 2*(i-1)+2] = dN[1,i]
|
||||
end
|
||||
add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B)
|
||||
end
|
||||
if haskey(element, "displacement load")
|
||||
b = element("displacement load", ip, time)
|
||||
add!(assembly.force_vector, gdofs, w*N'*b)
|
||||
end
|
||||
if haskey(element, "displacement traction force")
|
||||
T = element("displacement traction force", ip, time)
|
||||
L = w*T*N
|
||||
# dump(L)
|
||||
add!(assembly.force_vector, gdofs, vec(L))
|
||||
end
|
||||
end
|
||||
end
|
||||
@@ -6,7 +6,8 @@ module AsterReaderTests
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Test
|
||||
|
||||
using JuliaFEM: parse
|
||||
#using JuliaFEM: parse
|
||||
using JuliaFEM.Preprocess: aster_parse_nodes
|
||||
|
||||
function test_read_mesh()
|
||||
mesh = """
|
||||
@@ -44,7 +45,30 @@ mesh = """
|
||||
@test m["nsets"]["NALL"] == ["N1", "N2"]
|
||||
|
||||
end
|
||||
#test_read_mesh()
|
||||
|
||||
test_read_mesh()
|
||||
function test_parse_nodes()
|
||||
section = """
|
||||
N9 2.0 3.0 4.0
|
||||
COOR_3D
|
||||
N1 0.0 0.0 0.0
|
||||
N2 1.0 0.0 0.0
|
||||
N3 1.0 1.0 0.0
|
||||
N4 0.0 1.0 0.0
|
||||
N5 0.0 0.0 1.0
|
||||
N6 1.0 0.0 1.0
|
||||
N7 1.0 1.0 1.0
|
||||
N8 0.0 1.0 1.0
|
||||
FINSF
|
||||
absdflasdf
|
||||
N12 3.0 4.0 5.0 6.0
|
||||
N13 3.0 4.0 5.0
|
||||
"""
|
||||
nodes = aster_parse_nodes(section)
|
||||
@test nodes[1] == Float64[0.0, 0.0, 0.0]
|
||||
@test nodes[8] == Float64[0.0, 1.0, 1.0]
|
||||
@test length(nodes) == 8
|
||||
end
|
||||
test_parse_nodes()
|
||||
|
||||
end
|
||||
|
||||
+106
-2
@@ -4,7 +4,11 @@
|
||||
module ElasticityTests
|
||||
|
||||
using JuliaFEM.Test
|
||||
using JuliaFEM.Core: Seg2, Quad4, PlaneStressElasticityProblem, solve!
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Core: Seg2, Quad4, Hex8,
|
||||
ElasticityProblem, PlaneStressElasticityProblem,
|
||||
solve!, get_connectivity, DirichletProblem
|
||||
|
||||
|
||||
function test_elasticity_volume_load()
|
||||
element = Quad4([1, 2, 3, 4])
|
||||
@@ -19,9 +23,19 @@ function test_elasticity_volume_load()
|
||||
free_dofs = [3, 4, 5, 6]
|
||||
solve!(problem, free_dofs, 0.0; max_iterations=10)
|
||||
disp = element("displacement", [1.0, 1.0], 0.0)
|
||||
# function get_previous_ip(element::Element, current_ip::IntegrationPoint)
|
||||
# end
|
||||
# ipdata = element("integration points", time) => IntegrationPoint[ip1, ip2, ..., ipN]
|
||||
# for some_ip in ipdata
|
||||
# if isapprox(some_ip.xi, ip.xi)
|
||||
# info("found")
|
||||
# last_value = some_ip("material parameter", time)
|
||||
# break
|
||||
# end
|
||||
# end
|
||||
#ip1 = last(element["integration points"])[1]
|
||||
#ip2 = last(element["integration points"])[2]
|
||||
#strain = ip1["gl strain"]
|
||||
# strain = ip1("gl strain")
|
||||
info("displacement at tip: $disp")
|
||||
#info("strain in first ip: $strain. ip coord = $(ip1.xi) and weight = $(ip1.weight)")
|
||||
# verified using Code Aster, verification/2015-10-22-plane-stress/cplan_grot_gdep_volume_force.resu
|
||||
@@ -29,6 +43,7 @@ function test_elasticity_volume_load()
|
||||
end
|
||||
#test_elasticity_volume_load()
|
||||
|
||||
|
||||
function test_elasticity_surface_load()
|
||||
N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
|
||||
|
||||
@@ -55,4 +70,93 @@ function test_elasticity_surface_load()
|
||||
@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
|
||||
end
|
||||
|
||||
|
||||
function test_continuum_elasticity_with_surface_load()
|
||||
nodes = JuliaFEM.Preprocess.aster_parse_nodes("""
|
||||
COOR_3D
|
||||
N1 0.0 0.0 0.0
|
||||
N2 1.0 0.0 0.0
|
||||
N3 1.0 1.0 0.0
|
||||
N4 0.0 1.0 0.0
|
||||
N5 0.0 0.0 1.0
|
||||
N6 1.0 0.0 1.0
|
||||
N7 1.0 1.0 1.0
|
||||
N8 0.0 1.0 1.0
|
||||
FINSF
|
||||
""")
|
||||
|
||||
function set_geometry!(element, nodes)
|
||||
element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)]
|
||||
end
|
||||
element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8])
|
||||
set_geometry!(element1, nodes)
|
||||
# element1["youngs modulus"] = 900.0
|
||||
# element1["poissons ratio"] = 0.25
|
||||
element1["youngs modulus"] = 9000.0
|
||||
element1["poissons ratio"] = 0.25
|
||||
element1["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:8])
|
||||
|
||||
element2 = Quad4([5, 6, 7, 8])
|
||||
set_geometry!(element2, nodes)
|
||||
element2["displacement traction force"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]
|
||||
element2["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:4])
|
||||
|
||||
problem = ElasticityProblem()
|
||||
push!(problem, element1)
|
||||
push!(problem, element2)
|
||||
|
||||
#=
|
||||
free_dofs = zeros(Bool, 8, 3)
|
||||
x = 1
|
||||
y = 2
|
||||
z = 3
|
||||
free_dofs[2, x] = true
|
||||
free_dofs[3, [x, y]] = true
|
||||
free_dofs[4, y] = true
|
||||
free_dofs[5, z] = true
|
||||
free_dofs[6, [x, z]] = true
|
||||
free_dofs[7, [x, y, z]] = true
|
||||
free_dofs[8, [y, z]] = true
|
||||
free_dofs = find(vec(free_dofs'))
|
||||
info("free dofs: $free_dofs")
|
||||
|
||||
info("initial force vector")
|
||||
ass = JuliaFEM.Core.assemble(problem, 0.0)
|
||||
info(reshape(full(ass.force_vector), 3, 8))
|
||||
info("initial stiffness matrix")
|
||||
dump(round(Int, full(ass.stiffness_matrix))[free_dofs, free_dofs])
|
||||
solve!(problem, free_dofs, 0.0; max_iterations=10)
|
||||
=#
|
||||
|
||||
dx = Quad4([1, 4, 8, 5])
|
||||
dx["displacement 1"] = 0.0
|
||||
dy = Quad4([1, 5, 6, 2])
|
||||
dy["displacement 2"] = 0.0
|
||||
dz = Quad4([1, 2, 3, 4])
|
||||
dz["displacement 3"] = 0.0
|
||||
bc = DirichletProblem("displacement", 3)
|
||||
for el in [dx, dy, dz]
|
||||
set_geometry!(el, nodes)
|
||||
push!(bc, el)
|
||||
end
|
||||
solver = JuliaFEM.Core.DirectSolver()
|
||||
push!(solver, problem)
|
||||
push!(solver, bc)
|
||||
solver.dump_matrices = true
|
||||
solver.name = "3d_hex8"
|
||||
solver(0.0)
|
||||
|
||||
disp = element1("displacement", [1.0, 1.0, 1.0], 0.0)
|
||||
info("displacement at tip: $disp")
|
||||
info("displacement on element: ")
|
||||
for (i, d) in enumerate(element1("displacement", 0.0))
|
||||
@printf "%d %f %f %f\n" [i;d]...
|
||||
end
|
||||
# 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])
|
||||
@test isapprox(disp, [2.80559539222183E-03, 2.80559539222183E-03, -1.13019918093242E-02])
|
||||
end
|
||||
#test_continuum_elasticity_with_surface_load()
|
||||
|
||||
end
|
||||
|
||||
@@ -0,0 +1,123 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
module LinearElasticityTests
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Test
|
||||
|
||||
using JuliaFEM.Core: Seg2, Quad4, Hex8, LinearElasticityProblem, get_connectivity,
|
||||
assemble, PlaneStressLinearElasticityProblem
|
||||
using JuliaFEM.Preprocess: aster_parse_nodes
|
||||
|
||||
|
||||
function test_plane_stress_linear_elasticity_with_surface_load()
|
||||
nodes = Dict{Int64, Vector{Float64}}(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [1.0, 0.0],
|
||||
3 => [1.0, 1.0],
|
||||
4 => [0.0, 1.0])
|
||||
|
||||
function set_geometry!(element, nodes)
|
||||
element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)]
|
||||
end
|
||||
element1 = Quad4([1, 2, 3, 4])
|
||||
set_geometry!(element1, nodes)
|
||||
element1["youngs modulus"] = 9000.0
|
||||
element1["poissons ratio"] = 0.25
|
||||
|
||||
element2 = Seg2([3, 4])
|
||||
set_geometry!(element2, nodes)
|
||||
element2["displacement traction force"] = Vector{Float64}[[0.0, -100.0] for i=1:2]
|
||||
|
||||
problem = PlaneStressLinearElasticityProblem()
|
||||
push!(problem, element1)
|
||||
push!(problem, element2)
|
||||
|
||||
free_dofs = Int64[3, 5, 6, 8]
|
||||
|
||||
info("initial force vector")
|
||||
ass = assemble(problem, 0.0)
|
||||
f = full(ass.force_vector)
|
||||
K = full(ass.stiffness_matrix)
|
||||
dump(reshape(f, 2, 4))
|
||||
info("initial stiffness matrix")
|
||||
dump(round(Int, K)[free_dofs, free_dofs])
|
||||
|
||||
u = zeros(2, 4)
|
||||
u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs]
|
||||
|
||||
info("result vector")
|
||||
dump(u)
|
||||
# verified using Code Aster.
|
||||
# 2015-10-22-plane-stress/cplan_linear_traction_force.*
|
||||
@test isapprox(u[:,3], [2.77777777777778E-03, -1.11111111111111E-02])
|
||||
end
|
||||
#test_plane_stress_linear_elasticity_with_surface_load()
|
||||
|
||||
|
||||
function test_continuum_elasticity_with_surface_load()
|
||||
nodes = JuliaFEM.Preprocess.aster_parse_nodes("""
|
||||
COOR_3D
|
||||
N1 0.0 0.0 0.0
|
||||
N2 1.0 0.0 0.0
|
||||
N3 1.0 1.0 0.0
|
||||
N4 0.0 1.0 0.0
|
||||
N5 0.0 0.0 1.0
|
||||
N6 1.0 0.0 1.0
|
||||
N7 1.0 1.0 1.0
|
||||
N8 0.0 1.0 1.0
|
||||
FINSF
|
||||
""")
|
||||
|
||||
function set_geometry!(element, nodes)
|
||||
element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)]
|
||||
end
|
||||
element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8])
|
||||
set_geometry!(element1, nodes)
|
||||
element1["youngs modulus"] = 9000.0
|
||||
element1["poissons ratio"] = 0.25
|
||||
|
||||
element2 = Quad4([5, 6, 7, 8])
|
||||
set_geometry!(element2, nodes)
|
||||
element2["displacement traction force"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]
|
||||
|
||||
problem = LinearElasticityProblem()
|
||||
push!(problem, element1)
|
||||
push!(problem, element2)
|
||||
|
||||
free_dofs = zeros(Bool, 8, 3)
|
||||
x = 1
|
||||
y = 2
|
||||
z = 3
|
||||
free_dofs[2, x] = true
|
||||
free_dofs[3, [x, y]] = true
|
||||
free_dofs[4, y] = true
|
||||
free_dofs[5, z] = true
|
||||
free_dofs[6, [x, z]] = true
|
||||
free_dofs[7, [x, y, z]] = true
|
||||
free_dofs[8, [y, z]] = true
|
||||
free_dofs = find(vec(free_dofs'))
|
||||
info("free dofs: $free_dofs")
|
||||
|
||||
info("initial force vector")
|
||||
ass = assemble(problem, 0.0)
|
||||
f = full(ass.force_vector)
|
||||
K = full(ass.stiffness_matrix)
|
||||
dump(reshape(f, 3, 8))
|
||||
info("initial stiffness matrix")
|
||||
dump(round(Int, K)[free_dofs, free_dofs])
|
||||
|
||||
u = zeros(3, 8)
|
||||
u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs]
|
||||
|
||||
info("result vector")
|
||||
dump(u)
|
||||
# verified using Code Aster.
|
||||
# 2015-12-12-continuum-elasticity/c3d_linear.*
|
||||
@test isapprox(u[:,7], [2.77777777777778E-02, 2.77777777777778E-02, -1.11111111111111E-01])
|
||||
end
|
||||
#test_continuum_elasticity_with_surface_load()
|
||||
|
||||
|
||||
end
|
||||
Reference in New Issue
Block a user