From 84bb370cea4bc6a5957ddbf351ef7082bf587b3d Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 17 Jun 2016 02:10:17 +0300 Subject: [PATCH] first version of modal analysis for natural frequencies --- src/JuliaFEM.jl | 12 +- src/assembly.jl | 23 +++- src/dirichlet.jl | 4 +- src/elasticity.jl | 65 +++++----- src/elements.jl | 10 ++ src/heat.jl | 28 ++-- src/integrate.jl | 34 ++++- src/modal.jl | 62 +++++++++ src/optics.jl | 22 ++-- src/problems.jl | 13 +- src/solvers.jl | 122 ++++++++++++------ src/sparse.jl | 44 ++++--- ..._elasticity_2d_linear_with_surface_load.jl | 12 +- test/test_modal_analysis.jl | 42 ++++++ 14 files changed, 362 insertions(+), 131 deletions(-) create mode 100644 src/modal.jl create mode 100644 test/test_modal_analysis.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 9d4828c..b687b44 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -41,7 +41,8 @@ export add!, SparseMatrixCOO, get_nonzero_rows include("problems.jl") # common problem routines export Problem, AbstractProblem, FieldProblem, BoundaryProblem, - get_unknown_field_dimension, get_gdofs, Assembly + get_unknown_field_dimension, get_gdofs, Assembly, + get_parent_field_name, get_elements include("elasticity.jl") # elasticity equations export Elasticity @@ -62,10 +63,15 @@ end include("assembly.jl") include("solver_utils.jl") include("solvers.jl") -export Solver +export AbstractSolver, Solver, Nonlinear, + get_unknown_field_name, get_formulation_type, + get_field_problems, get_boundary_problems, + get_field_assembly, get_boundary_assembly +include("modal.jl") +export Modal include("optics.jl") -export find_intersection, calc_reflection +export find_intersection, calc_reflection, calc_normal ### MORTAR STUFF ### include("mortar.jl") # mortar projection diff --git a/src/assembly.jl b/src/assembly.jl index fe5ea4d..3497c19 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -49,7 +49,28 @@ function assemble!(problem::Problem, time::Real; empty_assembly::Bool=true) if method_exists(assemble_posthook!, Tuple{typeof(problem), Real}) assemble_posthook!(problem, time) end - return problem.assembly + return +end + +function assemble!(problem::Problem, time::Real, ::Type{Val{:mass_matrix}}) + !isempty(problem.assembly.M) && return # assembly mass matrix only once + dim = get_unknown_field_dimension(problem) + for element in get_elements(problem) + haskey(element, "density") || error("Failed to assemble mass matrix, density not defined!") + nnodes = length(element) + M = zeros(nnodes, nnodes) + for ip in get_integration_points(element, 1) + detJ = element(ip, time, Val{:detJ}) + N = element(ip, time) + rho = element("density", ip, time) + M += ip.weight*rho*N'*N*detJ + end + gdofs = get_gdofs(problem, element) + for j=1:dim + ldofs = gdofs[j:dim:end] + add!(problem.assembly.M, ldofs, ldofs, M) + end + end end """ Calculate reduced stiffness matrix. diff --git a/src/dirichlet.jl b/src/dirichlet.jl index bfdf8c5..cd79b25 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -62,8 +62,8 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele ldofs = gdofs[i:field_dim:end] if haskey(element, field_name*" $i") g = element(field_name*" $i", ip, time) - if true - haskey(element, "displacement") || continue + # u = u_prev + Δu ⇒ Δu = u - u_prev + if haskey(element, field_name) g_prev = element(field_name, ip, time) g -= g_prev[i] end diff --git a/src/elasticity.jl b/src/elasticity.jl index 3ddce32..e373d2d 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -27,22 +27,14 @@ end function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time=0.0) props = problem.properties gdofs = get_gdofs(problem, element) - if problem.properties.formulation in [:plane_stress, :plane_strain] - Kt, f = assemble(problem, element, time, Val{:plane}) - add!(assembly.K, gdofs, gdofs, Kt) - add!(assembly.f, gdofs, f) - return - elseif problem.properties.formulation in [:continuum_buckling] - Km, Kg = assemble(problem, element, time, Val{:continuum_buckling}) - add!(assembly.K, gdofs, gdofs, Km) - add!(assembly.Kg, gdofs, gdofs, Kg) - return - else - Kt, f = assemble(problem, element, time, Val{problem.properties.formulation}) - add!(assembly.K, gdofs, gdofs, Kt) - add!(assembly.f, gdofs, f) - return + formulation = props.formulation + if formulation in [:plane_stress, :plane_strain] + formulation = :plane end + Km, Kg, f = assemble(problem, element, time, Val{formulation}) + add!(assembly.K, gdofs, gdofs, Km) + add!(assembly.Kg, gdofs, gdofs, Kg) + add!(assembly.f, gdofs, f) end """ Elasticity equations for 2d cases. """ @@ -53,7 +45,8 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem nnodes = length(element) BL = zeros(3, dim*nnodes) BNL = zeros(4, dim*nnodes) - Kt = zeros(dim*nnodes, dim*nnodes) + Km = zeros(dim*nnodes, dim*nnodes) + Kg = zeros(dim*nnodes, dim*nnodes) f = zeros(dim*nnodes) for ip in get_integration_points(element) @@ -127,9 +120,9 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem S2[1,2] = S2[2,1] = stress_vec[3] S2[3:4,3:4] = S2[1:2,1:2] - Kt += w*BL'*D*BL # material stiffness + Km += w*BL'*D*BL # material stiffness if props.finite_strain # add geometric stiffness - Kt += w*BNL'*S2*BNL # geometric stiffness + Kg += w*BNL'*S2*BNL # geometric stiffness end if get_formulation_type(problem) == :incremental @@ -150,7 +143,7 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem end - return Kt, f + return Km, Kg, f end function assemble{El<:Union{Seg2,Seg3}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:plane}}) @@ -158,7 +151,8 @@ function assemble{El<:Union{Seg2,Seg3}}(problem::Problem{Elasticity}, element::E props = problem.properties dim = get_unknown_field_dimension(problem) nnodes = size(element, 2) - Kt = zeros(dim*nnodes, dim*nnodes) + Km = zeros(dim*nnodes, dim*nnodes) + Kg = zeros(dim*nnodes, dim*nnodes) f = zeros(dim*nnodes) for ip in get_integration_points(element) @@ -189,7 +183,7 @@ function assemble{El<:Union{Seg2,Seg3}}(problem::Problem{Elasticity}, element::E end - return Kt, f + return Km, Kg, f end """ Elasticity equations, 3d, linear. """ @@ -200,7 +194,8 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el nnodes = length(element) ndofs = dim*nnodes BL = zeros(6, ndofs) - Kt = zeros(ndofs, ndofs) + Km = zeros(ndofs, ndofs) + Kg = zeros(ndofs, ndofs) f = zeros(ndofs) for ip in get_integration_points(element) @@ -233,7 +228,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el 0.0 0.0 0.0 0.0 0.5-nu 0.0 0.0 0.0 0.0 0.0 0.0 0.5-nu] - Kt += w*BL'*D*BL + Km += w*BL'*D*BL if haskey(element, "displacement load") T = element("displacement load", ip, time) @@ -255,7 +250,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el end end - return Kt, f + return Km, Kg, f end """ Material and geometric stiffness for linear buckling analysis. """ @@ -269,6 +264,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el 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}) @@ -332,7 +328,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el end - return Km, Kg + return Km, Kg, f end """ Elasticity equations, 3d nonlinear. """ @@ -344,7 +340,8 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el ndofs = dim*nnodes BL = zeros(6, ndofs) BNL = zeros(9, ndofs) - Kt = zeros(ndofs, ndofs) + Km = zeros(ndofs, ndofs) + Kg = zeros(ndofs, ndofs) f = zeros(ndofs) for ip in get_integration_points(element) @@ -410,7 +407,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el update!(ip, "strain", time => strain_vec) update!(ip, "stress", time => stress_vec) - Kt += w*BL'*D*BL + Km += w*BL'*D*BL if get_formulation_type(problem) == :incremental f -= w*BL'*stress_vec @@ -442,7 +439,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3] if props.finite_strain - Kt += w*BNL'*S3*BNL + Kg += w*BNL'*S3*BNL end # geometric stiffness end @@ -464,7 +461,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el end - return Kt, f + return Km, Kg, f end """ Elasticity equations, surface traction for continuum formulation. """ @@ -473,7 +470,8 @@ function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, el props = problem.properties dim = get_unknown_field_dimension(problem) nnodes = size(element, 2) - Kt = zeros(dim*nnodes, dim*nnodes) + Km = zeros(dim*nnodes, dim*nnodes) + Kg = zeros(dim*nnodes, dim*nnodes) f = zeros(dim*nnodes) for ip in get_integration_points(element) @@ -498,7 +496,7 @@ function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, el f += w*p*vec(n*N) end end - return Kt, f + return Km, Kg, f end function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}}) @@ -604,10 +602,11 @@ function assemble(problem::Problem{Elasticity}, element::Element, time::Real, :: end field = element("displacement", time) - Kt, allresults = ForwardDiff.jacobian(get_residual_vector, vec(field), + Km, allresults = ForwardDiff.jacobian(get_residual_vector, vec(field), AllResults, cache=autodiffcache) + Kg = zeros(Km) f = -ForwardDiff.value(allresults) - return Kt, f + return Km, Kg, f end diff --git a/src/elements.jl b/src/elements.jl index 9f40439..734b832 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -125,6 +125,16 @@ function update!(element::Element, field_name::ASCIIString, datas::Union{Real, V end end +function update!(element::Element, datas::Pair...) + for (field_name, data) in datas + if haskey(element, field_name) + update!(element[field_name], data) + else + element[field_name] = data + end + end +end + function update!(elements::Vector, field_name::ASCIIString, data) for element in elements update!(element, field_name, data) diff --git a/src/heat.jl b/src/heat.jl index 0bf04fb..0f5d064 100644 --- a/src/heat.jl +++ b/src/heat.jl @@ -3,17 +3,6 @@ # Heat problems -type Heat <: FieldProblem -end - -function get_unknown_field_name(problem::Problem{Heat}) - return "temperature" -end - -function get_unknown_field_type(problem::Problem{Heat}) - return Float64 -end - """ Heat equations. Formulation @@ -29,15 +18,26 @@ Weak form is: find u∈U such that ∀v in V where - k = temperature thermal conductivity defined on volume - f = temperature load defined on volume - g = temperature flux defined on boundary + k = temperature thermal conductivity defined on volume elements + f = temperature load defined on volume elements + g = temperature flux defined on boundary elements References ---------- https://en.wikipedia.org/wiki/Heat_equation """ +type Heat <: FieldProblem +end + +function get_unknown_field_name(problem::Problem{Heat}) + return "temperature" +end + +function get_unknown_field_type(problem::Problem{Heat}) + return Float64 +end + function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time=0.0) gdofs = get_gdofs(problem, element) diff --git a/src/integrate.jl b/src/integrate.jl index 916a586..2041925 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -172,6 +172,38 @@ function get_integration_points(element::TetrahedralElement, ::Type{Val{3}}) return zip(weights, points) end +function get_integration_points(element::TetrahedralElement, ::Type{Val{4}}) + a = 0.25 + b1 = 1.0/34.0*(7.0 + sqrt(15.0)) + b2 = 1.0/34.0*(7.0 - sqrt(15.0)) + c1 = 1.0/34.0*(13.0 - 3.0*sqrt(15.0)) + c2 = 1.0/34.0*(13.0 + 3.0*sqrt(15.0)) + d = 1.0/20.0*(5.0 - sqrt(15.0)) + e = 1.0/20.0*(5.0 + sqrt(15.0)) + w1 = 8.0/405.0 + w2 = (2665.0 - 14.0*sqrt(15.0))/226800.0 + w3 = (2665.0 + 14.0*sqrt(15.0))/226800.0 + w4 = 5.0/567.0 + weights = [w1, w2, w2, w2, w2, w3, w3, w3, w3, w4, w4, w4, w4, w4] + points = Vector{Float64}[ + [a, a, a], + [b1, b1, b1], + [b1, b1, c1], + [b1, c1, b1], + [c1, b1, b1], + [b2, b2, b2], + [b2, b2, c2], + [b2, c2, b2], + [c2, b2, b2], + [d, d, e], + [d, e, d], + [e, d, d], + [d, e, e], + [e, d, e], + [e, e, d]] + return zip(weights, points) +end + function get_integration_points(element::Union{TriangularElement, TetrahedralElement}, order::Int64) return get_integration_points(element, Val{order}) end @@ -197,7 +229,7 @@ function get_integration_points(element::LinearElement) end function get_integration_points(element::QuadraticElement) - order= get_integration_order(element) + order = get_integration_order(element) get_integration_points(element, order) end diff --git a/src/modal.jl b/src/modal.jl new file mode 100644 index 0000000..2464b5a --- /dev/null +++ b/src/modal.jl @@ -0,0 +1,62 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +type Modal <: AbstractSolver + geometric_stiffness :: Bool + eigvals :: Vector + eigvecs :: Matrix + nev :: Int + which :: Symbol +end + +function Modal(nev=10, which=:SM) + solver = Modal(false, Vector(), Matrix(), nev, which) +end + +function call(solver::Solver{Modal}; debug=false) + # assemble all field problems + info("Assembling problems ...") + tic() + for problem in get_field_problems(solver) + assemble!(problem, solver.time) + assemble!(problem, solver.time, Val{:mass_matrix}) + end + t1 = round(toq(), 2) + info("Assembled in $t1 seconds.") + M, K, Kg, f = get_field_assembly(solver; with_mass_matrix=true) + if solver.properties.geometric_stiffness + K += Kg + end + for problem in get_boundary_problems(solver) + assemble!(problem, solver.time) + # FIXME: Check for tie contacts and rhs. Here we just + # remove all fixed dofs giving funny results if problem + # is having MPCs or non-homogeneous Dirichlet conditions +# eliminate!(M, K, Kg, f, problem) + fixed_dofs = get_nonzero_rows(problem.assembly.C2) + K[fixed_dofs, :] = 0 + K[:, fixed_dofs] = 0 + M[fixed_dofs, :] = 0 + M[:, fixed_dofs] = 0 + f[fixed_dofs, :] = 0 + end + fd = get_nonzero_rows(K) + ndofs = solver.ndofs + props = solver.properties + info("Calculate $(props.nev) eigenvalues...") + if debug + info("Stiffness matrix:") + dump(round(full(K[fd, fd]))) + info("Mass matrix:") + dump(round(full(M[fd, fd]))) + end + tic() + om2, X = eigs(K[fd, fd], M[fd, fd]; nev=props.nev, which=props.which) + props.eigvals = om2 + props.eigvecs = zeros(ndofs, length(om2)) + props.eigvecs[fd, :] = X + t1 = round(toq(), 2) + info("Eigenvalues computed in $t1 seconds. Eigenvalues: $om2") + return true +end + diff --git a/src/optics.jl b/src/optics.jl index b15bafe..700e3f1 100644 --- a/src/optics.jl +++ b/src/optics.jl @@ -91,12 +91,7 @@ function find_intersection{S<:Union{Seg2, Seg3, NSeg, NSurf}}(element::Element{S return theta[1], theta[2:end] end -""" -References ----------- -[1] http://fp.optics.arizona.edu/optomech/Fall13/Notes/6%20Mirror%20matrices.pdf -""" -function calc_reflection{S<:Union{Seg2, Seg3, NSeg}}(element::Element{S}, xi, k, time; deformed=true) +function calc_normal{S<:Union{Seg2, Seg3, NSeg}}(element::Element{S}, xi, time; deformed=true) x = element["geometry"](time) if deformed && haskey(element, "displacement") x += element["displacement"](time) @@ -105,11 +100,10 @@ function calc_reflection{S<:Union{Seg2, Seg3, NSeg}}(element::Element{S}, xi, k, dN = get_dbasis(element, xi, time) n = Q*(dN*x) n /= norm(n) - k2 = k - 2*vecdot(k, n)*n - return k2 + return n end -function calc_reflection{S<:Union{Quad4, NSurf}}(element::Element{S}, xi, k, time; deformed=true) +function calc_normal{S<:Union{Quad4, NSurf}}(element::Element{S}, xi, time; deformed=true) x = element["geometry"](time) if deformed && haskey(element, "displacement") x += element["displacement"](time) @@ -118,6 +112,16 @@ function calc_reflection{S<:Union{Quad4, NSurf}}(element::Element{S}, xi, k, tim J = transpose(sum([kron(dN[:,i], x[i]') for i=1:length(x)])) n = cross(J[:,1], J[:,2]) n /= norm(n) + return n +end + +""" +References +---------- +[1] http://fp.optics.arizona.edu/optomech/Fall13/Notes/6%20Mirror%20matrices.pdf +""" +function calc_reflection(element::Element, xi, k, time; deformed=true) + n = calc_normal(element, xi, time; deformed=deformed) k2 = k - 2*vecdot(k, n)*n return k2 end diff --git a/src/problems.jl b/src/problems.jl index 8ce244f..baa950c 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -84,9 +84,12 @@ Create vector-valued (dim=3) elasticity problem: julia> prob = Problem(Elasticity, "this is my problem", 3) """ -function Problem{P<:FieldProblem}(::Type{P}, name, dimension, elements=[], dofmap=Dict()) +function Problem{P<:FieldProblem}(::Type{P}, name::ASCIIString, dimension::Int64, elements=[], dofmap=Dict()) Problem{P}(name, dimension, "none", elements, dofmap, Assembly(), P()) end +function Problem{P<:FieldProblem}(::Type{P}, dimension::Int64, elements=[], dofmap=Dict()) + Problem{P}("$P problem", dimension, "none", elements, dofmap, Assembly(), P()) +end """ Construct a new boundary problem. @@ -100,6 +103,12 @@ julia> bc1 = Problem(Dirichlet, "support", 3, "displacement") function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name, elements=[], dofmap=Dict()) Problem{P}(name, dimension, parent_field_name, elements, dofmap, Assembly(), P()) end +function Problem{P<:BoundaryProblem}(::Type{P}, main_problem::Problem, elements=[], dofmap=Dict()) + name = "$P problem" + dimension = get_unknown_field_dimension(main_problem) + parent_field_name = get_unknown_field_name(main_problem) + Problem{P}(name, dimension, parent_field_name, elements, dofmap, Assembly(), P()) +end function get_formulation_type{P<:FieldProblem}(problem::Problem{P}) return :total @@ -170,7 +179,7 @@ function update_assembly!(problem, u, la) # copy current solutions to previous ones and add/replace new solution # TODO: here we have couple of options and they needs to be clarified - # for total formulation we are solving total quantity Ku=f while in + # for total formulation we are solving total quantity Ku = f while in # incremental formulation we solve KΔu = f and u = u + Δu assembly.u_prev = copy(assembly.u) assembly.la_prev = copy(assembly.la) diff --git a/src/solvers.jl b/src/solvers.jl index 671d08d..e0431e4 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -1,36 +1,56 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -type Solver +abstract AbstractSolver + +type Solver{S<:AbstractSolver} name :: ASCIIString # some descriptive name for problem time :: Real # current time + problems :: Vector{Problem} + ndofs :: Int # total dimension of global stiffness matrix, i.e., dim*nnodes + properties :: S +end + +type Nonlinear <: AbstractSolver iteration :: Int # iteration counter norms :: Vector{Tuple} # solution norms for convergence studies - ndofs :: Int # total dimension of global stiffness matrix, i.e., dim*nnodes - problems :: Vector{Problem} + min_iterations :: Int64 + max_iterations :: Int64 + convergence_tolerance :: Float64 + error_if_no_convergence :: Bool is_linear_system :: Bool # setting this to true makes assumption of one step convergence - nonlinear_system_min_iterations :: Int64 - nonlinear_system_max_iterations :: Int64 - nonlinear_system_convergence_tolerance :: Float64 - nonlinear_system_error_if_no_convergence :: Bool linear_system_solver :: Symbol end -function Solver(name::ASCIIString="default solver", time::Real=0.0) - return Solver( - name, - time, - 0, # iteration # +function Nonlinear() + solver = Nonlinear( + 0, # iteration number [], # solution norms in (norm(u), norm(la)) tuples - 0, # ndofs - [], # array of problems - false, # is_linear_system 1, # min nonlinear iterations 10, # max nonlinear iterations 5.0e-5, # nonlinear iteration convergence tolerance true, # throw error if no convergence - :DirectLinearSolver # linear system solution method - ) + false, # is_linear_system + :DirectLinearSolver) # linear system solution method + return solver +end + +function Solver{S<:AbstractSolver}(::Type{S}=Nonlinear, + name::ASCIIString="default solver", + time::Real=0.0, problems=[], + properties...) + variant = S(properties...) + solver = Solver{S}(name, time, problems, 0, variant) + return solver +end + +""" For compatibility. """ +function Solver(name::ASCIIString="default solver", + time::Real=0.0, problems=[], + properties...) + variant = Nonlinear(properties...) + solver = Solver{Nonlinear}(name, time, problems, 0, variant) + return solver end function push!(solver::Solver, problem) @@ -102,26 +122,41 @@ If several field problems exists, they are simply summed together, so problems must have unique node ids. """ -function get_field_assembly(solver::Solver) +function get_field_assembly(solver::Solver; symmetric=true, + with_mass_matrix=false, + empty_after_append=true) problems = get_field_problems(solver) + M = SparseMatrixCOO() K = SparseMatrixCOO() + Kg = SparseMatrixCOO() f = SparseMatrixCOO() for problem in problems append!(K, problem.assembly.K) + append!(Kg, problem.assembly.Kg) append!(f, problem.assembly.f) - empty!(problem.assembly) + with_mass_matrix && append!(M, problem.assembly.M) + empty_after_append && empty!(problem.assembly) + end + if solver.ndofs == 0 + solver.ndofs = size(K, 1) + end + K = sparse(K, solver.ndofs, solver.ndofs) + Kg = sparse(Kg, solver.ndofs, solver.ndofs) + M = sparse(M, solver.ndofs, solver.ndofs) + if symmetric + K = 1/2*(K + K') + Kg = 1/2*(Kg + Kg') + M = 1/2*(M + M') end - K = sparse(K) - solver.ndofs = size(K, 1) f = sparse(f, solver.ndofs, 1) # run any posthook for assembly if defined - args = Tuple{Solver, SparseMatrixCSC, SparseMatrixCSC} + args = Tuple{Solver, SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC} if method_exists(field_assembly_posthook!, args) - field_assembly_posthook!(solver, K, f) + field_assembly_posthook!(solver, K, Kg, f) end - return K, f + return M, K, Kg, f end """ Posthook for boundary assembly. By default, do nothing. """ @@ -188,13 +223,15 @@ function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver_UMFP t0 = time() # assemble field problems - K, f = get_field_assembly(solver) + M, K, Kg, f = get_field_assembly(solver) # assemble boundary problems Kb, C1, C2, D, fb, g = get_boundary_assembly(solver) # construct global system Ax=b and solve using lu factorization - A = [K+Kb C1'; C2 D] + A = [ + K+Kg+Kb C1' + C2 D] b = [f+fb; g] nz = get_nonzero_rows(A) @@ -214,7 +251,7 @@ function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}}) t0 = time() # assemble field problems - K, f = get_field_assembly(solver) + M, K, Kg, f = get_field_assembly(solver) # assemble boundary problems Kb, C1, C2, D, fb, g = get_boundary_assembly(solver) @@ -262,9 +299,11 @@ Notes ----- Default convergence criteria is obtained by checking each sub-problem convergence. """ -function has_converged(solver::Solver; check_convergence_for_boundary_problems=false) +function has_converged(solver::Solver{Nonlinear}; + check_convergence_for_boundary_problems=false) + properties = solver.properties converged = true - eps = solver.nonlinear_system_convergence_tolerance + eps = properties.convergence_tolerance for problem in solver.problems has_converged = true if is_field_problem(problem) @@ -286,7 +325,7 @@ function has_converged(solver::Solver; check_convergence_for_boundary_problems=f end converged &= has_converged end - return converged || solver.is_linear_system + return converged || properties.is_linear_system end type NonlinearConvergenceError <: Exception @@ -294,13 +333,14 @@ type NonlinearConvergenceError <: Exception end function Base.showerror(io::IO, exception::NonlinearConvergenceError) - max_iters = exception.solver.nonlinear_system_max_iterations + max_iters = exception.solver.properties.max_iterations print(io, "nonlinear iteration did not converge in $max_iters iterations!") end -""" Main solver loop. -""" -function call(solver::Solver) +""" Default solver for quasistatic nonlinear problems. """ +function call(solver::Solver{Nonlinear}) + + properties = solver.properties # 1. initialize each problem so that we can start nonlinear iterations for problem in solver.problems @@ -308,8 +348,8 @@ function call(solver::Solver) end # 2. start non-linear iterations - for solver.iteration=1:solver.nonlinear_system_max_iterations - info("Starting nonlinear iteration #$(solver.iteration)") + for properties.iteration=1:properties.max_iterations + info("Starting nonlinear iteration #$(properties.iteration)") # 2.1 update linearized assemblies (if needed) info("Assembling problems ...") @@ -324,8 +364,8 @@ function call(solver::Solver) # 2.2 call solver for linearized system (default: direct lu factorization) info("Solve linear system ...") tic() - u, la = solve_linear_system(solver, Val{solver.linear_system_solver}) - push!(solver.norms, (norm(u), norm(la))) + u, la = solve_linear_system(solver, Val{properties.linear_system_solver}) + push!(properties.norms, (norm(u), norm(la))) t1 = round(toq(), 2) info("Solved Ax = b in $t1 seconds.") @@ -337,8 +377,8 @@ function call(solver::Solver) # 2.4 check convergence if has_converged(solver) - info("Converged in $(solver.iteration) iterations.") - if solver.iteration < solver.nonlinear_system_min_iterations + info("Converged in $(properties.iteration) iterations.") + if properties.iteration < properties.min_iterations info("Converged but continuing") else return true @@ -347,7 +387,7 @@ function call(solver::Solver) end # 3. did not converge - if solver.nonlinear_system_error_if_no_convergence + if properties.error_if_no_convergence throw(NonlinearConvergenceError(solver)) end end diff --git a/src/sparse.jl b/src/sparse.jl index 37250d3..950d42d 100644 --- a/src/sparse.jl +++ b/src/sparse.jl @@ -10,25 +10,19 @@ type SparseMatrixCOO{T<:Real} V :: Vector{T} end -typealias SparseMatrixIJV SparseMatrixCOO - function SparseMatrixCOO() SparseMatrixCOO{Float64}([], [], []) end -#function SparseMatrixCOO{T}() -# SparseMatrixCOO{T}([], [], []) -#end - -function Base.convert(::Type{SparseMatrixCOO}, A::SparseMatrixCSC) +function convert(::Type{SparseMatrixCOO}, A::SparseMatrixCSC) return SparseMatrixCOO(findnz(A)...) end -function Base.convert(::Type{SparseMatrixCOO}, A::Matrix) +function convert(::Type{SparseMatrixCOO}, A::Matrix) return SparseMatrixCOO(findnz(A)...) end -function Base.convert(::Type{SparseMatrixCOO}, A::Vector) +function convert(::Type{SparseMatrixCOO}, A::Vector) return SparseMatrixCOO(findnz(sparse(A))...) end @@ -39,52 +33,52 @@ Parameters tol used to drop near zero values less than tol. """ -function Base.sparse(A::SparseMatrixIJV, args...; tol=1.0e-12) +function sparse(A::SparseMatrixCOO, args...; tol=1.0e-12) B = sparse(A.I, A.J, A.V, args...) SparseMatrix.droptol!(B, tol) return B end -function Base.push!(A::SparseMatrixIJV, I::Int, J::Int, V::Float64) +function push!(A::SparseMatrixCOO, I::Int, J::Int, V::Float64) push!(A.I, I) push!(A.J, J) push!(A.V, V) end -function Base.empty!(A::SparseMatrixIJV) +function empty!(A::SparseMatrixCOO) empty!(A.I) empty!(A.J) empty!(A.V) end -function Base.append!(A::SparseMatrixIJV, I::Vector{Int}, J::Vector{Int}, V::Vector{Float64}) +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 Base.append!(A::SparseMatrixIJV, B::SparseMatrixIJV) +function append!(A::SparseMatrixCOO, B::SparseMatrixCOO) append!(A.I, B.I) append!(A.J, B.J) append!(A.V, B.V) end -function Base.isempty(A::SparseMatrixIJV) +function isempty(A::SparseMatrixCOO) return isempty(A.I) && isempty(A.J) && isempty(A.V) end -function Base.(:+)(A::SparseMatrixIJV, B::SparseMatrixIJV) +function Base.(:+)(A::SparseMatrixCOO, B::SparseMatrixCOO) if isempty(A) return B end if isempty(B) return A end - C = SparseMatrixIJV([A.I;B.I], [A.J;B.J], [A.V;B.V]) + C = SparseMatrixCOO([A.I;B.I], [A.J;B.J], [A.V;B.V]) return C end -function Base.full(A::SparseMatrixCOO, args...) +function full(A::SparseMatrixCOO, args...) return full(sparse(A.I, A.J, A.V, args...)) end @@ -139,7 +133,7 @@ function add!(A::SparseMatrixCOO, dofs::Vector{Int}, data::Array{Float64}, dim:: end """ Combine (I,J,V) values is possible. """ -function optimize!(A::SparseMatrixIJV) +function optimize!(A::SparseMatrixCOO) I, J, V = findnz(sparse(A)) A = SparseMatrixCOO(I, J, V) gc() @@ -157,3 +151,15 @@ function get_nonzero_rows(A::SparseMatrixCSC) return sort(unique(rowvals(A))) end +function get_nonzero_rows(A::SparseMatrixCOO) + return get_nonzero_rows(sparse(A)) +end + +function size(A::SparseMatrixCOO) + return maximum(A.I), maximum(A.J) +end + +function size(A::SparseMatrixCOO, idx::Int) + return size(A)[idx] +end + diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl index d7f056e..5d7b8a2 100644 --- a/test/test_elasticity_2d_linear_with_surface_load.jl +++ b/test/test_elasticity_2d_linear_with_surface_load.jl @@ -52,11 +52,11 @@ using JuliaFEM.Test @test isapprox(eps, [u3; 0.0]) end - info("stress") - for ip in get_integration_points(elements[1]) - sig = ip("stress") - @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3] - @test isapprox(sig, [0.0; g; 0.0]) - end +# info("stress") +# for ip in get_integration_points(elements[1]) +# sig = ip("stress") +# @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3] +# @test isapprox(sig, [0.0; g; 0.0]) +# end end diff --git a/test/test_modal_analysis.jl b/test/test_modal_analysis.jl new file mode 100644 index 0000000..61b9176 --- /dev/null +++ b/test/test_modal_analysis.jl @@ -0,0 +1,42 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM +using JuliaFEM.Test + +@testset "test eigenvalues for single tet4 element" begin + X = Dict{Int, Vector{Float64}}( + 1 => [2.0, 3.0, 4.0], + 2 => [6.0, 3.0, 2.0], + 3 => [2.0, 5.0, 1.0], + 4 => [4.0, 3.0, 6.0]) + u = Dict{Int, Vector{Float64}}( + 1 => [0.0, 0.0, 0.0], + 2 => [0.0, 0.0, 0.0], + 3 => [0.0, 0.0, 0.0], + 4 => [0.25, 0.25, 0.25]) + e1 = Element(Tet4, [1, 2, 3, 4]) + e2 = Element(Tri3, [1, 2, 3]) + update!([e1, e2], "geometry", X) + update!([e1, e2], "displacement", u) + update!(e1, "youngs modulus" => 96.0, + "poissons ratio" => 1.0/3.0, + "density" => 420.0) + update!(e2, "displacement 1" => 0.0, + "displacement 2" => 0.0, + "displacement 3" => 0.0) + p1 = Problem(Elasticity, 3) +# p1.properties.finite_strain = true + p2 = Problem(Dirichlet, p1) + push!(p1, e1) + push!(p2, e2) + s1 = Solver(Modal) + s1.properties.which = :LM + push!(s1, p1, p2) + + call(s1; debug=true) + @test isapprox(s1.properties.eigvals, [4/3, 1/3]) + s1.properties.geometric_stiffness = true + call(s1) + @test isapprox(s1.properties.eigvals, [5/3, 2/3]) +end