diff --git a/README.md b/README.md
index a33fda8..99f340d 100644
--- a/README.md
+++ b/README.md
@@ -25,4 +25,51 @@ Testing package can be done using `Pkg.test`, i.e.
 julia> Pkg.test("Mortar2D")
 ```
 
-Probably the easiest way to test the functionality of package is to use [JuliaBox](https://juliabox.com/).
+Probably the easiest way to test the functionality of package is to
+use [JuliaBox](https://juliabox.com/).
+
+## Usage example
+
+Let us calculate projection matrices D and M for the following problem:
+
+![](docs/src/figs/poisson_problem_discretized.png)
+
+Problem setup:
+
+```julia
+Xs = Dict(1 => [0.0, 1.0], 2 => [5/4, 1.0], 3 => [2.0, 1.0])
+Xm = Dict(4 => [0.0, 1.0], 5 => [1.0, 1.0], 6 => [2.0, 1.0])
+coords = merge(Xm , Xs)
+Es = Dict(1 => [1, 2], 2 => [2, 3])
+Em = Dict(3 => [4, 5], 4 => [5, 6])
+elements = merge(Es, Em)
+element_types = Dict(1 => :Seg2, 2 => :Seg2, 3 => :Seg2, 4 => :Seg2)
+slave_element_ids = [1, 2]
+master_element_ids = [3, 4]
+```
+
+Calculate projection matrices D and M
+
+```julia
+s, m, D, M = calculate_mortar_assembly(
+    elements, element_types, coords,
+    slave_element_ids, master_element_ids)
+```
+
+According to theory, the interface should transfer constant without any
+error. Let's test that:
+
+```julia
+u_m = ones(3)
+u_s = D[s,s] \ (M[s,m]*um)
+
+# output
+
+3-element Array{Float64,1}:
+ 1.0
+ 1.0
+ 1.0
+```
+
+The rest of the story can be read from the [documentation](https://juliafem.github.io/Mortar2D.jl/latest/).
+There's also brief review to the theory behind non-conforming finite element meshes.
diff --git a/docs/src/figs/poisson_problem_discretized.png b/docs/src/figs/poisson_problem_discretized.png
new file mode 100644
index 0000000..8d3d528
Binary files /dev/null and b/docs/src/figs/poisson_problem_discretized.png differ
diff --git a/docs/src/theory.md b/docs/src/theory.md
index 5df1d8e..3b73a53 100644
--- a/docs/src/theory.md
+++ b/docs/src/theory.md
@@ -151,8 +151,8 @@ and ``\boldsymbol{M}`` or ``\boldsymbol{P}=\boldsymbol{D}^{-1}\boldsymbol{M}``.
 Calculation projection matrix ``\boldsymbol{P}`` is implemented as function `calculate_mortar_assembly`:
 
 ```@example 0
-P = calculate_mortar_assembly(elements, element_types, coords,
-                              slave_element_ids, master_element_ids)
+s, m, D, M = calculate_mortar_assembly(elements, element_types, coords,
+                                       slave_element_ids, master_element_ids)
 ```
 
 This last command combines everything above to single command to calculate
diff --git a/src/Mortar2D.jl b/src/Mortar2D.jl
index 98c5ffa..8138d39 100644
--- a/src/Mortar2D.jl
+++ b/src/Mortar2D.jl
@@ -7,4 +7,5 @@ include("calculate_projections.jl")
 include("calculate_segments.jl")
 include("calculate_mortar_matrices.jl")
 include("calculate_mortar_assembly.jl")
+export calculate_mortar_assembly
 end
diff --git a/src/calculate_mortar_assembly.jl b/src/calculate_mortar_assembly.jl
index 029b1dd..69d52a2 100644
--- a/src/calculate_mortar_assembly.jl
+++ b/src/calculate_mortar_assembly.jl
@@ -8,14 +8,17 @@
                               slave_element_ids::Vector{Int},
                               master_element_ids::Vector{Int})
 
-Given data, calculate projection matrix `P`. This is the main function of
-package.
+Given data, calculate projection matrices `D` and `M`. This is the main
+function of package. Relation between matrices is ``D u_s = M u_m``, where
+``u_s`` is slave nodes and ``u_m`` master nodes.
 
-Matrix ``P`` is defined as a projection between slave and master surface,
-i.e.
-```math
-    D u_s = M u_m \\Rightarrow u_s = D^{-1} M u_m = P u_m.
+# Example
+
+```julia
+s, m, D, M = calculate_mortar_assembly(elements, element_types, coords,
+                                       slave_element_ids, master_element_ids)
 ```
+
 """
 function calculate_mortar_assembly(elements::Dict{Int, Vector{Int}},
                                    element_types::Dict{Int, Symbol},
@@ -65,11 +68,8 @@ function calculate_mortar_assembly(elements::Dict{Int, Vector{Int}},
         end
     end
 
-    D = sparse(D_I, D_J, D_V)
-    Df = cholfact(1/2*(D+D')[S,S])
-
-    M_ = sparse(M_I, M_J, M_V)
-    P =  Df \ M_[S,M]
-    SparseArrays.droptol!(P, 1.0e-12)
-    return P
+    # return global matrices + slave and master dofs
+    Dg = sparse(D_I, D_J, D_V)
+    Mg = sparse(M_I, M_J, M_V)
+    return S, M, Dg, Mg
 end
diff --git a/test/test_calculate_mortar_assembly.jl b/test/test_calculate_mortar_assembly.jl
index dd2ef97..ef2fedc 100644
--- a/test/test_calculate_mortar_assembly.jl
+++ b/test/test_calculate_mortar_assembly.jl
@@ -19,18 +19,21 @@ using Mortar2D: calculate_mortar_assembly
     elements = merge(Es, Em)
     coords = merge(Xs, Xm)
 
-    P = calculate_mortar_assembly(elements, element_types, coords,
-                                  slave_element_ids, master_element_ids)
+    s, m, D, M = calculate_mortar_assembly(elements,
+                                           element_types,
+                                           coords,
+                                           slave_element_ids,
+                                           master_element_ids)
 
     # test case 1, constant pressure load
     um = ones(nm)
-    us = P*um
+    us = D[s,s] \ (M[s,m]*um)
     us_expected = ones(ns)
     @test isapprox(us, us_expected)
 
     # test case 2, linear load
     um = linspace(0, 1, nm)
-    us = P*um
+    us = D[s,s] \ (M[s,m]*um)
     us_expected = linspace(0, 1, ns)
     @test isapprox(us, us_expected)
 end