mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
Use FEMBeam.jl (#200)
Use FEMBeam.jl to solve beam problems. Added an example, where natural frequencies of frequencies of 3d frame structure are calculated. Some minor modifications to Modal analysis is done to make Xdmf writing of 6 dof nodes work.
This commit is contained in:
Jukka Aho
@@ -14,3 +14,4 @@ Reexport
|
||||
HeatTransfer
|
||||
MortarContact2D
|
||||
MortarContact2DAD
|
||||
FEMBeam
|
||||
|
||||
@@ -0,0 +1,95 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# # Natural frequency analysis of 3d frame structure
|
||||
|
||||
# For general information about Euler-Bernoulli beam theory, see
|
||||
# [this](https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory)
|
||||
# wikipedia page.
|
||||
|
||||
# The model is a 3d frame, shown in picture.
|
||||
|
||||
# 
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Preprocess
|
||||
using FEMBase.Test
|
||||
using Logging
|
||||
Logging.configure(level=INFO)
|
||||
add_elements! = JuliaFEM.add_elements!
|
||||
|
||||
# Reading mesh
|
||||
|
||||
datadir = Pkg.dir("JuliaFEM", "examples", "3d_frame")
|
||||
mesh = aster_read_mesh(joinpath(datadir, "model.med"))
|
||||
println("Number of nodes in a model: ", length(mesh.nodes))
|
||||
|
||||
# Create beam elements. For 3d model, we need to define at least
|
||||
# [Young's modulus](https://en.wikipedia.org/wiki/Young%27s_modulus),
|
||||
# [shear modulus](https://en.wikipedia.org/wiki/Shear_modulus),
|
||||
# [density](https://en.wikipedia.org/wiki/Density)
|
||||
# cross-section area, moment of inertia in local coordinate
|
||||
# system and polar moment of inertia.
|
||||
|
||||
beam_elements = create_elements(mesh, "FRAME")
|
||||
info("Number of elements: ", length(beam_elements))
|
||||
update!(beam_elements, "youngs modulus", 210.0e6)
|
||||
update!(beam_elements, "shear modulus", 84.0e6)
|
||||
update!(beam_elements, "density", 7850.0e-3)
|
||||
update!(beam_elements, "cross-section area", 20.0e-2)
|
||||
update!(beam_elements, "torsional moment of inertia 1", 10.0e-5)
|
||||
update!(beam_elements, "torsional moment of inertia 2", 10.0e-5)
|
||||
update!(beam_elements, "polar moment of inertia", 30.0e-5)
|
||||
|
||||
# The direction of beam is defined in same way than in ABAQUS.
|
||||
# That is, we have a tangent direction and one normal direction.
|
||||
# The third direction is then cross product of tangent and normal.
|
||||
# Because the second area moment is same in both directions, we can
|
||||
# choose normal direction freely.
|
||||
|
||||
for element in beam_elements
|
||||
X1, X2 = element("geometry", 0.0)
|
||||
t = (X2-X1)/norm(X2-X1)
|
||||
I = eye(3)
|
||||
k = indmax([norm(cross(t, I[:,k])) for k in 1:3])
|
||||
n = cross(t, I[:,k])/norm(cross(t, I[:,k]))
|
||||
update!(element, "normal", n)
|
||||
end
|
||||
|
||||
# Create boundary conditions: fix all degrees of freedom for nodes in
|
||||
# a set FIXED. Here we first create elements of type `Poi1` for each
|
||||
# node j in set FIXED, update geometry field and then create new fields
|
||||
# `fixed displacmeent 1`, `fixed displacement 2`, and so on, where the
|
||||
# displacement / rotation is prescribed.
|
||||
|
||||
bc_elements = [Element(Poi1, [j]) for j in mesh.node_sets[:FIXED]]
|
||||
update!(bc_elements, "geometry", mesh.nodes)
|
||||
for i=1:3
|
||||
update!(bc_elements, "fixed displacement $i", 0.0)
|
||||
update!(bc_elements, "fixed rotation $i", 0.0)
|
||||
end
|
||||
|
||||
# Create a problem, containing beam elements and boundary conditions:
|
||||
|
||||
frame = Problem(Beam, "3d frame", 6)
|
||||
add_elements!(frame, beam_elements)
|
||||
add_elements!(frame, bc_elements)
|
||||
|
||||
# Perform modal analysis
|
||||
|
||||
step = Analysis(Modal)
|
||||
xdmf = Xdmf(joinpath(datadir, "3d_frame_results"); overwrite=true)
|
||||
add_results_writer!(step, xdmf)
|
||||
add_problems!(step, [frame])
|
||||
run!(step)
|
||||
close(xdmf.hdf)
|
||||
|
||||
# Each `Analysis` can have properties, e.g. time, maximum number of iterations,
|
||||
# convergence tolerance and so on. Eigenvalues of calculation are stored as a
|
||||
# properties of analysis:
|
||||
|
||||
freqs = sqrt.(step.properties.eigvals) / (2*pi)
|
||||
println("Natural frequencies [Hz]: $(round.(freqs, 2))")
|
||||
|
||||
# [](https://www.youtube.com/watch?v=GzktCqeASmo)
|
||||
|
||||
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 14 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 84 KiB |
@@ -44,6 +44,9 @@ export Dirichlet
|
||||
|
||||
export assemble!, postprocess!
|
||||
|
||||
# Structural elements: beams
|
||||
@reexport using FEMBeam
|
||||
|
||||
### Mortar methods ###
|
||||
|
||||
@reexport using MortarContact2D
|
||||
|
||||
+108
-145
@@ -32,125 +32,72 @@ function Modal(nev=10, which=:SM)
|
||||
false, [], true, true, false, false, 0.0)
|
||||
end
|
||||
|
||||
""" Eliminate Dirichlet boundary condition from matrices K, M. """
|
||||
function FEMBase.eliminate_boundary_conditions!(K_red::SparseMatrixCSC,
|
||||
M_red::SparseMatrixCSC,
|
||||
problem::Problem{Dirichlet}, ndim::Int)
|
||||
K = sparse(problem.assembly.K, ndim, ndim)
|
||||
C1 = sparse(problem.assembly.C1, ndim, ndim)
|
||||
C2 = sparse(problem.assembly.C2, ndim, ndim)
|
||||
D = sparse(problem.assembly.D, ndim, ndim)
|
||||
f = sparse(problem.assembly.f, ndim, 1)
|
||||
g = sparse(problem.assembly.g, ndim, 1)
|
||||
Kg = sparse(problem.assembly.Kg, ndim, ndim)
|
||||
fg = sparse(problem.assembly.fg, ndim, ndim)
|
||||
# only homogenenous boundary condition u=0 is implemented at the moment.
|
||||
@assert nnz(K) == 0
|
||||
@assert nnz(D) == 0
|
||||
@assert nnz(Kg) == 0
|
||||
@assert nnz(fg) == 0
|
||||
@assert nnz(f) == 0
|
||||
@assert nnz(g) == 0
|
||||
"""
|
||||
calc_projection(problem, ndim)
|
||||
|
||||
A helper function to calculate P = D^-1*M
|
||||
"""
|
||||
function calc_projection(problem::T) where
|
||||
{T<:Union{Problem{Mortar}, Problem{Mortar2D}}}
|
||||
|
||||
C1 = sparse(problem.assembly.C1)
|
||||
C2 = sparse(problem.assembly.C2)
|
||||
@assert C1 == C2
|
||||
@assert isdiag(C1)
|
||||
fixed_dofs = get_nonzero_rows(C1)
|
||||
P = ones(ndim)
|
||||
P[fixed_dofs] = 0.0
|
||||
Q = spdiagm(P)
|
||||
K_red[:,:] = Q' * K_red * Q
|
||||
M_red[:,:] = Q' * M_red * Q
|
||||
return
|
||||
end
|
||||
|
||||
""" Given data vector, return slave displacements. """
|
||||
function calc_projection(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
|
||||
|
||||
@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)
|
||||
s = get_nonzero_rows(C2)
|
||||
m = setdiff(get_nonzero_columns(C2), s)
|
||||
|
||||
# Construct matrix P = D^-1*M
|
||||
D_ = C2[S,S]
|
||||
M_ = -C2[S,M]
|
||||
P = nothing
|
||||
D = C2[s,s]
|
||||
M = -C2[s,m]
|
||||
|
||||
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_
|
||||
if !isdiag(D)
|
||||
warn("Mortar matrix D is not diagonal. This might take a long time.")
|
||||
P = ldltfact(1/2*(D + D')) \ M
|
||||
else
|
||||
P = D_ \ M_
|
||||
P = D \ M
|
||||
end
|
||||
info("Matrix P ready.")
|
||||
|
||||
return S, M, P
|
||||
return s, m, P
|
||||
end
|
||||
|
||||
function FEMBase.eliminate_boundary_conditions!(problem::P, K, M, f) where {P}
|
||||
isempty(problem.assembly.C2) && return nothing
|
||||
C1 = sparse(problem.assembly.C1)
|
||||
C2 = sparse(problem.assembly.C2)
|
||||
C1 == C2 || error("Cannot eliminate boundary condition $P: C1 != C2.")
|
||||
isdiag(C1) || error("Cannot eliminate boundary condition $P: C is not diagonal")
|
||||
info("Eliminating boundary condition $(problem.name) from global system.")
|
||||
fixed_dofs = get_nonzero_rows(C1)
|
||||
K[fixed_dofs,:] = K[:,fixed_dofs] = 0.0
|
||||
M[fixed_dofs,:] = M[:,fixed_dofs] = 0.0
|
||||
dropzeros!(K)
|
||||
dropzeros!(M)
|
||||
return nothing
|
||||
end
|
||||
|
||||
""" Eliminate mesh tie constraints from matrices K, M. """
|
||||
function FEMBase.eliminate_boundary_conditions!(K_red::SparseMatrixCSC,
|
||||
M_red::SparseMatrixCSC,
|
||||
problem::Union{Problem{Mortar}, Problem{Mortar2D}},
|
||||
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
|
||||
|
||||
@assert C1 == C2
|
||||
#@assert problem.properties.dual_basis == true
|
||||
@assert problem.properties.adjust == false
|
||||
"""
|
||||
eliminate_boundary_conditions!(problem, K, M, f)
|
||||
|
||||
Eliminate Mortar boundary condition from matrices K, M and force vector f.
|
||||
"""
|
||||
function FEMBase.eliminate_boundary_conditions!(problem::T, K, M, f) where
|
||||
{T <: Union{Problem{Mortar}, Problem{Mortar2D}}}
|
||||
info("Eliminating mesh tie constraint $(problem.name) using static condensation")
|
||||
|
||||
S = get_nonzero_rows(C2)
|
||||
M = setdiff(get_nonzero_columns(C2), S)
|
||||
info("# slave dofs = $(length(S)), # master dofs = $(length(M))")
|
||||
|
||||
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
|
||||
s, m, P = calc_projection(problem)
|
||||
ndim = size(K, 1)
|
||||
Id = ones(ndim)
|
||||
Id[s] = 0.0
|
||||
Q = spdiagm(Id)
|
||||
Q[M,S] += P'
|
||||
|
||||
K_red[:,:] = Q*K_red*Q'
|
||||
K_red[S,:] = 0.0
|
||||
K_red[:,S] = 0.0
|
||||
|
||||
M_red[:,:] = Q*M_red*Q'
|
||||
M_red[S,:] = 0.0
|
||||
M_red[:,S] = 0.0
|
||||
|
||||
return true
|
||||
Q[s,m] += P
|
||||
K[:,:] = Q'*K*Q
|
||||
M[:,:] = Q'*M*Q
|
||||
return nothing
|
||||
end
|
||||
|
||||
function solve!(solver::Solver{Modal}, time::Float64)
|
||||
function FEMBase.run!(solver::Solver{Modal})
|
||||
time = solver.properties.time
|
||||
problems = get_problems(solver)
|
||||
properties = solver.properties
|
||||
info(repeat("-", 80))
|
||||
@@ -169,26 +116,17 @@ function solve!(solver::Solver{Modal}, time::Float64)
|
||||
dim = size(K, 1)
|
||||
ndofs = size(K, 1)
|
||||
|
||||
nboundary_problems = length(get_boundary_problems(solver))
|
||||
|
||||
K_red = K
|
||||
M_red = M
|
||||
|
||||
@timeit "eliminate boundary conditions" begin
|
||||
if length(properties.P) > 0
|
||||
info("Using custom P to make transform K_red = P'*K*P and M_red = P'*M*P")
|
||||
for P in properties.P
|
||||
K_red = P'*K_red*P
|
||||
M_red = P'*M_red*P
|
||||
end
|
||||
elseif nboundary_problems != 0
|
||||
info("Eliminate boundary conditions from system.")
|
||||
for boundary_problem in get_boundary_problems(solver)
|
||||
eliminate_boundary_conditions!(K_red, M_red, boundary_problem, dim)
|
||||
end
|
||||
else
|
||||
info("No boundary Dirichlet boundary conditions found for system.")
|
||||
end
|
||||
for P in properties.P
|
||||
info("Using P to make transformation K_red = P'*K*P and M_red = P'*M*P")
|
||||
K_red[:,:] = P'*K_red*P
|
||||
M_red[:,:] = P'*M_red*P
|
||||
end
|
||||
|
||||
for problem in get_problems(solver)
|
||||
eliminate_boundary_conditions!(problem, K_red, M_red, f)
|
||||
end
|
||||
|
||||
# free up some memory before solution
|
||||
@@ -234,7 +172,6 @@ function solve!(solver::Solver{Modal}, time::Float64)
|
||||
M_red = full(M_red)
|
||||
end
|
||||
|
||||
|
||||
om2 = nothing
|
||||
X = nothing
|
||||
passed = false
|
||||
@@ -245,34 +182,39 @@ function solve!(solver::Solver{Modal}, time::Float64)
|
||||
end
|
||||
passed = true
|
||||
catch
|
||||
info("failed to calculate eigenvalues for problem. Maybe stiffness matrix is not positive definite, checking...")
|
||||
info("is K symmetric? ", issymmetric(K_red))
|
||||
info("is M symmetric? ", issymmetric(M_red))
|
||||
info("is K positive definite? ", isposdef(K_red))
|
||||
info("is M positive definite? ", isposdef(M_red))
|
||||
info("Probably the reason is that stiffness matrix is not positive definite and Cholesky factorization is failing.")
|
||||
info("To work around this problem, use arguments `sigma = <some small value>` when calling solver, i.e.")
|
||||
info("solver(; sigma=1.0e-9")
|
||||
info("Be aware that using sigma shifts eigenvalues up and a bit different results can be expected.")
|
||||
if size(K_red, 1) < 2000
|
||||
info("stiffness matrix is small, using dense eigenvalue solver to check eigenvalues ...")
|
||||
om2 = eigvals(full(K_red))
|
||||
info("squared eigenvalues om2 = $om2")
|
||||
end
|
||||
info("Failed to calculate eigenvalues for problem.")
|
||||
b1 = issymmetric(K_red)
|
||||
b2 = issymmetric(M_red)
|
||||
b3 = isposdef(K_red)
|
||||
b4 = isposdef(M_red)
|
||||
info("Is K symmetric? $b1")
|
||||
info("Is M symmetric? $b2")
|
||||
info("Is K positive definite? $b3")
|
||||
info("Is M positive definite? $b4")
|
||||
if properties.sigma != 0.0
|
||||
info("sigma is manually set and did not work, giving up, try increase sigma.")
|
||||
if !b3
|
||||
info("Stiffness matrix is not positive definite and Cholesky ",
|
||||
"factorization is failing. Model is not supported enough ",
|
||||
"with boundary conditions. To work around this problem, ",
|
||||
"use `problem.properties.sigma = <some small value>` ",
|
||||
"To add artificial stiffness to model. (Or add boundary ",
|
||||
"conditions.)")
|
||||
end
|
||||
rethrow()
|
||||
end
|
||||
end
|
||||
|
||||
if !passed
|
||||
sigma = 1.0e-9
|
||||
info("Calculation of eigenvalues failed, trying again using sigma value sigma=$sigma")
|
||||
sigma = problem.properties.sigma = 1.0e-9
|
||||
info("Calculation of eigenvalues failed. Stiffness matrix is not ",
|
||||
"positive definite and Cholesky factorization is failing. Trying ",
|
||||
"again by adjusting problem.properties.sigma to $sigma.")
|
||||
try
|
||||
om2, X = eigs(K_red + sigma*I, M_red; nev=props.nev, which=props.which)
|
||||
passed = true
|
||||
catch
|
||||
info("Failed to calculate eigenvalues with sigma=$sigma, manually set sigma to something larger and try again.")
|
||||
info("Failed to calculate eigenvalues with sigma value $sigma. ",
|
||||
"Manually set sigma to something larger and try again.")
|
||||
rethrow()
|
||||
end
|
||||
end
|
||||
@@ -287,11 +229,11 @@ function solve!(solver::Solver{Modal}, time::Float64)
|
||||
props.eigvecs[nz,i] = X[:,i]
|
||||
for problem in get_boundary_problems(solver)
|
||||
isa(problem, Problem{Mortar}) || continue
|
||||
S, M, P = calc_projection(problem, dim)
|
||||
s, m, P = calc_projection(problem)
|
||||
# FIXME: store projection to boundary problem, i.e.
|
||||
# update!(problem, "master-slave projection", time => P)
|
||||
# us = P*um
|
||||
props.eigvecs[S,i] = P*props.eigvecs[M,i]
|
||||
props.eigvecs[s,i] = P*props.eigvecs[m,i]
|
||||
end
|
||||
end
|
||||
|
||||
@@ -321,7 +263,11 @@ function update_xdmf!(solver::Solver{Modal})
|
||||
node_ids = keys(X_)
|
||||
@timeit "create node permutation" P = Dict(j=>i for (i, j) in enumerate(node_ids))
|
||||
nnodes = length(X_)
|
||||
ndofs = round(Int, size(solver.properties.eigvecs, 1)/nnodes)
|
||||
ndim = length(X_[first(node_ids)])
|
||||
info("Number of nodes: $nnodes. ",
|
||||
"Number of dofs/node: $ndofs. ",
|
||||
"Dimension of geometry: $ndim.")
|
||||
@timeit "create ncoords array" begin
|
||||
X = zeros(ndim, nnodes)
|
||||
for j in node_ids
|
||||
@@ -410,15 +356,26 @@ function update_xdmf!(solver::Solver{Modal})
|
||||
end
|
||||
|
||||
@timeit "store eigenmode" begin
|
||||
mode_ = solver.properties.eigvecs[:,j]
|
||||
mode_ = reshape(mode_, ndim, round(Int, length(mode_)/ndim))
|
||||
mode_ = reshape(solver.properties.eigvecs[:,j], ndofs, nnodes)
|
||||
@timeit "create mode array" begin
|
||||
mode = zeros(ndim, nnodes)
|
||||
for i in node_ids
|
||||
mode[:,P[i]] = mode_[:,i]
|
||||
mode = zeros(ndofs, nnodes)
|
||||
for nid in node_ids
|
||||
mode[:,P[nid]] = mode_[:,nid]
|
||||
end
|
||||
end
|
||||
field_type = ndim == 1 ? "Scalar" : "Vector"
|
||||
if ndofs == 1
|
||||
field_type = "Scalar"
|
||||
elseif ndofs == 2 # extend to 3d
|
||||
field_type = "Vector"
|
||||
mode = vcat(mode, zeros(1, nnodes))
|
||||
elseif ndofs == 3
|
||||
field_type = "Vector"
|
||||
elseif ndofs == 6 # has rotation dofs, drop them
|
||||
field_type = "Vector"
|
||||
mode = mode[1:3, :]
|
||||
else
|
||||
error("Number of dofs / node = $ndofs, I don't know how to store results to Xdmf!")
|
||||
end
|
||||
field_center = "Node"
|
||||
attribute = new_child(frame, "Attribute")
|
||||
set_attribute(attribute, "Name", unknown_field_name)
|
||||
@@ -434,3 +391,9 @@ function update_xdmf!(solver::Solver{Modal})
|
||||
|
||||
save!(xdmf)
|
||||
end
|
||||
|
||||
function solve!(solver::Solver{Modal}, time::Float64)
|
||||
info("solve!(analysis, time) is deprecated. Use run!(analysis) instead.")
|
||||
solver.properties.time = time
|
||||
run!(solver)
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user