mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
Improve documentation (#199)
Let's use Literate.jl to automatically generate usage examples. * Automatically generate documentation from other packages (first try to include each package's docs/src/index.md, but if that fails, then use README.md to introduce the package). * Add example how to calculate local element matrices. * Add example how to perform 2d contact analysis.
This commit is contained in:
Jukka Aho
@@ -3,6 +3,8 @@
|
||||
.ipynb_checkpoints
|
||||
docs/build/
|
||||
docs/site/
|
||||
docs/src/examples
|
||||
docs/src/packages
|
||||
*.swp
|
||||
*.lnk
|
||||
*.mess
|
||||
|
||||
+107
-31
@@ -3,24 +3,59 @@
|
||||
|
||||
using Documenter, JuliaFEM
|
||||
|
||||
if !haskey(Pkg.installed(), "Literate")
|
||||
Pkg.add("Literate")
|
||||
end
|
||||
using Literate
|
||||
|
||||
if haskey(ENV, "TRAVIS") && get(ENV, "INSTALL_PYPLOT", "false") == "true"
|
||||
info("inside TRAVIS, installing PyPlot + matplotlib")
|
||||
Pkg.add("PyPlot")
|
||||
run(`pip install matplotlib`)
|
||||
end
|
||||
|
||||
"""
|
||||
copy_docs(pkg_name)
|
||||
|
||||
Copy documentation of some package `pkg_name` to `docs/src/pkg_name`,
|
||||
where `pkg_name` is the name of package. Return true if success, false
|
||||
otherwise.
|
||||
|
||||
If package is undocumented, i.e. directory `docs/src` is missing,
|
||||
but there still exists `README.md`, copy that file to
|
||||
`docs/src/pkg_name/index.md`.
|
||||
"""
|
||||
function copy_docs(pkg_name)
|
||||
src = joinpath(Pkg.dir(pkg_name), "docs", "src")
|
||||
if !isdir(src)
|
||||
warn("Cannot copy documentation of package $pkg_name from $src: ",
|
||||
"No such directory exists. (Is the package in REQUIRE of JuliaFEM?)")
|
||||
return false
|
||||
|
||||
src_dir = Pkg.dir(pkg_name, "docs", "src")
|
||||
pkg_dir = Pkg.dir("JuliaFEM", "docs", "src", "packages")
|
||||
dst_dir = joinpath(pkg_dir, pkg_name)
|
||||
isdir(pkg_dir) || mkpath(pkg_dir)
|
||||
|
||||
# if can find pkg_name/docs/src =>
|
||||
# copy that to docs/src/packages/pkg_name
|
||||
if isdir(src_dir)
|
||||
isdir(dst_dir) || cp(src_dir, dst_dir)
|
||||
info("Copied documentation of package $pkg_name from $src_dir succesfully.")
|
||||
return true
|
||||
end
|
||||
dst = joinpath(Pkg.dir("JuliaFEM"), "docs", "src", pkg_name)
|
||||
cp(src, dst; remove_destination=true)
|
||||
info("Copied documentation of package $pkg_name from $src succesfully.")
|
||||
return true
|
||||
|
||||
# if can find pkg_name/README.md =>
|
||||
# copy that to docs/src/packages/pkg_name/index.md
|
||||
readme_file = Pkg.dir(pkg_name, "README.md")
|
||||
if isfile(readme_file)
|
||||
isdir(dst_dir) || mkpath(dst_dir)
|
||||
dst_file = joinpath(dst_dir, "index.md")
|
||||
isfile(dst_file) || cp(readme_file, dst_file)
|
||||
info("Copied README.md of package $pkg_name from $readme_file succesfully.")
|
||||
return true
|
||||
end
|
||||
|
||||
warn("Cannot copy documentation of package $pkg_name from $src_dir: ",
|
||||
"No such directory exists. (Is the package in REQUIRE of JuliaFEM?)")
|
||||
|
||||
return false
|
||||
|
||||
end
|
||||
|
||||
"""
|
||||
@@ -48,22 +83,19 @@ add_page!(PACKAGES, "Theory" => "MyPackage/theory.md")
|
||||
"""
|
||||
function add_page!(dst, src)
|
||||
file = isa(src, Pair) ? src.second : src
|
||||
if !isfile(joinpath(Pkg.dir("JuliaFEM"), "docs", "src", file))
|
||||
warn("Cannot add page $file: no such file")
|
||||
return false
|
||||
src_dir = Pkg.dir("JuliaFEM", "docs", "src")
|
||||
if isfile(joinpath(src_dir, file))
|
||||
push!(dst, src)
|
||||
return true
|
||||
end
|
||||
push!(dst, src)
|
||||
return true
|
||||
if isfile(joinpath(src_dir, "packages", file))
|
||||
push!(dst, joinpath("packages", src))
|
||||
return true
|
||||
end
|
||||
warn("Cannot add page $file: no such file")
|
||||
return false
|
||||
end
|
||||
|
||||
#=
|
||||
if haskey(ENV, "TRAVIS")
|
||||
println("inside TRAVIS, installing PyPlot + matplotlib")
|
||||
Pkg.add("PyPlot")
|
||||
run(`pip install matplotlib`)
|
||||
end
|
||||
=#
|
||||
|
||||
USER_GUIDE = []
|
||||
|
||||
# Developer's guide is published in FEMBase.jl
|
||||
@@ -82,18 +114,61 @@ if copy_docs("FEMBase")
|
||||
end
|
||||
|
||||
# Let's construct here some description for packages
|
||||
|
||||
PACKAGES = []
|
||||
copy_docs("FEMBase") && add_page!(PACKAGES, "FEMBase/index.md")
|
||||
copy_docs("FEMBasis") && add_page!(PACKAGES, "FEMBasis/index.md")
|
||||
copy_docs("FEMQuad") && add_page!(PACKAGES, "FEMQuad/index.md")
|
||||
copy_docs("FEMSparse") && add_page!(PACKAGES, "FEMSparse/index.md")
|
||||
copy_docs("Materials") && add_page!(PACKAGES, "Materials/index.md")
|
||||
copy_docs("AsterReader") && add_page!(PACKAGES, "AsterReader/index.md")
|
||||
copy_docs("AbaqusReader") && add_page!(PACKAGES, "AbaqusReader/index.md")
|
||||
copy_docs("LinearImplicitDynamics") && add_page!(PACKAGES, "LinearImplicitDynamics/index.md")
|
||||
copy_docs("HeatTransfer") && add_page!(PACKAGES, "HeatTransfer/index.md")
|
||||
copy_docs("PlaneElasticity") && add_page!(PACKAGES, "PlaneElasticity/index.md")
|
||||
copy_docs("FEMBeam") && add_page!(PACKAGES, "FEMBeam/index.md")
|
||||
copy_docs("FEMCoupling") && add_page!(PACKAGES, "FEMCoupling/index.md")
|
||||
copy_docs("FEMTruss") && add_page!(PACKAGES, "FEMTruss/index.md")
|
||||
copy_docs("Mortar2D") && add_page!(PACKAGES, "Mortar2D/index.md")
|
||||
copy_docs("Mortar3D") && add_page!(PACKAGES, "Mortar3D/index.md")
|
||||
copy_docs("MortarContact2D") && add_page!(PACKAGES, "MortarContact2D/index.md")
|
||||
copy_docs("MortarContact2DAD") && add_page!(PACKAGES, "MortarContact2DAD/index.md")
|
||||
copy_docs("OptoMechanics") && add_page!(PACKAGES, "OptoMechanics/index.md")
|
||||
copy_docs("Miniball") && add_page!(PACKAGES, "Miniball/index.md")
|
||||
copy_docs("ModelReduction") && add_page!(PACKAGES, "ModelReduction/index.md")
|
||||
copy_docs("NodeNumbering") && add_page!(PACKAGES, "NodeNumbering/index.md")
|
||||
#copy_docs("Xdmf") && add_page!(PACKAGES, "Xdmf/index.md")
|
||||
copy_docs("UMAT") && add_page!(PACKAGES, "UMAT/index.md")
|
||||
|
||||
if copy_docs("FEMQuad")
|
||||
add_page!(PACKAGES, "FEMQuad/index.md")
|
||||
# Generate examples using Literate.jl
|
||||
|
||||
function license_stripper(s)
|
||||
# return replace(s, r"# This file is a part of JuliaFEM.*?\n\n"sm, "")
|
||||
lines = split(s, '\n')
|
||||
function myfilt(line)
|
||||
contains(line, "# This file is a part of JuliaFEM.") && return false
|
||||
contains(line, "# License is MIT") && return false
|
||||
return true
|
||||
end
|
||||
new_s = join(filter(myfilt, lines), '\n')
|
||||
return new_s
|
||||
end
|
||||
|
||||
if copy_docs("FEMBasis")
|
||||
add_page!(PACKAGES, "FEMBasis/index.md")
|
||||
end
|
||||
|
||||
if copy_docs("HeatTransfer")
|
||||
add_page!(PACKAGES, "HeatTransfer/index.md")
|
||||
EXAMPLES = []
|
||||
docs_src = Pkg.dir("JuliaFEM", "docs", "src")
|
||||
ex_src = Pkg.dir("JuliaFEM", "examples")
|
||||
ex_dst = joinpath(docs_src, "examples")
|
||||
for ex_file in readdir(ex_src)
|
||||
endswith(ex_file, ".jl") || continue
|
||||
startswith(ex_file, "test_") && continue
|
||||
dst_file = Literate.markdown(joinpath(ex_src, ex_file), ex_dst;
|
||||
documenter=true, preprocess=license_stripper)
|
||||
ex_dir = joinpath(ex_src, first(splitext(ex_file)))
|
||||
dst_ex_dir = joinpath(ex_dst, first(splitext(ex_file)))
|
||||
if isdir(ex_dir) && !isdir(dst_ex_dir)
|
||||
info("Copy $ex_dir to $dst_ex_dir")
|
||||
cp(ex_dir, dst_ex_dir)
|
||||
end
|
||||
add_page!(EXAMPLES, relpath(dst_file, docs_src))
|
||||
end
|
||||
|
||||
# API documentation
|
||||
@@ -103,6 +178,7 @@ LIBRARY = ["api.md"]
|
||||
PAGES = []
|
||||
push!(PAGES, "Home" => "index.md")
|
||||
#push!(PAGES, "User's guide" => USER_GUIDE)
|
||||
push!(PAGES, "Examples" => EXAMPLES)
|
||||
push!(PAGES, "Developer's guide" => DEVELOPER_GUIDE)
|
||||
push!(PAGES, "Description of packages" => PACKAGES)
|
||||
#push!("API documentation" => LIBRARY)
|
||||
|
||||
@@ -0,0 +1,221 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# # 2D Hertz contact problem
|
||||
|
||||
# 
|
||||
|
||||
# In the example, a cylinder is pressed agains block with a force of 35 kN.
|
||||
# A similar example can be found from NAFEMS report FENET D3613 (advanced
|
||||
# finite element contact benchmarks).
|
||||
#
|
||||
# Solution for maximum pressure ``p_0`` and contact radius ``a`` is
|
||||
# ```math
|
||||
# p_{0} = \sqrt{\frac{FE}{2\pi R}}, \\
|
||||
# a = \sqrt{\frac{8FR}{\pi E}},
|
||||
# ```
|
||||
# where
|
||||
#
|
||||
# ```math
|
||||
# E = \frac{2E_{1}E_{2}}{E_{2}\left(1-\nu_{1}^{2}\right)+E_{1}\left(1-\nu_{2}^{2}\right)}.
|
||||
# ```
|
||||
#
|
||||
# Substituting values, one gets accurate solution to be ``p_0 = 3585 \;\mathrm{MPa}`` and
|
||||
# ``a = 6.21 \;\mathrm{mm}``.
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Preprocess
|
||||
using JuliaFEM.Postprocess
|
||||
using Logging
|
||||
Logging.configure(level=INFO)
|
||||
add_elements! = JuliaFEM.add_elements!
|
||||
|
||||
# Simulation starts by reading the mesh. Model is constructed and meshed using
|
||||
# SALOME, thus mesh format is .med. Mesh type is quite simple structure,
|
||||
# containing things like `mesh.nodes`, `mesh.elements` and so on. Keep on mind,
|
||||
# that Mesh contains only standard Julia types and we think it as a structure
|
||||
# helping us to construct elements needed in simulation. In principle, we don't
|
||||
# need to use `Mesh` in simulation anyway if we figure some other way to define
|
||||
# the geometry for elements.
|
||||
|
||||
datadir = Pkg.dir("JuliaFEM", "examples", "2d_hertz_contact")
|
||||
meshfile = joinpath(datadir, "hertz_2d_full.med")
|
||||
mesh = aster_read_mesh(meshfile)
|
||||
for (elset_name, element_ids) in mesh.element_sets
|
||||
nel = length(element_ids)
|
||||
println("Element set $elset_name contains $nel elements.")
|
||||
end
|
||||
for (nset_name, node_ids) in mesh.node_sets
|
||||
nno = length(node_ids)
|
||||
println("Node set $nset_name contains $nno nodes.")
|
||||
end
|
||||
nnodes = length(mesh.nodes)
|
||||
println("Total number of nodes in mesh: $nnodes")
|
||||
nelements = length(mesh.elements)
|
||||
println("Total number of elements in mesh: $nelements")
|
||||
|
||||
# Next, define two bodies. Technically, we could have only one problem and add
|
||||
# elements from both bodies to the same problem, but defining two different
|
||||
# problems is recommended for clarity. Plain strain assumption is used.
|
||||
# To make clear what is happening here: we first create a set of elements
|
||||
# (elements are in vector called `upper_elements`), then we define new
|
||||
# problem which type is `Elasticity`, give it some meaningful name (this time
|
||||
# `cylinder`), and last value 2 means that problems does have two degrees of
|
||||
# freedom per node.
|
||||
|
||||
upper_elements = create_elements(mesh, "CYLINDER")
|
||||
update!(upper_elements, "youngs modulus", 70.0e3)
|
||||
update!(upper_elements, "poissons ratio", 0.3)
|
||||
upper = Problem(Elasticity, "cylinder", 2)
|
||||
upper.properties.formulation = :plane_strain
|
||||
add_elements!(upper, upper_elements)
|
||||
|
||||
lower_elements = create_elements(mesh, "BLOCK")
|
||||
update!(lower_elements, "youngs modulus", 210.0e3)
|
||||
update!(lower_elements, "poissons ratio", 0.3)
|
||||
lower = Problem(Elasticity, "block", 2)
|
||||
lower.properties.formulation = :plane_strain
|
||||
add_elements!(lower, lower_elements)
|
||||
|
||||
# Next we define some boundary conditions: creating "boundary" problems goes
|
||||
# in the same way than defining "field" problems, the only difference is that
|
||||
# we add extra argument giving what field are we tring to fix. This time,
|
||||
# we have 2 dofs / node and we fix displacement in direction 2.
|
||||
|
||||
bc_fixed_elements = create_elements(mesh, "FIXED")
|
||||
update!(bc_fixed_elements, "displacement 2", 0.0)
|
||||
bc_fixed = Problem(Dirichlet, "fixed", 2, "displacement")
|
||||
add_elements!(bc_fixed, bc_fixed_elements)
|
||||
|
||||
# Defining symmetry boundary condition goes with the same idea
|
||||
|
||||
bc_sym_23_elements = create_elements(mesh, "SYM23")
|
||||
update!(bc_sym_23_elements, "displacement 1", 0.0)
|
||||
bc_sym_23 = Problem(Dirichlet, "symmetry line 23", 2, "displacement")
|
||||
add_elements!(bc_sym_23, bc_sym_23_elements)
|
||||
|
||||
# Next we define point load. To define that, we first need to find some node
|
||||
# near the top of cylinder, using function `find_nearest_node`. Then we create
|
||||
# a new problem, again of type Elasticity. Like told already, we don't need to
|
||||
# use `Mesh` if we have some other procedure to define the geometry of the
|
||||
# element (and it's connectivity, of course). So we can directly create an
|
||||
# element of type `Poi1`, meaning 1-node point element, update it's geometry
|
||||
# and apply 35.0e3 kN load in negative y-direction:
|
||||
|
||||
nid = find_nearest_node(mesh, [0.0, 100.0])
|
||||
load = Problem(Elasticity, "point load", 2)
|
||||
load.properties.formulation = :plane_strain
|
||||
load.elements = [Element(Poi1, [nid])]
|
||||
update!(load.elements, "geometry", mesh.nodes)
|
||||
update!(load.elements, "displacement traction force 2", -35.0e3)
|
||||
|
||||
# Next, we define another boudary problem, this time the type of problem is
|
||||
# Contact2D, which is a mortar contact formulation for two dimensions.
|
||||
# Elements are added using `add_slave_elements!` and `add_master_elements!`.
|
||||
# Problems, in general, can have some properties defined, like the formulation
|
||||
# in `Elasticity` (we also have `:plane_stress`). For contact, we need to swap
|
||||
# normal direction for meshes created by SALOME because in Code Aster, element
|
||||
# orientation is defined opposite to what is used in ABAQUS, and in JuliaFEM in
|
||||
# general we follow the same conventions what are used in ABAQUS.
|
||||
|
||||
contact = Problem(Contact2D, "contact", 2, "displacement")
|
||||
contact.properties.rotate_normals = true
|
||||
contact_slave_elements = create_elements(mesh, "BLOCK_TO_CYLINDER")
|
||||
contact_master_elements = create_elements(mesh, "CYLINDER_TO_BLOCK")
|
||||
add_master_elements!(contact, contact_master_elements)
|
||||
add_slave_elements!(contact, contact_slave_elements)
|
||||
|
||||
# After all problems are defined, we define some `Analysis`, which can be e.g.
|
||||
# static analysis, dynamic analysis, modal analysis, linear perturbation
|
||||
# analysis and so on. Here, the analysis type is `Nonlinear`, which is nonlinear
|
||||
# quasistatic analysis. In the same manner as we do `add_elements!` to add
|
||||
# elements to `Problem`, we use `add_problems!` to add problems to analysis.
|
||||
# Because we are not restricted to some particular input and output formats,
|
||||
# we "connect" a `ResultsWriter` to our analysis, this time we want to visualize
|
||||
# results using ParaView, thus we write our results to Xdmf format, which uses
|
||||
# well defined standards XML and HDF to store model data.
|
||||
|
||||
step = Analysis(Nonlinear)
|
||||
add_problems!(step, [upper, lower, bc_fixed, bc_sym_23, load, contact])
|
||||
xdmf = Xdmf("2d_hertz_results"; overwrite=true)
|
||||
## todo for 2d
|
||||
## for body in (upper, lower)
|
||||
## push!(body.postprocess_fields, "stress")
|
||||
## end
|
||||
add_results_writer!(step, xdmf)
|
||||
|
||||
# In last part, we run the analysis.
|
||||
|
||||
step()
|
||||
|
||||
# # Results
|
||||
|
||||
# Results are stored in `2d_hertz_results.xmf` and `2d_hertz_results.h5` for
|
||||
# visual inspection. We can also postprocess results programmatically because
|
||||
# we are inside a real scripting / programming environment all the time. For
|
||||
# example, we can integrate the resultant force in normal and tangential direction
|
||||
# in contact surface to validate our result.
|
||||
|
||||
Rn = 0.0
|
||||
Rt = 0.0
|
||||
time = 0.0
|
||||
for sel in contact_slave_elements
|
||||
for ip in get_integration_points(sel)
|
||||
w = ip.weight*sel(ip, time, Val{:detJ})
|
||||
n = sel("normal", ip, time)
|
||||
t = sel("tangent", ip, time)
|
||||
la = sel("lambda", ip, time)
|
||||
Rn += w*dot(n, la)
|
||||
Rt += w*dot(t, la)
|
||||
end
|
||||
end
|
||||
|
||||
println("2d hertz contact resultant forces: Rn = $Rn, Rt = $Rt")
|
||||
|
||||
using Base.Test
|
||||
@test isapprox(Rn, 35.0e3)
|
||||
@test isapprox(Rt, 0.0)
|
||||
|
||||
# Visualization of the results can be done using ParaView:
|
||||
# 
|
||||
|
||||
# For optimization loops, we want to programmatically find, for example, maximum
|
||||
# contact pressure. We can, for example, get all the values in nodes:
|
||||
|
||||
lambda = contact("lambda", time)
|
||||
normal = contact("normal", time)
|
||||
p0 = 0.0
|
||||
p0_acc = 3585.0
|
||||
for (nid, n) in normal
|
||||
lan = dot(n, lambda[nid])
|
||||
println("$nid => $lan")
|
||||
p0 = max(p0, lan)
|
||||
end
|
||||
p0 = round(p0, 2)
|
||||
rtol = round(norm(p0-p0_acc)/max(p0,p0_acc)*100, 2)
|
||||
println("Maximum contact pressure p0 = $p0, p0_acc = $p0_acc, rtol = $rtol %")
|
||||
|
||||
# To get rough approximation where does the contact open, we can find the element
|
||||
# from slave contact surface, where contact pressure is zero in the other node
|
||||
# and something nonzero in the other node.
|
||||
|
||||
a_rad = 0.0
|
||||
for element in contact_slave_elements
|
||||
la1, la2 = element("lambda", time)
|
||||
p1, p2 = norm(la1), norm(la2)
|
||||
a, b = isapprox(p1, 0.0), isapprox(p2, 0.0)
|
||||
if (a && !b) || (b && !a)
|
||||
X1, X2 = element("geometry", time)
|
||||
println("Contact opening element geometry: X1 = $X1, X2 = $X2")
|
||||
println("Contact opening element lambda: la1 = $la1, la2 = $la2")
|
||||
x11, y11 = X1
|
||||
x12, y12 = X2
|
||||
a_rad = 1/2*abs(x11+x12)
|
||||
break
|
||||
end
|
||||
end
|
||||
println("Contact radius: $a_rad")
|
||||
|
||||
# This example briefly described some of the core features of JuliaFEM.
|
||||
|
||||
close(xdmf.hdf) # src
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 25 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 138 KiB |
@@ -0,0 +1,19 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# # Generating local matrices for problems
|
||||
|
||||
using JuliaFEM
|
||||
|
||||
# Plane stress Quad4 element with linear material model:
|
||||
|
||||
X = Dict(1 => [0.0, 0.0], 2 => [2.0, 0.0], 3 => [2.0, 2.0], 4 => [0.0, 2.0])
|
||||
element = Element(Quad4, [1, 2, 3, 4])
|
||||
update!(element, "geometry", X)
|
||||
update!(element, "youngs modulus", 288.0)
|
||||
update!(element, "poissons ratio", 1/3)
|
||||
problem = Problem(Elasticity, "test problem", 2)
|
||||
problem.properties.formulation = :plane_stress
|
||||
add_elements!(problem, [element])
|
||||
assemble!(problem, 0.0)
|
||||
K = round.(full(problem.assembly.K), 5)
|
||||
@@ -1,175 +0,0 @@
|
||||
# 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
|
||||
|
||||
testdir = joinpath(Pkg.dir("JuliaFEM"), "test")
|
||||
datadir = first(splitext(basename(@__FILE__)))
|
||||
|
||||
@testset "hertz contact, full 2d model, linear elements, flat slave surface" begin
|
||||
meshfile = joinpath(testdir, datadir, "hertz_2d_full.med")
|
||||
mesh = aster_read_mesh(meshfile)
|
||||
|
||||
upper = Problem(Elasticity, "CYLINDER", 2)
|
||||
upper.properties.formulation = :plane_strain
|
||||
upper.elements = create_elements(mesh, "CYLINDER")
|
||||
update!(upper, "youngs modulus", 70.0e3)
|
||||
update!(upper, "poissons ratio", 0.3)
|
||||
|
||||
lower = Problem(Elasticity, "BLOCK", 2)
|
||||
lower.properties.formulation = :plane_strain
|
||||
lower.elements = create_elements(mesh, "BLOCK")
|
||||
update!(lower, "youngs modulus", 210.0e3)
|
||||
update!(lower, "poissons ratio", 0.3)
|
||||
|
||||
# support block to ground
|
||||
bc_fixed = Problem(Dirichlet, "fixed", 2, "displacement")
|
||||
bc_fixed.elements = create_elements(mesh, "FIXED")
|
||||
update!(bc_fixed, "displacement 2", 0.0)
|
||||
|
||||
# symmetry line
|
||||
bc_sym_23 = Problem(Dirichlet, "symmetry line 23", 2, "displacement")
|
||||
bc_sym_23.elements = create_elements(mesh, "SYM23")
|
||||
update!(bc_sym_23, "displacement 1", 0.0)
|
||||
|
||||
nid = find_nearest_node(mesh, [0.0, 100.0])
|
||||
#load = Problem(Dirichlet, "load", 2, "displacement")
|
||||
load = Problem(Elasticity, "point load", 2)
|
||||
load.properties.formulation = :plane_strain
|
||||
load.elements = [Element(Poi1, [nid])]
|
||||
#update!(load.elements, "displacement 2", -10.0)
|
||||
update!(load, "displacement traction force 2", -35.0e3)
|
||||
|
||||
contact = Problem(Contact2D, "contact between block and cylinder", 2, "displacement")
|
||||
contact.properties.rotate_normals = true
|
||||
contact_slave_elements = create_elements(mesh, "BLOCK_TO_CYLINDER")
|
||||
contact_master_elements = create_elements(mesh, "CYLINDER_TO_BLOCK")
|
||||
add_master_elements!(contact, contact_master_elements)
|
||||
add_slave_elements!(contact, contact_slave_elements)
|
||||
|
||||
solver = Solver(Nonlinear)
|
||||
push!(solver, upper, lower, bc_fixed, bc_sym_23, load, contact)
|
||||
solver()
|
||||
|
||||
node_ids, la = get_nodal_vector(contact_slave_elements, "lambda", 0.0)
|
||||
node_ids, n = get_nodal_vector(contact_slave_elements, "normal", 0.0)
|
||||
pres = [dot(ni, lai) for (ni, lai) in zip(n, la)]
|
||||
#@test isapprox(maximum(pres), 4060.010799583303)
|
||||
# 12 % error in maximum pressure
|
||||
# integrate pressure in normal and tangential direction
|
||||
Rn = 0.0
|
||||
Rt = 0.0
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
time = 0.0
|
||||
for sel in contact_slave_elements
|
||||
for ip in get_integration_points(sel)
|
||||
w = ip.weight*sel(ip, time, Val{:detJ})
|
||||
n = sel("normal", ip, time)
|
||||
t = Q'*n
|
||||
la = sel("lambda", ip, time)
|
||||
Rn += w*dot(n, la)
|
||||
Rt += w*dot(t, la)
|
||||
end
|
||||
end
|
||||
info("2d hertz: Rn = $Rn, Rt = $Rt")
|
||||
info("2d hertz: maximum pressure pmax = ", maximum(pres))
|
||||
@test isapprox(maximum(pres), 3585.0; rtol = 0.13)
|
||||
# under 0.15 % error in resultant force
|
||||
@test isapprox(Rn, 35.0e3; rtol=0.020)
|
||||
@test isapprox(Rt, 0.0; atol=200.0)
|
||||
end
|
||||
|
||||
function get_model()
|
||||
meshfile = joinpath(testdir, datadir, "block_2d.med")
|
||||
mesh = aster_read_mesh(meshfile)
|
||||
println(mesh.nodes[1])
|
||||
|
||||
upper = Problem(mesh, Elasticity, "UPPER", 2)
|
||||
lower = Problem(mesh, Elasticity, "LOWER", 2)
|
||||
|
||||
for body in [upper, lower]
|
||||
body.properties.formulation = :plane_stress
|
||||
update!(body, "youngs modulus", 288.0)
|
||||
update!(body, "poissons ratio", 1/3)
|
||||
end
|
||||
|
||||
load = Problem(mesh, Elasticity, "UPPER_TOP", 2)
|
||||
load.properties.formulation = :plane_stress
|
||||
update!(load, "displacement traction force 2", -28.8)
|
||||
bc1 = Problem(mesh, Dirichlet, "LOWER_BOTTOM", 2, "displacement")
|
||||
update!(bc1, "displacement 2", 0.0)
|
||||
bc2 = Problem(mesh, Dirichlet, "LOWER_LEFT", 2, "displacement")
|
||||
update!(bc2, "displacement 1", 0.0)
|
||||
bc3 = Problem(mesh, Dirichlet, "UPPER_LEFT", 2, "displacement")
|
||||
update!(bc3, "displacement 1", 0.0)
|
||||
|
||||
interface = Problem(Contact2D, "interface", 2, "displacement")
|
||||
interface.properties.rotate_normals = true
|
||||
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)
|
||||
|
||||
# in LOWER_LEFT we have node belonging also to contact interface
|
||||
# let's remove it from dirichlet bc
|
||||
create_node_set_from_element_set!(mesh, "LOWER_LEFT")
|
||||
nid = find_nearest_node(mesh, [0.0, 0.5]; node_set="LOWER_LEFT")
|
||||
coords = mesh.nodes[nid]
|
||||
info("nearest node to (0.0, 0.5) = $nid, coordinates = $coords")
|
||||
dofs = [2*(nid-1)+1, 2*(nid-1)+2]
|
||||
info("removing nid $nid, dofs $dofs from LOWER_LEFT")
|
||||
push!(bc2.assembly.removed_dofs, dofs...)
|
||||
|
||||
solver = Solver(Nonlinear)
|
||||
#push!(solver, upper, lower, load, bc1, interface)
|
||||
push!(solver, upper, lower, load, bc1, bc2, bc3, interface)
|
||||
return solver
|
||||
end
|
||||
|
||||
@testset "small sliding 2d patch test, linear Seg2 elements, standard basis" begin
|
||||
|
||||
solver = get_model()
|
||||
interface = solver["interface"]
|
||||
solver()
|
||||
|
||||
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)
|
||||
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 0.0)
|
||||
node_ids, lambda = get_nodal_vector(interface.elements, "lambda", 0.0)
|
||||
u2 = [u[2] for u in displacement]
|
||||
f2 = [f[2] for f in lambda]
|
||||
maxabsu2 = maximum(abs.(u2))
|
||||
stdabsu2 = std(abs.(u2))
|
||||
info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2")
|
||||
@test isapprox(stdabsu2, 0.0; atol=1.0e-12)
|
||||
maxabsf2 = maximum(abs.(f2))
|
||||
stdabsf2 = std(abs.(f2))
|
||||
info("max(abs(f2)) = $maxabsf2, std(abs(f2)) = $stdabsf2")
|
||||
@test isapprox(stdabsf2, 0.0; atol=1.0e-12)
|
||||
@test isapprox(mean(abs.(f2)), 28.8; atol=1.0e-12)
|
||||
end
|
||||
|
||||
@testset "small sliding 2d patch test, linear Seg2 elements, dual basis" begin
|
||||
|
||||
solver = get_model()
|
||||
interface = solver["interface"]
|
||||
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)
|
||||
node_ids, lambda = get_nodal_vector(interface.elements, "lambda", 0.0)
|
||||
u2 = [u[2] for u in displacement]
|
||||
f2 = [f[2] for f in lambda]
|
||||
maxabsu2 = maximum(abs.(u2))
|
||||
stdabsu2 = std(abs.(u2))
|
||||
info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2")
|
||||
@test isapprox(stdabsu2, 0.0; atol=1.0e-12)
|
||||
maxabsf2 = maximum(abs.(f2))
|
||||
stdabsf2 = std(abs.(f2))
|
||||
info("max(abs(f2)) = $maxabsf2, std(abs(f2)) = $stdabsf2")
|
||||
@test isapprox(stdabsf2, 0.0; atol=1.0e-12)
|
||||
@test isapprox(mean(abs.(f2)), 28.8; atol=1.0e-12)
|
||||
end
|
||||
Binary file not shown.
@@ -1,86 +0,0 @@
|
||||
# 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
|
||||
|
||||
pkgdir = Pkg.dir("JuliaFEM")
|
||||
datadir = joinpath(pkgdir, "test", first(splitext(basename(@__FILE__))))
|
||||
|
||||
# from fenet d3613 advanced finite element contact benchmarks
|
||||
# a = 6.21 mm, pmax = 3585 MPa
|
||||
# this is a very sparse mesh and for that reason pmax is not very accurate
|
||||
# (only 6 elements in -20 .. 20 mm contact zone, 3 elements in contact
|
||||
|
||||
meshfile = joinpath(datadir, "hertz_2d_full.med")
|
||||
mesh = aster_read_mesh(meshfile)
|
||||
|
||||
upper = Problem(Elasticity, "CYLINDER", 2)
|
||||
upper.properties.formulation = :plane_strain
|
||||
upper.elements = create_elements(mesh, "CYLINDER")
|
||||
update!(upper, "youngs modulus", 70.0e3)
|
||||
update!(upper, "poissons ratio", 0.3)
|
||||
|
||||
lower = Problem(Elasticity, "BLOCK", 2)
|
||||
lower.properties.formulation = :plane_strain
|
||||
lower.elements = create_elements(mesh, "BLOCK")
|
||||
update!(lower, "youngs modulus", 210.0e3)
|
||||
update!(lower, "poissons ratio", 0.3)
|
||||
|
||||
# support block to ground
|
||||
bc_fixed = Problem(Dirichlet, "fixed", 2, "displacement")
|
||||
bc_fixed.elements = create_elements(mesh, "FIXED")
|
||||
update!(bc_fixed, "displacement 2", 0.0)
|
||||
|
||||
# symmetry line
|
||||
bc_sym_23 = Problem(Dirichlet, "symmetry line 23", 2, "displacement")
|
||||
bc_sym_23.elements = create_elements(mesh, "SYM23")
|
||||
update!(bc_sym_23, "displacement 1", 0.0)
|
||||
|
||||
nid = find_nearest_node(mesh, [0.0, 100.0])
|
||||
#load = Problem(Dirichlet, "load", 2, "displacement")
|
||||
load = Problem(Elasticity, "point load", 2)
|
||||
load.properties.formulation = :plane_strain
|
||||
load.elements = [Element(Poi1, [nid])]
|
||||
#update!(load.elements, "displacement 2", -10.0)
|
||||
update!(load, "displacement traction force 2", -35.0e3)
|
||||
|
||||
contact = Problem(Contact2D, "contact between block and cylinder", 2, "displacement")
|
||||
contact.properties.rotate_normals = true
|
||||
contact_slave_elements = create_elements(mesh, "CYLINDER_TO_BLOCK")
|
||||
contact_master_elements = create_elements(mesh, "BLOCK_TO_CYLINDER")
|
||||
add_master_elements!(contact, contact_master_elements)
|
||||
add_slave_elements!(contact, contact_slave_elements)
|
||||
|
||||
solver = Solver(Nonlinear)
|
||||
push!(solver, upper, lower, bc_fixed, bc_sym_23, load, contact)
|
||||
solver()
|
||||
|
||||
node_ids, la = get_nodal_vector(contact_slave_elements, "lambda", 0.0)
|
||||
node_ids, n = get_nodal_vector(contact_slave_elements, "normal", 0.0)
|
||||
pres = [dot(ni, lai) for (ni, lai) in zip(n, la)]
|
||||
#@test isapprox(maximum(pres), 4060.010799583303)
|
||||
# 12 % error in maximum pressure
|
||||
# integrate pressure in normal and tangential direction
|
||||
Rn = 0.0
|
||||
Rt = 0.0
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
time = 0.0
|
||||
for sel in contact_slave_elements
|
||||
for ip in get_integration_points(sel)
|
||||
w = ip.weight*sel(ip, time, Val{:detJ})
|
||||
n = sel("normal", ip, time)
|
||||
t = Q'*n
|
||||
la = sel("lambda", ip, time)
|
||||
Rn += w*dot(n, la)
|
||||
Rt += w*dot(t, la)
|
||||
end
|
||||
end
|
||||
info("2d hertz: Rn = $Rn, Rt = $Rt")
|
||||
info("2d hertz: maximum pressure pmax = ", maximum(pres))
|
||||
@test isapprox(maximum(pres), 3585.0; rtol = 0.13)
|
||||
# under 0.15 % error in resultant force
|
||||
@test isapprox(Rn, 35.0e3; rtol=0.020)
|
||||
@test isapprox(Rt, 0.0; atol=200.0)
|
||||
Binary file not shown.
Reference in New Issue
Block a user