diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb
index d791da8..f5314b4 100644
--- a/notebooks/2015-08-29-developing-juliafem.ipynb
+++ b/notebooks/2015-08-29-developing-juliafem.ipynb
@@ -1499,7 +1499,7 @@
     "    # solve problem\n",
     "    nz = unique(rowvals(A))\n",
     "    x = zeros(b)\n",
-    "    x[nz] = lufact(Atot[nz,nz]) \\ full(b[nz])\n",
+    "    x[nz] = lufact(A[nz,nz]) \\ full(b[nz])\n",
     "\n",
     "    # get \"problem-wise\" solution vectors\n",
     "    x1 = x[1:length(b1)]\n",
diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index de5b55b..e462ed1 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -7,6 +7,14 @@ using Lexicon
 using Logging
 @Logging.configure(level=DEBUG)
 
+"""
+Simple linspace extension to multidimensional values. Contribute to julialang?
+"""
+function Base.linspace(X1, X2, n)
+    [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]
+end
+
+
 include("types.jl")  # type definitions
 
 include("interpolate.jl")  # interpolation routines
@@ -16,16 +24,18 @@ include("elements.jl")
 include("lagrange.jl") # Lagrange elements
 #include("hierarchical.jl") # P-elements
 
+include("equations.jl")
+include("problems.jl")
+include("solvers.jl")
 
-include("equations.jl") # formulations
-include("problems.jl") # problems
-
-
-
-include("math.jl") # basic mathematical operations -- obsolete ..?
-include("elasticity_solver.jl")
+#include("math.jl") # basic mathematical operations -- obsolete ..?
+# pre- and postprocess
 include("xdmf.jl")
 include("abaqus_reader.jl")
-include("interfaces.jl")
+#include("interfaces.jl")
+
+include("dirichlet.jl")
+include("heat.jl")
+#include("elasticity_solver.jl")
 
 end # module
