Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-28 04:35:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

fixed deprecation warnings for 0.5

Browse Source
This commit is contained in:
Jukka Aho
2016-11-13 13:24:08 +02:00
parent 57a900e3d5
commit 50187466f4
18 changed files with 196 additions and 123 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/abaqus.jl
+2 -2
View File
@@ -625,7 +625,7 @@ function process_output_request(model::Model, solver::Solver, output_request::Ab
info("SECTION PRINT output request, with data $data and options $options")
end
function call(model::Model)
function (model::Model)()
info("Starting JuliaFEM-ABAQUS solver.")
# 1. create field problems and add elements
@@ -750,7 +750,7 @@ function create_surface_elements(mesh::Mesh, surface_name::Symbol)
get_child_element(parent_element_type, parent_element_side,
parent_element_connectivity)
child_element = Element(JuliaFEM.(child_element_type), child_element_connectivity)
child_element = Element(getfield(JuliaFEM, child_element_type), child_element_connectivity)
push!(elements, child_element)
end
update!(elements, "geometry", mesh.nodes)
src/elements.jl
+37 -21
View File
@@ -44,19 +44,40 @@ function setindex!(element::Element, data, field_name)
element.fields[field_name] = Field(data)
end
function call(element::Element, field_name::AbstractString)
""" Return a Field object from element.
Examples
--------
>>> element = Element(Seg2, [1, 2])
>>> data = Dict(1 => 1.0, 2 => 2.0)
>>> update!(element, "my field", data)
>>> element("my field")
"""
function (element::Element)(field_name::String)
return element[field_name]
end
function call(element::Element, field_name::AbstractString, time::Float64)
""" Return a Field object from element and interpolate in time direction.
Examples
--------
>>> element = Element(Seg2, [1, 2])
>>> data1 = Dict(1 => 1.0, 2 => 2.0)
>>> data2 = Dict(1 => 2.0, 2 => 3.0)
>>> update!(element, "my field", 0.0 => data1, 1.0 => data2)
>>> element("my field", 0.5)
"""
function (element::Element)(field_name::String, time::Float64)
return element[field_name](time)
end
function last(element::Element, field_name::AbstractString)
function last(element::Element, field_name::String)
return last(element[field_name])
end
function call(element::Element, ip, time::Float64=0.0)
function (element::Element)(ip, time::Float64=0.0)
return get_basis(element, ip, time)
end
@@ -75,7 +96,7 @@ julia> el([0.0, 0.0], 0.0, 2)
0.0 0.25 0.0 0.25 0.0 0.25 0.0 0.25
"""
function call(element::Element, ip, time::Float64, dim::Int)
function (element::Element)(ip, time::Float64, dim::Int)
dim == 1 && return get_basis(element, ip, time)
Ni = get_basis(element, ip, time)
N = zeros(dim, length(element)*dim)
@@ -85,7 +106,7 @@ function call(element::Element, ip, time::Float64, dim::Int)
return N
end
function call(element::Element, ip, time::Float64, ::Type{Val{:Jacobian}})
function (element::Element)(ip, time::Float64, ::Type{Val{:Jacobian}})
X = element("geometry", time)
dN = get_dbasis(element, ip, time)
nbasis = length(element)
@@ -98,7 +119,7 @@ function call(element::Element, ip, time::Float64, ::Type{Val{:Jacobian}})
return J
end
function call(element::Element, ip, time::Float64, ::Type{Val{:detJ}})
function (element::Element)(ip, time::Float64, ::Type{Val{:detJ}})
J = element(ip, time, Val{:Jacobian})
n, m = size(J)
if n == m # volume element
@@ -112,46 +133,41 @@ function call(element::Element, ip, time::Float64, ::Type{Val{:detJ}})
end
end
function call(element::Element, ip, time::Float64, ::Type{Val{:Grad}})
function (element::Element)(ip, time::Float64, ::Type{Val{:Grad}})
J = element(ip, time, Val{:Jacobian})
return inv(J)*get_dbasis(element, ip, time)
end
function call(element::Element, field_name::AbstractString, ip, time::Float64, ::Type{Val{:Grad}})
function (element::Element)(field_name::String, ip, time::Float64, ::Type{Val{:Grad}})
return element(ip, time, Val{:Grad})*element[field_name](time)
end
function call(element::Element, field::Field, time::Float64)
function (element::Element)(field::Field, time::Float64)
return field(time)
end
function call(element::Element, field::DCTI, time::Float64)
function (element::Element)(field::DCTI, time::Float64)
return field.data
end
function call(element::Element, field_name::AbstractString, time::Float64)
field = element[field_name]
return element(field, time)
end
function call(element::Element, field_name::AbstractString, ip, time::Float64)
function (element::Element)(field_name::String, ip, time::Float64)
field = element[field_name]
return element(field, ip, time)
end
function call(element::Element, field::DCTI, ip, time::Float64)
function (element::Element)(field::DCTI, ip, time::Float64)
return field.data
end
function call(element::Element, field::DCTV, ip, time::Float64)
function (element::Element)(field::DCTV, ip, time::Float64)
return field(time).data
end
function call(element::Element, field::CVTV, ip, time::Float64)
function (element::Element)(field::CVTV, ip, time::Float64)
return field(ip, time)
end
function call(element::Element, field::Field, ip, time::Float64)
function (element::Element)(field::Field, ip, time::Float64)
field_ = field(time)
basis = element(ip, time)
n = length(element)
src/elements_lagrange.jl
+1 -1
View File
@@ -26,7 +26,7 @@ function get_dbasis(element::Element{Poi1}, ip, time)
return [0]
end
function call(element::Element{Poi1}, ip, time::Float64, ::Type{Val{:detJ}})
function (element::Element{Poi1})(ip, time::Float64, ::Type{Val{:detJ}})
return 1.0
end
src/fields.jl
+24 -28
View File
@@ -38,7 +38,7 @@ function Base.getindex{T}(increment::Increment{Vector{T}}, i::Int64)
return increment.data[i]
end
function Base.(:*)(d, increment::Increment)
function Base.:*(d, increment::Increment)
return d*increment.data
end
@@ -49,11 +49,11 @@ type Basis
dbasis :: Function
end
function call(basis::Basis, xi::Vector)
function (basis::Basis)(xi::Vector)
basis.basis(xi)
end
function call(basis::Basis, xi::Vector, ::Type{Val{:grad}})
function (basis::Basis)(xi::Vector, ::Type{Val{:grad}})
basis.dbasis(xi)
end
@@ -170,10 +170,6 @@ function push!(field::DVTV, data::Pair)
push!(field.data, data)
end
function getindex(field::DVTV, i::Int64)
return field.data[i]
end
function getindex(field::DVTI, i::Int64)
return field.data[i]
end
@@ -218,19 +214,19 @@ for op = (:+, :*, :/, :-)
@eval ($op)(k::Number, field::DCTI) = ($op)(field.data, k)
end
function Base.(:+)(f1::DVTI, f2::DVTI)
function Base.:+(f1::DVTI, f2::DVTI)
return DVTI(f1.data + f2.data)
end
function Base.(:-)(f1::DVTI, f2::DVTI)
function Base.:-(f1::DVTI, f2::DVTI)
return DVTI(f1.data - f2.data)
end
function Base.(:*){T<:Real}(c::T, field::DVTI)
function Base.:*{T<:Real}(c::T, field::DVTI)
return DVTI(c*field.data)
end
function Base.(:*)(N::Matrix, f::DCTI)
function Base.:*(N::Matrix, f::DCTI)
return f.data*N'
end
@@ -240,7 +236,7 @@ end
# 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)
function Base.:*(T::Vector, f::DVTI)
@assert length(T) == length(f)
return sum([T[i]*f[i] for i=1:length(f)])
end
@@ -311,19 +307,19 @@ end
### Accessing continuous fields
function call(field::CVTI, xi::Vector)
function (field::CVTI)(xi::Vector)
return field.data(xi)
end
function call(field::CVTV, xi, time::Float64)
function (field::CVTV)(xi, time::Float64)
return field.data(xi, time)
end
function call(field::CVTI, xi::Vector, ::Type{Val{:Grad}})
function (field::CVTI)(xi::Vector, ::Type{Val{:Grad}})
return field.data(xi, Val{:Grad})
end
function call(field::CCTV, time::Float64)
function (field::CCTV)(time::Float64)
return field.data(time)
end
@@ -334,21 +330,21 @@ end
### Interpolation
""" Interpolate time-invariant field in time direction. """
function call(field::DVTI, time::Float64)
function (field::DVTI)(time::Float64)
return field
end
function call(field::DCTI, time::Float64)
function (field::DCTI)(time::Float64)
return field
end
function call(field::CVTI, time::Float64)
function (field::CVTI)(time::Float64)
return field.data()
end
function call(field::CCTI, time::Float64)
function (field::CCTI)(time::Float64)
return field.data()
end
""" Interpolate constant time-variant field in time direction. """
function call(field::DCTV, time::Real)
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))
@@ -368,7 +364,7 @@ function call(field::DCTV, time::Real)
error("interpolate DCTV: unknown failure when interpolating $(field.data) for time $time")
end
function call(field::DVTV, time::Float64)
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))
@@ -389,17 +385,17 @@ function call(field::DVTV, time::Float64)
end
""" Interpolate constant field in spatial dimension. """
function call(basis::CVTI, field::DCTI, xi::Vector)
function (basis::CVTI)(field::DCTI, xi::Vector)
return field.data
end
""" Interpolate variable field in spatial dimension. """
function call(basis::CVTI, values::DVTI, xi::Vector)
function (basis::CVTI)(values::DVTI, xi::Vector)
N = basis(xi)
return sum([N[i]*values[i] for i=1:length(N)])
end
function call(basis::CVTI, geometry::DVTI, xi::Vector, ::Type{Val{:grad}})
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)])
@@ -408,18 +404,18 @@ function call(basis::CVTI, geometry::DVTI, xi::Vector, ::Type{Val{:grad}})
return grad
end
function call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}})
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 call(basis::CVTI, xi::Vector, time::Number)
function (basis::CVTI)(xi::Vector, time::Number)
basis(xi)
end
function Base.(:*)(grad::Matrix, field::DVTI)
function Base.:*(grad::Matrix, field::DVTI)
return sum([kron(grad[:,i], field[i]') for i=1:length(field)])'
end
src/postprocess_utils.jl
+10 -9
View File
@@ -12,7 +12,7 @@ using Formatting
import HDF5: h5read, h5write
function h5read{T<:DataFrame}(::Type{T}, filename, name::ByteString)
function h5read{T<:DataFrame}(::Type{T}, filename, name::String)
raw_data = h5read(filename, name)
index = raw_data["index"]
column_names = raw_data["column_names"]
@@ -25,7 +25,7 @@ function h5read{T<:DataFrame}(::Type{T}, filename, name::ByteString)
return DataFrame(data, column_names)
end
function h5write(filename, name::ByteString, data::DataFrame)
function h5write(filename, name::String, data::DataFrame)
column_names = DataFrames._names(data)
column_names = map(string, column_names)
index = convert(Vector, data[:,1])
@@ -216,13 +216,13 @@ function to_dataframe(u::Dict, abbreviation::Symbol)
return df
end
function call(problem::Problem, ::Type{DataFrame}, field_name::AbstractString,
function (problem::Problem)(::Type{DataFrame}, field_name::AbstractString,
abbreviation::Symbol, time::Float64=0.0)
u = problem(field_name, time)
return to_dataframe(u, abbreviation)
end
function call(solver::Solver, ::Type{DataFrame}, field_name::AbstractString,
function (solver::Solver)(::Type{DataFrame}, field_name::AbstractString,
abbreviation::Symbol, time::Float64=0.0)
fields = [problem(field_name, time) for problem in get_problems(solver)]
fields = filter(f -> f != nothing, fields)
@@ -249,9 +249,9 @@ function get_components(n, m)
end
end
#=
""" Return T in integration points. """
function call{T}(problem::Problem, ::Type{DataFrame}, element::Element, time::Float64,
::Type{Val{T}})
function call{T}(problem::Problem, ::Type{DataFrame}, element::Element, time::Float64, ::Type{Val{T}})
column_names = [:ELEMENT, :IP]
ips = get_integration_points(element)
field = Any[problem(element, ip, time, Val{T}) for ip in ips]
@@ -300,9 +300,10 @@ function call{T}(solver::Solver, ::Type{DataFrame}, time::Float64, ::Type{Val{T}
results = [tables...;]
return results
end
=#
""" Interpolate field from a set of elements. """
function call(problem::Problem, field_name::AbstractString, X::Vector, time::Float64=0.0; fillna=NaN)
function (problem::Problem)(field_name::AbstractString, X::Vector, time::Float64=0.0; fillna=NaN)
for element in get_elements(problem)
if inside(element, X, time)
xi = get_local_coordinates(element, X, time)
@@ -313,7 +314,7 @@ function call(problem::Problem, field_name::AbstractString, X::Vector, time::Flo
end
""" Interpolate field from a set of elements. """
function call(problem::Problem, field_name::AbstractString, X::Vector, time::Float64, ::Type{Val{:Grad}}; fillna=NaN)
function (problem::Problem)(field_name::AbstractString, X::Vector, time::Float64, ::Type{Val{:Grad}}; fillna=NaN)
for element in get_elements(problem)
if inside(element, X, time)
xi = get_local_coordinates(element, X, time)
@@ -323,7 +324,7 @@ function call(problem::Problem, field_name::AbstractString, X::Vector, time::Flo
return fillna
end
function call(solver::Solver, field_name::AbstractString, X::Vector, time::Float64; fillna=NaN)
function (solver::Solver)(field_name::AbstractString, X::Vector, time::Float64; fillna=NaN)
for problem in get_problems(solver)
for element in get_elements(problem)
if inside(element, X, time)
src/preprocess.jl
+1 -1
View File
@@ -93,7 +93,7 @@ end
function create_element(mesh::Mesh, id::Int)
connectivity = mesh.elements[id]
element_type = JuliaFEM.(mesh.element_types[id])
element_type = getfield(JuliaFEM, mesh.element_types[id])
element = Element(element_type, connectivity)
update!(element, "geometry", mesh.nodes)
element.id = id
src/preprocess_abaqus_reader.jl
-1
View File
@@ -12,7 +12,6 @@ element_has_type( ::Type{Val{:C3D8}}) = :Hex8
element_has_nodes(::Type{Val{:C3D10}}) = 10
element_has_type(::Type{Val{:C3D10}}) = :Tet10
element_has_nodes(::Type{Val{:C3D10}}) = 10
element_has_nodes(::Type{Val{:C3D20}}) = 20
element_has_nodes(::Type{Val{:C3D20E}}) = 20
src/preprocess_aster_reader.jl
+5 -5
View File
@@ -60,7 +60,7 @@ type MEDFile
data :: Dict
end
function MEDFile(fn)
function MEDFile(fn::String)
return MEDFile(h5read(fn, "/"))
end
@@ -116,7 +116,7 @@ function get_element_sets(med::MEDFile, mesh_name)
for elset in keys(elsets)
k = split(elset, '_')
elset_id = parse(Int, k[2])
elset_name = ascii(pointer(convert(Vector{UInt8}, elsets[elset]["GRO"]["NOM"][1])))
elset_name = ascii(unsafe_string(pointer(convert(Vector{UInt8}, elsets[elset]["GRO"]["NOM"][1]))))
es[elset_id] = Symbol(elset_name)
end
return es
@@ -260,7 +260,7 @@ type RMEDFile
data :: Dict
end
function RMEDFile(fn)
function RMEDFile(fn::String)
return RMEDFile(h5read(fn, "/"))
end
@@ -280,7 +280,7 @@ function aster_read_nodes(rmed::RMEDFile)
# INFO: quite safe assumption is that id is in node name, i.e. N1 => 1, N123 => 123
node_id(node_name) = parse(matchall(r"\d+", node_name)[1])
node_ids = map(node_id, node_names)
nodes = Dict([j => node_coords[:,j] for j in node_ids])
nodes = Dict(j => node_coords[:,j] for j in node_ids)
return nodes
end
@@ -308,7 +308,7 @@ function aster_read_data(rmed::RMEDFile, field_name; field_type=:NODE,
increment = chdata[first(keys(chdata))]
if field_type == :NODE
data = increment["NOE"]["MED_NO_PROFILE_INTERNAL"]["CO"]
results = Dict([j => data[j] for j in node_ids])
results = Dict(j => data[j] for j in node_ids)
else
error("Unable to read result of type $field_type: not implemented")
end
src/problems.jl
+1 -1
View File
@@ -301,7 +301,7 @@ function getindex(problem::Problem, field_name::AbstractString)
end
""" Return field calculated to nodal points for elements in problem p. """
function call(problem::Problem, field_name::AbstractString, time::Float64=0.0)
function (problem::Problem)(field_name::AbstractString, time::Float64=0.0)
#if haskey(problem, field_name)
# return problem[field_name](time)
#end
src/problems_contact_2d_autodiff.jl
+2 -2
View File
@@ -273,8 +273,8 @@ function assemble!(problem::Problem{Contact}, time::Float64,
b = calculate_interface(x)
A = sparse(A)
b = sparse(b)
SparseMatrix.droptol!(A, 1.0e-12)
SparseMatrix.droptol!(b, 1.0e-12)
SparseArrays.droptol!(A, 1.0e-9)
SparseArrays.droptol!(b, 1.0e-9)
ndofs = round(Int, length(x)/2)
K = A[1:ndofs,1:ndofs]
src/problems_elasticity.jl
+3 -3
View File
@@ -834,14 +834,14 @@ end
=#
function call(problem::Problem, element::Element, ip, time::Float64, ::Type{Val{:E}})
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 call(problem::Problem, element::Element, ip, time::Float64, ::Type{Val{:S}})
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})
@@ -857,7 +857,7 @@ function call(problem::Problem, element::Element, ip, time::Float64, ::Type{Val{
return S
end
function call(problem::Problem, element::Element, ip, time::Float64, ::Type{Val{:COORD}})
function (problem::Problem)(element::Element, ip, time::Float64, ::Type{Val{:COORD}})
haskey(element, "geometry") || return nothing
return element("geometry", ip, time)
end
src/problems_mortar_2d_autodiff.jl
+4 -4
View File
@@ -237,8 +237,8 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
A = sparse(A)
b = sparse(b)
SparseMatrix.droptol!(A, 1.0e-12)
SparseMatrix.droptol!(b, 1.0e-12)
SparseArrays.droptol!(A, 1.0e-12)
SparseArrays.droptol!(b, 1.0e-12)
K = A[1:ndofs,1:ndofs]
C1 = transpose(A[1:ndofs,ndofs+1:end])
@@ -507,8 +507,8 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}})
A = sparse(A)
b = sparse(b)
SparseMatrix.droptol!(A, 1.0e-12)
SparseMatrix.droptol!(b, 1.0e-12)
SparseArrays.droptol!(A, 1.0e-12)
SparseArrays.droptol!(b, 1.0e-12)
K = A[1:ndofs,1:ndofs]
C1 = transpose(A[1:ndofs,ndofs+1:end])
src/problems_mortar_3d.jl
+16 -3
View File
@@ -282,6 +282,11 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
nm = length(master_element)
X2 = master_element("geometry", time)
if norm(mean(X1) - mean(X2)) > problem.properties.distval
# elements are "far enough"
continue
end
# 3.1 project master nodes to auxiliary plane and create polygon clipping
M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
P = get_polygon_clip(S, M, n0)
@@ -365,9 +370,6 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
# 6. add contribution to contact virtual work
sdofs = get_gdofs(problem, slave_element)
mdofs = get_gdofs(problem, master_element)
if problem.properties.dual_basis
De[abs(De) .< 1.0e-11] = 0
end
for i=1:field_dim
lsdofs = sdofs[i:field_dim:end]
@@ -382,6 +384,17 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{
end # master elements done
end # slave elements done, contact virtual work ready
if problem.properties.dual_basis
tol = 1.0e-9
debug && info("Dual basis is used, dropping small values for C1 & C2, tol = $tol")
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
SparseArrays.droptol!(C1, tol)
SparseArrays.droptol!(C2, tol)
problem.assembly.C1 = C1
problem.assembly.C2 = C2
end
debug && info("area of interface: $area")
src/solvers.jl
+14 -6
View File
@@ -216,7 +216,7 @@ function create_projection(C::SparseMatrixCSC, g; S=nothing, tol=1.0e-12)
resize!(P, n, m)
resize!(h, n, 1)
P = speye(n) - P
SparseMatrix.droptol!(P, tol)
droptol!(P, tol)
return P, h
end
@@ -468,8 +468,16 @@ function filter_by_element_type(element_type, elements)
return filter(element -> is_element_type(element, element_type), elements)
end
function call(solver::Solver, field_name::AbstractString, time::Float64)
fields = [problem(field_name, time) for problem in get_problems(solver)]
function (solver::Solver)(field_name::AbstractString, time::Float64)
fields = []
for problem in get_problems(solver)
field = problem(field_name, time)
if length(field) == 0
warn("no field $field_name found for problem $(problem.name)")
else
push!(fields, field)
end
end
return merge(fields...)
end
@@ -666,7 +674,7 @@ function Base.showerror(io::IO, exception::NonlinearConvergenceError)
end
""" Default solver for quasistatic nonlinear problems. """
function call(solver::Solver{Nonlinear}; show_info=true)
function (solver::Solver{Nonlinear})(; show_info=true)
properties = solver.properties
@@ -747,7 +755,7 @@ function assemble!(solver::Solver{Linear}; show_info=true)
show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
end
function call(solver::Solver{Linear}; show_info=true)
function (solver::Solver{Linear})(; show_info=true)
t0 = Base.time()
show_info && info(repeat("-", 80))
show_info && info("Starting linear solver")
@@ -807,7 +815,7 @@ function assemble!(solver::Solver{Postprocessor}; show_info=true)
show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
end
function call(solver::Solver{Postprocessor}; show_info=true)
function (solver::Solver{Postprocessor})(; show_info=true)
t0 = Base.time()
show_info && info(repeat("-", 80))
show_info && info("Starting postprocessor")
src/solvers_modal.jl
+70 -29
View File
@@ -24,7 +24,9 @@ function Modal(nev=10, which=:SM)
end
""" Eliminate Dirichlet boundary condition from matrices K, M. """
function eliminate_boundary_conditions!(K_red, M_red, problem::Problem{Dirichlet}, ndim::Int)
function eliminate_boundary_conditions!(K_red::SparseMatrixCSC,
M_red::SparseMatrixCSC,
problem::Problem{Dirichlet}, ndim::Int)
K = sparse(problem.assembly.K, ndim, ndim)
C1 = sparse(problem.assembly.C1, ndim, ndim)
C2 = sparse(problem.assembly.C2, ndim, ndim)
@@ -64,7 +66,7 @@ function calc_projection(problem::Problem{Mortar}, ndim::Int)
@assert nnz(sparse(problem.assembly.g)) == 0
@assert C1 == C2
@assert problem.properties.dual_basis == true
#@assert problem.properties.dual_basis == true
@assert problem.properties.adjust == false
# determine master and slave dofs
@@ -96,13 +98,15 @@ function calc_projection(problem::Problem{Mortar}, ndim::Int)
# Construct matrix P = D^-1*M
D_ = C2[S,S]
M_ = -C2[S,M]
P = nothing
if !isdiag(D_)
info("D is not diagonal, is dual basis used?")
println(D_)
warn("D is not diagonal, is dual basis used? This might take a long time.")
P = ldltfact(1/2*(D_ + D_')) \ M_
else
P = D_ \ M_
end
@assert isdiag(D_)
P = D_ \ M_
info("Projection P ready.")
info("Matrix P ready.")
return S, M, P
end
@@ -124,7 +128,7 @@ function eliminate_boundary_conditions!(K_red::SparseMatrixCSC,
@assert nnz(sparse(problem.assembly.g)) == 0
@assert C1 == C2
@assert problem.properties.dual_basis == true
#@assert problem.properties.dual_basis == true
@assert problem.properties.adjust == false
info("Eliminating mesh tie constraint $(problem.name) using static condensation")
@@ -160,30 +164,64 @@ function eliminate_boundary_conditions!(K_red::SparseMatrixCSC,
# Construct matrix P = D^-1*M
D_ = C2[S,S]
M_ = -C2[S,M]
if !isdiag(D_)
info("D is not diagonal, is dual basis used?")
println(D_)
end
@assert isdiag(D_)
P = D_ \ M_
info("Projection P ready.")
K_red[N,M] += K_red[N,S]*P
K_red[M,N] += P'*K_red[S,N]
K_red[M,M] += P'*K_red[S,S]*P
P = nothing
if !isdiag(D_)
warn("D is not diagonal, is dual basis used? This might take a long time.")
P = ldltfact(1/2*(D_ + D_')) \ M_
else
P = D_ \ M_
end
#@assert isdiag(D_)
info("Matrix P ready.")
Id = ones(ndim)
#Id[S] = 0
#Id[M] = 0
Q = spdiagm(Id)
Q[M,S] += P'
info("Matrix Q ready.")
# testing
#K_red_orig = copy(K_red)
#K_red[N,M] += K_red[N,S]*P
#K_red[M,N] += P'*K_red[S,N]
#K_red[M,M] += P'*K_red[S,S]*P
#K_red[S,:] = 0.0
#K_red[:,S] = 0.0
info("K transform")
K_red[:,:] = Q*K_red*Q'
K_red[S,:] = 0.0
K_red[:,S] = 0.0
M_red[N,M] += M_red[N,S]*P
M_red[M,N] += P'*M_red[S,N]
M_red[M,M] += P'*M_red[S,S]*P
#K_res = K_red - K_red_2
#SparseArrays.droptol!(K_res, 1.0e-9)
info("K transform ready")
#info("Create matrices, M")
#info("Sum matricse, M")
#M_red_orig = copy(M_red)
#M_red[N,M] += M_red[N,S]*P
#M_red[M,N] += P'*M_red[S,N]
#M_red[M,M] += P'*M_red[S,S]*P
#M_red[S,:] = 0.0
#M_red[:,S] = 0.0
info("M transform")
M_red[:,:] = Q*M_red*Q'
M_red[S,:] = 0.0
M_red[:,S] = 0.0
info("M transform ready")
#M_res = M_red - M_red_2
#SparseArrays.droptol!(M_res, 1.0e-9)
#info("Diff")
#println(K_res)
#println(M_res)
return true
end
function call(solver::Solver{Modal}; show_info=true, debug=false,
function (solver::Solver{Modal})(; show_info=true, debug=false,
bc_invertible=false, P=nothing, symmetric=true,
empty_assemblies_before_solution=true, dense=false,
real_eigenvalues=true, positive_eigenvalues=true)
@@ -227,8 +265,8 @@ function call(solver::Solver{Modal}; show_info=true, debug=false,
gc()
end
SparseArrays.droptol!(K_red, 1.0e-12)
SparseArrays.droptol!(M_red, 1.0e-12)
SparseArrays.droptol!(K_red, 1.0e-9)
SparseArrays.droptol!(M_red, 1.0e-9)
nz = get_nonzero_rows(K_red)
K_red = K_red[nz,nz]
M_red = M_red[nz,nz]
@@ -254,8 +292,8 @@ function call(solver::Solver{Modal}; show_info=true, debug=false,
M_red = 1/2*(M_red + transpose(M_red))
end
info("is K symmetric? ", issym(K_red))
info("is M symmetric? ", issym(M_red))
info("is K symmetric? ", issymmetric(K_red))
info("is M symmetric? ", issymmetric(M_red))
info("is K positive definite? ", isposdef(K_red))
info("is M positive definite? ", isposdef(M_red))
@@ -362,13 +400,14 @@ function update_xdmf!(solver::Solver{Modal}; show_info=true)
# 2. save topology
nid_mapping = Dict([j => i for (i, j) in enumerate(node_ids)])
nid_mapping = Dict(j => i for (i, j) in enumerate(node_ids))
all_elements = get_all_elements(solver)
nelements = length(all_elements)
element_types = unique(map(get_element_type, all_elements))
xdmf_element_mapping = Dict(
"Poi1" => "Polyvertex",
"Seg2" => "Polyline",
"Tri3" => "Triangle",
"Quad4" => "Quadrilateral",
@@ -386,6 +425,7 @@ function update_xdmf!(solver::Solver{Modal}; show_info=true)
topology = []
for element_type in element_types
info("Xdmf save: element type $element_type")
elements = filter_by_element_type(element_type, all_elements)
sort!(elements, by=get_element_id)
#elements = elements[1:5]
@@ -412,7 +452,7 @@ function update_xdmf!(solver::Solver{Modal}; show_info=true)
set_attribute(topology_, "NumberOfElements", length(elements))
add_child(topology_, dataitem)
push!(topology, topology_)
break
#break
end
# save modes
@@ -442,6 +482,7 @@ function update_xdmf!(solver::Solver{Modal}; show_info=true)
unknown_field_name = ucfirst(unknown_field_name)
freqn = freqs[j]
path = "/Results/Frequency $freqn/Nodal Fields/$unknown_field_name"
info("Storing data to $path")
dataitem = new_dataitem(xdmf, path, mode)
attribute = new_child(frame, "Attribute")
set_attribute(attribute, "Name", unknown_field_name)
src/sparse.jl
+2 -2
View File
@@ -35,7 +35,7 @@ tol
"""
function sparse(A::SparseMatrixCOO, args...; tol=1.0e-12)
B = sparse(A.I, A.J, A.V, args...)
SparseMatrix.droptol!(B, tol)
SparseArrays.droptol!(B, tol)
return B
end
@@ -67,7 +67,7 @@ function isempty(A::SparseMatrixCOO)
return isempty(A.I) && isempty(A.J) && isempty(A.V)
end
function Base.(:+)(A::SparseMatrixCOO, B::SparseMatrixCOO)
function Base.:+(A::SparseMatrixCOO, B::SparseMatrixCOO)
if isempty(A)
return B
end
src/types.jl
+1 -1
View File
@@ -29,7 +29,7 @@ function haskey(point::Point, field_name)
return haskey(point.fields, field_name)
end
function call(point::Point, field_name, time=0.0)
function (point::Point)(field_name, time=0.0)
point.fields[field_name](time).data
end
test/test_elasticity_2d_linear_with_surface_load.jl
+3 -4
View File
@@ -47,7 +47,7 @@ using JuliaFEM.Testing
bc_sym_13.elements = create_elements(mesh, "BOTTOM")
update!(bc_sym_13, "displacement 2", 0.0)
solver = LinearSolver(block, traction, bc_sym_23, bc_sym_13)
solver = Solver(Linear, block, traction, bc_sym_23, bc_sym_13)
solver()
info("u = ", block.assembly.u)
@@ -58,7 +58,7 @@ using JuliaFEM.Testing
E = 288.0
nu = 1/3
u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1]
#=
# fetch nodal results X + u and join them into one table using DataFrames
X = solver(DataFrame, "geometry", :COOR)
u = solver(DataFrame, "displacement", :U)
@@ -90,11 +90,9 @@ using JuliaFEM.Testing
S1 = block(DataFrame, 0.0, Val{:S})
E1 = block(DataFrame, 0.0, Val{:E})
C1 = block(DataFrame, 0.0, Val{:COORD})
println(S1)
println(E1)
println(C1)
S = solver(DataFrame, 0.0, Val{:S})
println(S)
@@ -111,6 +109,7 @@ using JuliaFEM.Testing
# u = solver2("displacement", 0.0)[3]
# info("nlsolver u3 = $u, expected = $u3_expected")
# @test isapprox(u, u3_expected; rtol=1.0e-5)
=#
end
#= TODO: to other file