From dce7472cda1a86f3deb2161bdf3f793cb77c2032 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 29 Jan 2018 07:53:14 +0200 Subject: [PATCH] Make JuliaFEM to use Analysis type from FEMBase `Analysis` is basically doing same than `Solver` before, but has a slighly simpler structure and is more general. --- src/JuliaFEM.jl | 5 +- src/io.jl | 2 +- src/problems_contact_2d.jl | 10 +- src/problems_contact_3d.jl | 11 +- src/solvers.jl | 257 ++++++++---------- src/solvers_modal.jl | 64 +++-- test/test_contact_2d_finite_sliding.jl | 3 + ..._elasticity_2d_linear_with_surface_load.jl | 5 +- ...asticity_2d_nonlinear_with_surface_load.jl | 6 +- test/test_heat_3.jl | 2 +- test/test_heat_4.jl | 18 +- test/test_modal_analysis_elasticity.jl | 1 - test/test_problems_contact_2d_autodiff.jl | 7 +- test/test_problems_contact_3d.jl | 8 +- test/test_problems_mortar_3d.jl | 16 +- test/test_solvers_postprocess.jl | 3 +- 16 files changed, 205 insertions(+), 213 deletions(-) diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 6686244..62201a9 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -8,11 +8,12 @@ This is JuliaFEM -- Finite Element Package """ module JuliaFEM -using FEMBase +using Reexport + +@reexport using FEMBase import FEMBase: get_unknown_field_name, get_unknown_field_dimension, assemble!, update!, initialize! - # from other packages TimerOutputs.jl and Logging.jl using TimerOutputs export @timeit, print_timer diff --git a/src/io.jl b/src/io.jl index c6ebc92..3424943 100644 --- a/src/io.jl +++ b/src/io.jl @@ -4,7 +4,7 @@ using HDF5 using LightXML -type Xdmf +type Xdmf <: AbstractResultsWriter name :: String xml :: XMLElement hdf :: HDF5File diff --git a/src/problems_contact_2d.jl b/src/problems_contact_2d.jl index c7b7d1c..6fddeb7 100644 --- a/src/problems_contact_2d.jl +++ b/src/problems_contact_2d.jl @@ -207,7 +207,15 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::T is_stick = Dict{Int64, Int}() la = problem.assembly.la - ndofs = length(la) + # FIXME: for matrix operations, we need to know the dimensions of the + # final matrices + ndofs = 0 + ndofs = max(ndofs, size(problem.assembly.K, 2)) + ndofs = max(ndofs, size(problem.assembly.C1, 2)) + ndofs = max(ndofs, size(problem.assembly.C2, 2)) + ndofs = max(ndofs, size(problem.assembly.D, 2)) + ndofs = max(ndofs, size(problem.assembly.g, 2)) + ndofs = max(ndofs, size(problem.assembly.c, 2)) C1 = sparse(problem.assembly.C1, ndofs, ndofs) C2 = sparse(problem.assembly.C2, ndofs, ndofs) diff --git a/src/problems_contact_3d.jl b/src/problems_contact_3d.jl index 6ce6028..d9c08a7 100644 --- a/src/problems_contact_3d.jl +++ b/src/problems_contact_3d.jl @@ -501,7 +501,16 @@ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::T is_stick = Dict{Int64, Int}() la = problem.assembly.la - ndofs = length(la) + + # FIXME: for matrix operations, we need to know the dimensions of the + # final matrices + ndofs = 0 + ndofs = max(ndofs, size(problem.assembly.K, 2)) + ndofs = max(ndofs, size(problem.assembly.C1, 2)) + ndofs = max(ndofs, size(problem.assembly.C2, 2)) + ndofs = max(ndofs, size(problem.assembly.D, 2)) + ndofs = max(ndofs, size(problem.assembly.g, 2)) + ndofs = max(ndofs, size(problem.assembly.c, 2)) C1 = sparse(problem.assembly.C1, ndofs, ndofs) C2 = sparse(problem.assembly.C2, ndofs, ndofs) diff --git a/src/solvers.jl b/src/solvers.jl index 9670337..586017d 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -1,29 +1,16 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -abstract type AbstractSolver end - -type Solver{S<:AbstractSolver} - name :: AbstractString # some descriptive name for problem - time :: Float64 # current time - problems :: Vector{Problem} - norms :: Vector{Tuple} # solution norms for convergence studies - ndofs :: Int # number of degrees of freedom in problem - xdmf :: Nullable{Xdmf} # input/output handle - initialized :: Bool - u :: Vector{Float64} - la :: Vector{Float64} - alpha :: Float64 # generalized alpha time integration coefficient - fields :: Dict{String, AbstractField} - properties :: S -end - +const Solver = Analysis +const AbstractSolver = AbstractAnalysis +#= function Solver{S<:AbstractSolver}(::Type{S}, name="solver", properties...) variant = S(properties...) solver = Solver{S}(name, 0.0, [], [], 0, nothing, false, [], [], 0.0, Dict(), variant) return solver end +=# function Solver{S<:AbstractSolver}(::Type{S}, problems::Problem...) solver = Solver(S, "$(S)Solver") @@ -31,10 +18,6 @@ function Solver{S<:AbstractSolver}(::Type{S}, problems::Problem...) return solver end -function get_problems(solver::Solver) - return solver.problems -end - function push!(solver::Solver, problem::Problem) push!(solver.problems, problem) end @@ -88,19 +71,17 @@ function get_field_assembly(solver::Solver) append!(fg, problem.assembly.fg) end - if solver.ndofs == 0 - solver.ndofs = size(K, 1) - info("automatically determined problem dimension, ndofs = $(solver.ndofs)") - end + N = size(K, 1) - M = sparse(M, solver.ndofs, solver.ndofs) - K = sparse(K, solver.ndofs, solver.ndofs) + M = sparse(M, N, N) + K = sparse(K, N, N) if nnz(K) == 0 - warn("Field assembly seems to be empty. Check that elements are pushed to problem and formulation is correct.") + warn("Field assembly seems to be empty. Check that elements are ", + "pushed to problem and formulation is correct.") end - Kg = sparse(Kg, solver.ndofs, solver.ndofs) - f = sparse(f, solver.ndofs, 1) - fg = sparse(fg, solver.ndofs, 1) + Kg = sparse(Kg, N, N) + f = sparse(f, N, 1) + fg = sparse(fg, N, 1) return M, K, Kg, f, fg end @@ -149,26 +130,24 @@ Returns K, C1, C2, D, f, g :: SparseMatrixCSC """ -function get_boundary_assembly(solver::Solver) +function get_boundary_assembly(solver::Solver, N) check_for_overconstrained_dofs(solver) - ndofs = solver.ndofs - @assert ndofs != 0 - K = spzeros(ndofs, ndofs) - C1 = spzeros(ndofs, ndofs) - C2 = spzeros(ndofs, ndofs) - D = spzeros(ndofs, ndofs) - f = spzeros(ndofs, 1) - g = spzeros(ndofs, 1) + K = spzeros(N, N) + C1 = spzeros(N, N) + C2 = spzeros(N, N) + D = spzeros(N, N) + f = spzeros(N, 1) + g = spzeros(N, 1) for problem in get_boundary_problems(solver) assembly = problem.assembly - K_ = sparse(assembly.K, ndofs, ndofs) - C1_ = sparse(assembly.C1, ndofs, ndofs) - C2_ = sparse(assembly.C2, ndofs, ndofs) - D_ = sparse(assembly.D, ndofs, ndofs) - f_ = sparse(assembly.f, ndofs, 1) - g_ = sparse(assembly.g, ndofs, 1) + K_ = sparse(assembly.K, N, N) + C1_ = sparse(assembly.C1, N, N) + C2_ = sparse(assembly.C2, N, N) + D_ = sparse(assembly.D, N, N) + f_ = sparse(assembly.f, N, 1) + g_ = sparse(assembly.g, N, 1) for dof in assembly.removed_dofs info("$(problem.name): removing dof $dof from assembly") C1_[dof,:] = 0.0 @@ -248,16 +227,17 @@ function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{2}}) A = [K C1'; C2 D] b = [f; g] + ndofs = size(K, 2) nz1 = get_nonzero_rows(A) nz2 = get_nonzero_columns(A) nz1 == nz2 || return false - x = zeros(2*solver.ndofs) + x = zeros(2*ndofs) x[nz1] = lufact(A[nz1,nz2]) \ full(b[nz1]) - u[:] = x[1:solver.ndofs] - la[:] = x[solver.ndofs+1:end] + u[:] = x[1:ndofs] + la[:] = x[ndofs+1:end] return true end @@ -272,14 +252,15 @@ function solve!(solver::Solver, K, C1, C2, D, f, g, u, la, ::Type{Val{3}}) A = [K C1'; C2 D] b = [f; g] - nz = ones(2*solver.ndofs) + ndofs = size(K, 2) + nz = ones(2*ndofs) nz[get_nonzero_rows(A)] = 0.0 A += spdiagm(nz) x = lufact(A) \ full(b) - u[:] = x[1:solver.ndofs] - la[:] = x[solver.ndofs+1:end] + u[:] = x[1:ndofs] + la[:] = x[ndofs+1:end] return true end @@ -296,7 +277,8 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric # M2, K2, Kg2, f2, fg2, C12, C22, D2, g2 = get_boundary_assembly(solver) M, K, Kg, f, fg = get_field_assembly(solver) - Kb, C1, C2, D, fb, g = get_boundary_assembly(solver) + N = size(K, 2) + Kb, C1, C2, D, fb, g = get_boundary_assembly(solver, N) K = K + Kg + Kb f = f + fg + fb @@ -312,7 +294,8 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric end gc() end - + + #= if !haskey(solver, "fint") solver.fields["fint"] = field(solver.time => f) else @@ -329,8 +312,9 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric C1 = (1-alpha)*C1 f = (1-alpha)*f + alpha*fint.data[end-1].second end + =# - ndofs = solver.ndofs + ndofs = N u = zeros(ndofs) la = zeros(ndofs) is_solved = false @@ -346,15 +330,15 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric end t1 = round(Base.time()-t0, 2) norms = (norm(u), norm(la)) - push!(solver.norms, norms) + #push!(solver.norms, norms) - solver.u = u - solver.la = la + #solver.u = u + #solver.la = la info("Solved problems in $t1 seconds using solver $i.") info("Solution norms = $norms.") - return + return u, la end """ @@ -367,24 +351,25 @@ standard assembler for them. As a result, each problem.assembly is populated with global stiffness matrix, force vector, and, optionally, mass matrix. """ -function assemble!(solver::Solver; with_mass_matrix=false) +function assemble!(solver::Solver, time::Float64; with_mass_matrix=false) info("Assembling problems ...") for problem in get_problems(solver) timeit("assemble $(problem.name)") do empty!(problem.assembly) - assemble!(problem, solver.time) + assemble!(problem, time) end end if with_mass_matrix for problem in get_field_problems(solver) timeit("assemble $(problem.name) mass matrix") do - assemble!(problem, solver.time, Val{:mass_matrix}) + assemble!(problem, time, Val{:mass_matrix}) end end end + #= ndofs = 0 for problem in solver.problems Ks = size(problem.assembly.K, 2) @@ -392,7 +377,7 @@ function assemble!(solver::Solver; with_mass_matrix=false) ndofs = max(ndofs, Ks, Cs) end solver.ndofs = ndofs - + =# info("Assembly done!") end @@ -491,9 +476,9 @@ function (solver::Solver)(field_name::String, time::Float64) end """ Default update for solver. """ -function update!{S}(solver::Solver{S}) - u = solver.u - la = solver.la +function update!{S}(solver::Solver{S}, u, la, time) + #u = solver.u + #la = solver.la info("Updating problems ...") t0 = Base.time() @@ -504,7 +489,7 @@ function update!{S}(solver::Solver{S}) # update solution, first for assembly (u,la) ... update!(problem, assembly, u, la) # .. and then from assembly (u,la) to elements - update!(problem, assembly, elements, solver.time) + update!(problem, assembly, elements, time) end t1 = round(Base.time()-t0, 2) @@ -515,35 +500,42 @@ end functions to calculate secondary fields, i.e. contact pressure, stress, heat flux, reaction force etc. quantities. """ -function postprocess!(solver::Solver) +function postprocess!(solver::Solver, time) info("Running postprocess scripts for solver...") for problem in get_problems(solver) for field_name in problem.postprocess_fields field = Val{Symbol(field_name)} info("Running postprocess for problem $(problem.name), field $field_name") - postprocess!(problem, solver.time, field) + postprocess!(problem, time, field) end end end -""" Default xdmf update for solver. Loop all problems and write them individually +""" + write_results!(solver, time) + +Default xdmf update for solver. Loop all problems and write them individually to Xdmf file. By default write the main unknown field (displacement, temperature, ...) and any fields requested separately in `problem.postprocess_fields` vector (stress, strain, ...) """ -function update_xdmf!(solver::Solver) - if isnull(solver.xdmf) - info("update_xdmf: xdmf not attached to solver, not writing output to file.") - info("turn Xdmf writing on to solver by typing: solver.xdmf = Xdmf(\"results\")") +function write_results!(solver, time) + results_writers = get_results_writers(solver) + if length(results_writers) == 0 + info("Xdmf is not attached to solver, not writing output to a file.") + info("To write results to Xdmf file, attach Xdmf to Solver, i.e.") + info("add_results_writer!(solver, Xdmf(\"results\"))") return end - xdmf = get(solver.xdmf) - for problem in get_problems(solver) - fields = [get_unknown_field_name(problem); problem.postprocess_fields] - if is_boundary_problem(problem) - fields = [fields; get_parent_field_name(problem)] + # FIXME: result writer can be anything, not only Xdmf + for xdmf in results_writers + for problem in get_problems(solver) + fields = [get_unknown_field_name(problem); problem.postprocess_fields] + if is_boundary_problem(problem) + fields = [fields; get_parent_field_name(problem)] + end + update_xdmf!(xdmf, problem, time, fields) end - update_xdmf!(xdmf, problem, solver.time, fields) end end @@ -588,36 +580,42 @@ function has_converged(solver::Solver{Nonlinear}) end """ Default solver for quasistatic nonlinear problems. """ -function (solver::Solver{Nonlinear})() +function solve!(solver::Solver{Nonlinear}, time::Float64) + problems = get_problems(solver) properties = solver.properties # 1. initialize each problem so that we can start nonlinear iterations - initialize!(solver) + for problem in problems + initialize!(problem, time) + end # 2. start non-linear iterations for properties.iteration=1:properties.max_iterations info(repeat("-", 80)) info("Starting nonlinear iteration #$(properties.iteration)") - info("Increment time t=$(round(solver.time, 3))") + info("Increment time t=$(round(time, 3))") info(repeat("-", 80)) # 2.1 update assemblies - assemble!(solver) + for problem in problems + empty!(problem.assembly) + assemble!(problem, time) + end # 2.2 call solver for linearized system - solve!(solver) + u, la = solve!(solver) # 2.3 update solution back to elements - update!(solver) + update!(solver, u, la, time) # 2.4 check convergence if properties.iteration >= properties.min_iterations && has_converged(solver) info("Converged in $(properties.iteration) iterations.") # 2.4.1 run any postprocessing of problems - postprocess!(solver) + postprocess!(solver, time) # 2.4.2 update Xdmf output - update_xdmf!(solver) + write_results!(solver, time) return true end end @@ -628,21 +626,6 @@ function (solver::Solver{Nonlinear})() end end -""" Convenience function to call nonlinear solver. """ -function NonlinearSolver(problems...) - solver = Solver(Nonlinear, "default nonlinear solver") - if length(problems) != 0 - push!(solver, problems...) - end - return solver -end -function NonlinearSolver(name::AbstractString, problems::Problem...) - solver = NonlinearSolver(problems...) - solver.name = name - return solver -end - - ### Linear quasistatic solver """ Quasistatic solver for linear problems. @@ -657,51 +640,41 @@ Main differences in this solver, compared to nonlinear solver are: type Linear <: AbstractSolver end -function assemble!(solver::Solver{Linear}) - info("Assembling problems ...") - tic() - nproblems = 0 - ndofs = 0 - for problem in get_problems(solver) - if isempty(problem.assembly) - assemble!(problem, solver.time) - nproblems += 1 - else - info("$(problem.name) already assembled, skipping.") - end - ndofs = max(ndofs, size(problem.assembly.K, 2)) +function solve!(solver::Solver{Linear}, time::Float64) + problems = get_problems(solver) + N = 0 + @timeit "assemble problems" for problem in problems + isempty(problem.assembly) || continue + initialize!(problem, time) + assemble!(problem, time) end - solver.ndofs = ndofs - t1 = round(toq(), 2) - info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.") + @timeit "solve linear system" u, la = solve!(solver) + @timeit "update problems" update!(solver, u, la, time) end -function (solver::Solver{Linear})() - t0 = Base.time() - info(repeat("-", 80)) - info("Starting linear solver") - info("Increment time t=$(round(solver.time, 3))") - info(repeat("-", 80)) - @timeit "initialize solver" initialize!(solver) - @timeit "assemble problems" assemble!(solver) - @timeit "solve linear system" solve!(solver) - @timeit "update problems" update!(solver) - t1 = round(Base.time()-t0, 2) - info("Linear solver ready in $t1 seconds.") +# Convenience functions + +function LinearSolver(name::String="Linear solver") + return Solver(Linear, name) end -""" Convenience function to call linear solver. """ function LinearSolver(problems::Problem...) - solver = Solver(Linear, "default linear solver") - if length(problems) != 0 - push!(solver, problems...) - end - return solver -end -function LinearSolver(name::AbstractString, problems::Problem...) - solver = LinearSolver(problems...) - solver.name = name + solver = LinearSolver() + add_problems!(solver, collect(problems)) return solver end -### End of linear quasistatic solver +function NonlinearSolver(name::String="Nonlinear solver") + return Solver(Nonlinear, name) +end + +function NonlinearSolver(problems::Problem...) + solver = NonlinearSolver() + add_problems!(solver, collect(problems)) + return solver +end + +# will be deprecated +function (solver::Solver)(time::Float64=0.0) + solve!(solver, time) +end diff --git a/src/solvers_modal.jl b/src/solvers_modal.jl index 3b7b033..450274e 100644 --- a/src/solvers_modal.jl +++ b/src/solvers_modal.jl @@ -17,10 +17,18 @@ type Modal <: AbstractSolver eigvecs :: Matrix nev :: Int which :: Symbol + bc_invertible :: Bool + P :: Vector{SparseMatrixCSC} + symmetric :: Bool + empty_assemblies_before_solution :: Bool + dense :: Bool + info_matrices :: Bool + sigma :: Float64 end function Modal(nev=10, which=:SM) - solver = Modal(false, Vector(), Matrix(0,0), nev, which) + solver = Modal(false, [], Matrix{Float64}(0,0), nev, which, + false, [], true, true, false, false, 0.0) end """ Eliminate Dirichlet boundary condition from matrices K, M. """ @@ -140,31 +148,24 @@ function eliminate_boundary_conditions!(K_red::SparseMatrixCSC, return true end -""" -Parameters ----------- - -sigma - Shift stiffness matrix by adding diagonal term, i.e. K_shifted = K + sigma*I -""" -function (solver::Solver{Modal})(; bc_invertible=false, P=nothing, symmetric=true, - empty_assemblies_before_solution=true, dense=false, - info_matrices=false, sigma=0.0) +function solve!(solver::Solver{Modal}, time::Float64) + problems = get_problems(solver) + properties = solver.properties info(repeat("-", 80)) info("Starting natural frequency solver") - info("Increment time t=$(round(solver.time, 3))") + info("Increment time t=$(round(time, 3))") info(repeat("-", 80)) - initialize!(solver) @timeit "assemble matrices" begin - assemble!(solver; with_mass_matrix=true) + assemble!(solver, time; with_mass_matrix=true) M, K, Kg, f = get_field_assembly(solver) - if solver.properties.geometric_stiffness + if properties.geometric_stiffness K += Kg end end dim = size(K, 1) + ndofs = size(K, 1) nboundary_problems = length(get_boundary_problems(solver)) @@ -172,10 +173,12 @@ function (solver::Solver{Modal})(; bc_invertible=false, P=nothing, symmetric=tru M_red = M @timeit "eliminate boundary conditions" begin - if !(P == nothing) + if length(properties.P) > 0 info("Using custom P to make transform K_red = P'*K*P and M_red = P'*M*P") - K_red = P'*K_red*P - M_red = P'*M_red*P + for P in properties.P + K_red = P'*K_red*P + M_red = P'*M_red*P + end elseif nboundary_problems != 0 info("Eliminate boundary conditions from system.") for boundary_problem in get_boundary_problems(solver) @@ -187,7 +190,7 @@ function (solver::Solver{Modal})(; bc_invertible=false, P=nothing, symmetric=tru end # free up some memory before solution - if empty_assemblies_before_solution + if properties.empty_assemblies_before_solution for problem in get_field_problems(solver) empty!(problem.assembly) end @@ -200,30 +203,29 @@ function (solver::Solver{Modal})(; bc_invertible=false, P=nothing, symmetric=tru K_red = K_red[nz,nz] M_red = M_red[nz,nz] - if sigma != 0.0 + if properties.sigma != 0.0 info("Adding diagonal term $sigma to stiffness matrix") end - ndofs = solver.ndofs props = solver.properties info("Calculate $(props.nev) eigenvalues...") tic() - if symmetric + if properties.symmetric K_red = 1/2*(K_red + transpose(K_red)) M_red = 1/2*(M_red + transpose(M_red)) end - if info_matrices + if properties.info_matrices 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)) end - if dense + if properties.dense K_red = full(K_red) M_red = full(M_red) end @@ -253,7 +255,7 @@ function (solver::Solver{Modal})(; bc_invertible=false, P=nothing, symmetric=tru om2 = eigvals(full(K_red)) info("squared eigenvalues om2 = $om2") end - if sigma != 0.0 + if properties.sigma != 0.0 info("sigma is manually set and did not work, giving up, try increase sigma.") rethrow() end @@ -270,7 +272,7 @@ function (solver::Solver{Modal})(; bc_invertible=false, P=nothing, symmetric=tru rethrow() end end - + t1 = round(toq(), 2) info("Eigenvalues computed in $t1 seconds. Squared eigenvalues: $om2") @@ -297,17 +299,19 @@ end function update_xdmf!(solver::Solver{Modal}) - if isnull(solver.xdmf) - info("update_xdmf: xdmf not attached to solver, not writing file output.") + results_writers = get_results_writers(solver) + if length(results_writers) == 0 + info("Xdmf is not attached to solver, not writing output to a file.") + info("To write results to Xdmf file, attach Xdmf to Solver, i.e.") + info("add_results_writer!(solver, Xdmf(\"results\"))") return end - if maximum(abs.(imag(solver.properties.eigvals))) > 1.0e-9 info("Writing imaginary eigenvalues for Xdmf not supported.") return end - xdmf = get(solver.xdmf) + xdmf = first(results_writers) @timeit "fetch geometry" X_ = solver("geometry", solver.time) node_ids = keys(X_) diff --git a/test/test_contact_2d_finite_sliding.jl b/test/test_contact_2d_finite_sliding.jl index 4155a5d..9a31a18 100644 --- a/test/test_contact_2d_finite_sliding.jl +++ b/test/test_contact_2d_finite_sliding.jl @@ -45,6 +45,9 @@ using JuliaFEM.Testing contact_master_elements = create_elements(mesh, "UPPER_BOTTOM") update!(contact_slave_elements, "master elements", contact_master_elements) contact.elements = [contact_master_elements; contact_slave_elements] + nnodes = length(mesh.nodes) + contact.assembly.u = zeros(2*nnodes) + contact.assembly.la = zeros(2*nnodes) solver = NonlinearSolver(upper, lower, bc_upper, bc_lower, contact) solver() diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl index 6cd4cce..a61623b 100644 --- a/test/test_elasticity_2d_linear_with_surface_load.jl +++ b/test/test_elasticity_2d_linear_with_surface_load.jl @@ -38,8 +38,9 @@ using JuliaFEM update!(bc.elements[1], "displacement 2", 0.0) update!(bc.elements[2], "displacement 1", 0.0) - solver = Solver(Linear, block, traction, bc) - solver() + solver = Solver(Linear, "solve 2d linear elasticity problem") + add_problems!(solver, [block, traction, bc]) + solve!(solver, 0.0) f = 288.0 g = 576.0 diff --git a/test/test_elasticity_2d_nonlinear_with_surface_load.jl b/test/test_elasticity_2d_nonlinear_with_surface_load.jl index 4e51d1f..4a99c5f 100644 --- a/test/test_elasticity_2d_nonlinear_with_surface_load.jl +++ b/test/test_elasticity_2d_nonlinear_with_surface_load.jl @@ -32,9 +32,9 @@ using JuliaFEM.Testing update!(bc_elements_bottom, "displacement 2", 0.0) push!(bc_sym, bc_elements_left..., bc_elements_bottom...) - solver = NonlinearSolver("solve block problem") - push!(solver, block, bc_sym) - solver() + solver = Solver(Nonlinear, "solve block problem") + add_problems!(solver, [block, bc_sym]) + solve!(solver, 0.0) # from code aster u3_expected = [-4.92316106779943E-01, 7.96321884292103E-01] diff --git a/test/test_heat_3.jl b/test/test_heat_3.jl index d247d79..8c332c4 100644 --- a/test/test_heat_3.jl +++ b/test/test_heat_3.jl @@ -46,7 +46,7 @@ using JuliaFEM.Testing gradT(X) = [2*X[1] 4*X[2]] X = [0.5, 0.5] gradT1 = gradT(X) - gradT2 = field("temperature", X, solver.time, Val{:Grad}) + gradT2 = field("temperature", X, 0.0, Val{:Grad}) info("gradT1 = $gradT1, gradT2 = $gradT2") # [1.1666666666666625 1.5000000000000018] quite big difference ..? @test isapprox(gradT1, gradT2; rtol=25.0e-2) diff --git a/test/test_heat_4.jl b/test/test_heat_4.jl index 2267abf..5bc5c84 100644 --- a/test/test_heat_4.jl +++ b/test/test_heat_4.jl @@ -26,20 +26,16 @@ using DataFrames push!(bc, boundary_element) solver = Solver(Linear, problem, bc) - solver.time = 0.0 - empty!(problem.assembly) - solver() + solve!(solver, 0.0) @test isapprox(solver("temperature", 0.0)[3], 1.0) empty!(problem.assembly) solver() @test isapprox(solver("temperature", 0.0)[3], 1.0) - solver.time = 1.0 - empty!(problem.assembly) - solver() + solve!(solver, 1.0) @test isapprox(solver("temperature", 1.0)[3], 2.0) empty!(problem.assembly) @@ -69,24 +65,20 @@ end push!(bc, boundary_element) solver = Solver(Nonlinear, problem, bc) - solver.time = 0.0 - empty!(problem.assembly) - solver() + solve!(solver, 0.0) @test isapprox(solver("temperature", 0.0)[3], 1.0) empty!(problem.assembly) solver() @test isapprox(solver("temperature", 0.0)[3], 1.0) - solver.time = 1.0 - empty!(problem.assembly) - solver() + solve!(solver, 1.0) @test isapprox(solver("temperature", 1.0)[3], 2.0) empty!(problem.assembly) - solver() + solve!(solver, 1.0) @test isapprox(solver("temperature", 1.0)[3], 2.0) end diff --git a/test/test_modal_analysis_elasticity.jl b/test/test_modal_analysis_elasticity.jl index a9e071f..4c6a52c 100644 --- a/test/test_modal_analysis_elasticity.jl +++ b/test/test_modal_analysis_elasticity.jl @@ -144,7 +144,6 @@ numéro fréquence (HZ) norme d'erreur solver = Solver(Modal, body, fixed1, fixed2) solver.properties.nev = 5 - solver.xdmf = Xdmf() solver() freqs_jf = sqrt.(solver.properties.eigvals)/(2.0*pi) # with Tet4 elements diff --git a/test/test_problems_contact_2d_autodiff.jl b/test/test_problems_contact_2d_autodiff.jl index f90ba2b..599032d 100644 --- a/test/test_problems_contact_2d_autodiff.jl +++ b/test/test_problems_contact_2d_autodiff.jl @@ -11,7 +11,7 @@ datadir = first(splitext(basename(@__FILE__))) function get_model() meshfile = joinpath(datadir, "block_2d.med") mesh = aster_read_mesh(meshfile) - println(mesh.nodes[1]) + #error("mesh has $(length(mesh.nodes)) nodes") upper = Problem(mesh, Elasticity, "UPPER", 2) lower = Problem(mesh, Elasticity, "LOWER", 2) @@ -59,6 +59,8 @@ end solver = get_model() interface = solver["interface"] + interface.assembly.u = zeros(48) + interface.assembly.la = zeros(48) upper = solver["UPPER"] lower = solver["LOWER"] for body in [upper, lower] @@ -70,8 +72,7 @@ end for time in [0.0, 1/3, 2/3, 1.0] interface.properties.iteration = 1 - solver.time = time - solver() + solve!(solver, time) end node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 1.0) diff --git a/test/test_problems_contact_3d.jl b/test/test_problems_contact_3d.jl index 6388182..0adf723 100644 --- a/test/test_problems_contact_3d.jl +++ b/test/test_problems_contact_3d.jl @@ -69,7 +69,7 @@ end @testset "small sliding contact patch test, tet4 + standard basis" begin solver = get_model(tet4_meshfile) - solver.xdmf = Xdmf("contact_sl_lin_disp_results"; overwrite=true) + add_results_writer!(solver, Xdmf("contact_sl_lin_disp_results"; overwrite=true)) interface = solver["LOWER_TO_UPPER"] interface.properties.dual_basis = false solver() @@ -91,7 +91,7 @@ end @testset "small sliding contact patch test, tet4 + dual basis" begin solver = get_model(tet4_meshfile) - solver.xdmf = Xdmf("contact_dl_lin_disp_results"; overwrite=true) + add_results_writer!(solver, Xdmf("contact_dl_lin_disp_results"; overwrite=true)) interface = solver["LOWER_TO_UPPER"] interface.properties.dual_basis = true solver() @@ -106,7 +106,7 @@ end @testset "small sliding contact patch test, tet10 + standard basis" begin solver = get_model(tet10_meshfile) - solver.xdmf = Xdmf("contact_sl_quad_disp_results"; overwrite=true) + add_results_writer!(solver, Xdmf("contact_sl_quad_disp_results"; overwrite=true)) interface = solver["LOWER_TO_UPPER"] interface.properties.dual_basis = false solver() @@ -121,7 +121,7 @@ end @testset "small sliding contact patch test, tet10 + dual basis, alpha=0.2" begin solver = get_model(tet10_meshfile) - solver.xdmf = Xdmf("contact_dl_quad_disp_results"; overwrite=true) + add_results_writer!(solver, Xdmf("contact_dl_quad_disp_results"; overwrite=true)) interface = solver["LOWER_TO_UPPER"] interface.properties.dual_basis = true interface.properties.alpha = 0.2 diff --git a/test/test_problems_mortar_3d.jl b/test/test_problems_mortar_3d.jl index 3461a1c..73caaa2 100644 --- a/test/test_problems_mortar_3d.jl +++ b/test/test_problems_mortar_3d.jl @@ -38,7 +38,7 @@ tet10_meshfile = "test_problems_mortar_3d/tet10.inp" JuliaFEM.diagnose_interface(interface, 0.0) solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface) - solver.xdmf = Xdmf("sl_lin_temp_results") + add_results_writer!(solver, Xdmf("sl_lin_temp_results"; overwrite=true)) solver() @@ -92,7 +92,7 @@ end JuliaFEM.diagnose_interface(interface, 0.0) solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface) - solver.xdmf = Xdmf("dl_lin_temp_results") + add_results_writer!(solver, Xdmf("dl_lin_temp_results"; overwrite=true)) solver() @@ -138,7 +138,7 @@ end # JuliaFEM.diagnose_interface(interface, 0.0) solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface) - solver.xdmf = Xdmf("sl_quad_temp_results") + add_results_writer!(solver, Xdmf("sl_quad_temp_results"; overwrite=true)) solver() node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0) @@ -185,7 +185,7 @@ end interface.properties.alpha = 0.2 solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface) - solver.xdmf = Xdmf("dl_quad_temp_results") + add_results_writer!(solver, Xdmf("dl_quad_temp_results"; overwrite=true)) solver() node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0) @@ -257,7 +257,7 @@ end interface.properties.dual_basis = false solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface) - solver.xdmf = Xdmf("sl_lin_disp_results") + add_results_writer!(solver, Xdmf("sl_lin_disp_results"; overwrite=true)) solver() node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) @@ -319,7 +319,7 @@ end interface.properties.dual_basis = true solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface) - solver.xdmf = Xdmf("dl_lin_disp_results") + add_results_writer!(solver, Xdmf("dl_lin_disp_results"; overwrite=true)) solver() node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) @@ -382,7 +382,7 @@ end interface.properties.dual_basis = false solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface) - solver.xdmf = Xdmf("sl_quad_disp_results") + add_results_writer!(solver, Xdmf("sl_quad_disp_results"; overwrite=true)) solver() node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) @@ -447,7 +447,7 @@ end interface.properties.alpha = 0.2 solver = LinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface) - solver.xdmf = Xdmf("dl_quad_disp_results") + add_results_writer!(solver, Xdmf("dl_quad_disp_results"; overwrite=true)) solver() node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0) diff --git a/test/test_solvers_postprocess.jl b/test/test_solvers_postprocess.jl index ef204cc..a732680 100644 --- a/test/test_solvers_postprocess.jl +++ b/test/test_solvers_postprocess.jl @@ -68,7 +68,8 @@ datadir = first(splitext(basename(@__FILE__))) solver = NonlinearSolver(upper, lower, bc_upper, bc_lower, bc_sym13, bc_sym23, interface) - solver.xdmf = Xdmf("contact_two_blocks_postprocess"; overwrite=true) + xdmf = Xdmf("contact_two_blocks_postprocess"; overwrite=true) + add_results_writer!(solver, xdmf) solver() node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 0.0)