diff --git a/src/dirichlet.jl b/src/dirichlet.jl
new file mode 100644
index 0000000..5451760
--- /dev/null
+++ b/src/dirichlet.jl
@@ -0,0 +1,52 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# Dirichlet boundary conditions in weak form
+
+abstract DirichletEquation <: Equation
+
+get_unknown_field_name(eq::DirichletEquation) = symbol("reaction force")
+
+### Dirichlet problem + equations
+
+type DirichletProblem <: BoundaryProblem
+    equations :: Array{DirichletEquation, 1}
+end
+function DirichletProblem()
+    DirichletProblem([])
+end
+get_dimension(pr::Type{DirichletProblem}) = 1  # ..?
+get_equation(pr::Type{DirichletProblem}, el::Type{Seg2}) = DBC2D2
+
+"""
+Dirichlet boundary condition element for 2 node line segment
+"""
+type DBC2D2 <: DirichletEquation
+    element :: Seg2
+    integration_points :: Array{IntegrationPoint, 1}
+    global_dofs :: Array{Int64, 1}
+    fieldval :: Function
+end
+function DBC2D2(el::Seg2)
+    integration_points = [
+        IntegrationPoint([-sqrt(1/3)], 1.0),
+        IntegrationPoint([+sqrt(1/3)], 1.0)]
+    new_fieldset!(el, "reaction force")
+    fieldval(X, t) = 0.0
+    DBC2D2(el, integration_points, [], fieldval)
+end
+function get_lhs(eq::DBC2D2, ip, t)
+    el = get_element(eq)
+    h = get_basis(el)(ip.xi)
+    return h*h'
+end
+function get_rhs(eq::DBC2D2, ip, t)
+    el = get_element(eq)
+    h = get_basis(el, ip.xi)
+    f = eq.fieldval
+    X = interpolate(el, "geometry", ip.xi, t)
+    return h*f(X, t)
+end
+has_lhs(eq::DBC2D2) = true
+has_rhs(eq::DBC2D2) = true
+
diff --git a/src/heat.jl b/src/heat.jl
new file mode 100644
index 0000000..c45dd49
--- /dev/null
+++ b/src/heat.jl
@@ -0,0 +1,75 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# Heat problems
+
+abstract HeatProblem <: Problem
+abstract HeatEquation <: Equation
+
+get_unknown_field_name(eq::HeatEquation) = symbol("temperature")
+
+
+### Plane heat problem + equations ###
+
+type PlaneHeatProblem <: HeatProblem
+    equations :: Array{HeatEquation, 1}
+end
+
+""" Default constructor for problem takes no arguments. """
+function PlaneHeatProblem()
+    return PlaneHeatProblem([])
+end
+
+""" Return dimension of unknown field variable, temperature is scalar field. """
+get_dimension(pr::Type{PlaneHeatProblem}) = 1
+
+""" Map Lagrange element Quad4 to equation DC2D4 """
+get_equation(pr::Type{PlaneHeatProblem}, el::Type{Quad4}) = DC2D4
+
+""" Map Lagrange element Seg2 to equation DC2D2 """
+get_equation(pr::Type{PlaneHeatProblem}, el::Type{Seg2}) = DC2D2
+
+
+""" Diffusive heat transfer for 4-node bilinear element. """
+type DC2D4 <: HeatEquation
+    element :: Quad4
+    integration_points :: Array{IntegrationPoint, 1}
+    global_dofs :: Array{Int64, 1}
+end
+function DC2D4(el::Quad4)
+    integration_points = [
+        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]
+    new_fieldset!(el, "temperature")
+    DC2D4(el, integration_points, [])
+end
+function get_lhs(eq::DC2D4, ip, t)
+    el = get_element(eq)
+    dNdX = get_dbasisdX(el, ip.xi, t)
+    k = interpolate(el, "temperature thermal conductivity", ip.xi, t)
+    return dNdX*k*dNdX'
+end
+JuliaFEM.has_lhs(eq::DC2D4) = true
+
+""" Diffusive heat transfer for 2-node linear segment. """
+type DC2D2 <: HeatEquation
+    element :: Seg2
+    integration_points :: Array{IntegrationPoint, 1}
+    global_dofs :: Array{Int64, 1}
+end
+function DC2D2(el::Seg2)
+    integration_points = [IntegrationPoint([0.0], 1.0)]
+    new_fieldset!(el, "temperature")
+    DC2D2(el, integration_points, [])
+end
+function get_rhs(eq::DC2D2, ip, t)
+    el = get_element(eq)
+    h = get_basis(el, ip.xi)
+    f = interpolate(el, "temperature flux", ip.xi, t)
+    println("h = $h")
+    println("f = $f")
+    return h*f
+end
+JuliaFEM.has_rhs(eq::DC2D2) = true
diff --git a/src/heat_equations.jl b/src/heat_equations.jl
deleted file mode 100644
index c09d93a..0000000
--- a/src/heat_equations.jl
+++ /dev/null
@@ -1,58 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-abstract Heat <: Equation
-
-"""
-Diffusive heat transfer for 4-node bilinear element.
-"""
-type DC2D4 <: Heat
-    element :: Quad4
-    integration_points :: Array{IntegrationPoint, 1}
-    global_dofs :: Array{Int64, 1}
-end
-function DC2D4(el::Quad4)
-    integration_points = [
-        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
-        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
-        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),
-        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]
-    set_field(el, :temperature, zeros(2, 4))
-    DC2D4(el, integration_points, [])
-end
-function get_lhs(eq::DC2D4)
-    function get_lhs_(eq, ip)
-        el = get_element(eq)
-        xi = ip.xi
-        dNdX = get_dbasisdX(el, xi)'
-        hc = interpolate(el, :"temperature heat coefficient", xi)
-        return dNdX'*hc*dNdX
-    end
-    integrate(eq, get_lhs_)
-end
-
-
-"""
-Diffusive heat transfer for 2-node linear segment.
-"""
-type DC2D2 <: Heat
-    element :: Seg2
-    integration_points :: Array{IntegrationPoint, 1}
-    global_dofs :: Array{Int64, 1}
-end
-function DC2D2(el::Seg2)
-    integration_points = [IntegrationPoint([0.0], 1.0)]
-    set_field(el, :temperature, zeros(2, 1))
-    set_field(el, :"temperature flux", zeros(2, 1))
-    DC2D2(el, integration_points, [])
-end
-function get_rhs(eq::DC2D2)
-    function get_rhs_(eq, ip)
-        el = get_element(eq)
-        xi = ip.xi
-        N = get_basis(el, xi)
-        f = interpolate(el, :"temperature flux", xi)
-        return f*N
-    end
-    integrate(eq, get_rhs_)
-end
\ No newline at end of file
diff --git a/src/problems.jl b/src/problems.jl
index bae79d6..b9426ca 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -2,10 +2,18 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
 abstract Problem
+abstract BoundaryProblem <: Problem
+abstract FieldProblem <: Problem
 
 get_equations(pr::Problem) = pr.equations
