diff --git a/src/abaqus.jl b/src/abaqus.jl
index b2da1dc..b8294d3 100644
--- a/src/abaqus.jl
+++ b/src/abaqus.jl
@@ -537,23 +537,35 @@ function process_output_request(model::Model, solver::Solver, output_request::Ab
                                 ::Type{Val{:NODE}}, ::Type{Val{:PRINT}})
     data = output_request.data
     options = output_request.options
-    info("nodal output request with data $data and options $options")
-    code_mapping = Dict(:U => "displacement", :COORD => "geometry")
-    abbr_mapping = Dict(:U => :U, :COORD => :COOR) # ..?
+    code_mapping = Dict(
+        :COORD => "geometry",
+        :U => "displacement",
+        :CF => "concentrated force",
+        :RF => "reaction force")
+    abbr_mapping = Dict(:COORD => :COOR)
     for row in data
+        info(repeat("-", 80))
+        codes = join(row, ", ")
+        info("*NODE OUTPUT request, with fields $codes")
+        if length(options) != 0
+            info("Additional options: $options")
+        end
+        info(repeat("-", 80))
         tables = Any[]
         for code in row
             haskey(code_mapping, code) || continue
             for problem in model.problems
                 field_name = code_mapping[code]
-                abbr = abbr_mapping[code]
+                abbr = get(abbr_mapping, code, code)
                 table = problem(DataFrame, field_name, abbr, solver.time)
                 push!(tables, table)
             end
         end
         length(tables) != 0 || continue
         results = join(tables..., on=:id, kind=:outer)
+        println()
         println(results)
+        println()
     end
 end
 
@@ -590,6 +602,9 @@ function call(model::Model)
         all_problems = [model.problems; boundary_problems; step_problems]
         push!(solver, all_problems...)
         solver()
+        info(repeat("-", 80))
+        info("Simulation ready, processing output requests")
+        info(repeat("-", 80))
         # 3.3 postprocessing based on output requests
         for output_request in step.output_requests
             process_output_request(model, solver, output_request)
diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl
index d582de3..7a2e9b5 100644
--- a/test/test_elasticity_2d_linear_with_surface_load.jl
+++ b/test/test_elasticity_2d_linear_with_surface_load.jl
@@ -6,43 +6,42 @@ using JuliaFEM.Preprocess
 using JuliaFEM.Postprocess
 using JuliaFEM.Testing
 
-function JuliaFEM.get_model(::Type{Val{Symbol("test 2d linear elasticity with surface + volume load")}})
+#=
+- solve 2d plane stress problem with known solution
+- test postprocessing of nodal fields: (geometry, displacement
+  reaction force, concentrated force)
+=#
+@testset "test 2d linear elasticity with surface + volume load" begin
     meshfile = "/geometry/2d_block/BLOCK_1elem.med"
     mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*meshfile)
 
     # field problem
     block = Problem(Elasticity, "BLOCK", 2)
-    block.properties.store_fields = ["stress", "strain"]
     block.properties.formulation = :plane_stress
     block.properties.finite_strain = false
     block.properties.geometric_stiffness = false
-
     block.elements = create_elements(mesh, "BLOCK")
     update!(block.elements, "youngs modulus", 288.0)
     update!(block.elements, "poissons ratio", 1/3)
     update!(block.elements, "displacement load 2", 576.0)
 
-    traction = create_elements(mesh, "TOP")
+    # traction
+    traction = Problem(Elasticity, "BLOCK", 2)
+    traction.properties.formulation = :plane_stress
+    traction.properties.finite_strain = false
+    traction.properties.geometric_stiffness = false
+    traction.elements = create_elements(mesh, "TOP")
     update!(traction, "displacement traction force 2", 288.0)
-    push!(block, traction...)
 
     # boundary conditions
-    bc_sym = Problem(Dirichlet, "symmetry bc", 2, "displacement")
-    bc_elements_left = create_elements(mesh, "LEFT")
-    bc_elements_bottom = create_elements(mesh, "BOTTOM")
-    update!(bc_elements_left, "displacement 1", 0.0)
-    update!(bc_elements_bottom, "displacement 2", 0.0)
-    push!(bc_sym, bc_elements_left..., bc_elements_bottom...)
+    bc_sym_23 = Problem(Dirichlet, "symmetry bc 23", 2, "displacement")
+    bc_sym_23.elements = create_elements(mesh, "LEFT")
+    update!(bc_sym_23, "displacement 1", 0.0)
+    bc_sym_13 = Problem(Dirichlet, "symmetry bc 13", 2, "displacement")
+    bc_sym_13.elements = create_elements(mesh, "BOTTOM")
+    update!(bc_sym_13, "displacement 2", 0.0)
 
-    solver = LinearSolver("solve block problem")
-    push!(solver, block, bc_sym)
-    return solver
-end
-
-@testset "test 2d linear elasticity with surface + volume load" begin
-
-    solver = get_model("test 2d linear elasticity with surface + volume load")
-    block, bc_sym = solver.problems
+    solver = LinearSolver(block, traction, bc_sym_23, bc_sym_13)
     solver()
 
     f = 288.0
@@ -51,11 +50,13 @@ end
     nu = 1/3
     u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1]
 
-    results = block(DataFrame, "displacement", :U, 0.0)
+    # fetch nodal results X + u and join them into one table using DataFrames
+    X = block(DataFrame, "geometry", :COOR, 0.0)
+    u = block(DataFrame, "displacement", :U, 0.0)
+    results = join(X, u, on=:id, kind=:outer)
     println(results)
+
     u3 = results[:N3, [:U1, :U2]]
-    info("(u1,u2) at node 3")
-    info(u3)
     @test isapprox(u3, u3_expected)
 
 #=