mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-10 13:17:42 +00:00
calculate stress, interpolate stress to nodes using lsq fitting
This commit is contained in:
Jukka Aho
+6
-6
@@ -35,11 +35,11 @@ include("integrate.jl") # default integration points for elements
|
||||
export get_integration_points
|
||||
|
||||
include("sparse.jl")
|
||||
export add!
|
||||
export add!, SparseMatrixCOO, get_nonzero_rows
|
||||
|
||||
include("problems.jl") # common problem routines
|
||||
export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
|
||||
get_unknown_field_dimension, get_gdofs
|
||||
get_unknown_field_dimension, get_gdofs, Assembly
|
||||
|
||||
include("elasticity.jl") # elasticity equations
|
||||
export Elasticity
|
||||
@@ -83,10 +83,10 @@ export aster_create_elements, parse_aster_med_file
|
||||
end
|
||||
|
||||
module Postprocess
|
||||
include("xdmf.jl")
|
||||
export xdmf_new_temporal_collection, xdmf_new_grid,
|
||||
xdmf_new_mesh!, xdmf_new_nodal_field!,
|
||||
xdmf_save_model, xdmf_new_model, xdmf_dump
|
||||
include("postprocess_utils.jl")
|
||||
export calc_nodal_values!
|
||||
include("postprocess_xdmf.jl")
|
||||
export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save!
|
||||
end
|
||||
|
||||
""" JuliaFEM testing routines. """
|
||||
|
||||
+28
-7
@@ -12,7 +12,7 @@ type Elasticity <: FieldProblem
|
||||
end
|
||||
function Elasticity()
|
||||
# formulations: plane_stress, plane_strain, continuum
|
||||
return Elasticity(:continuum, true)
|
||||
return Elasticity(:continuum, false)
|
||||
end
|
||||
|
||||
function get_unknown_field_name(problem::Problem{Elasticity})
|
||||
@@ -93,7 +93,8 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
|
||||
cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]]
|
||||
|
||||
update!(ip, "strain", time => strain_vec)
|
||||
update!(ip, "stress", time => cauchy_stress)
|
||||
update!(ip, "cauchy stress", time => cauchy_stress)
|
||||
update!(ip, "pk2 stress", time => stress_vec)
|
||||
|
||||
# add contributions: material and geometric stiffness + internal forces
|
||||
fill!(BL, 0.0)
|
||||
@@ -225,11 +226,31 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
|
||||
0.0 0.0 0.0 0.0 0.5-nu 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.5-nu]
|
||||
|
||||
# PK2 stress tensor in voigt notation
|
||||
# strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; 2*strain[2,3]; 2*strain[1,3]; 2*strain[1,2]]
|
||||
# order 11, 22, 33, 12, 23, 13 is in many text books ..?
|
||||
strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; 2*strain[1,2]; 2*strain[2,3]; 2*strain[1,3]]
|
||||
stress_vec = D*strain_vec
|
||||
# calculate stress
|
||||
strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]]
|
||||
stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec)
|
||||
|
||||
stress = [
|
||||
stress_vec[1] stress_vec[4] stress_vec[6]
|
||||
stress_vec[4] stress_vec[2] stress_vec[5]
|
||||
stress_vec[6] stress_vec[5] stress_vec[3]]
|
||||
cauchy_stress = F'*stress*F/det(F)
|
||||
cauchy_stress_vec = [
|
||||
cauchy_stress[1,1];
|
||||
cauchy_stress[2,2];
|
||||
cauchy_stress[3,3];
|
||||
cauchy_stress[1,2];
|
||||
cauchy_stress[2,3];
|
||||
cauchy_stress[1,3]]
|
||||
|
||||
s = cauchy_stress - 1.0/3.0*trace(cauchy_stress)*I
|
||||
J2 = 1/2*trace(s*s')
|
||||
|
||||
# update values to integration point
|
||||
update!(ip, "strain", time => strain_vec)
|
||||
update!(ip, "stress", time => stress_vec)
|
||||
update!(ip, "cauchy stress", time => cauchy_stress_vec)
|
||||
update!(ip, "von mises stress", time => sqrt(3.0*J2))
|
||||
|
||||
# add contributions: material and geometric stiffness + internal forces
|
||||
fill!(BL, 0.0)
|
||||
|
||||
@@ -0,0 +1,162 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using LightXML
|
||||
using JuliaFEM
|
||||
|
||||
# element codes: http://www.paraview.org/pipermail/paraview/2013-July/028859.html
|
||||
# > from ./VTK/ThirdParty/xdmf2/vtkxdmf2/libsrc/XdmfTopology.h
|
||||
# >
|
||||
# > // Topologies
|
||||
# > #define XDMF_NOTOPOLOGY 0x0
|
||||
# > #define XDMF_POLYVERTEX 0x1
|
||||
# > #define XDMF_POLYLINE 0x2
|
||||
# > #define XDMF_POLYGON 0x3
|
||||
# > #define XDMF_TRI 0x4
|
||||
# > #define XDMF_QUAD 0x5
|
||||
# > #define XDMF_TET 0x6
|
||||
# > #define XDMF_PYRAMID 0x7
|
||||
# > #define XDMF_WEDGE 0x8
|
||||
# > #define XDMF_HEX 0x9
|
||||
# > #define XDMF_EDGE_3 0x0022
|
||||
# > #define XDMF_TRI_6 0x0024
|
||||
# > #define XDMF_QUAD_8 0x0025
|
||||
# > #define XDMF_QUAD_9 0x0023
|
||||
# > #define XDMF_TET_10 0x0026
|
||||
# > #define XDMF_PYRAMID_13 0x0027
|
||||
# > #define XDMF_WEDGE_15 0x0028
|
||||
# > #define XDMF_WEDGE_18 0x0029
|
||||
# > #define XDMF_HEX_20 0x0030
|
||||
# > #define XDMF_HEX_24 0x0031
|
||||
# > #define XDMF_HEX_27 0x0032
|
||||
# > #define XDMF_MIXED 0x0070
|
||||
# > #define XDMF_2DSMESH 0x0100
|
||||
# > #define XDMF_2DRECTMESH 0x0101
|
||||
# > #define XDMF_2DCORECTMESH 0x0102
|
||||
# > #define XDMF_3DSMESH 0x1100
|
||||
# > #define XDMF_3DRECTMESH 0x1101
|
||||
# > #define XDMF_3DCORECTMESH 0x1102
|
||||
|
||||
get_xdmf_element_code(element::Element{Tri3}) = 0x0004
|
||||
get_xdmf_element_code(element::Element{Quad4}) = 0x0005
|
||||
get_xdmf_element_code(element::Element{Tet4}) = 0x0006
|
||||
get_xdmf_element_code(element::Element{Hex8}) = 0x0009
|
||||
get_xdmf_element_code(element::Element{Tet10}) = 0x0026
|
||||
|
||||
type XDMF
|
||||
dimension :: Int
|
||||
use_hdf :: Bool
|
||||
xdoc :: XMLDocument
|
||||
domain :: XMLElement
|
||||
temporal_collection :: XMLElement
|
||||
current_grid
|
||||
permutation :: Vector{Int}
|
||||
end
|
||||
|
||||
function XDMF()
|
||||
xdoc = XMLDocument()
|
||||
xroot = create_root(xdoc, "Xdmf")
|
||||
set_attribute(xroot, "xmlns:xi", "http://www.w3.org/2001/XInclude")
|
||||
set_attribute(xroot, "Version", "2.1")
|
||||
domain = new_child(xroot, "Domain")
|
||||
temporal_collection = new_child(domain, "Grid")
|
||||
set_attribute(temporal_collection, "CollectionType", "Temporal")
|
||||
set_attribute(temporal_collection, "GridType", "Collection")
|
||||
set_attribute(temporal_collection, "Name", "Collection")
|
||||
return XDMF(3, false, xdoc, domain, temporal_collection, Union{}, [])
|
||||
end
|
||||
|
||||
function xdmf_new_result!(xdmf::XDMF, elements, time)
|
||||
grid = new_child(xdmf.temporal_collection, "Grid")
|
||||
set_attribute(grid, "Name", "Grid")
|
||||
time_ = new_child(grid, "Time")
|
||||
set_attribute(time_, "Value", time)
|
||||
xdmf.current_grid = grid
|
||||
|
||||
# 1. calculate permutation
|
||||
nids = Set()
|
||||
X = Dict{Int64, Vector{Float64}}()
|
||||
for element in elements
|
||||
conn = get_connectivity(element)
|
||||
push!(nids, conn...)
|
||||
X_el = element["geometry"](time)
|
||||
for (i, c) in enumerate(conn)
|
||||
X[c] = X_el[i]
|
||||
end
|
||||
end
|
||||
xdmf.permutation = sort(collect(nids))
|
||||
iperm = Dict{Int64, Int64}()
|
||||
for (i, j) in enumerate(xdmf.permutation)
|
||||
iperm[j] = i
|
||||
end
|
||||
|
||||
# 2. write nodes
|
||||
geometry = new_child(grid, "Geometry")
|
||||
set_attribute(geometry, "Type", xdmf.dimension == 3 ? "XYZ" : "XY")
|
||||
dataitem = new_child(geometry, "DataItem")
|
||||
set_attribute(dataitem, "DataType", "Float")
|
||||
set_attribute(dataitem, "Format", "XML")
|
||||
#set_attribute(dataitem, "Precision", 8)
|
||||
s = ASCIIString[]
|
||||
ndim = 0
|
||||
for i in xdmf.permutation
|
||||
ndim += length(X[i])
|
||||
push!(s, join(round(X[i], 5), " "))
|
||||
end
|
||||
set_attribute(dataitem, "Dimensions", ndim)
|
||||
add_text(dataitem, "\n"*join(s, "\n")*"\n")
|
||||
|
||||
# 3. write elements
|
||||
topology = new_child(grid, "Topology")
|
||||
set_attribute(topology, "TopologyType", "Mixed")
|
||||
set_attribute(topology, "NumberOfElements", length(elements))
|
||||
dataitem = new_child(topology, "DataItem")
|
||||
set_attribute(dataitem, "Format", "XML")
|
||||
set_attribute(dataitem, "DataType", "Int")
|
||||
# set_attribute(dataitem, "Precision", 8)
|
||||
s = ASCIIString[]
|
||||
eldim = 0
|
||||
for element in elements
|
||||
eltype = get_xdmf_element_code(element)
|
||||
# note: id numbers start from 0 in Xdmf
|
||||
conn = [iperm[j] for j in get_connectivity(element)] - 1
|
||||
data = [eltype; conn]
|
||||
eldim += length(data)
|
||||
push!(s, join(data, " "))
|
||||
end
|
||||
set_attribute(dataitem, "Dimensions", eldim)
|
||||
add_text(dataitem, "\n"*join(s, "\n")*"\n")
|
||||
end
|
||||
|
||||
function xdmf_save_field!(xdmf, elements, time, field_name; field_type="Scalar")
|
||||
f = Dict()
|
||||
for element in elements
|
||||
g = element[field_name](time)
|
||||
conn = get_connectivity(element)
|
||||
for (i, c) in enumerate(conn)
|
||||
f[c] = g[i]
|
||||
end
|
||||
end
|
||||
|
||||
attribute = new_child(xdmf.current_grid, "Attribute")
|
||||
set_attribute(attribute, "Center", "Node")
|
||||
set_attribute(attribute, "Name", ucfirst(field_name))
|
||||
set_attribute(attribute, "Type", field_type)
|
||||
dataitem = new_child(attribute, "DataItem")
|
||||
set_attribute(dataitem, "DataType", "Float")
|
||||
set_attribute(dataitem, "Format", "XML")
|
||||
#set_attribute(dataitem, "Precision", 8)
|
||||
s = ASCIIString[]
|
||||
dim = 0
|
||||
for i in xdmf.permutation
|
||||
push!(s, join(round(f[i], 5), " "))
|
||||
dim += length(f[i])
|
||||
end
|
||||
set_attribute(dataitem, "Dimensions", dim)
|
||||
add_text(dataitem, "\n"*join(s, "\n")*"\n")
|
||||
end
|
||||
|
||||
function xdmf_save!(xdmf, filename)
|
||||
save_file(xdmf.xdoc, filename)
|
||||
end
|
||||
|
||||
+24
-11
@@ -66,6 +66,7 @@ type Problem{P<:AbstractProblem}
|
||||
dimension :: Int # degrees of freedom per node
|
||||
parent_field_name :: ASCIIString # (optional) name of parent field e.g. "displacement"
|
||||
elements :: Vector{Element}
|
||||
dofmap :: Dict{Element, Vector{Int64}} # connects element local dofs to global dofs
|
||||
assembly :: Assembly
|
||||
properties :: P
|
||||
end
|
||||
@@ -79,8 +80,8 @@ Create vector-valued (dim=3) elasticity problem:
|
||||
julia> prob = Problem(Elasticity, "this is my problem", 3)
|
||||
|
||||
"""
|
||||
function Problem{P<:FieldProblem}(::Type{P}, name, dimension, elements=[])
|
||||
Problem{P}(name, dimension, "none", elements, Assembly(), P())
|
||||
function Problem{P<:FieldProblem}(::Type{P}, name, dimension, elements=[], dofmap=Dict())
|
||||
Problem{P}(name, dimension, "none", elements, dofmap, Assembly(), P())
|
||||
end
|
||||
|
||||
""" Construct a new boundary problem.
|
||||
@@ -92,8 +93,8 @@ Create Dirichlet boundary problem for vector-valued (dim=3) elasticity problem.
|
||||
julia> bc1 = Problem(Dirichlet, "support", 3, "displacement")
|
||||
|
||||
"""
|
||||
function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name, elements=[])
|
||||
Problem{P}(name, dimension, parent_field_name, elements, Assembly(), P())
|
||||
function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name, elements=[], dofmap=Dict())
|
||||
Problem{P}(name, dimension, parent_field_name, elements, dofmap, Assembly(), P())
|
||||
end
|
||||
|
||||
function get_formulation_type{P<:FieldProblem}(problem::Problem{P})
|
||||
@@ -116,7 +117,7 @@ function initialize!(problem::Problem, time::Real)
|
||||
field_name = get_unknown_field_name(problem)
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
for element in get_elements(problem)
|
||||
gdofs = get_gdofs(element, problem)
|
||||
gdofs = get_gdofs(problem, element)
|
||||
if haskey(element, field_name)
|
||||
# if field is found, copy last known solution to new time as initial guess
|
||||
if !isapprox(last(element[field_name]).time, time)
|
||||
@@ -134,7 +135,7 @@ function initialize!(problem::Problem, time::Real)
|
||||
#is_dirichlet_problem(problem) && return
|
||||
field_name = get_parent_field_name(problem)
|
||||
for element in get_elements(problem)
|
||||
gdofs = get_gdofs(element, problem)
|
||||
gdofs = get_gdofs(problem, element)
|
||||
if haskey(element, field_name)
|
||||
# if field is found, copy last known solution to new time as initial guess
|
||||
if !isapprox(last(element[field_name]).time, time)
|
||||
@@ -269,16 +270,28 @@ end
|
||||
|
||||
function get_gdofs(element::Element, dim::Int)
|
||||
conn = get_connectivity(element)
|
||||
gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
|
||||
if length(conn) == 0
|
||||
error("element connectivity not defined, cannot determine global dofs for element: $element")
|
||||
end
|
||||
gdofs = vec([dim*(i-1)+j for j=1:dim, i in conn])
|
||||
return gdofs
|
||||
end
|
||||
|
||||
function get_gdofs(element::Element, problem::Problem)
|
||||
return get_gdofs(element, problem.dimension)
|
||||
end
|
||||
""" Return global degrees of freedom for element.
|
||||
|
||||
Notes
|
||||
-----
|
||||
First look dofs from problem.dofmap, it not found, update dofmap from
|
||||
element.element connectivity using formula gdofs = [dim*(nid-1)+j for j=1:dim]
|
||||
1. look element dofs from problem.dofmap
|
||||
2. if not found, use element.connectivity to update dofmap and 1.
|
||||
"""
|
||||
function get_gdofs(problem::Problem, element::Element)
|
||||
return get_gdofs(element, problem.dimension)
|
||||
if !haskey(element, problem.dofmap)
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
problem.dofmap[element] = get_gdofs(element, dim)
|
||||
end
|
||||
return problem.dofmap[element]
|
||||
end
|
||||
|
||||
""" Find dofs corresponding to nodes. """
|
||||
|
||||
@@ -1,57 +1,6 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
""" Calculate nodal vector from set of elements.
|
||||
|
||||
For example element 1 with dofs [1, 2, 3, 4] has [1, 1, 1, 1] and
|
||||
element 2 with dofs [3, 4, 5, 6] has [2, 2, 2, 2] the result will
|
||||
be sparse matrix with values [1, 1, 3, 3, 2, 2].
|
||||
|
||||
Parameters
|
||||
----------
|
||||
field_name
|
||||
name of field, e.g. "geometry"
|
||||
field_dim
|
||||
degrees of freedom / node
|
||||
elements
|
||||
elements used to calculate vector
|
||||
vec_dim
|
||||
used to resize solution vector if given
|
||||
time
|
||||
"""
|
||||
function calculate_nodal_vector(field_name, field_dim, elements::Vector{Element},
|
||||
time, vec_dim=0)
|
||||
A = SparseMatrixCOO()
|
||||
b = SparseMatrixCOO()
|
||||
for element in elements
|
||||
haskey(element, field_name) || continue
|
||||
gdofs = get_gdofs(element, 1)
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
J = get_jacobian(element, ip, time)
|
||||
w = ip.weight*norm(J)
|
||||
f = element(field_name, ip, time)
|
||||
N = element(ip, time)
|
||||
add!(A, gdofs, gdofs, w*kron(N', N))
|
||||
for dim=1:field_dim
|
||||
add!(b, gdofs, w*f[dim]*N, dim)
|
||||
end
|
||||
end
|
||||
end
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
nz = sort(unique(rowvals(A)))
|
||||
x = zeros(size(b)...)
|
||||
x[nz, :] = A[nz,nz] \ b[nz, :]
|
||||
x = vec(transpose(x))
|
||||
if vec_dim != 0
|
||||
v = zeros(vec_dim)
|
||||
v[1:length(x)] = x
|
||||
return v
|
||||
else
|
||||
return x
|
||||
end
|
||||
end
|
||||
|
||||
function calculate_rotated_nodal_vector(field_name, field_dim, elements::Vector{Element},
|
||||
time, vec_dim=0)
|
||||
A = SparseMatrixCOO()
|
||||
|
||||
-173
@@ -1,173 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using LightXML
|
||||
using JuliaFEM
|
||||
|
||||
# element codes: http://www.paraview.org/pipermail/paraview/2013-July/028859.html
|
||||
# > from ./VTK/ThirdParty/xdmf2/vtkxdmf2/libsrc/XdmfTopology.h
|
||||
# >
|
||||
# > // Topologies
|
||||
# > #define XDMF_NOTOPOLOGY 0x0
|
||||
# > #define XDMF_POLYVERTEX 0x1
|
||||
# > #define XDMF_POLYLINE 0x2
|
||||
# > #define XDMF_POLYGON 0x3
|
||||
# > #define XDMF_TRI 0x4
|
||||
# > #define XDMF_QUAD 0x5
|
||||
# > #define XDMF_TET 0x6
|
||||
# > #define XDMF_PYRAMID 0x7
|
||||
# > #define XDMF_WEDGE 0x8
|
||||
# > #define XDMF_HEX 0x9
|
||||
# > #define XDMF_EDGE_3 0x0022
|
||||
# > #define XDMF_TRI_6 0x0024
|
||||
# > #define XDMF_QUAD_8 0x0025
|
||||
# > #define XDMF_QUAD_9 0x0023
|
||||
# > #define XDMF_TET_10 0x0026
|
||||
# > #define XDMF_PYRAMID_13 0x0027
|
||||
# > #define XDMF_WEDGE_15 0x0028
|
||||
# > #define XDMF_WEDGE_18 0x0029
|
||||
# > #define XDMF_HEX_20 0x0030
|
||||
# > #define XDMF_HEX_24 0x0031
|
||||
# > #define XDMF_HEX_27 0x0032
|
||||
# > #define XDMF_MIXED 0x0070
|
||||
# > #define XDMF_2DSMESH 0x0100
|
||||
# > #define XDMF_2DRECTMESH 0x0101
|
||||
# > #define XDMF_2DCORECTMESH 0x0102
|
||||
# > #define XDMF_3DSMESH 0x1100
|
||||
# > #define XDMF_3DRECTMESH 0x1101
|
||||
# > #define XDMF_3DCORECTMESH 0x1102
|
||||
|
||||
global eltypes = Dict{Symbol, Int}(
|
||||
:Tri3 => 0x4,
|
||||
:Quad4 => 0x5,
|
||||
:Tet4 => 0x6,
|
||||
:Hex8 => 0x9,
|
||||
:Tet10 => 0x0026)
|
||||
|
||||
function xdmf_new_model(xdmf_version="2.1")
|
||||
xdoc = XMLDocument()
|
||||
xroot = create_root(xdoc, "Xdmf")
|
||||
set_attribute(xroot, "xmlns:xi", "http://www.w3.org/2001/XInclude")
|
||||
set_attribute(xroot, "Version", xdmf_version)
|
||||
domain = new_child(xroot, "Domain")
|
||||
return xdoc, domain
|
||||
end
|
||||
|
||||
function xdmf_new_temporal_collection(model)
|
||||
temporal_collection = new_child(model, "Grid")
|
||||
set_attribute(temporal_collection, "CollectionType", "Temporal")
|
||||
set_attribute(temporal_collection, "GridType", "Collection")
|
||||
set_attribute(temporal_collection, "Name", "Collection")
|
||||
# geometry = new_child(temporal_collection, "Geometry")
|
||||
# set_attribute(geometry, "Type", "None")
|
||||
# topology = new_child(temporal_collection, "Topology")
|
||||
# set_attribute(topology, "Dimensions", "0")
|
||||
# set_attribute(topology, "Type", "NoTopology")
|
||||
return temporal_collection
|
||||
end
|
||||
|
||||
function xdmf_new_grid(temporal_collection; time=0)
|
||||
grid = new_child(temporal_collection, "Grid")
|
||||
set_attribute(grid, "Name", "Grid")
|
||||
time_ = new_child(grid, "Time")
|
||||
set_attribute(time_, "Value", time)
|
||||
return grid
|
||||
end
|
||||
|
||||
function xdmf_new_mesh!(grid, nodes, elements; datatype="XYZ")
|
||||
|
||||
# 1. write nodes
|
||||
geometry = new_child(grid, "Geometry")
|
||||
set_attribute(geometry, "Type", datatype)
|
||||
dataitem = new_child(geometry, "DataItem")
|
||||
set_attribute(dataitem, "DataType", "Float")
|
||||
ndim = sum([length(node) for node in nodes])
|
||||
info("XDFM: ndim = $ndim")
|
||||
set_attribute(dataitem, "Dimensions", "$ndim")
|
||||
set_attribute(dataitem, "Format", "XML")
|
||||
set_attribute(dataitem, "Precision", 8)
|
||||
s = join([join(node, " ") for node in round(nodes, 5)], "\n")
|
||||
add_text(dataitem, "\n"*s*"\n")
|
||||
|
||||
# 2. write elements
|
||||
topology = new_child(grid, "Topology")
|
||||
eldim = sum([length(element[2]) for element in elements]) + length(elements)
|
||||
set_attribute(topology, "TopologyType", "Mixed")
|
||||
set_attribute(topology, "NumberOfElements", length(elements))
|
||||
dataitem = new_child(topology, "DataItem")
|
||||
set_attribute(dataitem, "Format", "XML")
|
||||
set_attribute(dataitem, "DataType", "Int")
|
||||
set_attribute(dataitem, "Dimensions", "$eldim")
|
||||
# set_attribute(dataitem, "Precision", 8)
|
||||
# note: id numbers start from 0 in Xdmf
|
||||
s = join([join([eltypes[eltype]; connectivity-1], " ") for (eltype, connectivity) in elements], "\n")
|
||||
add_text(dataitem, "\n"*s*"\n")
|
||||
|
||||
end
|
||||
|
||||
""" Write Vector field to nodes. """
|
||||
function xdmf_new_nodal_field!(grid, name, data)
|
||||
attribute = new_child(grid, "Attribute")
|
||||
set_attribute(attribute, "Center", "Node")
|
||||
set_attribute(attribute, "Name", name)
|
||||
set_attribute(attribute, "Type", "Vector")
|
||||
dataitem = new_child(attribute, "DataItem")
|
||||
set_attribute(dataitem, "DataType", "Float")
|
||||
ndim = sum([length(d) for d in data])
|
||||
set_attribute(dataitem, "Dimensions", "$ndim")
|
||||
set_attribute(dataitem, "Format", "XML")
|
||||
set_attribute(dataitem, "Precision", 8)
|
||||
s = join([join(d, " ") for d in round(data, 5)], "\n")
|
||||
add_text(dataitem, "\n"*s*"\n")
|
||||
end
|
||||
|
||||
function xdmf_save_model(xdoc, filename)
|
||||
save_file(xdoc, filename)
|
||||
end
|
||||
|
||||
function xdmf_dump(all_elements, eltype, elsym, time=0.0, filename="/tmp/xdmf_result.xmf")
|
||||
info("$(length(all_elements)) elements.")
|
||||
xdoc, xmodel = xdmf_new_model()
|
||||
coll = xdmf_new_temporal_collection(xmodel)
|
||||
grid = xdmf_new_grid(coll; time=time)
|
||||
|
||||
Xg = Dict{Int64, Vector{Float64}}()
|
||||
ug = Dict{Int64, Vector{Float64}}()
|
||||
nids = Dict{Int64, Int64}()
|
||||
for element in all_elements
|
||||
conn = get_connectivity(element)
|
||||
for (i, c) in enumerate(conn)
|
||||
nids[c] = c
|
||||
end
|
||||
X = element("geometry", time)
|
||||
for (i, c) in enumerate(conn)
|
||||
Xg[c] = X[i]
|
||||
end
|
||||
haskey(element, "displacement") || continue
|
||||
u = element("displacement", time)
|
||||
for (i, c) in enumerate(conn)
|
||||
ug[c] = u[i]
|
||||
end
|
||||
end
|
||||
perm = sort(collect(keys(Xg)))
|
||||
nodes = Vector{Float64}[Xg[i] for i in perm]
|
||||
disp = Vector{Float64}[ug[i] for i in perm]
|
||||
nids = Int[nids[i] for i in perm]
|
||||
inids = Dict{Int64, Int64}()
|
||||
for (i, nid) in enumerate(nids)
|
||||
inids[nid] = i
|
||||
end
|
||||
elements = []
|
||||
for element in all_elements
|
||||
isa(element, eltype) || continue
|
||||
conn = get_connectivity(element)
|
||||
nconn = [inids[i] for i in conn]
|
||||
push!(elements, (elsym, nconn))
|
||||
end
|
||||
|
||||
xdmf_new_mesh!(grid, nodes, elements)
|
||||
xdmf_new_nodal_field!(grid, "displacement", disp)
|
||||
xdmf_save_model(xdoc, filename)
|
||||
info("model dumped to $filename")
|
||||
end
|
||||
|
||||
@@ -0,0 +1,22 @@
|
||||
# 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
|
||||
|
||||
@testset "geometry missing" begin
|
||||
el = Element(Quad4, [1, 2, 3, 4])
|
||||
pr = Problem(Elasticity, "problem", 2)
|
||||
# this throws KeyError: geometry not found.
|
||||
# it's descriptive enough to give hint to user
|
||||
# what went wrong
|
||||
@test_throws KeyError assemble!(pr, el)
|
||||
end
|
||||
|
||||
@testset "connectivity information missing" begin
|
||||
el = Element(Quad4)
|
||||
nodes = Vector{Float64}[[0,0],[1,0],[1,1],[0,1]]
|
||||
update!(el, "geometry", nodes)
|
||||
pr = Problem(Elasticity, "problem", 2)
|
||||
@test_throws Exception assemble!(pr, el)
|
||||
end
|
||||
@@ -0,0 +1,71 @@
|
||||
# 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.Test
|
||||
|
||||
@testset "2d nonlinear elasticity: test nonhomogeneous boundary conditions and stress calculation" begin
|
||||
|
||||
# field problem
|
||||
block = Problem(Elasticity, "BLOCK", 2)
|
||||
block.properties.formulation = :plane_stress
|
||||
|
||||
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])
|
||||
|
||||
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)
|
||||
|
||||
# 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)
|
||||
|
||||
solver = Solver("solve block problem")
|
||||
push!(solver, block, bc)
|
||||
call(solver)
|
||||
|
||||
# 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]
|
||||
|
||||
u3 = reshape(block.assembly.u, 2, 4)[:, 3]
|
||||
info("u3 = $u3")
|
||||
@test isapprox(u3, u3_expected, atol=1.0e-5)
|
||||
|
||||
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)
|
||||
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
|
||||
|
||||
@@ -0,0 +1,34 @@
|
||||
# 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.Postprocess
|
||||
using JuliaFEM.Test
|
||||
|
||||
@testset "extrapolate stress from gauss points to nodes" begin
|
||||
X = Dict{Int, Vector{Float64}}(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [6.0, 0.0],
|
||||
3 => [6.0, 6.0],
|
||||
4 => [0.0, 6.0],
|
||||
5 => [12.0, 0.0],
|
||||
6 => [12.0, 6.0])
|
||||
el1 = Element(Quad4, [1, 2, 3, 4])
|
||||
el2 = Element(Quad4, [2, 5, 6, 3])
|
||||
el1.id = 1
|
||||
el2.id = 2
|
||||
elements = [el1, el2]
|
||||
time = 0.0
|
||||
update!(elements, "geometry", X)
|
||||
update!(get_integration_points(el1), "stress", time => [1.0, 2.0, 3.0])
|
||||
update!(get_integration_points(el2), "stress", time => [2.0, 3.0, 4.0])
|
||||
field_name = "stress"
|
||||
field_dim = 3
|
||||
calc_nodal_values!(elements, field_name, field_dim, time)
|
||||
s1 = el1("stress", [0.0, 0.0], time)
|
||||
s2 = el2("stress", [0.0, 0.0], time)
|
||||
# visually checked, see blog post "Postprocessing stress"
|
||||
@test isapprox(s1, [1.125, 2.125, 3.125])
|
||||
@test isapprox(s2, [1.875, 2.875, 3.875])
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user