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

@@ -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

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
 

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

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