Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-18 17:47:29 +00:00
Code Issues Packages Projects Releases Wiki Activity

Make JuliaFEM to use Analysis type from FEMBase

Browse Source 
`Analysis` is basically doing same than `Solver` before, but has a
slighly simpler structure and is more general.
This commit is contained in:
Jukka Aho
2018-01-29 07:53:14 +02:00
parent d7bef419ed
commit dce7472cda
16 changed files with 205 additions and 213 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+3 -2
View File
@@ -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
src/io.jl
+1 -1
View File
@@ -4,7 +4,7 @@
using HDF5
using LightXML
type Xdmf
type Xdmf <: AbstractResultsWriter
name :: String
xml :: XMLElement
hdf :: HDF5File
src/problems_contact_2d.jl
+9 -1
View File
@@ -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)
src/problems_contact_3d.jl
+10 -1
View File
@@ -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)
src/solvers.jl
+115 -142
View File
@@ -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
src/solvers_modal.jl
+34 -30
View File
@@ -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_)
test/test_contact_2d_finite_sliding.jl
+3
View File
@@ -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()
test/test_elasticity_2d_linear_with_surface_load.jl
+3 -2
View File
@@ -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
test/test_elasticity_2d_nonlinear_with_surface_load.jl
+3 -3
View File
@@ -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]
test/test_heat_3.jl
+1 -1
View File
@@ -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)
test/test_heat_4.jl
+5 -13
View File
@@ -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
test/test_modal_analysis_elasticity.jl
-1
View File
@@ -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
test/test_problems_contact_2d_autodiff.jl
+4 -3
View File
@@ -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)
test/test_problems_contact_3d.jl
+4 -4
View File
@@ -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
test/test_problems_mortar_3d.jl
+8 -8
View File
@@ -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)
test/test_solvers_postprocess.jl
+2 -1
View File
@@ -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)