From 42254e73adfbf7ff359537f6ae6a829acbac7284 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 1 Jan 2016 20:34:10 +0200 Subject: [PATCH] its now possible to define several pre and postprocessors for assemblies with custom arguments; see test/test_mortar_2d_contact.jl for concrete example --- src/assembly.jl | 14 ++--- src/directsolver.jl | 8 ++- src/problems.jl | 21 +++++-- test/test_mortar_2d_contact.jl | 102 +++++++++++++++++++++++++++------ 4 files changed, 112 insertions(+), 33 deletions(-) diff --git a/src/assembly.jl b/src/assembly.jl index 18bea4d..236b64c 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -64,11 +64,11 @@ function assemble(problem::AllProblems, elrange::UnitRange{Int64}, time::Real, o end """ Run preprocess for assembly. """ -function preprocess_assembly! +function assemble_preprocess! end """ Run postprocess for assembly. """ -function postprocess_assembly! +function assemble_postprocess! end function assemble(problem::AllProblems, time::Real, nchunks=10) @@ -77,9 +77,8 @@ function assemble(problem::AllProblems, time::Real, nchunks=10) slices = [kk[j]+1:kk[j+1] for j=1:nchunks] assembly = new_assembly(typeof(problem)) - args = Tuple{typeof(assembly), typeof(problem), Real} - if method_exists(preprocess_assembly!, args) - preprocess_assembly!(assembly, problem, time) + for (preprocessor, args, kwargs) in problem.preprocessors + assemble_preprocess!(assembly, problem, time, Val{preprocessor}, args...; kwargs...) end for (j, elrange) in enumerate(slices) @@ -90,9 +89,8 @@ function assemble(problem::AllProblems, time::Real, nchunks=10) end end - args = Tuple{typeof(assembly), typeof(problem), Real} - if method_exists(postprocess_assembly!, args) - postprocess_assembly!(assembly, problem, time) + for (postprocessor, args, kwargs) in problem.postprocessors + assemble_postprocess!(assembly, problem, time, Val{postprocessor}, args...; kwargs...) end return assembly diff --git a/src/directsolver.jl b/src/directsolver.jl index f001b58..ba9be5f 100644 --- a/src/directsolver.jl +++ b/src/directsolver.jl @@ -34,6 +34,10 @@ function DirectSolver(name="DirectSolver") ) 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 = Vector{Symbol}([method]) end @@ -262,7 +266,7 @@ function call(solver::DirectSolver, time::Real=0.0) last(element["reaction force"]).data = local_sol # <-- replaced # FIXME: Quick and dirty, updating dirichlet boundary problem # is causing drifting and convergence issue - if typeof(boundary_problem) <: BoundaryProblem{DirichletProblem{StandardBasis}} + if typeof(boundary_problem) <: BoundaryProblem{DirichletProblem} # info("skipping dirichlet problem update") continue end @@ -288,7 +292,7 @@ function call(solver::DirectSolver, time::Real=0.0) if (norm(sol) < solver.nonlinear_convergence_tolerance) toc(timing, "solver") - info("converged! solver finished in ", time_elapsed(timing, "solver"), " seconds.") + info("converged in $iter iterations! solver finished in ", time_elapsed(timing, "solver"), " seconds.") return (iter, true) end diff --git a/src/problems.jl b/src/problems.jl index aae9925..5cd5945 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -7,6 +7,8 @@ type FieldProblem{T<:AbstractProblem} name :: ASCIIString dim :: Int elements :: Vector{Element} + preprocessors :: Vector{Tuple{Symbol,Any,Any}} + postprocessors :: Vector{Tuple{Symbol,Any,Any}} end type BoundaryProblem{T<:AbstractProblem} @@ -15,6 +17,8 @@ type BoundaryProblem{T<:AbstractProblem} parent_field_dim :: Int dim :: Int elements :: Vector{Element} + preprocessors :: Vector{Tuple{Symbol,Any,Any}} + postprocessors :: Vector{Tuple{Symbol,Any,Any}} end """ Construct new field problem. @@ -26,8 +30,8 @@ Create vector-valued (dim=3) elasticity problem: julia> prob = FieldProblem(ElasticityProblem, "this is my problem", 3) """ -function FieldProblem(problem_type::DataType, name::ASCIIString, dim::Int, elements=[]) - FieldProblem{problem_type}(name, dim, elements) +function FieldProblem(problem_type::DataType, name::ASCIIString, dim::Int, elements=[], preprocessors=[], postprocessors=[]) + FieldProblem{problem_type}(name, dim, elements, preprocessors, postprocessors) end @@ -40,8 +44,17 @@ Create Dirichlet boundary problem for vector-valued (dim=3) elasticity problem. julia> bc1 = FieldProblem(DirichletProblem, "support dy=0", "displacement", 3) """ -function BoundaryProblem(problem_type::DataType, name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=[]) - BoundaryProblem{problem_type}(name, parent_field_name, parent_field_dim, dim, elements) +function BoundaryProblem(problem_type::DataType, name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=[], preprocessors=[], postprocessors=[]) + BoundaryProblem{problem_type}(name, parent_field_name, parent_field_dim, dim, elements, preprocessors, postprocessors) +end + + +function add_postprocessor!(problem::Union{FieldProblem, BoundaryProblem}, postprocessor_name::Symbol, args...; kwargs...) + push!(problem.postprocessors, (postprocessor_name, args, kwargs)) +end + +function add_preprocessor!(problem::Union{FieldProblem, BoundaryProblem}, postprocessor_name::Symbol, args...; kwargs...) + push!(problem.preprocessors, (postprocessor_name, args, kwargs)) end typealias Problem FieldProblem diff --git a/test/test_mortar_2d_contact.jl b/test/test_mortar_2d_contact.jl index fe28fd5..4a24598 100644 --- a/test/test_mortar_2d_contact.jl +++ b/test/test_mortar_2d_contact.jl @@ -6,16 +6,77 @@ using JuliaFEM.Test using JuliaFEM.Core: Node, Seg2, Tri3, update!, calculate_normal_tangential_coordinates!, PlaneStressLinearElasticityProblem, DirichletProblem, MortarProblem, get_elements, DirectSolver, calculate_nodal_vector, FieldAssembly, - FieldProblem, set_linear_system_solver!, set_nonlinear_max_iterations! + FieldProblem, set_linear_system_solver!, set_nonlinear_max_iterations!, + get_elements, Element, BoundaryAssembly, BoundaryProblem, + get_integration_points, get_jacobian, get_connectivity, StandardBasis, + add_postprocessor!, add_preprocessor! -import JuliaFEM.Core: postprocess_assembly!, linear_system_solver_preprocess! +import JuliaFEM.Core: assemble_postprocess!, linear_system_solver_preprocess! -function postprocess_assembly!(assembly::FieldAssembly, problem::FieldProblem{PlaneStressLinearElasticityProblem}, time::Real) - f = full(assembly.force_vector) - info("force vector = $(f')") +function calculate_normal_tangential_coordinates(elements::Vector{Element}, time::Real) + P = SparseMatrixCOO() + field_dim = 2 + for element in elements + haskey(element, "normal-tangential coordinates") || continue + for ip in get_integration_points(element, Val{2}) + J = get_jacobian(element, ip, time) + w = ip.weight*norm(J) + nt = transpose(element("normal-tangential coordinates", ip, time)) + normal = nt[1,:] + tangent = nt[2,:] + for nid in get_connectivity(element) + ndofs = [2*(nid-1)+1, 2*(nid-1)+2] + add!(P, [2*(nid-1)+1], ndofs, normal) + add!(P, [2*(nid-1)+2], ndofs, tangent) + end + end + end + P = sparse(P) + for i=1:size(P,1) + n = norm(P[i,:]) + if n > 0.0 + P[i,:] = P[i,:] / n + end + end + return P +end + +function assemble_postprocess!(assembly, problem, time::Real, ::Type{Val{:remove_tangential_constraints}}) + # example how to use postprocessor to manipulate constraint matrix before summing assemblies together + info("postprocess mortar assembly: remove contraints in tangent direction on boundary.") + C1 = sparse(assembly.C1) + C2 = sparse(assembly.C2) + dim = size(C1, 1) + P = calculate_normal_tangential_coordinates(get_elements(problem), time) + info("projection matrix for normals: ") + dump(round(full(P), 3)) + C1 = P*C1 + C2 = P*C2 + for i=2:2:dim + C1[i,:] = 0 + C2[i,:] = 0 + end + assembly.C1 = C1 + assembly.C2 = C2 + info("postprocess mortar assembly: done.") +end + +function assemble_postprocess!(assembly, problem, time::Real, ::Type{Val{:remove_constraint_from_dofs}}, dofs) + # giving additional arguments and keywords is possible too + C1 = sparse(assembly.C1) + C2 = sparse(assembly.C2) + info("removing constraints from dofs: $dofs") + info(round(full(C1), 3)) + for d in dofs + C1[d,:] = 0 + C2[d,:] = 0 + end + assembly.C1 = C1 + assembly.C2 = C2 end function linear_system_solver_preprocess!(solver, iter, time, K, f, C1, C2, D, g, sol, la, ::Type{Val{:foo}}) + # example how to use preprocessor to dump matrices before solution info("stiffness matrix") dump(round(full(K), 3)) info("constraint matrix C1") @@ -26,7 +87,6 @@ function linear_system_solver_preprocess!(solver, iter, time, K, f, C1, C2, D, g dump(round(full(f)', 3)) info("constraint vector") dump(round(full(g)', 3)) - end macro debug(msg) @@ -35,13 +95,14 @@ macro debug(msg) end @testset "2d frictionless contact" begin + gap = [0.0, 0.0] nodes = Node[ [6.0, 6.0], - [7.0, 8.0], + [6.0, 8.0]+gap, [0.0, 0.0], [6.0, 0.0], - [7.0, 0.0], - [19.0, 0.0]] + [6.0, 0.0]+gap, + [18.0, 0.0]+gap] fel1 = Tri3([1, 3, 4]) fel2 = Tri3([2, 5, 6]) force = Seg2([3, 1]) @@ -49,34 +110,37 @@ end bnd2 = Seg2([5, 6]) sel = Seg2([1, 4]) mel = Seg2([5, 2]) - update!([fel1, fel2, force, sel, mel, bnd1, bnd2], "geometry", nodes) + update!([fel1, fel2, force, sel, mel], "geometry", nodes) + update!([bnd1, bnd2], "geometry", nodes) - prob = PlaneStressLinearElasticityProblem() + prob = FieldProblem(PlaneStressLinearElasticityProblem, "bodies", 2) push!(prob, fel1, fel2) push!(prob, force) update!([fel1, fel2], "youngs modulus", 90.0) update!([fel1, fel2], "poissons ratio", 0.25) update!([force], "displacement traction force 1", 6/sqrt(2)) - fel1["displacement nodal load"] = Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [18.0, 0.0]] - bc = DirichletProblem("displacement", 2) + bc = BoundaryProblem(DirichletProblem, "support", "displacement", 2) push!(bc, bnd1, bnd2) update!(get_elements(bc), "displacement", 0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]]) - cont = MortarProblem("displacement", 2) + cont = BoundaryProblem(MortarProblem, "contact", "displacement", 2) push!(cont, sel, mel) calculate_normal_tangential_coordinates!(sel, 0.0) + nt = sel("normal-tangential coordinates", [0.0], 0.0) + info("normal direction = $(nt)") sel["master elements"] = [mel] + # remove coefficients from node 4 (dofs 7-8) because this conflicts with dirichlet bc. + add_postprocessor!(cont, :remove_constraint_from_dofs, [7, 8]) @debug begin info("fel1.fields = $(fel1.fields)") - info("bnd1.fields = $(bnd1.fields)") end solver = DirectSolver() push!(solver, prob) push!(solver, bc) -# push!(solver, cont) + push!(solver, cont) set_linear_system_solver!(solver, :UMFPACK) set_nonlinear_max_iterations!(solver, 2) push!(solver.linear_system_solver_preprocessors, :foo) @@ -87,9 +151,9 @@ end u = calculate_nodal_vector("displacement", 2, get_elements(prob), time) info("solution vector") dump(reshape(round(u, 8), 2, 6)) - la = calculate_nodal_vector("reaction force", 2, get_elements(prob), time) - info("reaction force") - dump(reshape(round(la, 8), 2, 6)) +# la = calculate_nodal_vector("reaction force", 2, get_elements(prob), time) +# info("reaction force") +# dump(reshape(round(la, 8), 2, 6)) end end