diff --git a/tests/conftest.py b/tests/conftest.py
index c83ff625..255ce71b 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,238 +1,267 @@
 # -*- coding: utf-8 -*-
 # @file
 # @section LICENSE
 #
 # Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
 # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Affero General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU Affero General Public License for more details.
 #
 # You should have received a copy of the GNU Affero General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 from __future__ import division, print_function
 
+import subprocess
+import sys
+
 import numpy as np
-import tamaas as tm
 import pytest
 
 from numpy.linalg import norm
 
+try:
+    import tamaas as tm
+except ModuleNotFoundError as e:
+    print(e, file=sys.stderr)
+    print("Use 'scons test' to run tests", file=sys.stderr)
+    sys.exit(1)
+
+
+def get_libtamaas_file():
+    try:
+        ldd_out = bytes(subprocess.check_output(['ldd', tm._tamaas.__file__]))
+        for line in filter(lambda x: 'Tamaas' in x,
+                           ldd_out.decode().split('\n')):
+            path = line.split(' => ')[1]
+            return path.split(' ')[0]
+
+    except subprocess.CalledProcessError:
+        return "N/A"
+
+
+def pytest_report_header(config):
+    print('tamaas build type: {}\nmodule file: {}\nwrapper: {}\nlibTamaas: {}'
+          .format(tm.TamaasInfo.build_type,
+                  tm.__file__,
+                  tm._tamaas.__file__,
+                  get_libtamaas_file()))
+
 
 def profile(func, size, mode, amplitude):
     x = np.linspace(0, 1, size[0], endpoint=False, dtype=tm.dtype)
     y = np.linspace(0, 1, size[1], endpoint=False, dtype=tm.dtype)
     x, y = np.meshgrid(x, y, indexing='ij')
     surface = amplitude * func(2*np.pi*x*mode[0]) * func(2*np.pi*y*mode[1])
     return surface.copy()
 
 
 @pytest.fixture(scope="session")
 def tamaas_fixture():
     tm.initialize()
     yield None
     tm.finalize()
 
 
 class HertzFixture:
     def __init__(self, n, load):
         self.domain_size = 1
         self.n = n
         self.load = load
         self.curvature = 0.1
         self.radius = 1. / self.curvature
         self.e_star = 1.
         self.a = (3 * load / (4 * self.curvature * self.e_star))**(1. / 3.)
 
         self.x = np.linspace(-self.domain_size / 2.,
                              self.domain_size / 2.,
                              self.n, dtype=tm.dtype)
         self.y = self.x.copy()
         self.x, self.y = np.meshgrid(self.x, self.y)
 
         self._computeSurface()
         self._computePressure()
         self._computeDisplacement()
 
     def _computeDisplacement(self):
         r = np.sqrt(self.x**2 + self.y**2)
         self.displacement = np.zeros_like(r)
         contact = r < self.a
         self.displacement[contact] = self.surface[contact]
         self.displacement[~contact] = \
             (self.surface[~contact]
              + self.a / (np.pi * self.radius) * np.sqrt(r[~contact]**2
                                                         - self.a**2)
              + (r[~contact]**2 - 2 * self.a**2) / (np.pi * self.radius)
              * np.arccos(self.a / r[~contact]))
 
     def _computePressure(self):
         r = np.sqrt(self.x**2 + self.y**2)
         self.pressure = np.zeros_like(r)
         contact = np.where(r < self.a)
         self.pressure[contact] = \
             2 * self.e_star / (np.pi * self.radius) \
             * np.sqrt(self.a**2 - r[contact]**2)
 
     def _computeSurface(self):
         self.surface = -1. / (2 * self.radius) * (self.x**2 + self.y**2)
 
 
 @pytest.fixture(scope="module")
 def hertz(tamaas_fixture):
     return HertzFixture(1024, 0.00001)
 
 
 @pytest.fixture(scope="module")
 def hertz_coarse(tamaas_fixture):
     return HertzFixture(512, 0.0001)
 
 @pytest.fixture(scope="module",
                 params=[tm.PolonskyKeerRey.pressure])
 def pkr(hertz, request):
     model = tm.ModelFactory.createModel(tm.model_type.basic_2d,
                                         [hertz.domain_size, hertz.domain_size],
                                         [hertz.n, hertz.n])
     model.setElasticity(hertz.e_star, 0)
     solver = tm.PolonskyKeerRey(model, hertz.surface, 1e-12,
                                 request.param,
                                 request.param)
     solver.solve(hertz.load)
 
     return model, hertz
 
 
 class WestergaardFixture:
     def __init__(self, n, load):
         self.domain_size = 1.
         self.lamda = 1.
         self.delta = 0.1
         self.e_star = 1.
 
         self.n = n
 
         self.p_star = np.pi * self.e_star * self.delta / self.lamda
 
         self.load = load * self.p_star
 
         self.a = self.lamda / np.pi \
             * np.arcsin(np.sqrt(self.load / self.p_star))
 
         self.x = np.linspace(-self.domain_size / 2.,
                              self.domain_size / 2.,
                              self.n, endpoint=False, dtype=tm.dtype)
 
         self._computeSurface()
         self._computePressure()
         self._computeDisplacement()
 
     def _computeSurface(self):
         self.surface = self.delta * np.cos(2 * np.pi * self.x / self.lamda)
 
     def _computePressure(self):
         self.pressure = np.zeros_like(self.surface)
         contact = np.where(np.abs(self.x) < self.a)
         self.pressure[contact] = 2 * self.load \
             * (np.cos(np.pi * self.x[contact] / self.lamda)
                / np.sin(np.pi * self.a / self.lamda)**2) \
             * np.sqrt(np.sin(np.pi * self.a / self.lamda)**2 -
                       np.sin(np.pi * self.x[contact] / self.lamda)**2)
 
     def _computeDisplacement(self):
         psi = np.pi * np.abs(self.x) / self.lamda
         psi_a = np.pi * self.a / self.lamda
 
         with np.errstate(invalid='ignore'):  # get some warnings out of the way
             self.displacement = (np.cos(2*psi) + 2 * np.sin(psi)
                                  * np.sqrt(np.sin(psi)**2 - np.sin(psi_a)**2)
                                  - 2 * np.sin(psi_a)**2 *
                                  np.log((np.sin(psi) +
                                          np.sqrt(np.sin(psi)**2
                                                  - np.sin(psi_a)**2))
                                         / np.sin(psi_a)))
 
         contact = np.where(np.abs(self.x) < self.a)
         self.displacement[contact] = np.cos(2*psi[contact])
         self.displacement *= self.load * self.lamda / (np.pi * self.e_star
                                                        * np.sin(psi_a)**2)
 
 
 @pytest.fixture(scope="module")
 def westergaard(tamaas_fixture):
     return WestergaardFixture(19683, 0.1)
 
 
 class PatchWestergaard:
     def __init__(self, model_type, domain, modes, size):
         self.E = 3.
         self.nu = 0.
         self.e_star = self.E / (1 - self.nu**2)
         self.n = size
         self.domain = domain
         self.model = tm.ModelFactory.createModel(model_type, domain,
                                                  size)
         self.model.setElasticity(self.E, self.nu)
 
         self.pressure = profile(np.cos, size, modes, 1)
         self.solution = profile(np.cos, size, modes,
                                 1 / (np.pi * self.e_star * norm(modes)))
 
 
 @pytest.fixture(scope="module", params=[tm.model_type.basic_2d])
 def patch_westergaard(tamaas_fixture, request):
     return PatchWestergaard(request.param, [1., 1.], [3, 1], [6, 6])
 
 
 class UniformPlasticity:
     def __init__(self, model_type, domain, sizes):
         self.n = sizes
         self.domain = domain
         self.model = tm.ModelFactory.createModel(model_type, domain,
                                                  sizes)
         self.E_h = 0.1
         self.sigma_y = 0.01
         self.residual = tm.ModelFactory.createResidual(self.model,
                                                        sigma_y=self.sigma_y,
                                                        hardening=self.E_h)
         self.model.E = 1.
         self.model.nu = 0.
 
     def solution(self, p):
         E, nu = self.model.E, self.model.nu
         E_h, sigma_y = self.E_h, self.sigma_y
         mu = E / (2 * (1 + nu))
 
         strain = -1 / (mu + E_h) * (p * (3 * mu + E_h) / (2 * mu)
                                     - np.sign(p) * sigma_y)
 
         dep = (2 * mu * np.abs(strain) - sigma_y) / (3 * mu + E_h)
         plastic_strain = np.sign(strain) / 2 * dep * np.array([
             -1, -1, 2, 0, 0, 0,
         ], dtype=tm.dtype)
 
         stress = 2 * mu * (np.array([
             0, 0, strain, 0, 0, 0
         ], dtype=tm.dtype) - plastic_strain)
 
         return {
             "stress": stress,
             "plastic_strain": plastic_strain,
             "cumulated_plastic_strain": dep
         }, {
             "stress": mu,
             "plastic_strain": 1,
             "cumulated_plastic_strain": 1
         }
 
 
 @pytest.fixture(scope="module", params=[tm.model_type.volume_2d])
 def patch_isotropic_plasticity(tamaas_fixture, request):
     return UniformPlasticity(request.param, [1., 1., 1.], [4, 4, 4])