-get_dimension(pr::Type{Problem}) = nothing
-get_equation(pr::Type{Problem}, el::Type{Element}) = nothing
+
+function get_dimension(pr::Type{Problem})
+    throw("Unable to determine problem dimension for problem $pr")
+end
+
+function get_equation(pr::Type{Problem}, el::Type{Element})
+    throw("Could not find corresponding equation for element $el in problem $pr")
+end
 
 """
 Add new element to problem
@@ -50,3 +58,83 @@ function set_global_dofs!(pr::Problem)
         set_global_dofs!(eq, gconn)
     end
 end
+
+function get_connectivity(pr::Problem)
+    conn = Int[]
+    for eq in get_equations(pr)
+        el = get_element(eq)
+        append!(conn, get_connectivity(el))
+    end
+    conn = unique(conn)
+    return conn
+end
+
+"""
+Calculate global dofs for equations, maybe using some bandwidth
+minimizing or fill reducing algorithm
+"""
+function calculate_global_dofs(pr::Problem)
+    conn = get_connectivity(pr)
+    dim = get_dimension(typeof(pr))
+    ndofs = dim*length(conn)
+    Logging.debug("total dofs: $ndofs")
+
+    mconn = maximum(conn)
+    gdofs = reshape(collect(1:mconn), dim, mconn)
+    dofmap = Dict{Int64, Array{Int64, 1}}()
+    for (i, c) in enumerate(conn)
+        dofmap[c] = gdofs[:, i]
+    end
+    return dofmap
+end
+
+"""
+Assign global dofs for equations.
+"""
+function assign_global_dofs!(pr::Problem, dofmap)
+    for eq in get_equations(pr)
+        el = get_element(eq)
+        c = get_connectivity(el)
+        #gdofs = [dofmap[ci] for ci in c]
+        gdofs = Int64[]
+        for ci in c
+            append!(gdofs, dofmap[ci])
+        end
+        set_global_dofs!(eq, gdofs)
+    end
+end
+
+function get_lhs(pr::Problem, t::Float64)
+    I = Int64[]
+    J = Int64[]
+    V = Float64[]
+    dim = get_dimension(typeof(pr))
+    for eq in filter(has_lhs, get_equations(pr))
+        dofs = get_global_dofs(eq)
+        lhs = integrate_lhs(eq, t)
+        for (li, i) in enumerate(dofs)
+            for (lj, j) in enumerate(dofs)
+                push!(I, i)
+                push!(J, j)
+                push!(V, lhs[li, lj])
+            end
+        end
+    end
+    return I, J, V
+end
+
+function get_rhs(pr::Problem, t::Float64)
+    I = Int64[]
+    V = Float64[]
+    dim = get_dimension(typeof(pr))
+    for eq in filter(has_rhs, get_equations(pr))
+        dofs = get_global_dofs(eq)
+        rhs = integrate_rhs(eq, t)
+        for (li, i) in enumerate(dofs)
+            push!(I, i)
+            push!(V, rhs[li])
+        end
+    end
+    return I, V
+end
+
diff --git a/src/solvers.jl b/src/solvers.jl
new file mode 100644
index 0000000..daf11ef
--- /dev/null
+++ b/src/solvers.jl
@@ -0,0 +1,95 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# Solver stuff
+
+abstract Solver
+
+"""
+Add new problem to solver
+"""
+function add_problem!(solver::Solver, problem::Problem)
+    push!(solver.problems, problem)
+end
+
+"""
+Get all problems assigned to solver
+"""
+function get_problems(s::Solver)
+    return s.problems
+end
+
+## SimpleSolver -- tiny direct demo solver
+""" Simple solver for educational purposes. """
+type SimpleSolver <: Solver
+    problems
+end
+""" Default initializer. """
+function SimpleSolver()
+    SimpleSolver(Problem[])
+end
+
+"""
+Call solver to solve a set of problems.
+
+This is simple serial solver for demonstration purposes. It handles the most
+common situation, i.e., some main field problem and it's Dirichlet boundary.
+"""
+function call(solver::SimpleSolver, t)
+    problems = get_problems(solver)
+    problem1 = problems[1]
+    problem2 = problems[2]
+
+    # calculate order of degrees of freedom in global matrix
+    # and set the ordering to problems
+    dofmap = calculate_global_dofs(problem1)
+    assign_global_dofs!(problem1, dofmap)
+    assign_global_dofs!(problem2, dofmap)
+
+    # assemble problem 1
+    A1 = sparse(get_lhs(problem1, t)...)
+    b1 = sparsevec(get_rhs(problem1, t)..., size(A1, 1))
+
+    # assemble problem 2
+    A2 = sparse(get_lhs(problem2, t)...)
+    b2 = sparsevec(get_rhs(problem2, t)..., size(A2, 1))
+    
+    # make one monolithic assembly
+    A = [A1 A2; A2' zeros(A2)]
+    b = [b1; b2]
+    dump(full(A))
+    dump(full(b))
+
+    # solve problem
+    nz = unique(rowvals(A))
+    x = zeros(b)
+    x[nz] = lufact(A[nz,nz]) \ full(b[nz])
+
+    # get "problem-wise" solution vectors
+    x1 = x[1:length(b1)]
+    x2 = x[length(b1)+1:end]
+
+    # check residual
+    R1 = A1*x1 - b1
+    R2 = A2*x2 - b2
+    println("Residual norm: $(norm(R1+R2))")
+
+    # update field for elements in problem 1
+    for equation in get_equations(problem1)
+        gdofs = get_global_dofs(equation)
+        element = get_element(equation)
+        field_name = get_unknown_field_name(equation)  # field we are solving
+        field = Field(t, full(x1[gdofs])[:])
+        add_field!(element, field_name, field)
+    end
+
+    # update field for elements in problem 2
+    for equation in get_equations(problem2)
+        gdofs = get_global_dofs(equation)
+        element = get_element(equation)
+        field_name = get_unknown_field_name(equation)
+        field = Field(t, full(x2[gdofs]))
+        add_field!(element, field_name, field)
+    end
+end
+