From d429e5b2fcd44328222f69730f712b65577763fa Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 28 May 2016 21:39:29 +0300 Subject: [PATCH] calculate stress, interpolate stress to nodes using lsq fitting --- src/JuliaFEM.jl | 12 +- src/elasticity.jl | 35 +++- src/postprocess_xdmf.jl | 162 ++++++++++++++++ src/problems.jl | 35 ++-- src/utils.jl | 51 ------ src/xdmf.jl | 173 ------------------ test/test_common_failures.jl | 22 +++ ...y_2d_nonhomogeneous_boundary_conditions.jl | 71 +++++++ test/test_extrapolate_to_nodes.jl | 34 ++++ 9 files changed, 347 insertions(+), 248 deletions(-) create mode 100644 src/postprocess_xdmf.jl delete mode 100644 src/xdmf.jl create mode 100644 test/test_common_failures.jl create mode 100644 test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl create mode 100644 test/test_extrapolate_to_nodes.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 9683a93..d8fa744 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -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. """ diff --git a/src/elasticity.jl b/src/elasticity.jl index 7c84ca3..be9feab 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -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) diff --git a/src/postprocess_xdmf.jl b/src/postprocess_xdmf.jl new file mode 100644 index 0000000..a87a7eb --- /dev/null +++ b/src/postprocess_xdmf.jl @@ -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 + diff --git a/src/problems.jl b/src/problems.jl index 8cfd39a..6e7270f 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -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. """ diff --git a/src/utils.jl b/src/utils.jl index f382cc8..76a70a8 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -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() diff --git a/src/xdmf.jl b/src/xdmf.jl deleted file mode 100644 index f248900..0000000 --- a/src/xdmf.jl +++ /dev/null @@ -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 - diff --git a/test/test_common_failures.jl b/test/test_common_failures.jl new file mode 100644 index 0000000..af43a2c --- /dev/null +++ b/test/test_common_failures.jl @@ -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 diff --git a/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl b/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl new file mode 100644 index 0000000..2876cb9 --- /dev/null +++ b/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl @@ -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 + diff --git a/test/test_extrapolate_to_nodes.jl b/test/test_extrapolate_to_nodes.jl new file mode 100644 index 0000000..5c44a25 --- /dev/null +++ b/test/test_extrapolate_to_nodes.jl @@ -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 +