diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index a5b5ff4..7bf8a6c 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -44,7 +44,9 @@ export Node, AbstractElement, Element, update!, get_connectivity, get_basis, get_dbasis, inside, get_local_coordinates include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements -export get_reference_coordinates, get_interpolation_polynomial +export get_reference_coordinates, + get_interpolation_polynomial, + description export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, @@ -62,7 +64,7 @@ include("integrate.jl") # default integration points for elements export get_integration_points include("sparse.jl") -export add!, SparseMatrixCOO, SparseVectorCOO, get_nonzero_rows, get_nonzero_columns +export add!, SparseMatrixCOO, SparseVectorCOO, get_nonzero_rows, get_nonzero_columns, optimize!, resize_sparse, resize_sparsevec include("problems.jl") # common problem routines export Problem, AbstractProblem, FieldProblem, BoundaryProblem, diff --git a/src/abaqus.jl b/src/abaqus.jl index 3907abc..64638fe 100644 --- a/src/abaqus.jl +++ b/src/abaqus.jl @@ -703,7 +703,7 @@ end function abaqus_open_results(name) path = abaqus_input_file_path(name) result_file = "$path/$name.xmf" - return XDMF(result_file) + return Xdmf(result_file) end ### JuliaFEM-ABAQUS interface entry point @@ -759,3 +759,7 @@ function create_surface_elements(mesh::Mesh, surface_name::Symbol) return elements end +function create_surface_elements(mesh::Mesh, surface_name::String) + return create_surface_elements(mesh, Symbol(surface_name)) +end + diff --git a/src/assembly.jl b/src/assembly.jl index 692a2f6..ed78052 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -1,46 +1,6 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -function optimize!(assembly::Assembly) - optimize!(assembly.K) - optimize!(assembly.Kg) - optimize!(assembly.f) - optimize!(assembly.fg) - optimize!(assembly.C1) - optimize!(assembly.C2) - optimize!(assembly.D) - optimize!(assembly.g) - optimize!(assembly.c) -end - -function append!(assembly::Assembly, sub_assembly::Assembly) - append!(assembly.M, sub_assembly.M) - append!(assembly.K, sub_assembly.K) - append!(assembly.Kg, sub_assembly.Kg) - append!(assembly.f, sub_assembly.f) - append!(assembly.fg, sub_assembly.fg) - append!(assembly.C1, sub_assembly.C1) - append!(assembly.C2, sub_assembly.C2) - append!(assembly.D, sub_assembly.D) - append!(assembly.g, sub_assembly.g) - append!(assembly.c, sub_assembly.c) -end - -""" Calculate norm of assembly, i.e., norm of each block of matrix. """ -function norm(assembly::Assembly, p=2) - N1 = norm(assembly.M, p) - N2 = norm(assembly.K, p) - N3 = norm(assembly.Kg, p) - N4 = norm(assembly.f, p) - N5 = norm(assembly.fg, p) - N6 = norm(assembly.C1, p) - N7 = norm(assembly.C2, p) - N8 = norm(assembly.D, p) - N9 = norm(assembly.g, p) - N10 = norm(assembly.c, p) - return [N1, N2, N3, N4, N5, N6, N7, N8, N9, N10] -end - function isapprox(a1::Assembly, a2::Assembly) T = isapprox(a1.K, a2.K) T &= isapprox(a1.C1, a2.C1) diff --git a/src/elements.jl b/src/elements.jl index 30eddb5..241e791 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -29,6 +29,22 @@ function setindex!(element::Element, data::Field, field_name) element.fields[field_name] = data end +function get_element_type{E}(element::Element{E}) + return E +end + +function get_element_id{E}(element::Element{E}) + return element.id +end + +function is_element_type{E}(element::Element{E}, element_type) + return is(E, element_type) +end + +function filter_by_element_type(element_type, elements) + return filter(element -> is_element_type(element, element_type), elements) +end + function setindex!(element::Element, data::Function, field_name) if method_exists(data, Tuple{Element, Vector, Float64}) # create enclosure to pass element as argument @@ -100,7 +116,7 @@ julia> el([0.0, 0.0], 0.0, 2) """ function (element::Element)(ip, time::Float64, dim::Int) dim == 1 && return get_basis(element, ip, time) - Ni = get_basis(element, ip, time) + Ni = vec(get_basis(element, ip, time)) N = zeros(dim, length(element)*dim) for i=1:dim N[i,i:dim:end] += Ni @@ -144,14 +160,6 @@ function (element::Element)(field_name::String, ip, time::Float64, ::Type{Val{:G return element(ip, time, Val{:Grad})*element[field_name](time) end -function (element::Element)(field::Field, time::Float64) - return field(time) -end - -function (element::Element)(field::DCTI, time::Float64) - return field.data -end - function (element::Element)(field_name::String, ip, time::Float64) field = element[field_name] return element(field, ip, time) @@ -220,35 +228,8 @@ function update!{K,V}(element::Element, field_name, data::Pair{Float64, Dict{K, time, field_data = data element_data = V[field_data[i] for i in get_connectivity(element)] update!(element, field_name, time => element_data) - #if haskey(element, field_name) - # update!(element[field_name], data) - #else - # element[field_name] = Field(data) - #end end -function update!(element::Element, field_name::AbstractString, datas::Union{Real, Vector, Pair{Float64, Union{Float64, Real, Vector{Any}}}}...) - for data in datas - if haskey(element, field_name) - update!(element[field_name], data) - else - if length(data) != length(element) - update!(element, field_name, DCTI(data)) - else - element[field_name] = data - end - end - end -end - -#= -function update!(element::Element, field_name, data::Pair...) - for data in datas - update!(element, field_name, data) - end -end -=# - function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Any}}) if haskey(element, field_name) update!(element[field_name], data) @@ -353,10 +334,6 @@ function get_integration_points(element::Element, change_order::Int) return [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)] end -function get_gdofs(element::Element) - return get_gdofs(element, 1) -end - """ Return dual basis transformation matrix Ae. """ function get_dualbasis(element::Element, time::Float64, order=1) nnodes = length(element) diff --git a/src/elements_lagrange.jl b/src/elements_lagrange.jl index a9266ba..83a3269 100644 --- a/src/elements_lagrange.jl +++ b/src/elements_lagrange.jl @@ -1,23 +1,66 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +global const ELEMENT_DESCRIPTIONS = Dict( + "Poi1" => "1 node discrete point element", + "Seg2" => "2 node linear segment/line element", + "Seg3" => "3 node quadratic segment/line element", + "Tri3" => "3 node linear triangle element", + "Tri6" => "6 node quadratic triangle element", + "Tri7" => "7 node quadratic triangle element (has middle node)", + "Quad4" => "4 node linear quadrangle element", + "Quad8" => "8 node quadratic quadrangle element (Serendip)", + "Quad9" => "9 node quadratic quadrangle element", + "Tet4" => "4 node linear tetrahedral element", + "Tet10" => "10 node quadratic tetrahedral element", + "Wedge6" => "6 node linear prismatic element (wedge)", + "Wedge15" => "15 node quadratic prismatic element (wedge)", + "Hex8" => "8 node linear hexahedral element", + "Hex20" => "20 node biquadratic hexahedral element", + "Hex27" => "27 node quadratic hexahedral element") + +global const ELEMENT_SIZES = Dict( + "Poi1" => (0, 1), + "Seg2" => (1, 2), + "Seg3" => (1, 3), + "Tri3" => (2, 3), + "Tri6" => (2, 6), + "Tri7" => (2, 7), + "Quad4" => (2, 4), + "Quad8" => (2, 8), + "Quad9" => (2, 9), + "Tet4" => (3, 4), + "Tet10" => (3, 10), + "Wedge6" => (3, 6), + "Wedge15" => (3, 15), + "Hex8" => (3, 8), + "Hex20" => (3, 20), + "Hex27" => (3, 27)) + +""" Return description line of element. """ +function description{T}(element::Element{T}) + element_type = last(split("$T", '.')) + return get(ELEMENT_DESCRIPTIONS, element_type, "Unknown element description") +end + +""" Return size of element, i.e. tuple (n, m) where n is dimension of element +(0, 1, 2, 3) and m is number of nodes. """ +function size{T}(element::Element{T}) + element_type = last(split("$T", '.')) + return ELEMENT_SIZES[element_type] +end + +""" Return length of element, i.e. number of nodes. """ +function length{T}(element::Element{T}) + return size(element)[end] +end + ### 0d element + type Poi1 <: AbstractElement end -function description(::Type{Poi1}) - "1 node point" -end - -function size(element::Element{Poi1}) - return (0, 1) -end - -function length(element::Element{Poi1}) - return 1 -end - function get_basis(element::Element{Poi1}, ip, time) return [1] end @@ -47,18 +90,6 @@ end type Seg2 <: AbstractElement end -function description(::Type{Seg2}) - "2 node segment" -end - -function size(element::Element{Seg2}) - return (1, 2) -end - -function length(element::Element{Seg2}) - return 2 -end - function get_reference_coordinates(::Type{Seg2}) Vector{Float64}[ [-1.0], # N1 @@ -78,18 +109,6 @@ end type Seg3 <: AbstractElement end -function description(::Type{Seg3}) - "3 node segment" -end - -function size(element::Element{Seg3}) - return (1, 3) -end - -function length(element::Element{Seg3}) - return 3 -end - function get_reference_coordinates(::Type{Seg3}) Vector{Float64}[ [-1.0], # N1 @@ -110,18 +129,6 @@ end type Tri3 <: AbstractElement end -function description(::Type{Tri3}) - "3 node triangle" -end - -function size(element::Element{Tri3}) - return (2, 3) -end - -function length(element::Element{Tri3}) - return 3 -end - function get_reference_coordinates(::Type{Tri3}) Vector{Float64}[ [0.0, 0.0], # N1 @@ -147,18 +154,6 @@ end type Tri6 <: AbstractElement end -function description(::Type{Tri6}) - "6 node triangle" -end - -function size(element::Element{Tri6}) - return (2, 6) -end - -function length(element::Element{Tri6}) - return 6 -end - function get_reference_coordinates(::Type{Tri6}) Vector{Float64}[ [0.0, 0.0], # N1 @@ -187,18 +182,6 @@ end type Tri7 <: AbstractElement end -function description(::Type{Tri7}) - "7 node triangle" -end - -function size(element::Element{Tri7}) - return (2, 7) -end - -function length(element::Element{Tri7}) - return 7 -end - function get_reference_coordinates(::Type{Tri7}) Vector{Float64}[ [0.0, 0.0], # N1 @@ -228,18 +211,6 @@ end type Quad4 <: AbstractElement end -function description(::Type{Quad4}) - "4 node quadrangle" -end - -function size(element::Element{Quad4}) - return (2, 4) -end - -function length(element::Element{Quad4}) - return 4 -end - function get_reference_coordinates(::Type{Quad4}) Vector{Float64}[ [-1.0, -1.0], # N1 @@ -266,18 +237,6 @@ end type Quad8 <: AbstractElement end -function description(::Type{Quad8}) - "8 node Serendip quadrangle" -end - -function size(element::Element{Quad8}) - return (2, 8) -end - -function length(element::Element{Quad8}) - return 8 -end - function get_reference_coordinates(::Type{Quad8}) Vector{Float64}[ [-1.0, -1.0], # N1 @@ -308,18 +267,6 @@ end type Quad9 <: AbstractElement end -function description(::Type{Quad9}) - "9 node quadrangle" -end - -function size(element::Element{Quad9}) - return (2, 9) -end - -function length(element::Element{Quad9}) - return 9 -end - function get_reference_coordinates(::Type{Quad9}) Vector{Float64}[ [-1.0, -1.0], # N1 @@ -351,18 +298,6 @@ end type Tet4 <: AbstractElement end -function description(::Type{Tet4}) - "4 node tetrahedral element" -end - -function size(element::Element{Tet4}) - return (3, 4) -end - -function length(element::Element{Tet4}) - return 4 -end - function get_reference_coordinates(::Type{Tet4}) Vector{Float64}[ [0.0, 0.0, 0.0], # N1 @@ -390,18 +325,6 @@ end type Tet10 <: AbstractElement end -function description(::Type{Tet10}) - "10 node tetrahedral element" -end - -function size(element::Element{Tet10}) - return (3, 10) -end - -function length(element::Element{Tet10}) - return 10 -end - function get_reference_coordinates(::Type{Tet10}) Vector{Float64}[ [0.0, 0.0, 0.0], # N1 @@ -435,18 +358,6 @@ end type Wedge6 <: AbstractElement end -function description(::Type{Wedge6}) - "6 node prismatic element (wedge)" -end - -function size(element::Element{Wedge6}) - return (3, 6) -end - -function length(element::Element{Wedge6}) - return 6 -end - function get_reference_coordinates(::Type{Wedge6}) Vector{Float64}[ [0.0, 0.0, -1.0], # N1 @@ -476,18 +387,6 @@ end type Wedge15 <: AbstractElement end -function description(::Type{Wedge15}) - "15 node prismatic element (wedge)" -end - -function size(element::Element{Wedge15}) - return (3, 15) -end - -function length(element::Element{Wedge15}) - return 15 -end - function get_reference_coordinates(::Type{Wedge15}) Vector{Float64}[ [0.0, 0.0, -1.0], # N1 @@ -526,18 +425,6 @@ end type Hex8 <: AbstractElement end -function description(::Type{Hex8}) - "8 node hexahedral element" -end - -function size(element::Element{Hex8}) - return (3, 8) -end - -function length(element::Element{Hex8}) - return 8 -end - function get_reference_coordinates(::Type{Hex8}) Vector{Float64}[ [-1.0, -1.0, -1.0], # N1 @@ -569,18 +456,6 @@ end type Hex20 <: AbstractElement end -function description(::Type{Hex20}) - "20 node hexahedral element" -end - -function size(element::Element{Hex20}) - return (3, 20) -end - -function length(element::Element{Hex20}) - return 20 -end - function get_reference_coordinates(::Type{Hex20}) Vector{Float64}[ [-1.0, -1.0, -1.0], # N1 @@ -624,18 +499,6 @@ end type Hex27 <: AbstractElement end -function description(::Type{Hex27}) - "27 node hexahedral element" -end - -function size(element::Element{Hex27}) - return (3, 27) -end - -function length(element::Element{Hex27}) - return 27 -end - function get_reference_coordinates(::Type{Hex27}) Vector{Float64}[ [-1.0, -1.0, -1.0], # N1 diff --git a/src/elements_nurbs.jl b/src/elements_nurbs.jl index 69e4aa0..7124557 100644 --- a/src/elements_nurbs.jl +++ b/src/elements_nurbs.jl @@ -92,12 +92,12 @@ function get_basis(element::Element{NSolid}, xi::Vector, time) tu = element.properties.knots_u tv = element.properties.knots_v tw = element.properties.knots_w - w = element.properties.weights - nu = length(tu) - nv = length(tv) - nw = length(tw) + weights = element.properties.weights + nu = length(tu)-pu-1 + nv = length(tv)-pv-1 + nw = length(tw)-pw-1 u, v, w = xi - N = [w[i,j,k]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv)*NURBS(k,pw,w,tw) for i=1:nu, j=1:nv, k=1:nw] + N = vec([weights[i,j,k]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv)*NURBS(k,pw,w,tw) for i=1:nu, j=1:nv, k=1:nw])' return N / sum(N) end diff --git a/src/fields.jl b/src/fields.jl index 45c40d5..b3e7541 100644 --- a/src/fields.jl +++ b/src/fields.jl @@ -10,48 +10,11 @@ abstract Variable <: AbstractField abstract TimeVariant <: AbstractField abstract TimeInvariant <: AbstractField - type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{TimeVariant,TimeInvariant}} data end -typealias FieldSet Dict{AbstractString, Field} - -### Basic data structure for discrete field - -type Increment{T} - time :: Float64 - data :: T -end - -function convert{T}(::Type{Increment{T}}, data::Pair{Float64,T}) - return Increment{T}(data[1], data[2]) -end - -function convert{T}(::Type{Increment{Vector{Vector{T}}}}, data::Pair{Float64, Matrix{T}}) - time = data[1] - content = data[2] - return Increment(time, Vector{T}[content[:,i] for i=1:size(content,2)]) -end - -function getindex{T}(increment::Increment{Vector{T}}, i::Int64) - return increment.data[i] -end - -### Basic data structure for continuous field - -type Basis - basis :: Function - dbasis :: Function -end - -function (basis::Basis)(xi::Vector) - basis.basis(xi) -end - -function (basis::Basis)(xi::Vector, ::Type{Val{:grad}}) - basis.dbasis(xi) -end +typealias FieldSet Dict{String, Field} ### Different field combinations and other typealiases @@ -64,156 +27,158 @@ typealias CVTI Field{Continuous, Variable, TimeInvariant} # can be used to inter typealias CCTV Field{Continuous, Constant, TimeVariant} # can be used to interpolate in time typealias CVTV Field{Continuous, Variable, TimeVariant} -typealias ScalarIncrement{T} Increment{T} -typealias VectorIncrement{T} Increment{Vector{T}} -typealias TensorIncrement{T} Increment{Matrix{T}} - -typealias DiscreteField Union{DCTI, DVTI, DCTV, DVTV} -typealias ContinuousField Union{CCTI, CVTI, CCTV, CVTV} -typealias ConstantField Union{DCTI, DCTV, CCTI, CCTV} -typealias VariableField Union{DVTI, DVTV, CVTI, CVTV} -typealias TimeInvariantField Union{DCTI, DVTI, CCTI, CVTI} -typealias TimeVariantField Union{DCTV, DVTV, CCTV, CVTV} +# Discrete fields -### Convenient functions to create fields +""" Discrete, constant, time-invariant field. This is constant in both spatial +direction and time direction, i.e. df/dX = 0 and df/dt = 0. -#function Base.convert(::Type{Field}, data) -# return Field(data) -#end +This is the most basic type of field having no anything special functionality. + +Examples +-------- + +julia> f = DCTI() +julia> update!(f, 1.0) + +Multiplying by constant works: + +julia> 2*f +2.0 + +Interpolation in time direction gives the same constant: + +julia> f(1.0) +1.0 + +By default, when calling Field with scalar, DCTI is assumed, i.e. + +julia> Field(0.0) == DCTI(0.0) +true + +""" +function DCTI() + return DCTI(nothing) +end + +function Field() + return DCTI() +end function Field(data) return DCTI(data) end +function ==(x::DCTI, y::DCTI) + return ==(x.data, y.data) +end + +function ==(x::DCTI, y) + return ==(x.data, y) +end + +function isapprox(x::DCTI, y::DCTI) + isapprox(x.data, y.data) +end + +function isapprox(x::DCTI, y) + isapprox(x.data, y) +end + +function length(f::DCTI) + return 1 +end + +function Base.:*(c::Number, f::DCTI) + return c*f.data +end + +""" Kind of spatial interpolation of DCTI. """ +function Base.:*(N::Matrix, f::DCTI) + @assert length(N) == 1 + return N[1]*f.data +end + +function update!(field::DCTI, data) + field.data = data +end + +""" Interpolate time-invariant field in time direction. """ +function (field::DCTI)(time::Float64) + return field.data +end + +""" Discrete, variable, time-invariant field. This is constant in time direction, +but not in spatial direction, i.e. df/dt = 0 but df/dX != 0. The basic structure +of data is Vector, and it is implicitly assumed that length of field matches to +the number of shape functions, so that interpolation in spatial direction works. + +Examples +-------- +""" +function DVTI() + return DVTI([]) +end + +""" For vector data, DVTI is automatically created. + +julia> DVTI([1.0, 2.0]) == Field([1.0, 2.0]) +true + +""" function Field(data::Vector) return DVTI(data) end -function Field{T}(data::Pair{Float64, T}...) - return DCTV([Increment{T}(d[1], d[2]) for d in data]) -end -#= -function Field{T}(data::Pair{Float64, Vector{T}}...) - return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data]) -end +""" For dictionary data, DVTI is automatically created. -function Field{T}(data::Pair{Float64, Dict{Int64, T}}...) - return DVTV([Increment{Dict{Int64, T}}(d[1], d[2]) for d in data]) -end -=# - -function Field{T<:Union{Vector, Dict}}(data::Pair{Float64, T}...) - return DVTV([Increment{T}(d[1], d[2]) for d in data]) -end +Define e.g. nodal coordinates in dictionary +julia> X = Dict(1 => [1.0, 2.0], 2 => [3.0, 4.0]) +julia> Field(X) == DVTI(X) +""" function Field(data::Dict) return DVTI(data) end -function convert{T}(::Type{DCTV}, data::Pair{Real, Vector{T}}...) - return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data]) +function ==(x::DVTI, y::DVTI) + return ==(x.data, y.data) end -""" Create new discrete, constant, time variant field. +function isapprox(x::DVTI, y) + return isapprox(x.data, y) +end -Examples --------- -julia> t0 = 0.0; t1=1.0; y0 = 0.0; y1 = 1.0 -julia> f = DCTV(t0 => y0, t1 => y1) +""" Default slicing of field. + +julia> f = DVTI([1.0, 2.0]) +julia> f[1] +1.0 """ -#function convert{T,v<:Real}(::Type{DCTV}, data::Pair{v, T}...) -# return DCTV([Increment(d[1],d[2]) for d in data]) -#end -function DCTV(data::Pair...) - return DCTV([Increment(d[1],d[2]) for d in data]) -end - -function Field(func::Function) - if method_exists(func, Tuple{}) - return CCTI(func) - elseif method_exists(func, Tuple{Float64}) - return CCTV(func) - elseif method_exists(func, Tuple{Vector}) - return CVTI(func) - elseif method_exists(func, Tuple{Vector, Number}) - return CVTV(func) - else - error("no proper definition found for function: check methods.") - end -end - -function CVTI(basis::Function, dbasis::Function) - return CVTI(Basis(basis, dbasis)) -end - -function Field(basis::Function, dbasis::Function) - return CVTI(basis, dbasis) -end - -### Accessing and manipulating discrete fields - -function getindex(field::DVTV, i::Int64) - return field.data[i] -end - -function push!(field::DCTV, data::Pair) - push!(field.data, data) -end - -function push!(field::DVTV, data::Pair) - push!(field.data, data) -end - function getindex(field::DVTI, i::Int64) return field.data[i] end +""" Multi-slicing of field. + +julia> f = DVTI([1.0, 2.0, 3.0]) +julia> f[[1, 3]] +[1.0, 3.0] + +""" function getindex(field::DVTI, I::Array{Int64, 1}) return [field.data[i] for i in I] end -function getindex(field::DCTV, i::Int64) - return field.data[i] -end - -function getindex(field::Field, i::Int64) - return field.data[i] -end - function length(field::DVTI) return length(field.data) end -function length(field::DCTI) +function start(field::DVTI) return 1 end -function length(field::DVTV) - return length(field.data) -end - -function length(field::DCTV) - return length(field.data) -end - -function first(field::Union{DCTV, DVTV}) - return field[1] -end - -function isapprox(f1::DCTI, f2::DCTI) - isapprox(f1.data, f2.data) -end - -for op = (:+, :*, :/, :-) - @eval ($op)(increment::Increment, field::DCTI) = ($op)(increment.data, field.data) - @eval ($op)(field::DCTI, increment::Increment) = ($op)(increment.data, field.data) - @eval ($op)(field1::DCTI, field2::DCTI) = ($op)(field1.data, field2.data) - @eval ($op)(field::DCTI, k::Number) = ($op)(field.data, k) - @eval ($op)(k::Number, field::DCTI) = ($op)(field.data, k) -end - function Base.:+(f1::DVTI, f2::DVTI) return DVTI(f1.data + f2.data) end @@ -222,49 +187,55 @@ function Base.:-(f1::DVTI, f2::DVTI) return DVTI(f1.data - f2.data) end -function Base.:*{T<:Real}(c::T, field::DVTI) - return DVTI(c*field.data) +function update!(field::DVTI, data::Union{Vector, Dict}) + field.data = data end -function Base.:*(N::Matrix, f::DCTI) - return f.data*N' +""" Take scalar product of DVTI and constant T. """ +function Base.:*(T::Number, field::DVTI) + return DVTI(T*field.data) end - - -# Multiply DVTI field with another vector T. Vector length -# must match to the field length and this can be used mainly -# for interpolation purposes, i.e., u = ∑ Nᵢuᵢ +""" Take dot product of DVTI field and vector T. Vector length must match to the +field length and this can be used mainly for interpolation purposes, i.e., u = ∑ Nᵢuᵢ. +""" function Base.:*(T::Vector, f::DVTI) @assert length(T) <= length(f) return sum([T[i]*f[i] for i=1:length(T)]) end +""" Take outer product of DVTI field and matrix T. """ +function Base.:*(T::Matrix, f::DVTI) + n, m = size(T) + return sum([kron(T[:,i], f[i]') for i=1:m])' +end + function vec(field::DVTI) return [field.data...;] end -function vec(field::DCTV) - error("trying to vectorize $field does not make sense") +""" Interpolate time-invariant field in time direction. """ +function (field::DVTI)(time::Float64) + return field end -function endof(field::Field) - return endof(field.data) -end +""" Create a similar DVTI field from vector data. -#function Base.similar{T}(field::DVTI, data::Vector{T}) -# return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment))) -#end +julia> f1 = DVTI(Vector[[1.0, 2.0], [3.0, 4.0]]) +julia> f2 = similar(f1, [2.0, 3.0, 4.0, 5.0]) +julia> f2 == DVTI(Vector[[2.0, 3.0], [4.0, 5.0]]) +true -function similar{T}(field::DVTI, data::Vector{T}) - n = length(field.data) - data = reshape(data, round(Int, length(data)/n), n) - newdata = Vector[data[:,i] for i=1:n] - return typeof(field)(newdata) -end - -function start(::DVTI) - return 1 +""" +function similar(field::DVTI, data::Vector) + n = length(field) + m = length(data) + dim = round(Int, m/n) + @assert dim*n == m + new_data = reshape(data, dim, n) + new_field = DVTI() + new_field.data = [new_data[:,i] for i=1:n] + return new_field end function next(f::DVTI, state) @@ -275,6 +246,114 @@ function done(f::DVTI, s) return s > length(f.data) end +""" Simple time frame / increment to contain both time and data. """ +type Increment{T} + time :: Float64 + data :: T +end + +""" Discrete, constant, time variant field. This is constant in spatial +direction but non-constant in time direction, i.e. df/dX = 0 but df/dt != 0. + +Examples +-------- +julia> t0 = 0.0; t1=1.0; y0 = 0.0; y1 = 1.0 +julia> f = DCTV(t0 => y0, t1 => y1) + +""" +function DCTV(data::Pair...) + return DCTV([Increment(d[1],d[2]) for d in data]) +end + +function Field{T}(data::Pair{Float64, T}...) + return DCTV([Increment{T}(d[1], d[2]) for d in data]) +end + +function getindex(field::DCTV, i::Int64) + return field.data[i] +end + +function length(field::DCTV) + return length(field.data) +end + +function first(field::DCTV) + return field[1] +end + +""" Interpolate constant time-variant field in time direction. """ +function (field::DCTV)(time::Number) + time < first(field).time && return DCTI(first(field).data) + time > last(field).time && return DCTI(last(field).data) + for i=reverse(1:length(field)) + isapprox(field[i].time, time) && return DCTI(field[i].data) + end + for i=reverse(2:length(field)) + t0 = field[i-1].time + t1 = field[i].time + if t0 < time < t1 + y0 = field[i-1].data + y1 = field[i].data + dt = t1-t0 + new_data = y0*(1-(time-t0)/dt) + y1*(1-(t1-time)/dt) + return DCTI(new_data) + end + end +end + +function endof(field::DCTV) + return endof(field.data) +end + +""" Discrete, variable, time variant fields. """ +function DVTV() + return DVTV(Increment[]) +end + +function DVTV{T<:Union{Vector, Dict}}(data::Pair{Float64, T}...) + return DVTV([Increment{T}(d[1], d[2]) for d in data]) +end + +function Field{T<:Union{Vector, Dict}}(data::Pair{Float64, T}...) + return DVTV([Increment{T}(d[1], d[2]) for d in data]) +end + +function length(field::DVTV) + return length(field.data) +end + +function getindex(field::DVTV, i::Int64) + return field.data[i] +end + +function first(field::DVTV) + return field[1] +end + +function endof(field::DVTV) + return endof(field.data) +end + +""" Interpolate discrete, variable, time-variant field in time direction. """ +function (field::DVTV)(time::Float64) + time < first(field).time && return DVTI(first(field).data) + time > last(field).time && return DVTI(last(field).data) + for i=reverse(1:length(field)) + isapprox(field[i].time, time) && return DVTI(field[i].data) + end + for i=reverse(2:length(field)) + t0 = field[i-1].time + t1 = field[i].time + if t0 < time < t1 + y0 = field[i-1].data + y1 = field[i].data + dt = t1-t0 + new_data = y0*(1-(time-t0)/dt) + y1*(1-(t1-time)/dt) + return DVTI(new_data) + end + end +end + """ Update time-dependent fields with new values. Examples @@ -300,146 +379,44 @@ function update!{T}(field::Union{DCTV, DVTV}, val::Pair{Float64, T}) end end -function update!{T}(field::Union{DCTI, DVTI}, val::T) - field.data = val +### Basic data structure for continuous field + +type Basis + basis :: Function + dbasis :: Function +end + +### Convenient functions to create fields + +function Field(func::Function) + if method_exists(func, Tuple{}) + return CCTI(func) + elseif method_exists(func, Tuple{Float64}) + return CCTV(func) + elseif method_exists(func, Tuple{Vector}) + return CVTI(func) + elseif method_exists(func, Tuple{Vector, Float64}) + return CVTV(func) + else + error("no proper definition found for function: check methods.") + end end ### Accessing continuous fields -function (field::CVTI)(xi::Vector) +function (field::CCTI)(xi::Vector, time::Number) + return field.data() +end + +function (field::CVTI)(xi::Vector, time::Number) return field.data(xi) end -function (field::CVTV)(xi, time::Float64) - return field.data(xi, time) -end - -function (field::CVTI)(xi::Vector, ::Type{Val{:Grad}}) - return field.data(xi, Val{:Grad}) -end - -function (field::CCTV)(time::Float64) +function (field::CCTV)(xi::Vector, time::Number) return field.data(time) end -function convert(::Type{Basis}, field::CVTI) - return field.data +function (field::CVTV)(xi::Vector, time::Number) + return field.data(xi, time) end -### Interpolation - -""" Interpolate time-invariant field in time direction. """ -function (field::DVTI)(time::Float64) - return field -end -function (field::DCTI)(time::Float64) - return field.data -end -function (field::CVTI)(time::Float64) - return field.data() -end -function (field::CCTI)(time::Float64) - return field.data() -end - -""" Interpolate constant time-variant field in time direction. """ -function (field::DCTV)(time::Real) - time < first(field).time && return DCTI(first(field).data) - time > last(field).time && return DCTI(last(field).data) - for i=reverse(1:length(field)) - isapprox(field[i].time, time) && return DCTI(field[i].data) - end - for i=reverse(2:length(field)) - t0 = field[i-1].time - t1 = field[i].time - if t0 < time < t1 - y0 = field[i-1].data - y1 = field[i].data - dt = t1-t0 - new_data = y0*(1-(time-t0)/dt) + y1*(1-(t1-time)/dt) - return DCTI(new_data) - end - end - error("interpolate DCTV: unknown failure when interpolating $(field.data) for time $time") -end - -function (field::DVTV)(time::Float64) - time < first(field).time && return DVTI(first(field).data) - time > last(field).time && return DVTI(last(field).data) - for i=reverse(1:length(field)) - isapprox(field[i].time, time) && return DVTI(field[i].data) - end - for i=reverse(2:length(field)) - t0 = field[i-1].time - t1 = field[i].time - if t0 < time < t1 - y0 = field[i-1].data - y1 = field[i].data - dt = t1-t0 - new_data = y0*(1-(time-t0)/dt) + y1*(1-(t1-time)/dt) - return DVTI(new_data) - end - end - error("interpolate DVTV: unknown failure when interpolating $(field.data) for time $time") -end - -""" Interpolate constant field in spatial dimension. """ -function (basis::CVTI)(field::DCTI, xi::Vector) - return field.data -end - -""" Interpolate variable field in spatial dimension. """ -function (basis::CVTI)(values::DVTI, xi::Vector) - N = basis(xi) - return sum([N[i]*values[i] for i=1:length(N)]) -end - -function (basis::CVTI)(geometry::DVTI, xi::Vector, ::Type{Val{:grad}}) - dbasis = basis(xi, Val{:grad}) -# J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)]) - J = sum([kron(dbasis[:,i], geometry[i]') for i=1:length(geometry)]) - invJ = isa(J, Vector) ? inv(J[1]) : inv(J) - grad = invJ * dbasis - return grad -end - -function (basis::CVTI)(geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}}) - grad = basis(geometry, xi, Val{:grad}) -# gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])' - gradf = sum([kron(grad[:,i], values[i]') for i=1:length(values)])' - return length(gradf) == 1 ? gradf[1] : gradf -end - -function (basis::CVTI)(xi::Vector, time::Number) - basis(xi) -end - -function Base.:*(grad::Matrix, field::DVTI) - n, m = size(grad) - return sum([kron(grad[:,i], field[i]') for i=1:m])' -end - -function DVTV(data::Pair{Float64, Vector}...) - return DVTV([Increment(d[1], d[2]) for d in data]) -end - -function start(f::DVTV) - return start(f.data) -end - -function next(f::DVTV, state) - return next(f.data, state) -end - -function done(f::DVTV, state) - return done(f.data, state) -end - -""" Return time vector from time variable field. """ -function keys(field::DVTV) - return Float64[increment.time for increment in field] -end - -function setindex!(field::Field, val, idx::Int64) - field.data[idx] = val -end diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl index 3e76dc2..d4c4016 100644 --- a/src/postprocess_utils.jl +++ b/src/postprocess_utils.jl @@ -26,8 +26,8 @@ function calc_nodal_values!(elements::Vector, field_name, field_dim, time; add!(A, gdofs, gdofs, w*kron(N', N)) end end - nz = get_nonzero_rows(A) A = sparse(A) + nz = get_nonzero_rows(A) A = 1/2*(A + A') F = ldltfact(A[nz,nz]) end diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl index 69bdf5c..07c324f 100644 --- a/src/preprocess_aster_reader.jl +++ b/src/preprocess_aster_reader.jl @@ -33,27 +33,6 @@ function aster_parse_nodes(section; strip_characters=true) return nodes end -function parse(mesh, ::Type{Val{:CODE_ASTER_MAIL}}) - model = Dict() - header = nothing - data = [] - for line in split(mesh, '\n') - length(line) != 0 || continue - info("line: $line") - if is_aster_mail_keyword(strip(line)) - header = parse_aster_header(line) - empty!(data) - continue - end - if line == "FINSF" - info(data) - header = nothing - process_aster_section!(model, join(data, ""), header, Val{header[1]}) - end - end - return model -end - """ Code Aster binary file (.med). """ type MEDFile @@ -161,7 +140,7 @@ Returns Dict containing fields "nodes" and "connectivity". """ -function parse_aster_med_file(fn, mesh_name=nothing; debug=false) +function parse_aster_med_file(fn, mesh_name=nothing) med = MEDFile(fn) mesh_names = get_mesh_names(med::MEDFile) all_meshes = join(mesh_names, ", ") @@ -173,12 +152,10 @@ function parse_aster_med_file(fn, mesh_name=nothing; debug=false) end nsets = get_node_sets(med, mesh_name) elsets = get_element_sets(med, mesh_name) - if debug - elset_names = join(values(elsets), ", ") - info("Code Aster .med reader: found $(length(elsets)) element sets: $elset_names") - nset_names = join(values(nsets), ", ") - info("Code ASter .med reader: found $(length(nsets)) node sets: $nset_names") - end + elset_names = join(values(elsets), ", ") + debug("Code Aster .med reader: found $(length(elsets)) element sets: $elset_names") + nset_names = join(values(nsets), ", ") + debug("Code ASter .med reader: found $(length(nsets)) node sets: $nset_names") nodes = get_nodes(med, nsets, mesh_name) conn = get_connectivity(med, elsets, mesh_name) result = Dict("nodes" => nodes, "connectivity" => conn) diff --git a/src/problems.jl b/src/problems.jl index 28ab40c..1e2a5f3 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -81,10 +81,6 @@ function isempty(assembly::Assembly) return T end -function get_dofs(assembly::Assembly) - return sort(unique(assembly.K.J)) -end - type Problem{P<:AbstractProblem} name :: AbstractString # descriptive name for problem dimension :: Int # degrees of freedom per node @@ -359,10 +355,6 @@ function push!(problem::Problem, elements_::Vector...) end end -function get_connectivity(problem::Problem) - return union([get_connectivity(element) for element in get_elements(problem)]...) -end - function get_gdofs(element::Element, dim::Int) conn = get_connectivity(element) if length(conn) == 0 @@ -372,10 +364,6 @@ function get_gdofs(element::Element, dim::Int) return gdofs end -function get_dofs(problem::Problem) - return get_dofs(problem.assembly) -end - function empty!(problem::Problem) empty!(problem.assembly) end @@ -396,33 +384,3 @@ function get_gdofs(problem::Problem, element::Element) end return problem.dofmap[element] end - -""" Find dofs corresponding to nodes. """ -function find_dofs_by_nodes(problem::Problem, nodes) - dim = get_unknown_field_dimension(problem) - return find_dofs_by_nodes(dim, nodes) -end -function find_dofs_by_nodes(dim::Int, nodes) - dofs = Int64[] - for node in nodes - for j=1:dim - push!(dofs, dim*(node-1)+j) - end - end - return dofs -end - -""" Find nodes corresponding to dofs. """ -function find_nodes_by_dofs(problem::Problem, dofs) - dim = get_unknown_field_dimension(problem) - return find_nodes_by_dofs(dim, dofs) -end -function find_nodes_by_dofs(dim, dofs) - nodes = Int64[] - for dof in dofs - j = Int(ceil(dof/dim)) - j in nodes && continue - push!(nodes, j) - end - return nodes -end diff --git a/src/problems_contact_2d.jl b/src/problems_contact_2d.jl index 92d4efb..4f9f9b9 100644 --- a/src/problems_contact_2d.jl +++ b/src/problems_contact_2d.jl @@ -287,273 +287,3 @@ function assemble!(problem::Problem{Contact}, time::Float64, problem.assembly.g = g end - - -""" -Frictionless 2d small sliding contact without forwarddiff. - -true/false flags: finite_sliding, friction, use_forwarddiff -""" -function _assemble!(problem::Problem{Contact}, time::Float64, - ::Type{Val{1}}, ::Type{Val{false}}, - ::Type{Val{false}}, ::Type{Val{false}}; debug=false) - - props = problem.properties - field_dim = get_unknown_field_dimension(problem) - field_name = get_parent_field_name(problem) - slave_elements = get_slave_elements(problem) - - # 1. calculate nodal normals and tangents for slave element nodes j ∈ S - normals, tangents = calculate_normals(slave_elements, time, Val{1}; - rotate_normals=props.rotate_normals) - update!(slave_elements, "normal", time => normals) - update!(slave_elements, "tangent", time => tangents) - - Rn = 0.0 - - # 2. loop all slave elements - for slave_element in slave_elements - - nsl = length(slave_element) - X1 = slave_element("geometry", time) - u1 = slave_element("displacement", time) - la1 = slave_element("reaction force", time) - n1 = slave_element("normal", time) - t1 = slave_element("tangent", time) - x1 = X1 + u1 - Q1_ = [n1[1] t1[1]] - Q2_ = [n1[2] t1[2]] - Z = zeros(2, 2) - Q2 = [Q1_ Z; Z Q2_] - contact_area = 0.0 - contact_error = 0.0 - - if "element area" in props.store_fields - element_area = 0.0 - for ip in get_integration_points(slave_element) - detJ = slave_element(ip, time, Val{:detJ}) - w = ip.weight*detJ - element_area += w - end - update!(slave_element, "element area", time => element_area) - end - - # 3. loop all master elements - for master_element in slave_element("master elements", time) - - nm = length(master_element) - X2 = master_element("geometry", time) - u2 = master_element("displacement", time) - x2 = X2 + u2 - - if norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) > props.distval - continue - end - - if norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) > props.distval - continue - end - - # 3.1 calculate segmentation - xi1a = project_from_master_to_slave(slave_element, X2[1], time) - xi1b = project_from_master_to_slave(slave_element, X2[2], time) - xi1 = clamp([xi1a; xi1b], -1.0, 1.0) - l = 1/2*abs(xi1[2]-xi1[1]) - isapprox(l, 0.0) && continue # no contribution in this master element - - # 3.2. bi-orthogonal basis - De = zeros(nsl, nsl) - Me = zeros(nsl, nsl) - Ae = zeros(nsl, nsl) - if props.dual_basis - for ip in get_integration_points(slave_element, 3) - detJ = slave_element(ip, time, Val{:detJ}) - w = ip.weight*detJ*l - xi = ip.coords[1] - xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) - N1 = vec(get_basis(slave_element, xi_s, time)) - De += w*diagm(N1) - Me += w*N1*N1' - end - Ae = De*inv(Me) - else - Ae = eye(nsl) - end - - # 3.3. loop integration points of one integration segment and calculate - # local mortar matrices - fill!(De, 0.0) - fill!(Me, 0.0) - ge = zeros(field_dim*nsl) - for ip in get_integration_points(slave_element, 3) - detJ = slave_element(ip, time, Val{:detJ}) - w = ip.weight*detJ*l - xi = ip.coords[1] - xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) - N1 = vec(get_basis(slave_element, xi_s, time)) - Phi = Ae*N1 - - # project gauss point from slave element to master element in direction n_s - X_s = N1*X1 # coordinate in gauss point - n_s = N1*n1 # normal direction in gauss point - t_s = N1*t1 # tangent condition in gauss point - n_s /= norm(n_s) - t_s /= norm(t_s) - xi_m = project_from_slave_to_master(master_element, X_s, n_s, time) - N2 = vec(get_basis(master_element, xi_m, time)) - X_m = N2*X2 - - u_s = N1*u1 - u_m = N2*u2 - x_s = X_s + u_s - x_m = X_m + u_m - la_s = Phi*la1 - ge += w*vec((x_m-x_s)*Phi') - - # virtual work - De += w*Phi*N1' - Me += w*Phi*N2' - - contact_area += w - contact_error += 1/2*w*dot(n_s, x_s-x_m)^2 - end - - sdofs = get_gdofs(problem, slave_element) - mdofs = get_gdofs(problem, master_element) - - # add contribution to contact virtual work - D2 = zeros(field_dim*nsl, field_dim*nsl) - M2 = zeros(field_dim*nsl, field_dim*nsl) - for i=1:field_dim - D2[i:field_dim:end, i:field_dim:end] += De - M2[i:field_dim:end, i:field_dim:end] += Me - end - - add!(problem.assembly.C1, sdofs, sdofs, D2) - add!(problem.assembly.C1, sdofs, mdofs, -M2) - add!(problem.assembly.C2, sdofs, sdofs, Q2'*D2) - add!(problem.assembly.C2, sdofs, mdofs, -Q2'*M2) - ge = -D2*vec(x1)+M2*vec(x2) - add!(problem.assembly.g, sdofs, Q2'*ge) - ce = vec(la1) + ge - add!(problem.assembly.c, sdofs, Q2'*ce) - - end # master elements done - - if "contact area" in props.store_fields - update!(slave_element, "contact area", time => contact_area) - end - - if "contact error" in props.store_fields - update!(slave_element, "contact error", time => contact_error) - end - - end # slave elements done, contact virtual work ready - - S = sort(collect(keys(normals))) # slave element nodes - weighted_gap = Dict{Int64, Vector{Float64}}() - contact_pressure = Dict{Int64, Vector{Float64}}() - complementarity_condition = Dict{Int64, Vector{Float64}}() - is_active = Dict{Int64, Int}() - is_inactive = Dict{Int64, Int}() - is_slip = Dict{Int64, Int}() - is_stick = Dict{Int64, Int}() - - la = problem.assembly.la - ndofs = length(la) - - C1 = sparse(problem.assembly.C1) - C2 = sparse(problem.assembly.C2, ndofs, ndofs) - D = spzeros(ndofs, ndofs) - c = full(problem.assembly.c, ndofs, 1) - g = full(problem.assembly.g, ndofs, 1) - - # active / inactive node detection - for j in S - dofs = [2*(j-1)+1, 2*(j-1)+2] - weighted_gap[j] = g[dofs] - - if length(la) != 0 - p = dot(normals[j], la[dofs]) - t = dot(tangents[j], la[dofs]) - contact_pressure[j] = [p, t] - else - contact_pressure[j] = [0.0, 0.0] - end - - #complementarity_condition[j] = contact_pressure[j] - weighted_gap[j] - complementarity_condition[j] = c[dofs] - if complementarity_condition[j][1] < 0 - is_inactive[j] = 1 - is_active[j] = 0 - is_slip[j] = 0 - is_stick[j] = 0 - else - is_inactive[j] = 0 - is_active[j] = 1 - is_slip[j] = 1 - is_stick[j] = 0 - end - end - - if "weighted gap" in props.store_fields - update!(slave_elements, "weighted gap", time => weighted_gap) - end - if "contact pressure" in props.store_fields - update!(slave_elements, "contact pressure", time => contact_pressure) - end - if "complementarity condition" in props.store_fields - update!(slave_elements, "complementarity condition", time => complementarity_condition) - end - if "active nodes" in props.store_fields - update!(slave_elements, "active nodes", time => is_active) - end - if "inactive nodes" in props.store_fields - update!(slave_elements, "inactive nodes", time => is_inactive) - end - if "stick nodes" in props.store_fields - update!(slave_elements, "stick nodes", time => is_stick) - end - if "slip nodes" in props.store_fields - update!(slave_elements, "slip nodes", time => is_slip) - end - -# info("# | active | inactive | stick | slip | gap | pres | comp") -# for j in S -# str1 = "$j | $(is_active[j]) | $(is_inactive[j]) | $(is_stick[j]) | $(is_slip[j]) | " -# str2 = "$(round(weighted_gap[j], 3)) | $(round(contact_pressure[j], 3)) | $(round(complementarity_condition[j], 3))" -# info(str1 * str2) -# end - - # solve variational inequality - - # constitutive modelling in tangent direction, frictionless contact - for j in S - dofs = [2*(j-1)+1, 2*(j-1)+2] - if (is_active[j] == 1) && (is_slip[j] == 1) -# info("$j is in active/slip, removing tangential constraint $(dofs[2])") - C2[dofs[2],:] = 0.0 - g[dofs[2]] = 0.0 - D[dofs[2], dofs] = tangents[j] - end - end - - # remove inactive nodes from assembly - for j in S - dofs = [2*(j-1)+1, 2*(j-1)+2] - if is_inactive[j] == 1 -# info("$j is inactive, removing dofs $dofs") - C1[dofs,:] = 0.0 - C2[dofs,:] = 0.0 - D[dofs,:] = 0.0 - g[dofs,:] = 0.0 - end - end - - problem.assembly.C1 = C1 - problem.assembly.C2 = C2 - problem.assembly.D = D - problem.assembly.g = g - -end - diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 1cced9a..b01184a 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -304,151 +304,6 @@ function assemble{El<:Elasticity2DSurfaceElements}(problem::Problem{Elasticity}, return Km, Kg, f end -""" Elasticity equations, 3d, linear. """ -function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}}) - - props = problem.properties - dim = get_unknown_field_dimension(problem) - nnodes = length(element) - ndofs = dim*nnodes - BL = zeros(6, ndofs) - Km = zeros(ndofs, ndofs) - Kg = zeros(ndofs, ndofs) - f = zeros(ndofs) - - for ip in get_integration_points(element) - detJ = element(ip, time, Val{:detJ}) - w = ip.weight*detJ - N = element(ip, time) - dN = element(ip, time, Val{:Grad}) - - fill!(BL, 0.0) - for i=1:nnodes - BL[1, 3*(i-1)+1] = dN[1,i] - BL[2, 3*(i-1)+2] = dN[2,i] - BL[3, 3*(i-1)+3] = dN[3,i] - BL[4, 3*(i-1)+1] = dN[2,i] - BL[4, 3*(i-1)+2] = dN[1,i] - BL[5, 3*(i-1)+2] = dN[3,i] - BL[5, 3*(i-1)+3] = dN[2,i] - BL[6, 3*(i-1)+1] = dN[3,i] - BL[6, 3*(i-1)+3] = dN[1,i] - end - - E = element("youngs modulus", ip, time) - nu = element("poissons ratio", ip, time) - - D = E/((1.0+nu)*(1.0-2.0*nu)) * [ - 1.0-nu nu nu 0.0 0.0 0.0 - nu 1.0-nu nu 0.0 0.0 0.0 - nu nu 1.0-nu 0.0 0.0 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 0.0 - 0.0 0.0 0.0 0.0 0.0 0.5-nu] - - Km += w*BL'*D*BL - - if haskey(element, "displacement load") - T = element("displacement load", ip, time) - f += w*vec(T*N) - end - for i=1:dim - if haskey(element, "displacement load $i") - b = element("displacement load $i", ip, time) - f[i:dim:end] += w*vec(b*N) - end - end - - end - - if get_formulation_type(problem) == :incremental - if haskey(element, "displacement") - u = vec(element["displacement"](time)) - f -= Kt*u - end - end - - return Km, Kg, f -end - -""" Material and geometric stiffness for linear buckling analysis. """ -function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_buckling}}) - - props = problem.properties - dim = get_unknown_field_dimension(problem) - nnodes = length(element) - ndofs = dim*nnodes - BL = zeros(6, ndofs) - BNL = zeros(9, ndofs) - Km = zeros(ndofs, ndofs) - Kg = zeros(ndofs, ndofs) - f = zeros(ndofs) - - for ip in get_integration_points(element) - detJ = element(ip, time, Val{:detJ}) - w = ip.weight*detJ - N = element(ip, time) - dN = element(ip, time, Val{:Grad}) - - gradu = element("displacement", ip, time, Val{:Grad}) - strain = 1/2*(gradu' + gradu) - - fill!(BL, 0.0) - for i=1:nnodes - BL[1, 3*(i-1)+1] = dN[1,i] - BL[2, 3*(i-1)+2] = dN[2,i] - BL[3, 3*(i-1)+3] = dN[3,i] - BL[4, 3*(i-1)+1] = dN[2,i] - BL[4, 3*(i-1)+2] = dN[1,i] - BL[5, 3*(i-1)+2] = dN[3,i] - BL[5, 3*(i-1)+3] = dN[2,i] - BL[6, 3*(i-1)+1] = dN[3,i] - BL[6, 3*(i-1)+3] = dN[1,i] - end - - fill!(BNL, 0.0) - for i=1:size(dN, 2) - BNL[1, 3*(i-1)+1] = dN[1,i] - BNL[2, 3*(i-1)+1] = dN[2,i] - BNL[3, 3*(i-1)+1] = dN[3,i] - BNL[4, 3*(i-1)+2] = dN[1,i] - BNL[5, 3*(i-1)+2] = dN[2,i] - BNL[6, 3*(i-1)+2] = dN[3,i] - BNL[7, 3*(i-1)+3] = dN[1,i] - BNL[8, 3*(i-1)+3] = dN[2,i] - BNL[9, 3*(i-1)+3] = dN[3,i] - end - - E = element("youngs modulus", ip, time) - nu = element("poissons ratio", ip, time) - D = E/((1.0+nu)*(1.0-2.0*nu)) * [ - 1.0-nu nu nu 0.0 0.0 0.0 - nu 1.0-nu nu 0.0 0.0 0.0 - nu nu 1.0-nu 0.0 0.0 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 0.0 - 0.0 0.0 0.0 0.0 0.0 0.5-nu] - - 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) - - S3 = zeros(3*dim, 3*dim) - S3[1,1] = stress_vec[1] - S3[2,2] = stress_vec[2] - S3[3,3] = stress_vec[3] - S3[1,2] = S3[2,1] = stress_vec[4] - S3[2,3] = S3[3,2] = stress_vec[5] - S3[1,3] = S3[3,1] = stress_vec[6] - S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3] - - Km += w*BL'*D*BL - Kg += w*BNL'*S3*BNL - - end - - return Km, Kg, f -end - """ Elasticity equations, 3d nonlinear. """ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}}) props = problem.properties @@ -679,185 +534,3 @@ function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, end return Km, Kg, f end - -function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}}) - return assemble(problem, element, time, Val{:continuum}) -end - -""" Elasticity equations using ForwardDiff -""" -function assemble(problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:forwarddiff}}) - - dim = get_unknown_field_dimension(problem) - nnodes = size(element, 2) - - function get_residual_vector(u::Vector) - u = reshape(u, dim, nnodes) - u = Field([u[:,i] for i=1:nnodes]) - r = zeros(dim, nnodes) - - for ip in get_integration_points(element) - - JT = transpose(get_jacobian(element, ip, time)) - n, m = size(JT) - if n == m - w = ip.weight*det(JT) - elseif m == 1 - w = ip.weight*norm(JT) - elseif m == 2 - w = ip.weight*norm(cross(JT[:,1], JT[:,2])) - else - error("jacobian $JT") - end - - # calculate internal forces - if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") - grad = element(ip, time, Val{:grad}) - gradu = grad*u - - # kinematics - F = I + gradu - E = 1/2*(F'*F - I) - - # material - young = element("youngs modulus", ip, time) - poisson = element("poissons ratio", ip, time) - mu = young/(2*(1+poisson)) - lambda = young*poisson/((1+poisson)*(1-2*poisson)) - if problem.properties.formulation == :plane_stress - lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for plane stress - end - - # stress - S = lambda*trace(E)*I + 2*mu*E - - r += w*F*S*grad - end - - # calculate external forces - volume load - if haskey(element, "displacement load") - basis = element(ip, time) - b = element("displacement load", ip, time) - r -= w*b*basis - end - - # external forces - surface traction force - if haskey(element, "displacement traction force") - basis = element(ip, time) - T = element("displacement traction force", ip, time) - r -= w*T*basis - end - - end - - return vec(r) - - end - - field = element("displacement", time) - Km, allresults = ForwardDiff.jacobian(get_residual_vector, vec(field), - AllResults, cache=autodiffcache) - Kg = zeros(Km) - f = -ForwardDiff.value(allresults) - return Km, Kg, f -end - - -############################### -# Plastic material # -############################### -#= -abstract PlaneStressLinearElasticPlasticProblem <: LinearElasticityProblem - -function PlaneStressLinearElasticPlasticProblem(name="plane stress linear elasticity", dim::Int=2, elements=[]) - return Problem{PlaneStressLinearElasticPlasticProblem}(name, dim, elements) -end - -""" Elasticity equations, plane stress. """ -function assemble!{E<:CG, P<:PlaneStressLinearElasticPlasticProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real) - - gdofs = get_gdofs(element, problem.dim) - ndim, nnodes = size(E) - B = zeros(3, 2*nnodes) - for ip in get_integration_points(element) - w = ip.weight - J = get_jacobian(element, ip, time) - N = element(ip, time) - if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") - nu = element("poissons ratio", ip, time) - E_ = element("youngs modulus", ip, time) - C = E_/(1.0 - nu^2) .* [ - 1.0 nu 0.0 - nu 1.0 0.0 - 0.0 0.0 (1.0-nu)/2.0] - dN = element(ip, time, Val{:grad}) - fill!(B, 0.0) - for i=1:size(dN, 2) - B[1, 2*(i-1)+1] = dN[1,i] - B[2, 2*(i-1)+2] = dN[2,i] - B[3, 2*(i-1)+1] = dN[2,i] - B[3, 2*(i-1)+2] = dN[1,i] - end - add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B*det(J)) - end - if haskey(element, "displacement load") - b = element("displacement load", ip, time) - add!(assembly.force_vector, gdofs, w*N'*b*det(J)) - end - if haskey(element, "displacement traction force") - T = element("displacement traction force", ip, time) - L = w*T*N*norm(J) - add!(assembly.force_vector, gdofs, vec(L)) - end - end -end - - - - - -include("elasticplastic.jl") - -# Elasticity problems -abstract ElasticityProblem <: AbstractProblem -abstract PlaneStressElasticityProblem <: ElasticityProblem - -function get_unknown_field_name{P<:ElasticityProblem}(::Type{P}) - return "displacement" -end - -function get_unknown_field_type{P<:ElasticityProblem}(::Type{P}) - return Vector{Float64} -end - - - -=# - -function (problem::Problem)(element::Element, ip, time::Float64, ::Type{Val{:E}}) - haskey(element, "displacement") || return nothing - gradu = element("displacement", ip, time, Val{:Grad}) - eps = 0.5*(gradu + gradu') - return eps -end - -function (problem::Problem)(element::Element, ip, time::Float64, ::Type{Val{:S}}) - haskey(element, "displacement") || return nothing - props = problem.properties - eps = problem(element, ip, time, Val{:E}) - eps == nothing && return nothing - E = element("youngs modulus", ip, time) - nu = element("poissons ratio", ip, time) - mu = E/(2.0*(1.0+nu)) - la = E*nu/((1.0+nu)*(1.0-2.0*nu)) - if props.formulation in [:plane_stress, :plane_strain] - la = 2.0*la*mu/(la+2.0*mu) - end - S = la*trace(eps)*I + 2.0*mu*eps - return S -end - -function (problem::Problem)(element::Element, ip, time::Float64, ::Type{Val{:COORD}}) - haskey(element, "geometry") || return nothing - return element("geometry", ip, time) -end diff --git a/src/problems_elasticplastic.jl b/src/problems_elasticplastic.jl deleted file mode 100644 index cc1b092..0000000 --- a/src/problems_elasticplastic.jl +++ /dev/null @@ -1,237 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -# Elasticity problems -abstract ElasticPlasticProblem <: AbstractProblem -abstract PlaneStressElasticPlasticProblem <: ElasticPlasticProblem - -function get_unknown_field_name{P<:ElasticPlasticProblem}(::Type{P}) - return "displacement" -end - -function get_unknown_field_type{P<:ElasticPlasticProblem}(::Type{P}) - return Vector{Float64} -end - -# 3D Elasticity problems -function ElasticPlasticProblem(dim::Int=3, elements=[]) - return Problem{ElasticPlasticProblem}("elasticplastic problem", dim, elements) -end - -# 2D Plane stress elasticity problems -function PlaneStressElasticPlasticProblem(dim::Int=2, elements=[]) - return Problem{PlaneStressElasticPlasticProblem}("plane stress elasticplastic problem", dim, elements) -end - - -function get_residual_vector{P<:PlaneStressElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing) - r = zeros(Float64, problem.dim, length(element)) - - J = get_jacobian(element, ip, time) - - # internal forces - if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") - if !haskey(element, "integration points") - if P == PlaneStressElasticPlasticProblem - last_stress = zeros(2,2) - last_strain = zeros(2,2) - else - last_stress = zeros(3,3) - last_strain = zeros(3,3) - end - else - for each_ip in element("integration points", time) - if isapprox(each_ip.xi, ip.xi) - last_stress = ip("stress", time) - last_strain = ip("stress", time) - break - end - end - end - u = element("displacement", time, variation) - grad = element(ip, time, Val{:grad}) - gradu = grad*u - - # deformation gradient - F = I + gradu - E = 1/2*(F'*F - I) - - #E = 1/2*(gradu + gradu') # finite strain (total) - - # material - young = element("youngs modulus", ip, time) - poisson = element("poissons ratio", ip, time) - stress_y = element("yield stress", time).data - dstrain = E - last_strain - de_v = [dstrain[1,1], dstrain[2,2], dstrain[1,2]] - material_model = element("material model", time) - s = last_stress - de = copy(ForwardDiff.get_value(dstrain)) - - if P == PlaneStressElasticPlasticProblem - C = stiffnessTensorPlaneStress(young, poisson) - s_v = [s[1,1], s[2,2], s[1,2]] - de_ = [de[1,1], de[2,2], de[1,2]] - problem_stress_type = :PlaneStressElasticPlasticProblem - else - C = stiffnessTensor(young, poisson) - s_v = [s[1,1], s[2,2], s[3,3], s[2,3], s[1,3], s[1,2]] - de_ = [de[1,1], de[2,2], de[3,3], de[2,3], de[1,3], de[1,2]] - problem_stress_type = :ElasticPlasticProblem - end - dep = zeros(3) - stress_inc, dep = calculate_stress(de_, - s_v, - C, - stress_y, - Val{:vonMises}, - Val{problem_stress_type}) - info("%% ", dep) - s_v += C * (de_v - dep) - info("--: ", ForwardDiff.get_value(s_v)) - # stress - if P == PlaneStressElasticPlasticProblem - S = [s_v[1] s_v[3]; - s_v[3] s_v[2]] - else - S = [s_v[1] s_v[6] s_v[5]; - s_v[6] s_v[2] s_v[4]; - s_v[5] s_v[4] s_v[3]] - end - r += F*S*grad*det(J) - end - - # external forces - volume load - if haskey(element, "displacement load") - basis = element(ip, time) - b = element("displacement load", ip, time) - r -= b*basis*det(J) - end - - # external forces - surface traction force - if haskey(element, "displacement traction force") - basis = element(ip, time) - T = element("displacement traction force", ip, time) - JT = transpose(J) - s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2]) - r -= T*basis*norm(s) - end - - return vec(r) -end - - - -#= -function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing) - r = zeros(Float64, problem.dim, length(element)) - - J = get_jacobian(element, ip, time) - - info("_____________________") - # internal forces - if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") - - if !haskey(element, "integration points") - if P == PlaneStressElasticPlasticProblem - last_stress = zeros(2,2) - last_strain = zeros(2,2) - else - last_stress = zeros(3,3) - last_strain = zeros(3,3) - end - else - for each_ip in element("integration points", time) - if isapprox(each_ip.xi, ip.xi) - last_stress = ip("stress", time) - last_strain = ip("stress", time) - break - end - end - end - u = element("displacement", time, variation) - grad = element(ip, time, Val{:grad}) - gradu = grad*u - - # deformation gradient - F = I + gradu - - # material - young = element("youngs modulus", ip, time) - poisson = element("poissons ratio", ip, time) - mu = young/(2*(1+poisson)) - lambda = young*poisson/((1+poisson)*(1-2*poisson)) - if P == PlaneStressElasticityProblem - lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d problems - end - - # strain - E = 1/2*(F'*F - I) - #E = 1/2*(gradu + gradu') # finite strain (total) - - young = element("youngs modulus", ip, time) - poisson = element("poissons ratio", ip, time) - stress_y = element("yield stress", time).data - dstrain = E - last_strain - material_model = element("material model", time) - s = last_stress - de = ForwardDiff.get_value(dstrain) - - if P == PlaneStressElasticPlasticProblem - C = stiffnessTensorPlaneStress(young, poisson) - s_v = [s[1,1], s[2,2], s[1,2]] - de_ = [de[1,1], de[2,2], de[1,2]] - problem_stress_type = :PlaneStressElasticPlasticProblem - else - C = stiffnessTensor(young, poisson) - s_v = [s[1,1], s[2,2], s[3,3], s[2,3], s[1,3], s[1,2]] - de_ = [de[1,1], de[2,2], de[3,3], de[2,3], de[1,3], de[1,2]] - problem_stress_type = :ElasticPlasticProblem - end - - stress_inc, lambda = plastic_multiplier = calculate_stress(de_, - s_v, - C, - stress_y, - Val{:vonMises}, - Val{problem_stress_type}) - - # dep = lambda * dfds(s) - # upate_material_parameters!(...) - s_new = s_v + stress_inc - #S = [s_v[1] s_v[6] s_v[5]; - # s_v[6] s_v[2] s_v[4]; -# s_v[5] s_v[4] s_v[3]] - S = [s_new[1] s_new[3]; - s_new[3] s_new[2]] - # S = C * (E - dep) - - - info("Stress: ", vec(ForwardDiff.get_value(S))) - # stress - #S = lambda*trace(E)*I + 2*mu*E - - r += F*S*grad*det(J) - - end - - - # external forces - volume load - if haskey(element, "displacement load") - basis = element(ip, time) - b = element("displacement load", ip, time) - r -= b*basis*det(J) - end - - # external forces - surface traction force - if haskey(element, "displacement traction force") - basis = element(ip, time) - T = element("displacement traction force", ip, time) - JT = transpose(J) - s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2]) - r -= T*basis*norm(s) - end - - return vec(r) -end -=# #fff diff --git a/src/problems_mortar_3d.jl b/src/problems_mortar_3d.jl index 21ed2c8..08d6d3f 100644 --- a/src/problems_mortar_3d.jl +++ b/src/problems_mortar_3d.jl @@ -34,15 +34,15 @@ function vertex_inside_polygon(q, P; atol=1.0e-3) cosa = dot(A,B)/c isapprox(cosa, 1.0; atol=atol) && return false isapprox(cosa, -1.0; atol=atol) && return true - try - angle += acos(cosa) - catch - info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.") - info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))") - info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))") - info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))") - rethrow() - end + #try + angle += acos(cosa) + #catch + # info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.") + # info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))") + # info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))") + # info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))") + # rethrow() + #end end return isapprox(angle, 2*pi; atol=atol) end @@ -68,26 +68,8 @@ function get_cells(P, C; allow_quads=false) if N == 4 && allow_quads return Vector[P] end - #V = sum([cross(P[i], P[mod(i,N)+1]) for i=1:N]) - #A = 1/2*abs(dot(n, V)) - #info("A = $A") cells = Vector[Vector[C, P[i], P[mod(i,N)+1]] for i=1:N] return cells - - maxa = 0.0 - maxj = 0 - for i=1:N - A = P[i] - C - B = P[mod(i,N)+1] - C - theta = acos(dot(A,B)/(norm(A)*norm(B))) - if theta > maxa - maxa = theta - maxj = i - end - end - info("max angle $(maxa/pi*180) at index $maxj, N=$N") - indices = mod(collect(maxj:maxj+N), N) - info("indices = $indices") end """ Test does P contain q. """ @@ -175,7 +157,7 @@ function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector, return theta[1:2], theta[3] end end - + #= info("failed to project vertex from auxiliary plane back to surface") info("element type: $E") info("element connectivity: $(get_connectivity(element))") @@ -199,7 +181,7 @@ function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector, info("dtheta = $(dtheta)") theta -= dtheta end - + =# throw(error("project_point_to_surface: did not converge in $max_iterations iterations!")) end @@ -272,8 +254,10 @@ function split_quadratic_element(element::Element{Tri6}, time::Float64) new_element = Element(Tri3, connectivity[elmap]) X = element("geometry", time) update!(new_element, "geometry", time => X[elmap]) - u = element("displacement", time) - update!(new_element, "displacement", time => u[elmap]) + if haskey(element, "displacement") + u = element("displacement", time) + update!(new_element, "displacement", time => u[elmap]) + end #n = element("normal", time) #update!(new_element, "normal", time => n[elmap]) if haskey(element, "master elements") @@ -372,6 +356,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ check_orientation!(P, n0) N_P = length(P) P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P]) + if first_slave_element debug("Polygon clip info for first slave element:") debug("S = $S") @@ -380,11 +365,15 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ debug("N_P = $N_P") debug("P_area = $P_area") end + if isapprox(P_area, 0.0) info("Polygon P has zero area: $P_area") continue end + C0 = calculate_centroid(P) + + #= if isnan(C0[1]) info("C0 = $C0") info("P = $P") @@ -393,6 +382,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ info("n0 = $n0") error("Calculation of centroid of polygon clip P failed.") end + =# De = zeros(nsl, nsl) Me = zeros(nsl, nm) @@ -412,7 +402,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ Me = zeros(nnodes, nnodes) for ip in get_integration_points(virtual_element, 3) x_gauss = nothing - try + #try x_gauss = virtual_element("geometry", ip, time) xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time) detJ = virtual_element(ip, time, Val{:detJ}) @@ -420,17 +410,17 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ N1 = vec(get_basis(slave_element, xi_s, time)) De += w*diagm(vec(N1)) Me += w*N1*N1' - catch - info("Failed to construct bi-orthogonal basis: cannot project vertex from auxiliary plane back to sufface.") - info("x_gauss = $x_gauss") - info("cell = $cell") - info("C0 = $C0") - info("P = $P") - info("S = $S") - info("M = $M") - info("n0 = $n0") - rethrow() - end + #catch + # info("Failed to construct bi-orthogonal basis: cannot project vertex from auxiliary plane back to sufface.") + # info("x_gauss = $x_gauss") + # info("cell = $cell") + # info("C0 = $C0") + # info("P = $P") + # info("S = $S") + # info("M = $M") + # info("n0 = $n0") + # rethrow() + #end end Ae = De*inv(Me) else @@ -449,6 +439,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ # project gauss point from auxiliary plane to master and slave element #x_gauss = N*x_cell x_gauss = virtual_element("geometry", ip, time) + #= if isnan(x_gauss[1]) info("is nan") info("x_gauss = $x_gauss") @@ -460,25 +451,26 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{ info("n0 = $n0") error("nan, unable to continue") end + =# xi_s = nothing xi_m = nothing alpha = nothing - try + #try xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time) xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time) - catch - info("projecting vertex back to surface has failed.") - info("x_gauss = $x_gauss") - info("cell = $cell") - info("C0 = $C0") - info("P = $P") - info("S = $S") - info("M = $M") - info("n0 = $n0") - rethrow() - end + #catch + # info("projecting vertex back to surface has failed.") + # info("x_gauss = $x_gauss") + # info("cell = $cell") + # info("C0 = $C0") + # info("P = $P") + # info("S = $S") + # info("M = $M") + # info("n0 = $n0") + # rethrow() + #end # add contributions N1 = vec(get_basis(slave_element, xi_s, time)) diff --git a/src/solvers.jl b/src/solvers.jl index 79b5709..b3aecb2 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -55,24 +55,6 @@ is_boundary_problem{P<:BoundaryProblem}(problem::Problem{P}) = true get_field_problems(solver::Solver) = filter(is_field_problem, get_problems(solver)) get_boundary_problems(solver::Solver) = filter(is_boundary_problem, get_problems(solver)) -""" -Posthook for field assembly. By default, do nothing. -This can be used to make some modifications for assembly -after all elements are assembled. - -Examples --------- -function field_assembly_posthook!(solver::Solver, - K::SparseMatrixCSC, - Kg::SparseMatrixCSC, - f::SparseMatrixCSC, - fg::SpareMatrixCSC) - info("doing stuff, size(K) = ", size(K)) -end -""" -function field_assembly_posthook! -end - """Return one combined field assembly for a set of field problems. Parameters @@ -120,12 +102,6 @@ function get_field_assembly(solver::Solver; show_info=true) f = sparse(f, solver.ndofs, 1) fg = sparse(fg, solver.ndofs, 1) - # run any posthook for assembly if defined - args = Tuple{Solver, SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC} - if method_exists(field_assembly_posthook!, args) - field_assembly_posthook!(solver, K, Kg, fg, fg) - end - return M, K, Kg, f, fg end @@ -191,9 +167,11 @@ function get_boundary_assembly(solver::Solver) g_ = sparse(assembly.g, ndofs, 1) for dof in assembly.removed_dofs info("$(problem.name): removing dof $dof from assembly") - C1_[:,dof] = 0.0 + C1_[dof,:] = 0.0 C2_[dof,:] = 0.0 end + SparseArrays.dropzeros!(C1_) + SparseArrays.dropzeros!(C2_) already_constrained = get_nonzero_rows(C2) new_constraints = get_nonzero_rows(C2_) @@ -203,8 +181,6 @@ function get_boundary_assembly(solver::Solver) warn("already constrained = $already_constrained") warn("new constraints = $new_constraints") overconstrained_dofs = sort(overconstrained_dofs) - overconstrained_nodes = find_nodes_by_dofs(problem, overconstrained_dofs) - warn("in overconstrained nodes $overconstrained_nodes") error("overconstrained dofs, not solving problem.") end @@ -224,10 +200,9 @@ Solve linear system using LDLt factorization (SuiteSparse). This version requires that final system is symmetric and positive definite, so boundary conditions are first eliminated before solution. """ -function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; debug=false) +function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}}) nnz(D) == 0 || return false - C1 == C2 || return false A = get_nonzero_rows(K) B = get_nonzero_rows(C2) @@ -235,28 +210,14 @@ function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; debug B == B2 || return false I = setdiff(A, B) - if debug - info("# A = $(length(A))") - info("# B = $(length(B))") - info("# I = $(length(I))") - end + debug("# A = $(length(A))") + debug("# B = $(length(B))") + debug("# I = $(length(I))") if length(B) == 0 warn("No rows in C2, forget to set Dirichlet boundary conditions to model?") else - # solver boundary dofs (usually a trivial solution Iu = g - try - u[B] = lufact(C2[B,B2]) \ full(g[B]) - catch - info("solver #1 failed to solve boundary dofs (you should not see this message).") - info("# A = $(length(A))") - info("# B = $(length(B))") - info("# B2 = $(length(B2))") - info("# I = $(length(I))") - info("B = $B") - info("B2 = $B2") - rethrow() - end + u[B] = lufact(C2[B,B2]) \ full(g[B]) end # solve interior domain using LDLt factorization @@ -287,12 +248,9 @@ function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{2}}) end """ Default linear system solver for solver. """ -function solve!(solver::Solver; empty_assemblies_before_solution=true, -show_info=true, symmetric=true, optimize=false, fill_D_diagonal=false) +function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric=true) - if show_info - info("Solving problems ...") - end + info("Solving problems ...") t0 = Base.time() # assemble field & boundary problems @@ -310,21 +268,10 @@ show_info=true, symmetric=true, optimize=false, fill_D_diagonal=false) M = 1/2*(M + M') end - if fill_D_diagonal - nz = ones(solver.ndofs) - nz[get_nonzero_rows(C2)] = 0.0 - nz[get_nonzero_rows(D)] = 0.0 - D += spdiagm(nz) - end - - # free up some memory before solution by either emptying field assemblies - # or combining values with same indices in sparse COO matrices. Small - # boundary problems are untouched. - for problem in get_field_problems(solver) - if empty_assemblies_before_solution + if empty_assemblies_before_solution + # free up some memory before solution by emptying field assemblies from problems + for problem in get_field_problems(solver) empty!(problem.assembly) - elseif optimize - optimize!(problem.assembly) end gc() end @@ -332,13 +279,17 @@ show_info=true, symmetric=true, optimize=false, fill_D_diagonal=false) ndofs = solver.ndofs u = zeros(ndofs) la = zeros(ndofs) - status = false + is_solved = false i = 0 for i in [1, 2] - status = solve!(solver, K, C1, C2, D, f, g, u, la, Val{i}) - status && break + is_solved = solve!(solver, K, C1, C2, D, f, g, u, la, Val{i}) + if is_solved + break + end + end + if !is_solved + error("Failed to solve linear system!") end - status || error("Failed to solve linear system!") t1 = round(Base.time()-t0, 2) norms = (norm(u), norm(la)) push!(solver.norms, norms) @@ -346,10 +297,8 @@ show_info=true, symmetric=true, optimize=false, fill_D_diagonal=false) solver.u = u solver.la = la - if show_info - info("Solved problems in $t1 seconds using solver $i.") - info("Solution norms = $norms.") - end + info("Solved problems in $t1 seconds using solver $i.") + info("Solution norms = $norms.") return end @@ -459,22 +408,6 @@ function get_all_elements(solver::Solver) return [elements...;] end -function get_element_type{E}(element::Element{E}) - return E -end - -function get_element_id{E}(element::Element{E}) - return element.id -end - -function is_element_type{E}(element::Element{E}, element_type) - return is(E, element_type) -end - -function filter_by_element_type(element_type, elements) - return filter(element -> is_element_type(element, element_type), elements) -end - function (solver::Solver)(field_name::AbstractString, time::Float64) fields = [] for problem in get_problems(solver) @@ -644,47 +577,27 @@ Notes ----- Default convergence criteria is obtained by checking each sub-problem convergence. """ -function has_converged(solver::Solver{Nonlinear}; show_info=false, - check_convergence_for_boundary_problems=false) +function has_converged(solver::Solver{Nonlinear}) properties = solver.properties converged = true eps = properties.convergence_tolerance - for problem in solver.problems - has_converged = true - if is_field_problem(problem) - has_converged = problem.assembly.u_norm_change < eps - if isapprox(norm(problem.assembly.u), 0.0) - # trivial solution - has_converged = true - end - show_info && info("Details for problem $(problem.name)") - show_info && info("Norm: $(norm(problem.assembly.u))") - show_info && info("Norm change: $(problem.assembly.u_norm_change)") - show_info && info("Has converged? $(has_converged)") - end - if is_boundary_problem(problem) && check_convergence_for_boundary_problems - has_converged = problem.assembly.la_norm_change/norm(problem.assembly.la) < eps - show_info && info("Details for problem $(problem.name)") - show_info && info("Norm: $(norm(problem.assembly.la))") - show_info && info("Norm change: $(problem.assembly.la_norm_change)") - show_info && info("Has converged? $(has_converged)") + for problem in get_field_problems(solver) + has_converged = problem.assembly.u_norm_change < eps + if isapprox(norm(problem.assembly.u), 0.0) + # trivial solution + has_converged = true end + debug("Details for problem $(problem.name)") + debug("Norm: $(norm(problem.assembly.u))") + debug("Norm change: $(problem.assembly.u_norm_change)") + debug("Has converged? $(has_converged)") converged &= has_converged end return converged end -type NonlinearConvergenceError <: Exception - solver :: Solver -end - -function Base.showerror(io::IO, exception::NonlinearConvergenceError) - max_iters = exception.solver.properties.max_iterations - print(io, "nonlinear iteration did not converge in $max_iters iterations!") -end - """ Default solver for quasistatic nonlinear problems. """ -function (solver::Solver{Nonlinear})(; show_info=true) +function (solver::Solver{Nonlinear})() properties = solver.properties @@ -693,10 +606,10 @@ function (solver::Solver{Nonlinear})(; show_info=true) # 2. start non-linear iterations for properties.iteration=1:properties.max_iterations - show_info && info(repeat("-", 80)) - show_info && info("Starting nonlinear iteration #$(properties.iteration)") - show_info && info("Increment time t=$(round(solver.time, 3))") - show_info && info(repeat("-", 80)) + info(repeat("-", 80)) + info("Starting nonlinear iteration #$(properties.iteration)") + info("Increment time t=$(round(solver.time, 3))") + info(repeat("-", 80)) # 2.1 update linearized assemblies assemble!(solver) @@ -714,7 +627,9 @@ function (solver::Solver{Nonlinear})(; show_info=true) end # 3. did not converge - properties.error_if_no_convergence && throw(NonlinearConvergenceError(solver)) + if properties.error_if_no_convergence + error("nonlinear iteration did not converge in $(properties.iteration) iterations!") + end end """ Convenience function to call nonlinear solver. """ @@ -851,9 +766,3 @@ function Postprocessor(problems::Problem...) end return solver end - -function Postprocessor(name::AbstractString, problems::Problem...) - solver = Postprocessor(problems...) - solver.name = name - return solver -end diff --git a/src/sparse.jl b/src/sparse.jl index 309f218..38ffaa0 100644 --- a/src/sparse.jl +++ b/src/sparse.jl @@ -32,10 +32,6 @@ function convert(::Type{SparseMatrixCOO}, A::Matrix) return SparseMatrixCOO(findnz(A)...) end -function convert(::Type{SparseMatrixCOO}, A::Vector) - return SparseMatrixCOO(findnz(sparse(A))...) -end - """ Convert from COO format to CSC. Parameters @@ -61,12 +57,6 @@ function empty!(A::SparseMatrixCOO) empty!(A.V) end -function append!(A::SparseMatrixCOO, I::Vector{Int}, J::Vector{Int}, V::Vector{Float64}) - append!(A.I, I) - append!(A.J, J) - append!(A.V, V) -end - function append!(A::SparseMatrixCOO, B::SparseMatrixCOO) append!(A.I, B.I) append!(A.J, B.J) @@ -77,17 +67,6 @@ function isempty(A::SparseMatrixCOO) return isempty(A.I) && isempty(A.J) && isempty(A.V) end -function Base.:+(A::SparseMatrixCOO, B::SparseMatrixCOO) - if isempty(A) - return B - end - if isempty(B) - return A - end - C = SparseMatrixCOO([A.I;B.I], [A.J;B.J], [A.V;B.V]) - return C -end - function full(A::SparseMatrixCOO, args...) return full(sparse(A.I, A.J, A.V, args...)) end @@ -156,6 +135,7 @@ function optimize!(A::SparseMatrixCOO) A.I = I A.J = J A.V = V + return end """ Find all nonzero rows from sparse matrix. @@ -173,20 +153,6 @@ function get_nonzero_columns(A::SparseMatrixCSC) return get_nonzero_rows(transpose(A)) end -function get_nonzero_rows(A::Union{SparseMatrixCOO, Matrix}) - return get_nonzero_rows(sparse(A)) -end - -function get_nonzero_columns(A::Union{SparseMatrixCOO, Matrix}) - return get_nonzero_columns(sparse(A)) -end - -function get_nonzeros(C::Union{SparseMatrixCSC, Matrix}) - nz1 = get_nonzero_rows(C) - nz2 = get_nonzero_columns(C) - return (nz1, nz2) -end - function size(A::SparseMatrixCOO) isempty(A) && return (0, 0) return maximum(A.I), maximum(A.J) @@ -206,20 +172,6 @@ function resize_sparsevec(b, n) return sparsevec(findnz(b)..., n) end -""" Matrix norm. Automatically convert to dense when asking for 2-norm for small matrices. """ -function norm(A::SparseMatrixCOO, p=Inf; maxdim=1000) - dim = size(A, 1) - if p == 2 && dim > maxdim - warn("Assembly norm: dim = $dim > $maxdim and p=$p, not making dense matrices for operation.") - return 0.0 - end - if p == 2 - return norm(full(A), p) - else - return norm(sparse(A), p) - end -end - """ Approximative comparison of two matricse A and B. """ function isapprox(A::SparseMatrixCOO, B::SparseMatrixCOO) A2 = sparse(A) diff --git a/src/types.jl b/src/types.jl index edcbb61..f0190e1 100644 --- a/src/types.jl +++ b/src/types.jl @@ -57,9 +57,3 @@ typealias IP Point{IntegrationPoint} function IP(id, weight, coords) return IP(id, weight, coords, Dict(), IntegrationPoint()) end - -function convert(::Type{IP}, data::Tuple{Float64, Vector{Float64}}) - weight, coords = data - return IP(-1, weight, coords) -end - diff --git a/test/test_elements.jl b/test/test_elements.jl index f2e10fa..4d4e2a1 100644 --- a/test/test_elements.jl +++ b/test/test_elements.jl @@ -124,3 +124,10 @@ end @test isa(lst, Vector) end +@testset "extend basis" begin + el = Element(Quad4, [1, 2, 3, 4]) + expected = [ + 0.25 0.00 0.25 0.00 0.25 0.00 0.25 0.00 + 0.00 0.25 0.00 0.25 0.00 0.25 0.00 0.25] + @test isapprox(el([0.0, 0.0], 0.0, 2), expected) +end diff --git a/test/test_elements_nurbs.jl b/test/test_elements_nurbs.jl new file mode 100644 index 0000000..328ee2a --- /dev/null +++ b/test/test_elements_nurbs.jl @@ -0,0 +1,28 @@ +# 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.Testing + +@testset "NSeg interpolate" begin + element = Element(NSeg, [1, 2]) + @test element([0.0], 0.0) == [0.5 0.5] + @test size(element) == (1, 2) + @test is_nurbs(element) + element2 = Element(Seg2, [1, 2]) + @test !is_nurbs(element2) +end + +@testset "NSurf interpolate" begin + element = Element(NSurf, [1, 2, 3, 4]) + @test element([0.0, 0.0], 0.0) == [0.25 0.25 0.25 0.25] + @test size(element) == (2, 4) + @test is_nurbs(element) +end + +@testset "NSolid interpolate" begin + element = Element(NSolid, [1, 2, 3, 4, 5, 6, 7, 8]) + @test element([0.0, 0.0, 0.0], 0.0) == [0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125] + @test size(element) == (3, 8) + @test is_nurbs(element) +end diff --git a/test/test_fields.jl b/test/test_fields.jl index 667a616..b845843 100644 --- a/test/test_fields.jl +++ b/test/test_fields.jl @@ -3,13 +3,83 @@ using JuliaFEM using JuliaFEM.Testing -using Logging -Logging.configure(level=DEBUG) -@testset "create and manipulate fields" begin +@testset "discrete, constant, time invariant field" begin + @test isa(DCTI(), DCTI) + @test DCTI(0.0).data == 0.0 + @test isa(Field(0.0), DCTI) + @test isa(Field(), DCTI) + f = DCTI() + update!(f, 1.0) + @test f.data == 1.0 + @test DCTI(1) == 1 + @test length(DCTI(1)) == 1 + @test f == DCTI(1.0) + @test isapprox(f, DCTI(1.0)) + @test isapprox(f, 1.0) + @test 2*f == 2.0 # multiply by constant + @test f(1.0) == 1.0 # time interpolation + @test isapprox([2.0]''*f, 2.0) # wanted behavior? +end -@testset "updating time dependent fields" begin +@testset "discrete, variable, time invariant field" begin + @test isa(DVTI(), DVTI) + @test DVTI([1.0, 2.0]).data == [1.0, 2.0] + @test isa(Field([1.0, 2.0]), DVTI) + + f = DVTI() + update!(f, [2.0, 3.0]) + @test isapprox(f.data, [2.0, 3.0]) + @test length(f) == 2 + + # slicing + @test isapprox(f[1], 2.0) + @test isapprox(f[[1, 2]], [2.0, 3.0]) + + # boolean comparison and multiplying by a constant + @test f == DVTI([2.0, 3.0]) + @test isapprox(2*f, [4.0, 6.0]) + + f3 = 2*f + @test isa(f3, DVTI) + @test f3+f == 3*f + @test f3-f == f + + # spatial interpolation + N = [1.0, 2.0] + @test isapprox(N*f, 8.0) + + # time interpolation + @test isapprox(f(1.0), [2.0, 3.0]) + + # spatial interpolation of vector valued variable field + f2 = DVTI(Vector[[1.0, 2.0], [3.0, 4.0]]) + @test isapprox(f2[1], [1.0, 2.0]) + @test isapprox(f2[2], [3.0, 4.0]) + @test length(f2) == 2 + @test isapprox(N*f2, [1.0, 2.0] + [6.0, 8.0]) + + # iteration of DVTI field + s = zeros(2) + for j in f2 + s += j + end + @test isapprox(s, [4.0, 6.0]) + + @test vec(f2) == [1.0, 2.0, 3.0, 4.0] + @test isapprox([1.0 2.0]*f, [8.0]'') + + new_data = [2.0, 3.0, 4.0, 5.0] + f4 = similar(f2, new_data) + @test isa(f4, DVTI) + @test isapprox(f4.data[1], [2.0, 3.0]) + @test isapprox(f4.data[2], [4.0, 5.0]) +end + +@testset "discrete, constant, time-variant field" begin + @test isa(DCTV(), DCTV) f = Field(0.0 => 1.0) + @test isa(f, DCTV) @test last(f).time == 0.0 @test last(f).data == 1.0 update!(f, 0.0 => 2.0) @@ -20,27 +90,81 @@ Logging.configure(level=DEBUG) @test last(f).time == 1.0 @test last(f).data == 3.0 @test length(f) == 2 + + @testset "interpolation in time direction" begin + @test isa(f(0.0), DCTI) # converts to time-invariant after time interpolation + @test isapprox(f(-1.0), 2.0) + @test isapprox(f(0.0), 2.0) + @test isapprox(f(0.5), 2.5) + @test isapprox(f(1.0), 3.0) + @test isapprox(f(2.0), 3.0) + end + + # create several time steps at once + f = DCTV(0.0 => 1.0, 1.0 => 2.0) + @test isapprox(f(0.5), 1.5) + end -@testset "updating time invariant fields" begin - f = Field(1.0) - @test f.data == 1.0 - update!(f, 2.0) - @test f.data == 2.0 +@testset "discrete, variable, time-variant field" begin + @test isa(DVTV(), DVTV) + f = Field(0.0 => [1.0, 2.0]) + @test isa(f, DVTV) + @test last(f).time == 0.0 + @test last(f).data == [1.0, 2.0] + update!(f, 0.0 => [2.0, 3.0]) + @test last(f).time == 0.0 + @test last(f).data == [2.0, 3.0] + @test length(f) == 1 + update!(f, 1.0 => [3.0, 4.0]) + @test last(f).time == 1.0 + @test last(f).data == [3.0, 4.0] + @test length(f) == 2 + + @testset "interpolation in time direction" begin + @test isa(f(0.0), DVTI) # converts to time-invariant after time interpolation + @test isapprox(f(-1.0), [2.0, 3.0]) + @test isapprox(f(0.0), [2.0, 3.0]) + @test isapprox(f(0.5), [2.5, 3.5]) + @test isapprox(f(1.0), [3.0, 4.0]) + @test isapprox(f(2.0), [3.0, 4.0]) + end + + # create several time steps at once + f = DVTV(0.0 => [1.0, 2.0], 1.0 => [2.0, 3.0]) + @test isapprox(f(0.5), [1.5, 2.5]) end -@testset "field defined using function" begin - g(xi, t) = xi[1]*t - f = Field(g) - v = f([1.0], 2.0) - @test isapprox(v, 2.0) +@testset "continuous, constant, time-invariant field" begin + f = Field(() -> 2.0) + @test isapprox(f([1.0], 2.0), 2.0) + +end + +@testset "continuous, constant, time variant field" begin + f = Field((time::Float64) -> 2.0*time) + @test isapprox(f([1.0], 2.0), 4.0) + +end + +@testset "continuous, variable, time invariant field" begin + f = Field((xi::Vector) -> sum(xi)) + @test isapprox(f([1.0, 2.0], 2.0), 3.0) +end + +@testset "continuous, variable, time variant field" begin + f = Field((xi::Vector, t::Float64) -> xi[1]*t) + @test isapprox(f([1.0], 2.0), 2.0) +end + +@testset "unknown function argument for continuous field" begin + @test_throws ErrorException Field((a, b, c) -> a*b*c) end @testset "dictionary fields" begin f1 = Dict{Int64, Vector{Float64}}(1 => [0.0, 0.0], 2 => [0.0, 0.0]) f2 = Dict{Int64, Vector{Float64}}(1 => [1.0, 1.0], 2 => [1.0, 1.0]) f = Field(0.0 => f1, 1.0 => f2) - debug("field = $f") @test isa(f, DVTV) @test isapprox(f(0.0)[1], [0.0, 0.0]) @test isapprox(f(1.0)[2], [1.0, 1.0]) @@ -58,5 +182,3 @@ end f = Field(f1) @test isa(f, DVTI) end - -end diff --git a/test/test_lagrange.jl b/test/test_lagrange.jl index 72f4e73..8fcb931 100644 --- a/test/test_lagrange.jl +++ b/test/test_lagrange.jl @@ -3,7 +3,6 @@ using JuliaFEM using JuliaFEM.Testing -using JuliaFEM: description ALL_ELEMENTS = [ Seg2, Seg3, @@ -14,6 +13,16 @@ ALL_ELEMENTS = [ Hex8, Hex20, Hex27 ] +info("basic data for elements implemented so far:") +for element_type in [Poi1; ALL_ELEMENTS] + element = Element(element_type, Int[]) + element_length = length(element) + element_size = size(element) + element_description = description(element) + info("Element $element_type, description = $element_description, length = $element_length, size = $element_size") +end + + ALL_ELEMENTS_NODES = [ [1,2], [1,2,3], [1,2,3], [1,2,3,4,5,6], [1,2,3,4,5,6,7], @@ -78,17 +87,3 @@ end @test length(el) == length(vec) end end - -DESC = ["2 node segment", "3 node segment", "3 node triangle", - "6 node triangle", "7 node triangle", "4 node quadrangle", - "8 node Serendip quadrangle", "9 node quadrangle", - "4 node tetrahedral element", "10 node tetrahedral element", - "6 node prismatic element (wedge)", - "8 node hexahedral element", "20 node hexahedral element", - "27 node hexahedral element"] - -@testset "element description" begin - for (T, res) in zip(ALL_ELEMENTS, DESC) - @test description(Type(T)) == res - end -end diff --git a/test/test_mortar_3d_mesh_tie.jl b/test/test_problems_mortar_3d.jl similarity index 97% rename from test/test_mortar_3d_mesh_tie.jl rename to test/test_problems_mortar_3d.jl index 5425ada..0e36c2d 100644 --- a/test/test_mortar_3d_mesh_tie.jl +++ b/test/test_problems_mortar_3d.jl @@ -5,6 +5,7 @@ using JuliaFEM using JuliaFEM.Preprocess using JuliaFEM.Postprocess using JuliaFEM.Testing +using JuliaFEM.Abaqus: create_surface_elements @testset "test that interface transfers constant field without error" begin meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_3d.med" diff --git a/test/test_solver.jl b/test/test_solvers.jl similarity index 58% rename from test/test_solver.jl rename to test/test_solvers.jl index 29e3fc1..99e275f 100644 --- a/test/test_solver.jl +++ b/test/test_solvers.jl @@ -47,3 +47,32 @@ using JuliaFEM.Testing info("Temperature at point X = $X is T = $T") @test isapprox(T, 100.0) end + +@testset "problem not found from solver" begin + s = Solver(Linear, "demo solver") + @test_throws KeyError getindex(s, "not_found") +end + +@testset "automatic determination of problem dimension if not spesified" begin + s = Solver(Linear, "demo solver") + p = Problem(Elasticity, "demo problem", 2) + push!(s, p) + get_field_assembly(s) + @test s.ndofs == 0 + add!(p.assembly.K, [4], [4], [4.0]'') + get_field_assembly(s) + @test s.ndofs == 4 +end + +@testset "test for error when overdetermined system and requesting boundary assembly" begin + s = Solver(Linear, "demo solver") + @test_throws AssertionError get_boundary_assembly(s) # ndofs = 0 + p1 = Problem(Dirichlet, "bc1", 2, "displacement") + p2 = Problem(Dirichlet, "bc2", 2, "displacement") + # third dofs constrained + add!(p1.assembly.C2, [3], [3], [1.0]'') + add!(p2.assembly.C2, [3], [4], [1.0]'') + s.ndofs = 4 + push!(s, p1, p2) + @test_throws ErrorException get_boundary_assembly(s) +end diff --git a/test/test_sparse.jl b/test/test_sparse.jl index bdf4515..edd0ecf 100644 --- a/test/test_sparse.jl +++ b/test/test_sparse.jl @@ -17,3 +17,31 @@ end add!(b, sparse(b2)) @test isapprox(full(b), full(b2)) end + +@testset "Failure to add data to sparse vector due dimensino mismatch" begin + b = SparseVectorCOO() + @test_throws ErrorException add!(b, [1, 2], [1.0, 2.0, 3.0]) +end + +@testset "Test combining of SparseMatrixCOO" begin + k = convert(Matrix{Float64}, reshape(collect(1:9), 3, 3)) + dofs1 = [1, 2, 3] + dofs2 = [2, 3, 4] + A = SparseMatrixCOO() + add!(A, dofs1, dofs1, k) + add!(A, dofs2, dofs2, k) + A1 = full(A) + optimize!(A) + A2 = full(A) + @test isapprox(A1, A2) +end + +@testset "resize of sparse matrix and sparse vector" begin + A = sparse(rand(3, 3)) + B = resize_sparse(A, 4, 4) + @test size(B) == (4, 4) + a = sparse(rand(3)) + b = resize_sparsevec(a, 4) + @test size(b) == (4, ) +end + diff --git a/test/testdata/test_problems_mortar_3d_tet10.inp b/test/testdata/test_problems_mortar_3d_tet10.inp new file mode 100755 index 0000000..6579ac9 --- /dev/null +++ b/test/testdata/test_problems_mortar_3d_tet10.inp @@ -0,0 +1,1354 @@ +**NSET COUNT = 607 +*NODE +1, 0.74534, 0.75046, 0.82163 +2, 0.62267, 0.75023, 0.91082 +3, 0.62267, 0.87523, 0.91082 +4, 0.74767, 0.75023, 0.91082 +5, 0.87500, 0.50000, 0.87500 +6, 0.74767, 0.62523, 0.91082 +7, 0.87267, 0.62523, 0.91082 +8, 0.87267, 0.62523, 0.78582 +9, 0.87267, 0.75023, 0.91082 +10, 0.87267, 0.75023, 0.78582 +11, 0.87500, 0.62500, 0.62500 +12, 0.87500, 0.75000, 0.62500 +13, 0.74767, 0.75023, 0.66082 +14, 0.66838, 0.54400, 0.72860 +15, 0.87500, 0.50000, 0.62500 +16, 0.70919, 0.52200, 0.61430 +17, 0.83419, 0.52200, 0.73930 +18, 0.70919, 0.64700, 0.61430 +19, 0.62500, 0.87500, 0.62500 +20, 0.50000, 0.87500, 0.62500 +21, 0.62267, 0.75023, 0.66082 +22, 0.62267, 0.87523, 0.78582 +23, 0.75000, 0.87500, 0.62500 +24, 0.74767, 0.87523, 0.78582 +25, 0.87500, 0.87500, 0.75000 +26, 0.87500, 0.87500, 0.87500 +27, 0.74767, 0.87523, 0.91082 +28, 0.58419, 0.52200, 0.86430 +29, 0.58419, 0.64700, 0.86430 +30, 0.70919, 0.52200, 0.86430 +31, 0.58419, 0.39700, 0.86430 +32, 0.70919, 0.39700, 0.86430 +33, 0.87500, 0.37500, 0.87500 +34, 0.87500, 0.25000, 0.87500 +35, 0.83419, 0.39700, 0.73930 +36, 0.58419, 0.39700, 0.61430 +37, 0.70919, 0.39700, 0.61430 +38, 0.58419, 0.52200, 0.61430 +39, 0.87500, 0.25000, 0.62500 +40, 0.87500, 0.37500, 0.62500 +41, 0.70686, 0.64723, 0.77512 +42, 0.87500, 0.87500, 0.62500 +43, 0.58419, 0.64700, 0.61430 +44, 0.75000, 0.12500, 0.62500 +45, 0.87500, 0.12500, 0.62500 +46, 0.87500, 0.12500, 0.87500 +47, 0.87500, 0.12500, 0.75000 +48, 0.37836, 0.85595, 0.74983 +49, 0.43918, 0.80298, 0.87491 +50, 0.56185, 0.80321, 0.78573 +51, 0.52337, 0.69998, 0.73922 +52, 0.43918, 0.67798, 0.87491 +53, 0.31418, 0.80298, 0.87491 +54, 0.31418, 0.67798, 0.87491 +55, 0.12500, 0.62500, 0.87500 +56, 0.12500, 0.75000, 0.87500 +57, 0.18918, 0.80298, 0.74991 +58, 0.12500, 0.50000, 0.87500 +59, 0.18918, 0.67798, 0.74991 +60, 0.43918, 0.67798, 0.62491 +61, 0.43918, 0.80298, 0.62491 +62, 0.43918, 0.92798, 0.74991 +63, 0.43918, 0.92798, 0.62491 +64, 0.31418, 0.92798, 0.62491 +65, 0.31418, 0.80298, 0.62491 +66, 0.31418, 0.67798, 0.62491 +67, 0.12500, 0.62500, 0.62500 +68, 0.12500, 0.50000, 0.62500 +69, 0.12500, 0.75000, 0.62500 +70, 0.12500, 0.87500, 0.62500 +71, 0.31418, 0.92798, 0.74991 +72, 0.31418, 0.92798, 0.87491 +73, 0.25000, 0.87500, 0.87500 +74, 0.43918, 0.92798, 0.87491 +75, 0.12500, 0.87500, 0.75000 +76, 0.62345, 0.16661, 0.75334 +77, 0.50000, 0.12500, 0.87500 +78, 0.56172, 0.08330, 0.87667 +79, 0.56172, 0.08330, 0.75167 +80, 0.56172, 0.20830, 0.87667 +81, 0.68672, 0.08330, 0.75167 +82, 0.68672, 0.08330, 0.87667 +83, 0.68672, 0.20830, 0.87667 +84, 0.64591, 0.35531, 0.74097 +85, 0.81172, 0.20830, 0.75167 +86, 0.68672, 0.20830, 0.62667 +87, 0.56172, 0.20830, 0.62667 +88, 0.68672, 0.08330, 0.62667 +89, 0.56172, 0.08330, 0.62667 +90, 0.22841, 0.22025, 0.65164 +91, 0.23920, 0.23513, 0.57582 +92, 0.36420, 0.11013, 0.57582 +93, 0.36420, 0.23513, 0.57582 +94, 0.41565, 0.41697, 0.69068 +95, 0.33283, 0.45849, 0.84534 +96, 0.39700, 0.63646, 0.72025 +97, 0.45783, 0.45849, 0.84534 +98, 0.20783, 0.45849, 0.72034 +99, 0.33283, 0.33349, 0.84534 +100, 0.45783, 0.33349, 0.84534 +101, 0.33283, 0.33349, 0.59534 +102, 0.45783, 0.33349, 0.59534 +103, 0.33283, 0.45849, 0.59534 +104, 0.51955, 0.29179, 0.72201 +105, 0.54202, 0.48049, 0.70964 +106, 0.45783, 0.45849, 0.59534 +107, 0.37500, 0.12500, 0.87500 +108, 0.25000, 0.12500, 0.87500 +109, 0.23920, 0.23513, 0.82582 +110, 0.36420, 0.11013, 0.70082 +111, 0.23920, 0.11013, 0.70082 +112, 0.12500, 0.12500, 0.75000 +113, 0.12500, 0.25000, 0.87500 +114, 0.11420, 0.23513, 0.70082 +115, 0.12500, 0.37500, 0.87500 +116, 0.11420, 0.36013, 0.70082 +117, 0.32203, 0.31861, 0.67116 +118, 0.36420, 0.23513, 0.82582 +119, 0.42593, 0.19343, 0.70249 +120, 0.12500, 0.12500, 0.87500 +121, 0.23920, 0.11013, 0.57582 +122, 0.12500, 0.12500, 0.62500 +123, 0.11420, 0.23513, 0.57582 +124, 0.11420, 0.36013, 0.57582 +125, 0.23920, 0.36013, 0.57582 +126, 0.12500, 0.87500, 0.87500 +127, 1.00000, 0.00000, 0.75000 +128, 1.00000, 0.00000, 1.00000 +129, 0.75000, 0.00000, 1.00000 +130, 1.00000, 0.00000, 0.87500 +131, 0.87500, 0.00000, 1.00000 +132, 0.87500, 0.00000, 0.87500 +133, 0.75000, 0.00000, 0.75000 +134, 0.75000, 0.00000, 0.87500 +135, 0.87500, 0.00000, 0.75000 +136, 1.00000, 0.00000, 0.50000 +137, 1.00000, 0.00000, 0.62500 +138, 0.87500, 0.00000, 0.62500 +139, 0.75000, 0.00000, 0.50000 +140, 0.75000, 0.00000, 0.62500 +141, 0.87500, 0.00000, 0.50000 +142, 0.50000, 0.00000, 1.00000 +143, 0.62500, 0.00000, 1.00000 +144, 0.62500, 0.00000, 0.87500 +145, 0.50000, 0.00000, 0.75000 +146, 0.50000, 0.00000, 0.87500 +147, 0.62500, 0.00000, 0.75000 +148, 0.62500, 0.00000, 0.62500 +149, 0.50000, 0.00000, 0.50000 +150, 0.50000, 0.00000, 0.62500 +151, 0.62500, 0.00000, 0.50000 +152, 0.25000, 0.00000, 1.00000 +153, 0.37500, 0.00000, 1.00000 +154, 0.37500, 0.00000, 0.87500 +155, 0.25000, 0.00000, 0.75000 +156, 0.25000, 0.00000, 0.87500 +157, 0.37500, 0.00000, 0.75000 +158, 0.37500, 0.00000, 0.62500 +159, 0.25000, 0.00000, 0.50000 +160, 0.25000, 0.00000, 0.62500 +161, 0.37500, 0.00000, 0.50000 +162, 0.00000, 0.00000, 1.00000 +163, 0.12500, 0.00000, 1.00000 +164, 0.12500, 0.00000, 0.87500 +165, 0.00000, 0.00000, 0.75000 +166, 0.00000, 0.00000, 0.87500 +167, 0.12500, 0.00000, 0.75000 +168, 0.12500, 0.00000, 0.62500 +169, 0.00000, 0.00000, 0.50000 +170, 0.00000, 0.00000, 0.62500 +171, 0.12500, 0.00000, 0.50000 +172, 0.75000, 1.00000, 1.00000 +173, 1.00000, 1.00000, 1.00000 +174, 1.00000, 1.00000, 0.75000 +175, 0.87500, 1.00000, 1.00000 +176, 1.00000, 1.00000, 0.87500 +177, 0.87500, 1.00000, 0.87500 +178, 0.75000, 1.00000, 0.75000 +179, 0.75000, 1.00000, 0.87500 +180, 0.87500, 1.00000, 0.75000 +181, 1.00000, 1.00000, 0.50000 +182, 1.00000, 1.00000, 0.62500 +183, 0.87500, 1.00000, 0.62500 +184, 0.75000, 1.00000, 0.50000 +185, 0.75000, 1.00000, 0.62500 +186, 0.87500, 1.00000, 0.50000 +187, 0.50000, 1.00000, 1.00000 +188, 0.62500, 1.00000, 1.00000 +189, 0.62500, 1.00000, 0.87500 +190, 0.50000, 1.00000, 0.75000 +191, 0.50000, 1.00000, 0.87500 +192, 0.62500, 1.00000, 0.75000 +193, 0.62500, 1.00000, 0.62500 +194, 0.50000, 1.00000, 0.50000 +195, 0.50000, 1.00000, 0.62500 +196, 0.62500, 1.00000, 0.50000 +197, 0.25000, 1.00000, 1.00000 +198, 0.37500, 1.00000, 1.00000 +199, 0.37500, 1.00000, 0.87500 +200, 0.25000, 1.00000, 0.75000 +201, 0.25000, 1.00000, 0.87500 +202, 0.37500, 1.00000, 0.75000 +203, 0.37500, 1.00000, 0.62500 +204, 0.25000, 1.00000, 0.50000 +205, 0.25000, 1.00000, 0.62500 +206, 0.37500, 1.00000, 0.50000 +207, 0.00000, 1.00000, 1.00000 +208, 0.12500, 1.00000, 1.00000 +209, 0.12500, 1.00000, 0.87500 +210, 0.00000, 1.00000, 0.75000 +211, 0.00000, 1.00000, 0.87500 +212, 0.12500, 1.00000, 0.75000 +213, 0.12500, 1.00000, 0.62500 +214, 0.00000, 1.00000, 0.50000 +215, 0.00000, 1.00000, 0.62500 +216, 0.12500, 1.00000, 0.50000 +217, 1.00000, 0.75000, 1.00000 +218, 1.00000, 0.87500, 1.00000 +219, 1.00000, 0.87500, 0.87500 +220, 1.00000, 0.75000, 0.75000 +221, 1.00000, 0.75000, 0.87500 +222, 1.00000, 0.87500, 0.75000 +223, 1.00000, 0.87500, 0.62500 +224, 1.00000, 0.75000, 0.50000 +225, 1.00000, 0.75000, 0.62500 +226, 1.00000, 0.87500, 0.50000 +227, 1.00000, 0.50000, 1.00000 +228, 1.00000, 0.62500, 1.00000 +229, 1.00000, 0.62500, 0.87500 +230, 1.00000, 0.50000, 0.75000 +231, 1.00000, 0.50000, 0.87500 +232, 1.00000, 0.62500, 0.75000 +233, 1.00000, 0.62500, 0.62500 +234, 1.00000, 0.50000, 0.50000 +235, 1.00000, 0.50000, 0.62500 +236, 1.00000, 0.62500, 0.50000 +237, 1.00000, 0.25000, 1.00000 +238, 1.00000, 0.37500, 1.00000 +239, 1.00000, 0.37500, 0.87500 +240, 1.00000, 0.25000, 0.75000 +241, 1.00000, 0.25000, 0.87500 +242, 1.00000, 0.37500, 0.75000 +243, 1.00000, 0.37500, 0.62500 +244, 1.00000, 0.25000, 0.50000 +245, 1.00000, 0.25000, 0.62500 +246, 1.00000, 0.37500, 0.50000 +247, 1.00000, 0.12500, 1.00000 +248, 1.00000, 0.12500, 0.87500 +249, 1.00000, 0.12500, 0.75000 +250, 1.00000, 0.12500, 0.62500 +251, 1.00000, 0.12500, 0.50000 +252, 0.00000, 0.75000, 1.00000 +253, 0.00000, 0.87500, 1.00000 +254, 0.00000, 0.87500, 0.87500 +255, 0.00000, 0.75000, 0.75000 +256, 0.00000, 0.75000, 0.87500 +257, 0.00000, 0.87500, 0.75000 +258, 0.00000, 0.87500, 0.62500 +259, 0.00000, 0.75000, 0.50000 +260, 0.00000, 0.75000, 0.62500 +261, 0.00000, 0.87500, 0.50000 +262, 0.00000, 0.50000, 1.00000 +263, 0.00000, 0.62500, 1.00000 +264, 0.00000, 0.62500, 0.87500 +265, 0.00000, 0.50000, 0.75000 +266, 0.00000, 0.50000, 0.87500 +267, 0.00000, 0.62500, 0.75000 +268, 0.00000, 0.62500, 0.62500 +269, 0.00000, 0.50000, 0.50000 +270, 0.00000, 0.50000, 0.62500 +271, 0.00000, 0.62500, 0.50000 +272, 0.00000, 0.25000, 1.00000 +273, 0.00000, 0.37500, 1.00000 +274, 0.00000, 0.37500, 0.87500 +275, 0.00000, 0.25000, 0.75000 +276, 0.00000, 0.25000, 0.87500 +277, 0.00000, 0.37500, 0.75000 +278, 0.00000, 0.37500, 0.62500 +279, 0.00000, 0.25000, 0.50000 +280, 0.00000, 0.25000, 0.62500 +281, 0.00000, 0.37500, 0.50000 +282, 0.00000, 0.12500, 1.00000 +283, 0.00000, 0.12500, 0.87500 +284, 0.00000, 0.12500, 0.75000 +285, 0.00000, 0.12500, 0.62500 +286, 0.00000, 0.12500, 0.50000 +287, 0.87500, 0.87500, 1.00000 +288, 0.75000, 0.75000, 1.00000 +289, 0.87500, 0.75000, 1.00000 +290, 0.75000, 0.87500, 1.00000 +291, 0.62500, 0.87500, 1.00000 +292, 0.50000, 0.75000, 1.00000 +293, 0.62500, 0.75000, 1.00000 +294, 0.50000, 0.87500, 1.00000 +295, 0.37500, 0.87500, 1.00000 +296, 0.25000, 0.75000, 1.00000 +297, 0.37500, 0.75000, 1.00000 +298, 0.25000, 0.87500, 1.00000 +299, 0.12500, 0.87500, 1.00000 +300, 0.12500, 0.75000, 1.00000 +301, 0.87500, 0.62500, 1.00000 +302, 0.75000, 0.50000, 1.00000 +303, 0.87500, 0.50000, 1.00000 +304, 0.75000, 0.62500, 1.00000 +305, 0.62500, 0.62500, 1.00000 +306, 0.50000, 0.50000, 1.00000 +307, 0.62500, 0.50000, 1.00000 +308, 0.50000, 0.62500, 1.00000 +309, 0.37500, 0.62500, 1.00000 +310, 0.25000, 0.50000, 1.00000 +311, 0.37500, 0.50000, 1.00000 +312, 0.25000, 0.62500, 1.00000 +313, 0.12500, 0.62500, 1.00000 +314, 0.12500, 0.50000, 1.00000 +315, 0.87500, 0.37500, 1.00000 +316, 0.75000, 0.25000, 1.00000 +317, 0.87500, 0.25000, 1.00000 +318, 0.75000, 0.37500, 1.00000 +319, 0.62500, 0.37500, 1.00000 +320, 0.50000, 0.25000, 1.00000 +321, 0.62500, 0.25000, 1.00000 +322, 0.50000, 0.37500, 1.00000 +323, 0.37500, 0.37500, 1.00000 +324, 0.25000, 0.25000, 1.00000 +325, 0.37500, 0.25000, 1.00000 +326, 0.25000, 0.37500, 1.00000 +327, 0.12500, 0.37500, 1.00000 +328, 0.12500, 0.25000, 1.00000 +329, 0.87500, 0.12500, 1.00000 +330, 0.75000, 0.12500, 1.00000 +331, 0.62500, 0.12500, 1.00000 +332, 0.50000, 0.12500, 1.00000 +333, 0.37500, 0.12500, 1.00000 +334, 0.25000, 0.12500, 1.00000 +335, 0.12500, 0.12500, 1.00000 +336, 0.87500, 0.87500, 0.50000 +337, 0.75000, 0.75000, 0.50000 +338, 0.87500, 0.75000, 0.50000 +339, 0.75000, 0.87500, 0.50000 +340, 0.62500, 0.87500, 0.50000 +341, 0.50000, 0.75000, 0.50000 +342, 0.62500, 0.75000, 0.50000 +343, 0.50000, 0.87500, 0.50000 +344, 0.37500, 0.87500, 0.50000 +345, 0.25000, 0.75000, 0.50000 +346, 0.37500, 0.75000, 0.50000 +347, 0.25000, 0.87500, 0.50000 +348, 0.12500, 0.87500, 0.50000 +349, 0.12500, 0.75000, 0.50000 +350, 0.87500, 0.62500, 0.50000 +351, 0.75000, 0.50000, 0.50000 +352, 0.87500, 0.50000, 0.50000 +353, 0.75000, 0.62500, 0.50000 +354, 0.62500, 0.62500, 0.50000 +355, 0.50000, 0.50000, 0.50000 +356, 0.62500, 0.50000, 0.50000 +357, 0.50000, 0.62500, 0.50000 +358, 0.37500, 0.62500, 0.50000 +359, 0.25000, 0.50000, 0.50000 +360, 0.37500, 0.50000, 0.50000 +361, 0.25000, 0.62500, 0.50000 +362, 0.12500, 0.62500, 0.50000 +363, 0.12500, 0.50000, 0.50000 +364, 0.87500, 0.37500, 0.50000 +365, 0.75000, 0.25000, 0.50000 +366, 0.87500, 0.25000, 0.50000 +367, 0.75000, 0.37500, 0.50000 +368, 0.62500, 0.37500, 0.50000 +369, 0.50000, 0.25000, 0.50000 +370, 0.62500, 0.25000, 0.50000 +371, 0.50000, 0.37500, 0.50000 +372, 0.37500, 0.37500, 0.50000 +373, 0.25000, 0.25000, 0.50000 +374, 0.37500, 0.25000, 0.50000 +375, 0.25000, 0.37500, 0.50000 +376, 0.12500, 0.37500, 0.50000 +377, 0.12500, 0.25000, 0.50000 +378, 0.87500, 0.12500, 0.50000 +379, 0.75000, 0.12500, 0.50000 +380, 0.62500, 0.12500, 0.50000 +381, 0.50000, 0.12500, 0.50000 +382, 0.37500, 0.12500, 0.50000 +383, 0.25000, 0.12500, 0.50000 +384, 0.12500, 0.12500, 0.50000 +385, 0.83333, 0.33333, 0.12500 +386, 0.83333, 0.50000, 0.12500 +387, 0.83333, 0.83333, 0.25000 +388, 0.83333, 0.16667, 0.12500 +389, 0.66667, 0.16667, 0.12500 +390, 0.16667, 0.50000, 0.12500 +391, 0.16667, 0.66667, 0.12500 +392, 0.83333, 0.33333, 0.37500 +393, 0.83333, 0.50000, 0.37500 +394, 0.83333, 0.66667, 0.12500 +395, 0.33333, 0.16667, 0.12500 +396, 0.16667, 0.33333, 0.37500 +397, 0.16667, 0.50000, 0.37500 +398, 0.33333, 0.83333, 0.12500 +399, 0.50000, 0.83333, 0.12500 +400, 0.16667, 0.83333, 0.12500 +401, 0.83333, 0.16667, 0.37500 +402, 0.83333, 0.16667, 0.25000 +403, 0.50000, 0.83333, 0.37500 +404, 0.66667, 0.83333, 0.37500 +405, 0.16667, 0.33333, 0.12500 +406, 0.16667, 0.16667, 0.12500 +407, 0.16667, 0.16667, 0.25000 +408, 0.16667, 0.16667, 0.37500 +409, 0.83333, 0.66667, 0.37500 +410, 0.33333, 0.83333, 0.37500 +411, 0.50000, 0.50000, 0.25000 +412, 0.41667, 0.41667, 0.37500 +413, 0.25000, 0.41667, 0.25000 +414, 0.25000, 0.58333, 0.25000 +415, 0.33333, 0.16667, 0.37500 +416, 0.41667, 0.25000, 0.25000 +417, 0.50000, 0.16667, 0.37500 +418, 0.58333, 0.25000, 0.25000 +419, 0.66667, 0.16667, 0.37500 +420, 0.58333, 0.41667, 0.37500 +421, 0.41667, 0.58333, 0.37500 +422, 0.16667, 0.66667, 0.37500 +423, 0.16667, 0.83333, 0.25000 +424, 0.41667, 0.58333, 0.12500 +425, 0.41667, 0.75000, 0.25000 +426, 0.58333, 0.75000, 0.25000 +427, 0.58333, 0.58333, 0.37500 +428, 0.75000, 0.58333, 0.25000 +429, 0.75000, 0.41667, 0.25000 +430, 0.58333, 0.41667, 0.12500 +431, 0.50000, 0.16667, 0.12500 +432, 0.41667, 0.41667, 0.12500 +433, 0.58333, 0.58333, 0.12500 +434, 0.66667, 0.83333, 0.12500 +435, 0.16667, 0.83333, 0.37500 +436, 0.83333, 0.83333, 0.37500 +437, 0.83333, 0.83333, 0.12500 +438, 1.00000, 0.00000, 0.25000 +439, 1.00000, 0.00000, 0.50000 +440, 0.66667, 0.00000, 0.50000 +441, 1.00000, 0.00000, 0.37500 +442, 0.83333, 0.00000, 0.50000 +443, 0.83333, 0.00000, 0.37500 +444, 0.66667, 0.00000, 0.25000 +445, 0.66667, 0.00000, 0.37500 +446, 0.83333, 0.00000, 0.25000 +447, 1.00000, 0.00000, 0.00000 +448, 1.00000, 0.00000, 0.12500 +449, 0.83333, 0.00000, 0.12500 +450, 0.66667, 0.00000, 0.00000 +451, 0.66667, 0.00000, 0.12500 +452, 0.83333, 0.00000, 0.00000 +453, 0.33333, 0.00000, 0.50000 +454, 0.50000, 0.00000, 0.50000 +455, 0.50000, 0.00000, 0.37500 +456, 0.33333, 0.00000, 0.25000 +457, 0.33333, 0.00000, 0.37500 +458, 0.50000, 0.00000, 0.25000 +459, 0.50000, 0.00000, 0.12500 +460, 0.33333, 0.00000, 0.00000 +461, 0.33333, 0.00000, 0.12500 +462, 0.50000, 0.00000, 0.00000 +463, 0.00000, 0.00000, 0.50000 +464, 0.16667, 0.00000, 0.50000 +465, 0.16667, 0.00000, 0.37500 +466, 0.00000, 0.00000, 0.25000 +467, 0.00000, 0.00000, 0.37500 +468, 0.16667, 0.00000, 0.25000 +469, 0.16667, 0.00000, 0.12500 +470, 0.00000, 0.00000, 0.00000 +471, 0.00000, 0.00000, 0.12500 +472, 0.16667, 0.00000, 0.00000 +473, 0.66667, 1.00000, 0.50000 +474, 1.00000, 1.00000, 0.50000 +475, 1.00000, 1.00000, 0.25000 +476, 0.83333, 1.00000, 0.50000 +477, 1.00000, 1.00000, 0.37500 +478, 0.83333, 1.00000, 0.37500 +479, 0.66667, 1.00000, 0.25000 +480, 0.66667, 1.00000, 0.37500 +481, 0.83333, 1.00000, 0.25000 +482, 1.00000, 1.00000, 0.00000 +483, 1.00000, 1.00000, 0.12500 +484, 0.83333, 1.00000, 0.12500 +485, 0.66667, 1.00000, 0.00000 +486, 0.66667, 1.00000, 0.12500 +487, 0.83333, 1.00000, 0.00000 +488, 0.33333, 1.00000, 0.50000 +489, 0.50000, 1.00000, 0.50000 +490, 0.50000, 1.00000, 0.37500 +491, 0.33333, 1.00000, 0.25000 +492, 0.33333, 1.00000, 0.37500 +493, 0.50000, 1.00000, 0.25000 +494, 0.50000, 1.00000, 0.12500 +495, 0.33333, 1.00000, 0.00000 +496, 0.33333, 1.00000, 0.12500 +497, 0.50000, 1.00000, 0.00000 +498, 0.00000, 1.00000, 0.50000 +499, 0.16667, 1.00000, 0.50000 +500, 0.16667, 1.00000, 0.37500 +501, 0.00000, 1.00000, 0.25000 +502, 0.00000, 1.00000, 0.37500 +503, 0.16667, 1.00000, 0.25000 +504, 0.16667, 1.00000, 0.12500 +505, 0.00000, 1.00000, 0.00000 +506, 0.00000, 1.00000, 0.12500 +507, 0.16667, 1.00000, 0.00000 +508, 1.00000, 0.66667, 0.50000 +509, 1.00000, 0.83333, 0.50000 +510, 1.00000, 0.83333, 0.37500 +511, 1.00000, 0.66667, 0.25000 +512, 1.00000, 0.66667, 0.37500 +513, 1.00000, 0.83333, 0.25000 +514, 1.00000, 0.83333, 0.12500 +515, 1.00000, 0.66667, 0.00000 +516, 1.00000, 0.66667, 0.12500 +517, 1.00000, 0.83333, 0.00000 +518, 1.00000, 0.33333, 0.50000 +519, 1.00000, 0.50000, 0.50000 +520, 1.00000, 0.50000, 0.37500 +521, 1.00000, 0.33333, 0.25000 +522, 1.00000, 0.33333, 0.37500 +523, 1.00000, 0.50000, 0.25000 +524, 1.00000, 0.50000, 0.12500 +525, 1.00000, 0.33333, 0.00000 +526, 1.00000, 0.33333, 0.12500 +527, 1.00000, 0.50000, 0.00000 +528, 1.00000, 0.16667, 0.50000 +529, 1.00000, 0.16667, 0.37500 +530, 1.00000, 0.16667, 0.25000 +531, 1.00000, 0.16667, 0.12500 +532, 1.00000, 0.16667, 0.00000 +533, 0.00000, 0.66667, 0.50000 +534, 0.00000, 0.83333, 0.50000 +535, 0.00000, 0.83333, 0.37500 +536, 0.00000, 0.66667, 0.25000 +537, 0.00000, 0.66667, 0.37500 +538, 0.00000, 0.83333, 0.25000 +539, 0.00000, 0.83333, 0.12500 +540, 0.00000, 0.66667, 0.00000 +541, 0.00000, 0.66667, 0.12500 +542, 0.00000, 0.83333, 0.00000 +543, 0.00000, 0.33333, 0.50000 +544, 0.00000, 0.50000, 0.50000 +545, 0.00000, 0.50000, 0.37500 +546, 0.00000, 0.33333, 0.25000 +547, 0.00000, 0.33333, 0.37500 +548, 0.00000, 0.50000, 0.25000 +549, 0.00000, 0.50000, 0.12500 +550, 0.00000, 0.33333, 0.00000 +551, 0.00000, 0.33333, 0.12500 +552, 0.00000, 0.50000, 0.00000 +553, 0.00000, 0.16667, 0.50000 +554, 0.00000, 0.16667, 0.37500 +555, 0.00000, 0.16667, 0.25000 +556, 0.00000, 0.16667, 0.12500 +557, 0.00000, 0.16667, 0.00000 +558, 0.83333, 0.83333, 0.50000 +559, 0.66667, 0.66667, 0.50000 +560, 0.83333, 0.66667, 0.50000 +561, 0.66667, 0.83333, 0.50000 +562, 0.50000, 0.83333, 0.50000 +563, 0.33333, 0.66667, 0.50000 +564, 0.50000, 0.66667, 0.50000 +565, 0.33333, 0.83333, 0.50000 +566, 0.16667, 0.83333, 0.50000 +567, 0.16667, 0.66667, 0.50000 +568, 0.83333, 0.50000, 0.50000 +569, 0.66667, 0.33333, 0.50000 +570, 0.83333, 0.33333, 0.50000 +571, 0.66667, 0.50000, 0.50000 +572, 0.50000, 0.50000, 0.50000 +573, 0.33333, 0.33333, 0.50000 +574, 0.50000, 0.33333, 0.50000 +575, 0.33333, 0.50000, 0.50000 +576, 0.16667, 0.50000, 0.50000 +577, 0.16667, 0.33333, 0.50000 +578, 0.83333, 0.16667, 0.50000 +579, 0.66667, 0.16667, 0.50000 +580, 0.50000, 0.16667, 0.50000 +581, 0.33333, 0.16667, 0.50000 +582, 0.16667, 0.16667, 0.50000 +583, 0.83333, 0.83333, 0.00000 +584, 0.66667, 0.66667, 0.00000 +585, 0.83333, 0.66667, 0.00000 +586, 0.66667, 0.83333, 0.00000 +587, 0.50000, 0.83333, 0.00000 +588, 0.33333, 0.66667, 0.00000 +589, 0.50000, 0.66667, 0.00000 +590, 0.33333, 0.83333, 0.00000 +591, 0.16667, 0.83333, 0.00000 +592, 0.16667, 0.66667, 0.00000 +593, 0.83333, 0.50000, 0.00000 +594, 0.66667, 0.33333, 0.00000 +595, 0.83333, 0.33333, 0.00000 +596, 0.66667, 0.50000, 0.00000 +597, 0.50000, 0.50000, 0.00000 +598, 0.33333, 0.33333, 0.00000 +599, 0.50000, 0.33333, 0.00000 +600, 0.33333, 0.50000, 0.00000 +601, 0.16667, 0.50000, 0.00000 +602, 0.16667, 0.33333, 0.00000 +603, 0.83333, 0.16667, 0.00000 +604, 0.66667, 0.16667, 0.00000 +605, 0.50000, 0.16667, 0.00000 +606, 0.33333, 0.16667, 0.00000 +607, 0.16667, 0.16667, 0.00000 +** +**ELSET COUNT = 177 +**HWCOLOR COMP 54 0 +*ELEMENT, TYPE=C3D10, ELSET=UPPER + 1, 292, 187, 288, 1, 294, 291, 293, + 2, 3, 4 + 2, 302, 227, 230, 1, 303, 231, 5, + 6, 7, 8 + 3, 302, 288, 227, 1, 304, 301, 303, + 6, 4, 7 + 4, 227, 288, 217, 1, 301, 289, 228, + 7, 4, 9 + 5, 227, 217, 220, 1, 228, 221, 229, + 7, 9, 10 + 6, 227, 220, 230, 1, 229, 232, 231, + 7, 10, 8 + 7, 337, 230, 220, 1, 11, 232, 12, + 13, 8, 10 + 8, 351, 230, 337, 14, 15, 11, 353, + 16, 17, 18 + 9, 341, 337, 190, 1, 342, 19, 20, + 21, 13, 22 + 10, 337, 178, 190, 1, 23, 192, 19, + 13, 24, 22 + 11, 337, 220, 178, 1, 12, 25, 23, + 13, 10, 24 + 12, 217, 178, 220, 1, 26, 25, 221, + 9, 24, 10 + 13, 217, 172, 178, 1, 287, 179, 26, + 9, 27, 24 + 14, 288, 172, 217, 1, 290, 287, 289, + 4, 27, 9 + 15, 288, 187, 172, 1, 291, 188, 290, + 4, 3, 27 + 16, 187, 178, 172, 1, 189, 179, 188, + 3, 24, 27 + 17, 306, 292, 302, 14, 308, 305, 307, + 28, 29, 30 + 18, 320, 306, 316, 14, 322, 319, 321, + 31, 28, 32 + 19, 316, 306, 302, 14, 319, 307, 318, + 32, 28, 30 + 20, 316, 302, 240, 14, 318, 33, 34, + 32, 30, 35 + 21, 369, 365, 355, 14, 370, 368, 371, + 36, 37, 38 + 22, 365, 351, 355, 14, 367, 356, 368, + 37, 16, 38 + 23, 365, 240, 351, 14, 39, 40, 367, + 37, 35, 16 + 24, 351, 240, 230, 14, 40, 242, 15, + 16, 35, 17 + 25, 302, 230, 240, 14, 5, 242, 33, + 30, 17, 35 + 26, 302, 1, 230, 14, 6, 8, 5, + 30, 41, 17 + 27, 302, 292, 288, 1, 305, 293, 304, + 6, 2, 4 + 28, 302, 292, 1, 14, 305, 2, 6, + 30, 29, 41 + 29, 337, 224, 184, 178, 338, 336, 339, + 23, 42, 185 + 30, 355, 351, 341, 14, 356, 354, 357, + 38, 16, 43 + 31, 351, 337, 341, 14, 353, 342, 354, + 16, 18, 43 + 32, 337, 230, 1, 14, 11, 8, 13, + 18, 17, 41 + 33, 365, 136, 244, 133, 378, 251, 366, + 44, 138, 45 + 34, 341, 337, 1, 14, 342, 13, 21, + 43, 18, 41 + 35, 129, 240, 127, 133, 46, 249, 132, + 134, 47, 135 + 36, 292, 1, 14, 48, 2, 41, 29, + 49, 50, 51 + 37, 306, 292, 14, 48, 308, 29, 28, + 52, 49, 51 + 38, 306, 296, 292, 48, 309, 297, 308, + 52, 53, 49 + 39, 310, 296, 306, 48, 312, 309, 311, + 54, 53, 52 + 40, 310, 255, 296, 48, 55, 56, 312, + 54, 57, 53 + 41, 310, 265, 255, 48, 58, 267, 55, + 54, 59, 57 + 42, 355, 14, 341, 48, 38, 43, 357, + 60, 51, 61 + 43, 341, 14, 1, 48, 43, 41, 21, + 61, 51, 50 + 44, 341, 1, 190, 48, 21, 22, 20, + 61, 50, 62 + 45, 341, 190, 194, 48, 20, 195, 343, + 61, 62, 63 + 46, 341, 194, 204, 48, 343, 206, 344, + 61, 63, 64 + 47, 345, 341, 204, 48, 346, 344, 347, + 65, 61, 64 + 48, 355, 341, 345, 48, 357, 346, 358, + 60, 61, 65 + 49, 359, 355, 345, 48, 360, 358, 361, + 66, 60, 65 + 50, 359, 345, 265, 48, 361, 67, 68, + 66, 65, 59 + 51, 345, 255, 265, 48, 69, 267, 67, + 65, 57, 59 + 52, 259, 345, 214, 255, 349, 348, 261, + 260, 69, 258 + 53, 345, 204, 255, 48, 347, 70, 69, + 65, 64, 57 + 54, 204, 194, 200, 48, 206, 203, 205, + 64, 63, 71 + 55, 194, 190, 200, 48, 195, 202, 203, + 63, 62, 71 + 56, 197, 200, 190, 48, 201, 202, 199, + 72, 71, 62 + 57, 296, 200, 197, 48, 73, 201, 298, + 53, 71, 72 + 58, 296, 197, 292, 48, 298, 295, 297, + 53, 72, 49 + 59, 292, 197, 187, 48, 295, 198, 294, + 49, 72, 74 + 60, 292, 187, 1, 48, 294, 3, 2, + 49, 74, 50 + 61, 197, 190, 187, 48, 199, 191, 198, + 72, 62, 74 + 62, 296, 255, 200, 48, 56, 75, 73, + 53, 57, 71 + 63, 142, 145, 320, 76, 146, 77, 332, + 78, 79, 80 + 64, 142, 133, 145, 76, 144, 147, 146, + 78, 81, 79 + 65, 142, 129, 133, 76, 143, 134, 144, + 78, 82, 81 + 66, 142, 320, 129, 76, 332, 331, 143, + 78, 80, 82 + 67, 320, 316, 129, 76, 321, 330, 331, + 80, 83, 82 + 68, 320, 14, 316, 76, 31, 32, 321, + 80, 84, 83 + 69, 316, 14, 240, 76, 32, 35, 34, + 83, 84, 85 + 70, 129, 240, 133, 76, 46, 47, 134, + 82, 85, 81 + 71, 129, 316, 240, 76, 330, 34, 46, + 82, 83, 85 + 72, 365, 133, 240, 76, 44, 47, 39, + 86, 81, 85 + 73, 365, 240, 14, 76, 39, 35, 37, + 86, 85, 84 + 74, 369, 365, 14, 76, 370, 37, 36, + 87, 86, 84 + 75, 369, 139, 365, 76, 380, 379, 370, + 87, 88, 86 + 76, 149, 139, 369, 76, 151, 380, 381, + 89, 88, 87 + 77, 149, 145, 139, 76, 150, 148, 151, + 89, 79, 88 + 78, 373, 149, 369, 90, 382, 381, 374, + 91, 92, 93 + 79, 139, 145, 133, 76, 148, 147, 140, + 88, 79, 81 + 80, 139, 133, 365, 76, 140, 44, 379, + 88, 81, 86 + 81, 310, 48, 306, 94, 54, 52, 311, + 95, 96, 97 + 82, 310, 265, 48, 94, 58, 59, 54, + 95, 98, 96 + 83, 324, 310, 320, 94, 326, 323, 325, + 99, 95, 100 + 84, 373, 369, 359, 94, 374, 372, 375, + 101, 102, 103 + 85, 369, 76, 14, 94, 87, 84, 36, + 102, 104, 105 + 86, 320, 14, 76, 94, 31, 84, 80, + 100, 105, 104 + 87, 320, 306, 14, 94, 322, 28, 31, + 100, 97, 105 + 88, 320, 310, 306, 94, 323, 311, 322, + 100, 95, 97 + 89, 306, 48, 14, 94, 52, 51, 28, + 97, 96, 105 + 90, 355, 14, 48, 94, 38, 51, 60, + 106, 105, 96 + 91, 369, 14, 355, 94, 36, 38, 371, + 102, 105, 106 + 92, 369, 355, 359, 94, 371, 360, 372, + 102, 106, 103 + 93, 359, 355, 48, 94, 360, 60, 66, + 103, 106, 96 + 94, 359, 48, 265, 94, 66, 59, 68, + 103, 96, 98 + 95, 324, 145, 155, 90, 107, 157, 108, + 109, 110, 111 + 96, 324, 155, 275, 90, 108, 112, 113, + 109, 111, 114 + 97, 324, 275, 265, 90, 113, 277, 115, + 109, 114, 116 + 98, 324, 265, 310, 94, 115, 58, 326, + 99, 98, 95 + 99, 224, 181, 184, 178, 226, 186, 336, + 42, 183, 185 + 100, 324, 94, 320, 90, 99, 100, 325, + 109, 117, 118 + 101, 324, 320, 145, 90, 325, 77, 107, + 109, 118, 110 + 102, 320, 76, 145, 90, 80, 79, 77, + 118, 119, 110 + 103, 320, 94, 76, 90, 100, 104, 80, + 118, 117, 119 + 104, 373, 369, 94, 90, 374, 102, 101, + 91, 93, 117 + 105, 149, 145, 76, 90, 150, 79, 89, + 92, 110, 119 + 106, 162, 275, 152, 155, 283, 120, 163, + 164, 112, 156 + 107, 149, 76, 369, 90, 89, 87, 381, + 92, 119, 93 + 108, 159, 149, 373, 90, 161, 382, 383, + 121, 92, 91 + 109, 159, 155, 149, 90, 160, 158, 161, + 121, 111, 92 + 110, 159, 275, 155, 90, 122, 112, 160, + 121, 114, 111 + 111, 279, 275, 159, 90, 280, 122, 384, + 123, 114, 121 + 112, 279, 159, 373, 90, 384, 383, 377, + 123, 121, 91 + 113, 279, 373, 269, 90, 377, 376, 281, + 123, 91, 124 + 114, 279, 269, 275, 90, 281, 278, 280, + 123, 124, 114 + 115, 269, 265, 275, 90, 270, 277, 278, + 124, 116, 114 + 116, 269, 359, 265, 90, 363, 68, 270, + 124, 125, 116 + 117, 373, 359, 269, 90, 375, 363, 376, + 91, 125, 124 + 118, 373, 94, 359, 90, 101, 103, 375, + 91, 117, 125 + 119, 359, 94, 265, 90, 103, 98, 68, + 125, 117, 116 + 120, 324, 265, 94, 90, 115, 98, 99, + 109, 116, 117 + 121, 369, 76, 94, 90, 87, 104, 102, + 93, 119, 117 + 122, 149, 155, 145, 90, 158, 157, 150, + 92, 111, 110 + 123, 217, 174, 220, 178, 219, 222, 221, + 26, 180, 25 + 124, 345, 204, 214, 255, 347, 216, 348, + 69, 70, 258 + 125, 337, 184, 190, 178, 339, 193, 19, + 23, 185, 192 + 126, 337, 220, 224, 178, 12, 225, 338, + 23, 25, 42 + 127, 217, 172, 174, 178, 287, 177, 219, + 26, 179, 180 + 128, 204, 210, 255, 200, 213, 257, 70, + 205, 212, 75 + 129, 234, 224, 337, 230, 236, 338, 350, + 235, 233, 11 + 130, 181, 220, 174, 178, 223, 222, 182, + 183, 25, 180 + 131, 337, 230, 224, 220, 11, 233, 338, + 12, 232, 225 + 132, 341, 337, 194, 190, 342, 340, 343, + 20, 19, 195 + 133, 237, 240, 302, 230, 241, 33, 315, + 239, 242, 5 + 134, 359, 345, 259, 265, 361, 349, 362, + 68, 67, 268 + 135, 237, 302, 227, 230, 315, 303, 238, + 239, 5, 231 + 136, 351, 240, 234, 230, 40, 243, 352, + 15, 242, 235 + 137, 351, 234, 337, 230, 352, 350, 353, + 15, 235, 11 + 138, 244, 127, 240, 133, 250, 249, 245, + 45, 135, 47 + 139, 136, 127, 244, 133, 137, 250, 251, + 138, 135, 45 + 140, 337, 184, 194, 190, 339, 196, 340, + 19, 193, 195 + 141, 244, 234, 351, 240, 246, 352, 364, + 245, 243, 40 + 142, 365, 244, 351, 240, 366, 364, 367, + 39, 245, 40 + 143, 152, 155, 324, 145, 156, 108, 334, + 154, 157, 107 + 144, 316, 237, 128, 240, 317, 247, 329, + 34, 241, 248 + 145, 129, 316, 128, 240, 330, 329, 131, + 46, 34, 248 + 146, 252, 210, 207, 200, 254, 211, 253, + 126, 212, 209 + 147, 269, 359, 259, 265, 363, 362, 271, + 270, 68, 268 + 148, 324, 262, 310, 265, 327, 314, 326, + 115, 266, 58 + 149, 310, 252, 296, 255, 313, 300, 312, + 55, 256, 56 + 150, 324, 320, 142, 145, 325, 332, 333, + 107, 77, 146 + 151, 162, 165, 275, 155, 166, 284, 283, + 164, 167, 112 + 152, 279, 169, 159, 165, 286, 171, 384, + 285, 170, 168 + 153, 214, 255, 204, 210, 258, 70, 216, + 215, 257, 213 + 154, 217, 172, 173, 174, 287, 175, 218, + 219, 177, 176 + 155, 272, 262, 324, 265, 273, 327, 328, + 274, 266, 115 + 156, 316, 302, 237, 240, 318, 315, 317, + 34, 33, 241 + 157, 129, 128, 127, 240, 131, 130, 132, + 46, 248, 249 + 158, 139, 136, 365, 133, 141, 378, 379, + 140, 138, 44 + 159, 159, 275, 165, 155, 122, 284, 168, + 160, 112, 167 + 160, 152, 272, 324, 275, 335, 328, 334, + 120, 276, 113 + 161, 259, 265, 345, 255, 268, 67, 349, + 260, 267, 69 + 162, 262, 252, 310, 255, 263, 313, 314, + 264, 256, 55 + 163, 252, 296, 255, 200, 300, 56, 256, + 126, 73, 75 + 164, 204, 200, 255, 48, 205, 75, 70, + 64, 71, 57 + 165, 365, 244, 240, 133, 366, 245, 39, + 44, 45, 47 + 166, 279, 165, 159, 275, 285, 168, 384, + 280, 284, 122 + 167, 224, 220, 181, 178, 225, 223, 226, + 42, 25, 183 + 168, 296, 207, 197, 200, 299, 208, 298, + 73, 209, 201 + 169, 272, 324, 275, 265, 328, 113, 276, + 274, 115, 277 + 170, 252, 207, 296, 200, 253, 299, 300, + 126, 209, 73 + 171, 252, 255, 210, 200, 256, 257, 254, + 126, 75, 212 + 172, 162, 272, 152, 275, 282, 335, 163, + 283, 276, 120 + 173, 262, 310, 265, 255, 314, 58, 266, + 264, 55, 267 + 174, 152, 275, 324, 155, 120, 113, 334, + 156, 112, 108 + 175, 152, 324, 142, 145, 334, 333, 153, + 154, 107, 146 + 176, 187, 190, 178, 1, 191, 192, 189, + 3, 22, 24 + 177, 187, 190, 1, 48, 191, 22, 3, + 74, 62, 50 +** +**ELSET COUNT = 92 +**HWCOLOR COMP 1 0 +*ELEMENT, TYPE=C3D10, ELSET=LOWER + 354, 594, 525, 584, 521, 595, 593, 596, + 385, 526, 386 + 355, 482, 511, 475, 479, 514, 513, 483, + 484, 387, 481 + 356, 594, 447, 438, 444, 603, 448, 388, + 389, 449, 446 + 357, 540, 546, 588, 536, 549, 390, 592, + 541, 548, 391 + 358, 569, 518, 521, 511, 570, 522, 392, + 393, 520, 523 + 359, 584, 521, 515, 511, 386, 524, 585, + 394, 523, 516 + 360, 460, 450, 598, 456, 462, 605, 606, + 461, 459, 395 + 361, 543, 573, 546, 536, 577, 396, 547, + 545, 397, 548 + 362, 588, 584, 495, 491, 589, 587, 590, + 398, 399, 496 + 363, 588, 495, 501, 491, 590, 504, 400, + 398, 496, 503 + 364, 588, 495, 505, 501, 590, 507, 591, + 400, 504, 506 + 365, 588, 505, 536, 501, 591, 539, 391, + 400, 506, 538 + 366, 450, 447, 594, 444, 452, 603, 604, + 451, 449, 389 + 367, 440, 521, 438, 444, 401, 530, 443, + 445, 402, 446 + 368, 563, 488, 559, 479, 565, 562, 564, + 403, 490, 404 + 369, 550, 598, 546, 456, 602, 405, 551, + 406, 395, 407 + 370, 463, 546, 543, 456, 554, 547, 553, + 465, 407, 408 + 371, 440, 569, 439, 521, 579, 578, 442, + 401, 392, 529 + 372, 569, 559, 518, 511, 571, 568, 570, + 393, 409, 520 + 373, 463, 543, 453, 456, 553, 582, 464, + 465, 408, 457 + 374, 550, 466, 460, 456, 556, 469, 607, + 406, 468, 461 + 375, 563, 491, 488, 479, 410, 492, 565, + 403, 493, 490 + 376, 573, 546, 536, 411, 396, 548, 397, + 412, 413, 414 + 377, 573, 456, 546, 411, 415, 407, 396, + 412, 416, 413 + 378, 573, 444, 456, 411, 417, 458, 415, + 412, 418, 416 + 379, 573, 569, 444, 411, 574, 419, 417, + 412, 420, 418 + 380, 573, 563, 569, 411, 575, 572, 574, + 412, 421, 420 + 381, 573, 536, 563, 411, 397, 422, 575, + 412, 414, 421 + 382, 588, 491, 536, 411, 398, 423, 391, + 424, 425, 414 + 383, 563, 536, 491, 411, 422, 423, 410, + 421, 414, 425 + 384, 563, 491, 479, 411, 410, 493, 403, + 421, 425, 426 + 385, 563, 479, 559, 411, 403, 404, 564, + 421, 426, 427 + 386, 569, 563, 559, 411, 572, 564, 571, + 420, 421, 427 + 387, 569, 559, 511, 411, 571, 409, 393, + 420, 427, 428 + 388, 569, 511, 521, 411, 393, 523, 392, + 420, 428, 429 + 389, 594, 444, 521, 411, 389, 402, 385, + 430, 418, 429 + 390, 569, 521, 444, 411, 392, 402, 419, + 420, 429, 418 + 391, 594, 456, 444, 411, 431, 458, 389, + 430, 416, 418 + 392, 598, 456, 594, 411, 395, 431, 599, + 432, 416, 430 + 393, 598, 594, 588, 411, 599, 597, 600, + 432, 430, 424 + 394, 598, 588, 546, 411, 600, 390, 405, + 432, 424, 413 + 395, 598, 546, 456, 411, 405, 407, 395, + 432, 413, 416 + 396, 588, 536, 546, 411, 391, 548, 390, + 424, 414, 413 + 397, 594, 584, 588, 411, 596, 589, 597, + 430, 433, 424 + 398, 594, 521, 584, 411, 385, 386, 596, + 430, 429, 433 + 399, 584, 521, 511, 411, 386, 523, 394, + 433, 429, 428 + 400, 559, 479, 511, 411, 404, 387, 409, + 427, 426, 428 + 401, 584, 511, 479, 411, 394, 387, 434, + 433, 428, 426 + 402, 584, 479, 491, 411, 434, 493, 399, + 433, 426, 425 + 403, 588, 584, 491, 411, 589, 399, 398, + 424, 433, 425 + 404, 533, 498, 563, 491, 534, 566, 567, + 435, 500, 410 + 405, 559, 508, 511, 479, 560, 512, 409, + 404, 436, 387 + 406, 515, 511, 482, 479, 516, 514, 517, + 437, 387, 484 + 407, 594, 438, 525, 521, 388, 531, 595, + 385, 530, 526 + 408, 594, 447, 525, 438, 603, 532, 595, + 388, 448, 531 + 409, 598, 450, 594, 456, 605, 604, 599, + 395, 459, 431 + 410, 508, 473, 475, 479, 558, 478, 510, + 436, 480, 481 + 411, 515, 482, 485, 479, 517, 487, 583, + 437, 484, 486 + 412, 508, 473, 474, 475, 558, 476, 509, + 510, 478, 477 + 413, 559, 473, 508, 479, 561, 558, 560, + 404, 480, 436 + 414, 518, 559, 508, 511, 568, 560, 519, + 520, 409, 512 + 415, 563, 498, 488, 491, 566, 499, 565, + 410, 500, 492 + 416, 559, 488, 473, 479, 562, 489, 561, + 404, 490, 480 + 417, 584, 515, 485, 479, 585, 583, 586, + 434, 437, 486 + 418, 588, 501, 536, 491, 400, 538, 391, + 398, 503, 423 + 419, 594, 438, 521, 444, 388, 530, 385, + 389, 446, 402 + 420, 533, 501, 498, 491, 535, 502, 534, + 435, 503, 500 + 421, 543, 533, 573, 536, 544, 576, 577, + 545, 537, 397 + 422, 569, 518, 439, 521, 570, 528, 578, + 392, 522, 529 + 423, 533, 563, 536, 491, 567, 422, 537, + 435, 410, 423 + 424, 440, 569, 521, 444, 579, 392, 401, + 445, 419, 402 + 425, 450, 594, 456, 444, 604, 431, 459, + 451, 389, 458 + 426, 573, 569, 440, 444, 574, 579, 580, + 417, 419, 445 + 427, 550, 546, 466, 456, 551, 555, 556, + 406, 407, 468 + 428, 440, 439, 438, 521, 442, 441, 443, + 401, 529, 530 + 429, 508, 475, 511, 479, 510, 513, 512, + 436, 481, 387 + 430, 453, 573, 440, 444, 581, 580, 454, + 455, 417, 445 + 431, 453, 456, 573, 444, 457, 415, 581, + 455, 458, 417 + 432, 550, 460, 598, 456, 607, 606, 602, + 406, 461, 395 + 433, 453, 543, 573, 456, 582, 577, 581, + 457, 408, 415 + 434, 543, 546, 573, 456, 547, 396, 577, + 408, 407, 415 + 435, 540, 588, 505, 536, 592, 591, 542, + 541, 391, 539 + 436, 573, 533, 563, 536, 576, 567, 575, + 397, 537, 422 + 437, 550, 598, 540, 546, 602, 601, 552, + 551, 405, 549 + 438, 533, 536, 501, 491, 537, 538, 535, + 435, 423, 503 + 439, 550, 470, 460, 466, 557, 472, 607, + 556, 471, 469 + 440, 598, 588, 540, 546, 600, 592, 601, + 405, 390, 549 + 441, 584, 485, 491, 479, 586, 494, 399, + 434, 486, 493 + 442, 584, 485, 495, 491, 586, 497, 587, + 399, 494, 496 + 443, 584, 511, 515, 479, 394, 516, 585, + 434, 387, 437 + 444, 525, 515, 584, 521, 527, 585, 593, + 526, 524, 386 + 445, 463, 466, 546, 456, 467, 555, 554, + 465, 468, 407 +** +**ELSET COUNT = 18 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_TO_UPPER +373,S1 +433,S1 +430,S1 +426,S1 +371,S1 +422,S1 +421,S1 +436,S1 +380,S1 +386,S1 +372,S1 +414,S1 +404,S1 +415,S1 +368,S1 +416,S1 +413,S1 +412,S1 +** +**ELSET COUNT = 18 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_BOTTOM +439,S1 +432,S1 +360,S1 +409,S1 +366,S1 +408,S1 +437,S1 +440,S1 +393,S1 +397,S1 +354,S1 +444,S1 +435,S1 +364,S1 +362,S1 +442,S1 +417,S1 +411,S1 +** +**ELSET COUNT = 12 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM13 +428,S1 +445,S2 +360,S2 +431,S2 +366,S2 +439,S3 +374,S3 +356,S3 +373,S4 +425,S4 +430,S4 +367,S4 +** +**ELSET COUNT = 12 +*SURFACE,TYPE=ELEMENT,NAME=LOWER_SYM23 +427,S1 +445,S1 +370,S1 +438,S1 +420,S1 +439,S2 +357,S2 +421,S2 +365,S3 +437,S4 +361,S4 +435,S4 +** +**ELSET COUNT = 32 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_TOP +172,S1 +160,S1 +175,S1 +150,S1 +66,S1 +67,S1 +145,S1 +144,S1 +155,S1 +148,S1 +83,S1 +88,S1 +18,S1 +19,S1 +156,S1 +135,S1 +162,S1 +149,S1 +39,S1 +38,S1 +17,S1 +27,S1 +3,S1 +4,S1 +170,S1 +168,S1 +58,S1 +59,S1 +1,S1 +15,S1 +14,S1 +154,S1 +** +**ELSET COUNT = 32 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_TO_LOWER +152,S1 +112,S1 +108,S1 +78,S1 +76,S1 +75,S1 +158,S1 +33,S1 +113,S1 +117,S1 +84,S1 +92,S1 +21,S1 +22,S1 +142,S1 +141,S1 +147,S1 +134,S1 +49,S1 +48,S1 +30,S1 +31,S1 +137,S1 +129,S1 +52,S1 +124,S1 +47,S1 +46,S1 +132,S1 +140,S1 +29,S1 +99,S1 +** +**ELSET COUNT = 16 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM13 +109,S1 +122,S1 +77,S1 +79,S1 +64,S1 +65,S1 +157,S1 +151,S2 +143,S2 +158,S2 +139,S2 +152,S3 +159,S4 +106,S4 +175,S4 +35,S4 +** +**ELSET COUNT = 16 +*SURFACE,TYPE=ELEMENT,NAME=UPPER_SYM23 +151,S1 +114,S1 +115,S1 +171,S1 +146,S1 +152,S2 +166,S2 +172,S2 +155,S2 +161,S2 +162,S2 +153,S2 +169,S4 +147,S4 +173,S4 +52,S4 +** +**Property Definitions +** +*SOLID SECTION, ELSET=UPPER, MATERIAL=Def_Material +*SOLID SECTION, ELSET=LOWER, MATERIAL=Def_Material +** +**Material Definitions +** +**Material:Def_Material +*MATERIAL,NAME=Def_Material +*ELASTIC,TYPE=ISO +2.08000e+005,3.00000e-001 +*DENSITY +7.80000e-009, +*SPECIFIC HEAT +5.00000e-001 +*CONDUCTIVITY +4.98100e-002 +** \ No newline at end of file