diff --git a/src/assembly.jl b/src/assembly.jl index fa875f8..0d21ec7 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -20,29 +20,22 @@ function optimize!(assembly::Assembly) end function append!(assembly::Assembly, sub_assembly::Assembly) - append!(assembly.mass_matrix, sub_assembly.mass_matrix) - append!(assembly.stiffness_matrix, sub_assembly.stiffness_matrix) - append!(assembly.force_vector, sub_assembly.force_vector) -end - -function append!(assembly::BoundaryAssembly, sub_assembly::BoundaryAssembly) + append!(assembly.M, sub_assembly.M) + append!(assembly.K, sub_assembly.K) + append!(assembly.f, sub_assembly.f) append!(assembly.C1, sub_assembly.C1) append!(assembly.C2, sub_assembly.C2) append!(assembly.D, sub_assembly.D) append!(assembly.g, sub_assembly.g) end -function assemble!(problem::Union{FieldProblem, BoundaryProblem}, time::Float64; empty_assembly::Bool=true) +function assemble!(problem::Problem, time::Float64; empty_assembly::Bool=true) !problem.assembly.changed && return empty_assembly && empty!(problem.assembly) for element in get_elements(problem) assemble!(problem.assembly, problem, element, time) end problem.assembly.changed = true -end - -function assemble(problem::Union{FieldProblem, BoundaryProblem}, time::Real) - assemble!(problem, time; empty_assembly=true) return problem.assembly end diff --git a/src/directsolver.jl b/src/directsolver.jl index 090f318..0dbdd71 100644 --- a/src/directsolver.jl +++ b/src/directsolver.jl @@ -1,100 +1,7 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -## Direct solver - -using JuliaFEM - -type DirectSolver - name :: ASCIIString - field_problems :: Vector{Problem} - boundary_problems :: Vector{BoundaryProblem} - parallel :: Bool - solve_residual :: Bool - nonlinear_max_iterations :: Int64 - nonlinear_convergence_tolerance :: Float64 - linear_system_solver_preprocessors :: Vector{Tuple{Symbol,Any,Any}} - linear_system_solvers :: Vector{Tuple{Symbol,Any,Any}} - linear_system_solver_postprocessors :: Vector{Tuple{Symbol,Any,Any}} -end - -""" Default initializer. """ -function DirectSolver(name="DirectSolver") - DirectSolver( - name, - [], # field problems - [], # boundary problems - false, # parallel run? - true, # solve residual or total quantity - 10, # nonlinear problem max iterations - 5.0e-6, # nonlinear convergence tolerance - [], # default solution preprocessors - [(:UMFPACK, (), [])], # linear system solver: CHOLMOD, UMFPACK - [], # default solution postprocessors - ) -end - -function set_name!(solver::DirectSolver, name::ASCIIString) - solver.name = name -end - -function set_linear_system_solver!(solver::DirectSolver, method::Symbol) - solver.linear_system_solvers = [(method, (), [])] -end - -function set_nonlinear_max_iterations!(solver::DirectSolver, max_iterations::Int) - solver.nonlinear_max_iterations = max_iterations -end - -function push!(solver::DirectSolver, problem::FieldProblem) - push!(solver.field_problems, problem) -end - -function push!(solver::DirectSolver, problem::BoundaryProblem) - push!(solver.boundary_problems, problem) -end - -function add_linear_system_solver_preprocessor!(solver::DirectSolver, preprocessor_name::Symbol, args...; kwargs...) - push!(solver.linear_system_solver_preprocessors, (preprocessor_name, args, kwargs)) -end - -function add_linear_system_solver_postprocessor!(solver::DirectSolver, postprocessor_name::Symbol, args...; kwargs...) - push!(solver.linear_system_solver_postprocessors, (postprocessor_name, args, kwargs)) -end - -function tic(timing, what::ASCIIString) - timing[what * " start"] = time() -end - -function toc(timing, what::ASCIIString) - timing[what * " finish"] = time() -end - -function time_elapsed(timing, what::ASCIIString) - return timing[what * " finish"] - timing[what * " start"] -end - -""" -Linear system solver for problem - - Ku + C₁'λ = f - C₂u + Dλ = g - -""" -function linear_system_solver_solve!(solver, iter, time, K, f, C1, C2, D, g, sol, la, ::Type{Val{:UMFPACK}}) - t0 = Base.time() - dim = size(K, 1) - A = [K C1'; C2 D] - b = [f; g] - nz1 = sort(unique(rowvals(A))) - nz2 = sort(unique(rowvals(A'))) - u = zeros(length(b)) - u[nz1] = lufact(A[nz1,nz2]) \ full(b[nz1]) - sol[:] = u[1:dim] - la[:] = u[dim+1:end] - info("UMFPACK: solved in ", Base.time()-t0, " seconds. norm = ", norm(sol)) -end - +#= """ Solution preprocessor: dump matrices to disk before solution. Examples @@ -117,232 +24,6 @@ end function linear_system_solver_postprocess! end - -""" Initialize unknown field ready for nonlinear iterations, i.e., - take last known value and set it as a initial quess for next - time increment. -""" -function initialize!(problem::FieldProblem, time::Real) - field_name = get_unknown_field_name(problem) - field_dim = get_unknown_field_dimension(problem) - for element in get_elements(problem) - gdofs = get_gdofs(element, problem) - if haskey(element, field_name) - if !isapprox(last(element[field_name]).time, time) - last_data = copy(last(element[field_name]).data) - push!(element[field_name], time => last_data) - end - else # if field not found at all, initialize new zero field. - data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] - element[field_name] = (time => data) - end - end -end - -function initialize!(problem::BoundaryProblem, time::Real; initialize_primary_field=false) - field_name = problem.parent_field_name - field_dim = problem.parent_field_dim - - for element in get_elements(problem) - gdofs = get_gdofs(element, problem) - data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] - - # add new field "reaction force" for boundary element if not found - if haskey(element, "reaction force") - if !isapprox(last(element["reaction force"]).time, time) - push!(element["reaction force"], time => data) - end - else - element["reaction force"] = (time => data) - end - - if initialize_primary_field - # add new primary field for boundary element if not found - if haskey(element, field_name) - if !isapprox(last(element[field_name]).time, time) - last_data = copy(last(element[field_name]).data) - push!(element[field_name], time => last_data) - end - else - data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] - element[field_name] = (time => data) - end - end - end -end - -function update!(problem::FieldProblem, solution::Vector, ::Type{Val{:elements}}) - field_name = get_unknown_field_name(problem) - field_dim = get_unknown_field_dimension(problem) - for element in get_elements(problem) - gdofs = get_gdofs(element, problem) - local_sol = solution[gdofs] - local_sol = reshape(local_sol, field_dim, length(element)) - local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)] - last(element[field_name]).data = local_sol - end -end - -function update!(problem::BoundaryProblem, solution::Vector, ::Type{Val{:elements}}) - field_name = problem.parent_field_name - field_dim = problem.parent_field_dim - for element in get_elements(problem) - gdofs = get_gdofs(element, field_dim) - local_sol = solution[gdofs] - local_sol = reshape(local_sol, field_dim, length(element)) - local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)] - last(element["reaction force"]).data = local_sol - end -end +=# -""" Call solver to solve a set of problems. """ -function call(solver::DirectSolver, time::Real=0.0) - info("Starting solver $(solver.name)") - info("# of field problems: $(length(solver.field_problems))") - info("# of boundary problems: $(length(solver.boundary_problems))") - (length(solver.field_problems) != 0) || error("no field problems defined for solver, use push!(solver, problem, ...) to define field problems.") - - timing = Dict{ASCIIString, Float64}() - tic(timing, "solver") - - # check that all problems are "same kind" - field_name = get_unknown_field_name(solver.field_problems[1]) - field_dim = get_unknown_field_dimension(solver.field_problems[1]) - for field_problem in solver.field_problems - get_unknown_field_name(field_problem) == field_name || error("several different fields not supported yet") - get_unknown_field_dimension(field_problem) == field_dim || error("several different field dimensions not supported yet") - end - - tic(timing, "initialization") - for field_problem in solver.field_problems - initialize!(field_problem, time) - end - for boundary_problem in solver.boundary_problems - initialize!(boundary_problem, time) - end - toc(timing, "initialization") - - dim = nothing - sol = nothing - last_sol = nothing - la = nothing - last_la = nothing - - for iter=1:solver.nonlinear_max_iterations - info("Starting nonlinear iteration $iter") - tic(timing, "non-linear iteration") - - tic(timing, "field assembly") - info("Assembling field problems...") - field_assembly = FieldAssembly() - for (i, problem) in enumerate(solver.field_problems) - info("Assembling body $i: $(problem.name)") - append!(field_assembly, assemble(problem, time)) - end - K = sparse(field_assembly.stiffness_matrix) - dim = size(K, 1) - f = sparse(field_assembly.force_vector, dim, 1) - field_assembly = nothing - gc() - toc(timing, "field assembly") - - tic(timing, "boundary assembly") - info("Assembling boundary problems...") - boundary_assembly = BoundaryAssembly() - for (i, problem) in enumerate(solver.boundary_problems) - info("Assembling boundary $i: $(problem.name)") - append!(boundary_assembly, assemble(problem, time)) - end - - C1 = sparse(boundary_assembly.C1, dim, dim) - C2 = sparse(boundary_assembly.C2, dim, dim) - D = sparse(boundary_assembly.D, dim, dim) - g = sparse(boundary_assembly.g, dim, 1) - boundary_assembly = nothing - gc() - toc(timing, "boundary assembly") - - if iter == 1 - # initialize vectors in first iteration - sol = zeros(dim) - la = zeros(dim) - last_sol = zeros(dim) - last_la = zeros(dim) - end - - tic(timing, "preprocess solution") - # NOTE: sol and la are vectors from previous solution - for (preprocessor, args, kwargs) in solver.linear_system_solver_preprocessors - linear_system_solver_preprocess!(solver, iter, time, K, f, C1, C2, D, g, sol, la, Val{preprocessor}, args...; kwargs...) - end - toc(timing, "preprocess solution") - - gc() - tic(timing, "solution of system") - info("Solving linear system Ax=b") - for (linear_solver, args, kwargs) in solver.linear_system_solvers - last_sol = copy(sol) - last_la = copy(la) - sol = fill!(sol, 0.0) - la = fill!(la, 0.0) - linear_system_solver_solve!(solver, iter, time, K, f, C1, C2, D, g, sol, la, Val{linear_solver}, args...; kwargs...) - # if solved only difference, add to last known solution, i.e. x(i+1) = x(i) + Δx - solver.solve_residual && (sol += last_sol) - end - toc(timing, "solution of system") - gc() - - tic(timing, "postprocess solution") - for (postprocessor, args, kwargs) in solver.linear_system_solver_postprocessors - linear_system_solver_postprocess!(solver, iter, time, K, f, C1, C2, D, g, sol, la, Val{postprocessor}, args...; kwargs...) - end - toc(timing, "postprocess solution") - - tic(timing, "update element data") - for problem in solver.field_problems - update!(problem, sol, Val{:elements}) - end - for problem in solver.boundary_problems - update!(problem, la, Val{:elements}) - end - toc(timing, "update element data") - - toc(timing, "non-linear iteration") - - if false - info("timing info for iteration:") - info("boundary assembly : ", time_elapsed(timing, "boundary assembly")) - info("field assembly : ", time_elapsed(timing, "field assembly")) - info("preprocess of solution : ", time_elapsed(timing, "preprocess solution")) - info("solve linearized problem : ", time_elapsed(timing, "solution of system")) - info("update element data : ", time_elapsed(timing, "update element data")) - info("non-linear iteration : ", time_elapsed(timing, "non-linear iteration")) - end - - # check convergence - function is_converged(solver, sol, last_sol, la, last_la) - if solver.solve_residual - if norm(sol) < solver.nonlinear_convergence_tolerance - return true - end - else - if abs(norm(sol) - norm(last_sol)) < solver.nonlinear_convergence_tolerance - return true - end - end - return false - end - - if is_converged(solver, sol, last_sol, la, last_la) - toc(timing, "solver") - info("converged in $iter iterations! solver finished in ", time_elapsed(timing, "solver"), " seconds.") - return (iter, true) - end - - end - - info("Warning: did not coverge in $(solver.nonlinear_max_iterations) iterations!") - return (solver.nonlinear_max_iterations, false) - -end diff --git a/src/dirichlet.jl b/src/dirichlet.jl index e8906a0..4207e72 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -1,12 +1,7 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -abstract DirichletProblem{T} <: AbstractProblem -abstract StandardBasis -abstract DualBasis -global const BiorthogonalBasis = DualBasis - -type Dirichlet <: AbstractProblem +type Dirichlet <: BoundaryProblem dual_basis :: Bool end @@ -14,12 +9,12 @@ function Dirichlet() Dirichlet(true) end -function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{Dirichlet}, element::Element, time::Real) +function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Element, time::Real) @assert problem.properties.dual_basis # get dimension and name of PARENT field - field_dim = problem.parent_field_dim + field_dim = problem.dimension field_name = problem.parent_field_name gdofs = get_gdofs(element, field_dim) @@ -79,8 +74,9 @@ function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{Dirichle end end +#= -function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{DirichletProblem}, element::Element, time::Real) +function assemble!(assembly::Assembly, problem::Problem{DirichletProblem}, element::Element, time::Real) # get dimension and name of PARENT field field_dim = problem.parent_field_dim @@ -127,4 +123,4 @@ function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{Dirichle end end end - +=# diff --git a/src/elasticity.jl b/src/elasticity.jl index cd510b4..d037abc 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -1,6 +1,215 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +""" Concrete Elasticity type. """ +type Elasticity <: FieldProblem + # these are found from problem.properties for type Problem{Elasticity} + formulation :: Symbol + nonlinear_geometry :: Bool +end +function Elasticity() + Elasticity( + :continuum, # formulations: :plane_stress, :continuum + false, # geometrically nonlinear analysis + ) +end + +# in case of experimenting new things; +# 1. import JuliaFEM.Core: assemble! +# 2. copy/paste assemble! code to notebook +# 3. change to last argument, i.e. ::Type{Val{:plane_stress}} to ::Type{Val{:my_formulation}} +# 4. when running code: set problem.properties.formulation = :my_formulation +# 5. let multiple dispatch do the magic for you + +function get_unknown_field_name(::Type{Elasticity}) + return "displacement" +end + +function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real) + return assemble!(assembly, problem, element, time, Val{problem.properties.formulation}) +end + +""" Elasticity equations, plane stress formulation. """ +function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:plane_stress}}) + + gdofs = get_gdofs(problem, element) + ndim, nnodes = size(element) + B = zeros(3, 2*nnodes) + for ip in get_integration_points(element) + w = ip.weight + J = get_jacobian(element, ip, time) + N = element(ip, time) + if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") + nu = element("poissons ratio", ip, time) + E_ = element("youngs modulus", ip, time) + C = E_/(1.0 - nu^2) .* [ + 1.0 nu 0.0 + nu 1.0 0.0 + 0.0 0.0 (1.0-nu)/2.0] + dN = element(ip, time, Val{:grad}) + fill!(B, 0.0) + for i=1:size(dN, 2) + B[1, 2*(i-1)+1] = dN[1,i] + B[2, 2*(i-1)+2] = dN[2,i] + B[3, 2*(i-1)+1] = dN[2,i] + B[3, 2*(i-1)+2] = dN[1,i] + end + Kt = w*B'*C*B*det(J) + add!(assembly.K, gdofs, gdofs, Kt) + end + if haskey(element, "displacement load") + b = element("displacement load", ip, time) + add!(assembly.f, gdofs, w*N'*b*det(J)) + end + if haskey(element, "displacement traction force") + T = element("displacement traction force", ip, time) + L = w*T*N*norm(J) + add!(assembly.f, gdofs, vec(L)) + end + for dim in 1:get_unknown_field_dimension(problem) + if haskey(element, "displacement traction force $dim") + T = element("displacement traction force $dim", ip, time) + ldofs = gdofs[dim:problem.dim:end] + L = w*T*N*norm(J) + add!(assembly.f, ldofs, vec(L)) + end + end + if haskey(element, "displacement traction force N") + # surface pressure + p = zeros(2) + p[1] = element("displacement traction force N", ip, time) + R = element("normal-tangential coordinates", ip, time) + T = R'*p + L = w*T*N*norm(J) + add!(assembly.f, gdofs, vec(L)) + end + end +end + + +""" Elasticity equations, continuum formulation. """ +function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:continuum}}) + + gdofs = get_gdofs(problem, element) + ndim, nnodes = size(element) + B = zeros(6, 3*nnodes) + for ip in get_integration_points(element) + w = ip.weight + J = get_jacobian(element, ip, time) + N = element(ip, time) + if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") + v = element("poissons ratio", ip, time) + E_ = element("youngs modulus", ip, time) + a = 1 - v + b = 1 - 2*v + c = 1 + v + C = E_/(b*c) .* [ + a v v 0 0 0 + v a v 0 0 0 + v v a 0 0 0 + 0 0 0 b 0 0 + 0 0 0 0 b 0 + 0 0 0 0 0 b] + dN = element(ip, time, Val{:grad}) + fill!(B, 0.0) + for i=1:size(dN, 2) + B[1, 3*(i-1)+1] = dN[1,i] + B[2, 3*(i-1)+2] = dN[2,i] + B[3, 3*(i-1)+3] = dN[3,i] + B[4, 3*(i-1)+1] = dN[2,i] + B[4, 3*(i-1)+2] = dN[1,i] + B[5, 3*(i-1)+2] = dN[3,i] + B[5, 3*(i-1)+3] = dN[2,i] + B[6, 3*(i-1)+1] = dN[3,i] + B[6, 3*(i-1)+3] = dN[1,i] + end + # L = b * B' + # D = 0.5 * (L' + L) + # F = ... + # E = 0.5 * (F'*F - I) + # de = E - E_last + # S = vonMisesStress(de, stress) + # K = B' * S * J * w + Kt = w*B'*C*B*det(J) + add!(assembly.K, gdofs, gdofs, Kt) + end + if haskey(element, "displacement load") + b = element("displacement load", ip, time) + add!(assembly.f, gdofs, w*N'*b*det(J)) + end + if haskey(element, "displacement traction force") + T = element("displacement traction force", ip, time) + JT = transpose(J) + L = w*T*N*norm(cross(JT[:,1], JT[:,2])) + add!(assembly.f, gdofs, vec(L)) + end + for dim in 1:problem.dim + if haskey(element, "displacement traction force $dim") + T = element("displacement traction force $dim", ip, time) + ldofs = gdofs[dim:problem.dim:end] + JT = transpose(J) + L = w*T*N*norm(cross(JT[:,1], JT[:,2])) + add!(assembly.f, ldofs, vec(L)) + end + end + end +end + + +############################### +# Plastic material # +############################### +#= +include("vonmises.jl") +abstract PlaneStressLinearElasticPlasticProblem <: LinearElasticityProblem + +function PlaneStressLinearElasticPlasticProblem(name="plane stress linear elasticity", dim::Int=2, elements=[]) + return Problem{PlaneStressLinearElasticPlasticProblem}(name, dim, elements) +end + +""" Elasticity equations, plane stress. """ +function assemble!{E<:CG, P<:PlaneStressLinearElasticPlasticProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real) + + gdofs = get_gdofs(element, problem.dim) + ndim, nnodes = size(E) + B = zeros(3, 2*nnodes) + for ip in get_integration_points(element) + w = ip.weight + J = get_jacobian(element, ip, time) + N = element(ip, time) + if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") + nu = element("poissons ratio", ip, time) + E_ = element("youngs modulus", ip, time) + C = E_/(1.0 - nu^2) .* [ + 1.0 nu 0.0 + nu 1.0 0.0 + 0.0 0.0 (1.0-nu)/2.0] + dN = element(ip, time, Val{:grad}) + fill!(B, 0.0) + for i=1:size(dN, 2) + B[1, 2*(i-1)+1] = dN[1,i] + B[2, 2*(i-1)+2] = dN[2,i] + B[3, 2*(i-1)+1] = dN[2,i] + B[3, 2*(i-1)+2] = dN[1,i] + end + add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B*det(J)) + end + if haskey(element, "displacement load") + b = element("displacement load", ip, time) + add!(assembly.force_vector, gdofs, w*N'*b*det(J)) + end + if haskey(element, "displacement traction force") + T = element("displacement traction force", ip, time) + L = w*T*N*norm(J) + add!(assembly.force_vector, gdofs, vec(L)) + end + end +end + + + + + include("elasticplastic.jl") # Elasticity problems @@ -15,19 +224,6 @@ function get_unknown_field_type{P<:ElasticityProblem}(::Type{P}) return Vector{Float64} end -# 3D Elasticity problems -function ElasticityProblem(dim::Int=3, elements=[]) - return Problem{ElasticityProblem}("elasticity problem", dim, elements) -end - -# 2D Plane stress elasticity problems -function PlaneStressElasticityProblem(dim::Int=2, elements=[]) - return Problem{PlaneStressElasticityProblem}("plane stress elasticity problem", dim, elements) -end -function PlaneStressElasticityProblem(problem_name::ASCIIString, dim::Int=2, elements=[]) - return Problem{PlaneStressElasticityProblem}(problem_name, dim, elements) -end - """ Elasticity equations. @@ -78,7 +274,7 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element: poisson = element("poissons ratio", ip, time) mu = young/(2*(1+poisson)) lambda = young*poisson/((1+poisson)*(1-2*poisson)) - if P == PlaneStressElasticityProblem + if problem.properties.formulation == :plane_stress lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d problems end @@ -108,3 +304,5 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element: return vec(r) end + +=# diff --git a/src/equations.jl b/src/equations.jl index 2034eac..603e4e4 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -47,14 +47,12 @@ function get_gdofs(element::Element, dim::Int) return gdofs end -function get_gdofs(element::Element, problem::FieldProblem) - dim = problem.dim - return get_gdofs(element, dim) +function get_gdofs(element::Element, problem::Problem) + return get_gdofs(element, problem.dimension) end -function get_gdofs(element::Element, problem::BoundaryProblem) - dim = problem.parent_field_dim - return get_gdofs(element, dim) +function get_gdofs(problem::Problem, element::Element) + return get_gdofs(element, problem.dimension) end """ Assemble element. """ diff --git a/src/linear_elasticity.jl b/src/linear_elasticity.jl index c68eadd..96e2ef9 100644 --- a/src/linear_elasticity.jl +++ b/src/linear_elasticity.jl @@ -3,215 +3,4 @@ # Linear elasticity -""" Concrete Elasticity type. """ -type Elasticity <: AbstractProblem - plane_stress :: Bool - nonlinear_geometry :: Bool -end -function Elasticity() - Elasticity(false, false) -end -function get_unknown_field_name(::Type{Elasticity}) - return "displacement" -end - -function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real) - - # assemble plane stress problem - if problem.properties.plane_stress - return assemble!(assembly, problem, element, time, Val{:plane_stress}) - end - -end - -""" Elasticity equations, plane stress. """ -function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:plane_stress}}) - - gdofs = get_gdofs(element, problem.dim) - ndim, nnodes = size(element) - B = zeros(3, 2*nnodes) - for ip in get_integration_points(element) - w = ip.weight - J = get_jacobian(element, ip, time) - N = element(ip, time) - if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") - nu = element("poissons ratio", ip, time) - E_ = element("youngs modulus", ip, time) - C = E_/(1.0 - nu^2) .* [ - 1.0 nu 0.0 - nu 1.0 0.0 - 0.0 0.0 (1.0-nu)/2.0] - dN = element(ip, time, Val{:grad}) - fill!(B, 0.0) - for i=1:size(dN, 2) - B[1, 2*(i-1)+1] = dN[1,i] - B[2, 2*(i-1)+2] = dN[2,i] - B[3, 2*(i-1)+1] = dN[2,i] - B[3, 2*(i-1)+2] = dN[1,i] - end - Kt = w*B'*C*B*det(J) - add!(assembly.stiffness_matrix, gdofs, gdofs, Kt) - end - if haskey(element, "displacement load") - b = element("displacement load", ip, time) - add!(assembly.force_vector, gdofs, w*N'*b*det(J)) - end - if haskey(element, "displacement traction force") - T = element("displacement traction force", ip, time) - L = w*T*N*norm(J) - add!(assembly.force_vector, gdofs, vec(L)) - end - for dim in 1:problem.dim - if haskey(element, "displacement traction force $dim") - T = element("displacement traction force $dim", ip, time) - ldofs = gdofs[dim:problem.dim:end] - L = w*T*N*norm(J) - add!(assembly.force_vector, ldofs, vec(L)) - end - end - if haskey(element, "displacement traction force N") - # surface pressure - p = zeros(2) - p[1] = element("displacement traction force N", ip, time) - R = element("normal-tangential coordinates", ip, time) - T = R'*p - L = w*T*N*norm(J) - add!(assembly.force_vector, gdofs, vec(L)) - end - end -end - - - - -abstract LinearElasticityProblem <: ElasticityProblem - -function LinearElasticityProblem(name="linear elasticity", dim::Int=3, elements=[]) - return Problem{LinearElasticityProblem}(name, dim, elements) -end - -""" Elasticity equations, general 3D case. """ -function assemble!{E<:CG, P<:LinearElasticityProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real) - - gdofs = get_gdofs(element, problem.dim) - ndim, nnodes = size(E) - B = zeros(6, 3*nnodes) - for ip in get_integration_points(element) - w = ip.weight - J = get_jacobian(element, ip, time) - N = element(ip, time) - if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") - v = element("poissons ratio", ip, time) - E_ = element("youngs modulus", ip, time) - a = 1 - v - b = 1 - 2*v - c = 1 + v - C = E_/(b*c) .* [ - a v v 0 0 0 - v a v 0 0 0 - v v a 0 0 0 - 0 0 0 b 0 0 - 0 0 0 0 b 0 - 0 0 0 0 0 b] - dN = element(ip, time, Val{:grad}) - fill!(B, 0.0) - for i=1:size(dN, 2) - B[1, 3*(i-1)+1] = dN[1,i] - B[2, 3*(i-1)+2] = dN[2,i] - B[3, 3*(i-1)+3] = dN[3,i] - B[4, 3*(i-1)+1] = dN[2,i] - B[4, 3*(i-1)+2] = dN[1,i] - B[5, 3*(i-1)+2] = dN[3,i] - B[5, 3*(i-1)+3] = dN[2,i] - B[6, 3*(i-1)+1] = dN[3,i] - B[6, 3*(i-1)+3] = dN[1,i] - end - # L = b * B' - # D = 0.5 * (L' + L) - # F = ... - # E = 0.5 * (F'*F - I) - # de = E - E_last - # S = vonMisesStress(de, stress) - # K = B' * S * J * w - Kt = w*B'*C*B*det(J) - add!(assembly.stiffness_matrix, gdofs, gdofs, Kt) - # solve residual, i.e. K du = K \ -(Ku(prev) - F) - # in first iteration u(prev) typically 0 -> no effect - # but if geometrical or material nonlinearities iterations are needed -# if haskey(element, "displacement") -# u_prev = vec(element("displacement", ip, time)) -# add!(assembly.force_vector, gdofs, -Kt*u_prev) -# end - end - if haskey(element, "displacement load") - b = element("displacement load", ip, time) - add!(assembly.force_vector, gdofs, w*N'*b*det(J)) - end - if haskey(element, "displacement traction force") - T = element("displacement traction force", ip, time) - JT = transpose(J) - L = w*T*N*norm(cross(JT[:,1], JT[:,2])) - add!(assembly.force_vector, gdofs, vec(L)) - end - for dim in 1:problem.dim - if haskey(element, "displacement traction force $dim") - T = element("displacement traction force $dim", ip, time) - ldofs = gdofs[dim:problem.dim:end] - JT = transpose(J) - L = w*T*N*norm(cross(JT[:,1], JT[:,2])) - add!(assembly.force_vector, ldofs, vec(L)) - end - end - end -end - - -############################### -# Plastic material # -############################### -include("vonmises.jl") -abstract PlaneStressLinearElasticPlasticProblem <: LinearElasticityProblem - -function PlaneStressLinearElasticPlasticProblem(name="plane stress linear elasticity", dim::Int=2, elements=[]) - return Problem{PlaneStressLinearElasticPlasticProblem}(name, dim, elements) -end - -""" Elasticity equations, plane stress. """ -function assemble!{E<:CG, P<:PlaneStressLinearElasticPlasticProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real) - - gdofs = get_gdofs(element, problem.dim) - ndim, nnodes = size(E) - B = zeros(3, 2*nnodes) - for ip in get_integration_points(element) - w = ip.weight - J = get_jacobian(element, ip, time) - N = element(ip, time) - if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") - nu = element("poissons ratio", ip, time) - E_ = element("youngs modulus", ip, time) - C = E_/(1.0 - nu^2) .* [ - 1.0 nu 0.0 - nu 1.0 0.0 - 0.0 0.0 (1.0-nu)/2.0] - dN = element(ip, time, Val{:grad}) - fill!(B, 0.0) - for i=1:size(dN, 2) - B[1, 2*(i-1)+1] = dN[1,i] - B[2, 2*(i-1)+2] = dN[2,i] - B[3, 2*(i-1)+1] = dN[2,i] - B[3, 2*(i-1)+2] = dN[1,i] - end - add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B*det(J)) - end - if haskey(element, "displacement load") - b = element("displacement load", ip, time) - add!(assembly.force_vector, gdofs, w*N'*b*det(J)) - end - if haskey(element, "displacement traction force") - T = element("displacement traction force", ip, time) - L = w*T*N*norm(J) - add!(assembly.force_vector, gdofs, vec(L)) - end - end -end diff --git a/src/mortar.jl b/src/mortar.jl index 4a59d1f..c1ca1a2 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -638,31 +638,23 @@ function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix error("project_point_to_surface: did not converge in $max_iterations iterations!") end - ### Mortar problem - -""" -Parameters ----------- -node_csys - coordinate system in node, normal + tangent + "binormal" - in 3d 3x3 matrix, in 2d 2x2 matrix, respectively -""" -abstract MortarProblem{T} <: AbstractProblem - +# abstract MortarProblem{T} <: AbstractProblem # Mortar assembly 2d +type Mortar <: BoundaryProblem +end + typealias MortarElements2D Union{Seg2, Seg3} -function assemble!{E<:MortarElements2D}(assembly::BoundaryAssembly, - problem::BoundaryProblem{MortarProblem}, +function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar}, slave_element::Element{E}, time::Real) # slave element must have a set of master elements haskey(slave_element, "master elements") || return # get dimension and name of PARENT field - field_dim = problem.parent_field_dim + field_dim = problem.dimension field_name = problem.parent_field_name slave_dofs = get_gdofs(slave_element, field_dim) @@ -760,8 +752,9 @@ function find_master_elements(slave_element::Element, time::Real) return master_elements end -function assemble!{E<:MortarElements3D}(assembly::BoundaryAssembly, problem::BoundaryProblem{MortarProblem}, slave_element::Element{E}, time::Real) - field_dim = problem.parent_field_dim +function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::Problem{Mortar}, + slave_element::Element{E}, time::Real) + field_dim = problem.dimension field_name = problem.parent_field_name slave_dofs = get_gdofs(slave_element, field_dim) # info("Slave dofs: $slave_dofs") diff --git a/src/problems.jl b/src/problems.jl index 3d228c2..a1c8bb1 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -2,106 +2,51 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md abstract AbstractProblem +abstract FieldProblem <: AbstractProblem +abstract BoundaryProblem <: AbstractProblem +abstract MixedProblem <: AbstractProblem -function get_formulation_type{P<:AbstractProblem}(::Type{P}) - return :total -end +""" +General linearized problem to solve + K*u + C1'*la = f + C2*u + D*la = g +""" +type Assembly + # for field assembly + M :: SparseMatrixCOO # mass matrix + K :: SparseMatrixCOO # stiffness matrix + f :: SparseMatrixCOO # force vector + # for boundary assembly + C1 :: SparseMatrixCOO + C2 :: SparseMatrixCOO + D :: SparseMatrixCOO + g :: SparseMatrixCOO -type FieldAssembly - mass_matrix :: SparseMatrixCOO - stiffness_matrix :: SparseMatrixCOO - force_vector :: SparseMatrixCOO - solution :: Vector{Float64} - previous_solution :: Vector{Float64} - solution_norm_change :: Real - prehooks :: Vector{Tuple{Symbol,Any,Any}} - posthooks :: Vector{Tuple{Symbol,Any,Any}} + solution :: Vector{Float64} # full solution vector when solving problem Ax = b + previous_solution :: Vector{Float64} # previous solution vector + solution_norm_change :: Real # for convergence studies + prehooks :: Vector{Tuple{Symbol,Any,Any}} # assign possible prehooks before assembly + posthooks :: Vector{Tuple{Symbol,Any,Any}} # assign possible posthooks after assembly changed :: Bool # flag to control is reassembly needed end -function FieldAssembly() - return FieldAssembly( - SparseMatrixCOO(), - SparseMatrixCOO(), - SparseMatrixCOO(), - [], [], Inf, [], [], true) -end - -function Base.empty!(assembly::FieldAssembly) - empty!(assembly.mass_matrix) - empty!(assembly.stiffness_matrix) - empty!(assembly.force_vector) - assembly.changed = true -end - -typealias Assembly FieldAssembly - -""" Construct a new field problem. - -Examples --------- -Create vector-valued (dim=3) elasticity problem: - -julia> prob = FieldProblem(ElasticityProblem, "this is my problem", 3) - -""" -type FieldProblem{T} - name :: ASCIIString - dim :: Int - elements :: Vector{Element} - assembly :: FieldAssembly - properties :: T -end -function FieldProblem(problem::DataType, name::ASCIIString, dim::Int, - elements=[]) - FieldProblem{problem}(name, dim, elements, FieldAssembly(), problem()) -end - - -function update!{P}(problem::FieldProblem{P}, solution::Vector{Float64}) - # resize & fill with zeros solution vector if length mismatch with current solution - if length(solution) != length(problem.assembly.solution) - resize!(problem.assembly.solution, length(solution)) - fill!(problem.assembly.solution, 0.0) - end - problem.assembly.previous_solution = copy(problem.assembly.solution) - if get_formulation_type(P) == :incremental - problem.assembly.solution += solution - else - problem.assembly.solution = solution - end - problem.assembly.solution_norm_change = norm(problem.assembly.solution - problem.assembly.previous_solution) -end - - -""" -Interface matrices C₁, C₂ & D, g for general problem type -Au + C₁'λ = f -C₂u + Dλ = g -""" -type BoundaryAssembly - C1 :: SparseMatrixCOO - C2 :: SparseMatrixCOO - D :: SparseMatrixCOO - g :: SparseMatrixCOO - solution :: Vector{Float64} - previous_solution :: Vector{Float64} - solution_norm_change :: Real - prehooks :: Vector{Tuple{Symbol,Any,Any}} - posthooks :: Vector{Tuple{Symbol,Any,Any}} - changed :: Bool # flag to control is reassembly needed -end - -function BoundaryAssembly() - return BoundaryAssembly( +function Assembly() + return Assembly( SparseMatrixCOO(), SparseMatrixCOO(), SparseMatrixCOO(), SparseMatrixCOO(), - [], [], Inf, [], [], true) + SparseMatrixCOO(), + SparseMatrixCOO(), + SparseMatrixCOO(), + [], [], Inf, + [], [], true) end -function Base.empty!(assembly::BoundaryAssembly) +function Base.empty!(assembly::Assembly) + empty!(assembly.M) + empty!(assembly.K) + empty!(assembly.f) empty!(assembly.C1) empty!(assembly.C2) empty!(assembly.D) @@ -109,6 +54,27 @@ function Base.empty!(assembly::BoundaryAssembly) assembly.changed = true end +type Problem{P<:AbstractProblem} + name :: ASCIIString # descriptive name for problem + dimension :: Int # degrees of freedom per node + parent_field_name :: ASCIIString # (optional) name of parent field e.g. "displacement" + elements :: Vector{Element} + assembly :: Assembly + properties :: P +end + +""" Construct a new field problem. + +Examples +-------- +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=[]) + Problem{P}(name, dimension, "none", elements, Assembly(), P()) +end """ Construct a new boundary problem. @@ -116,38 +82,41 @@ Examples -------- Create Dirichlet boundary problem for vector-valued (dim=3) elasticity problem. -julia> bc1 = FieldProblem(DirichletProblem, "support", "displacement", 3) +julia> bc1 = Problem(Dirichlet, "support", 3, "displacement") """ -type BoundaryProblem{T} - name :: ASCIIString - parent_field_name :: ASCIIString - parent_field_dim :: Int - elements :: Vector{Element} - assembly :: BoundaryAssembly - properties :: T -end -function BoundaryProblem(problem::DataType, - name::ASCIIString, - parent_field_name::ASCIIString, - parent_field_dim::Int, - elements=[]) - BoundaryProblem{problem}(name, parent_field_name, parent_field_dim, - elements, BoundaryAssembly(), problem()) +function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name, elements=[]) + Problem{P}(name, dimension, parent_field_name, elements, Assembly(), P()) end -function update!{P}(problem::BoundaryProblem{P}, solution::Vector{Float64}) - if length(solution) != length(problem.assembly.solution) - resize!(problem.assembly.solution, length(solution)) - fill!(problem.assembly.solution, 0.0) +function get_formulation_type{P<:FieldProblem}(problem::Problem{P}) + return :total +end + +function get_formulation_type{P<:BoundaryProblem}(problem::Problem{P}) + return :total +end + +function get_assembly(problem) + return problem.assembly +end + +""" Update problem solution vector. +""" +function update!(problem::Problem, solution::Vector{Float64}) + assembly = get_assembly(problem) + # resize & fill with zeros solution vector if length mismatch with current solution + if length(solution) != length(assembly.solution) + resize!(assembly.solution, length(solution)) + fill!(assembly.solution, 0.0) end - problem.assembly.previous_solution = copy(problem.assembly.solution) - if get_formulation_type(P) == :incremental - problem.assembly.solution += solution + assembly.previous_solution = copy(assembly.solution) + if get_formulation_type(problem) == :incremental + assembly.solution += solution else - problem.assembly.solution = solution + assembly.solution = solution end - problem.assembly.solution_norm_change = norm(problem.assembly.solution - problem.assembly.previous_solution) + assembly.solution_norm_change = norm(assembly.solution - assembly.previous_solution) end #= @@ -161,17 +130,13 @@ end =# -typealias Problem FieldProblem - -typealias AllProblems Union{FieldProblem, BoundaryProblem} - -function get_elements(problem::AllProblems) +function get_elements(problem) return problem.elements end """ Return the dimension of the unknown field of this problem. """ function get_unknown_field_dimension(problem::Problem) - return problem.dim + return problem.dimension end """ Return the name of the unknown field of this problem. """ @@ -179,10 +144,12 @@ function get_unknown_field_name{P}(problem::Problem{P}) return get_unknown_field_name(P) end -function Base.push!(problem::AllProblems, element::Element) +function push!(problem::Problem, element) push!(problem.elements, element) end + + # TODO: better place for utility functions? """ Calculate "nodal" vector from set of elements. diff --git a/src/solvers.jl b/src/solvers.jl index 677de75..59ec8cf 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -54,19 +54,19 @@ end type LinearSolver name :: ASCIIString field_problems :: Vector{Problem} - boundary_problems :: Vector{BoundaryProblem} + boundary_problems :: Vector{Problem} end function LinearSolver(name="LinearSolver") LinearSolver(name, [], []) end -function push!(solver::LinearSolver, problem::Problem) +function push!{P<:FieldProblem}(solver::LinearSolver, problem::Problem{P}) length(solver.field_problems) == 0 || error("Only one field problem allowed for LinearSolver") push!(solver.field_problems, problem) end -function push!(solver::LinearSolver, problem::BoundaryProblem) +function push!{P<:BoundaryProblem}(solver::LinearSolver, problem::Problem{P}) length(solver.boundary_problems) == 0 || error("Only one boundary problem allowed for LinearSolver") push!(solver.boundary_problems, problem) end @@ -141,7 +141,7 @@ type Solver name :: ASCIIString time :: Real iteration :: Int - problems :: Vector{Union{FieldProblem, BoundaryProblem}} + problems :: Vector{Problem} is_linear_system :: Bool nonlinear_system_max_iterations :: Int64 nonlinear_system_convergence_tolerance :: Float64 @@ -161,27 +161,38 @@ function Solver(name::ASCIIString="default solver", time::Real=0.0) ) end -function push!(solver::Solver, problem::Union{FieldProblem, BoundaryProblem}) +function push!(solver::Solver, problem) push!(solver.problems, problem) end # one-liner helpers to identify problem types function is_field_problem(problem) - return typeof(problem) <: FieldProblem + return false +end +function is_field_problem{P<:FieldProblem}(problem::Problem{P}) + return true end function is_boundary_problem(problem) - return typeof(problem) <: BoundaryProblem + return false +end +function is_boundary_problem{P<:BoundaryProblem}(problem::Problem{P}) + return true end function is_dirichlet_problem(problem) - return typeof(problem) <: Union{BoundaryProblem{DirichletProblem}, BoundaryProblem{DirichletProblem{DualBasis}}} + return false +end +function is_dirichlet_problem{P<:Problem{Dirichlet}}(problem::P) + return true end -function is_mortar_problem(problem) - return typeof(problem) <: BoundaryProblem{MortarProblem} +#= +function is_mortar_problem{P<:Problem{Mortar}}(problem::P) + return true end +=# function get_field_problems(solver::Solver) filter(is_field_problem, solver.problems) @@ -218,12 +229,12 @@ problems must have unique node ids. function get_field_assembly(solver::Solver) return get_field_assembly(get_field_problems(solver)) end -function get_field_assembly(problems::Vector{Union{BoundaryProblem, FieldProblem}}) +function get_field_assembly(problems::Vector{Problem}) K = SparseMatrixCOO() f = SparseMatrixCOO() for problem in problems - append!(K, problem.assembly.stiffness_matrix) - append!(f, problem.assembly.force_vector) + append!(K, problem.assembly.K) + append!(f, problem.assembly.f) end return K, f end @@ -238,7 +249,7 @@ C1, C2, D, g :: SparseMatrixCOO function get_boundary_assembly(solver::Solver) return get_boundary_assembly(get_boundary_problems(solver)) end -function get_boundary_assembly(problems::Vector{Union{BoundaryProblem, FieldProblem}}) +function get_boundary_assembly(problems::Vector{Problem}) C1 = SparseMatrixCOO() C2 = SparseMatrixCOO() D = SparseMatrixCOO() @@ -283,8 +294,8 @@ function solve_linear_system!(solver::Solver, ::Type{Val{:DirectLinearSolver}}) u = x[1:dim] la = x[dim+1:end] for problem in solver.problems - typeof(problem) <: FieldProblem && update!(problem, u) - typeof(problem) <: BoundaryProblem && update!(problem, la) + is_field_problem(problem) && update!(problem, u) + is_boundary_problem(problem) && update!(problem, la) end info("UMFPACK: solved in ", time()-t0, " seconds. norm = ", norm(u)) @@ -317,6 +328,82 @@ function Base.showerror(io::IO, exception::NonlinearConvergenceError) print(io, "nonlinear iteration did not converge in $max_iters iterations!") end +""" Initialize unknown field ready for nonlinear iterations, i.e., + take last known value and set it as a initial quess for next + time increment. +""" +function initialize!{P<:FieldProblem}(problem::Problem{P}, time::Real) + field_name = get_unknown_field_name(problem) + field_dim = get_unknown_field_dimension(problem) + for element in get_elements(problem) + gdofs = get_gdofs(element, problem) + if haskey(element, field_name) + if !isapprox(last(element[field_name]).time, time) + last_data = copy(last(element[field_name]).data) + push!(element[field_name], time => last_data) + end + else # if field not found at all, initialize new zero field. + data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] + element[field_name] = (time => data) + end + end +end + +function initialize!{P<:BoundaryProblem}(problem::Problem{P}, time::Real; initialize_primary_field=false) + field_name = problem.parent_field_name + field_dim = problem.dimension + + for element in get_elements(problem) + gdofs = get_gdofs(element, problem) + data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] + + # add new field "reaction force" for boundary element if not found + if haskey(element, "reaction force") + if !isapprox(last(element["reaction force"]).time, time) + push!(element["reaction force"], time => data) + end + else + element["reaction force"] = (time => data) + end + + if initialize_primary_field + # add new primary field for boundary element if not found + if haskey(element, field_name) + if !isapprox(last(element[field_name]).time, time) + last_data = copy(last(element[field_name]).data) + push!(element[field_name], time => last_data) + end + else + data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] + element[field_name] = (time => data) + end + end + end +end + +function update!{P<:FieldProblem}(problem::Problem{P}, solution::Vector, ::Type{Val{:elements}}) + field_name = get_unknown_field_name(problem) + field_dim = get_unknown_field_dimension(problem) + for element in get_elements(problem) + gdofs = get_gdofs(element, problem) + local_sol = solution[gdofs] + local_sol = reshape(local_sol, field_dim, length(element)) + local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)] + last(element[field_name]).data = local_sol + end +end + +function update!{P<:BoundaryProblem}(problem::Problem{P}, solution::Vector, ::Type{Val{:elements}}) + field_dim = get_unknown_field_dimension(problem) + for element in get_elements(problem) + gdofs = get_gdofs(element, field_dim) + local_sol = solution[gdofs] + local_sol = reshape(local_sol, field_dim, length(element)) + local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)] + last(element["reaction force"]).data = local_sol + end +end + """ Main solver loop. """ function call(solver::Solver) diff --git a/test/test_elasticity_surface_load.jl b/test/test_elasticity_surface_load.jl index 907f495..3f27743 100644 --- a/test/test_elasticity_surface_load.jl +++ b/test/test_elasticity_surface_load.jl @@ -2,8 +2,7 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md using JuliaFEM.Test -using JuliaFEM.Core: Node, update!, Quad4, Seg2, assemble, BoundaryProblem, - Problem, Elasticity, Solver, Dirichlet +using JuliaFEM.Core: Node, update!, Quad4, Seg2, Problem, Elasticity, Solver, Dirichlet @testset "test 2d linear elasticity with surface load." begin @@ -26,8 +25,9 @@ using JuliaFEM.Core: Node, update!, Quad4, Seg2, assemble, BoundaryProblem, update!(element1, "poissons ratio", nu) # update!(element2, "displacement traction force", [0.0, f]) update!(element2, "displacement traction force", Vector{Float64}[[0.0, f], [0.0, f]]) + # type, name, dimension elasticity_problem = Problem(Elasticity, "block", 2) - elasticity_problem.properties.plane_stress = true + elasticity_problem.properties.formulation = :plane_stress push!(elasticity_problem, element1, element2) # dirichlet boundary condition, symmetry @@ -36,8 +36,8 @@ using JuliaFEM.Core: Node, update!, Quad4, Seg2, assemble, BoundaryProblem, update!([sym13, sym23], "geometry", nodes) update!(sym13, "displacement 2", 0.0) update!(sym23, "displacement 1", 0.0) - # name, unknown field, unknown field dimension - boundary_problem = BoundaryProblem(Dirichlet, "symmetry boundaries", "displacement", 2) + # type, name, dimension, unknown_field_name + boundary_problem = Problem(Dirichlet, "symmetry boundaries", 2, "displacement") push!(boundary_problem, sym13, sym23) solver = Solver("solve block problem")