diff --git a/SConstruct b/SConstruct
index f714fcf..787cc53 100644
--- a/SConstruct
+++ b/SConstruct
@@ -1,315 +1,318 @@
 # ------------------------------------------------------------------------------
 # Imports
 # ------------------------------------------------------------------------------
 from __future__ import print_function
 import os
 from os.path import join, abspath, basename
 from version import write_info_file
 # ------------------------------------------------------------------------------
 
 EnsurePythonVersion(2, 7)
 EnsureSConsVersion(2, 4)
 
 
 # ------------------------------------------------------------------------------
 # Helper functions
 # ------------------------------------------------------------------------------
 def detect_fftw(env):
     """Detect fftw on clusters"""
     fftw_include = ""
     fftw_library = ""
 
     # If FFTW is provided by module system (on clusters)
     if 'FFTW_ROOT' in env['ENV']:
         fftw_include = join(env['ENV']['FFTW_ROOT'], 'include')
         fftw_library = join(env['ENV']['FFTW_ROOT'], 'lib')
 
     # Setting up FFTW
     env['FFTW_LIBRARY_WISH'] = ['main', 'omp']
     env['FFTW_INCLUDE_DIR'] = fftw_include
     env['FFTW_LIBRARY_DIR'] = fftw_library
     env.Tool(fftw)
 
 
 # ------------------------------------------------------------------------------
 def detect_cuda(env):
     """Detect cuda on clusters"""
     if 'CUDA_ROOT' in env['ENV']:
         env['CUDA_TOOLKIT_PATH'] = env['ENV']['CUDA_ROOT']
     else:
         env['CUDA_TOOLKIT_PATH'] = '/opt/cuda'
 
     env['CUDA_COMPONENTS'] = ['cufft']
     env['CUDA_ARCH_FLAG'] = '-arch=sm_60'
 
     colors = env['COLOR_DICT']
     if not env['verbose']:
         env['NVCCCOMSTR'] = env['SHCXXCOMSTR']
         env['SHLINKCOMSTR'] = u'{0}[Linking (cuda)] {1}$TARGET{2}'.format(colors['purple'],
                                                                           colors['blue'],
                                                                           colors['end'])
     env.AppendUnique(CXXFLAGS="-expt-extended-lambda")  # experimental lambda support
     env.AppendUnique(CXXFLAGS="-expt-relaxed-constexpr")  # experimental lambda support
     if env['build_type'] == 'debug':
         env.AppendUnique(CXXFLAGS="-G")
     env.Tool('nvcc')
 
 
 # ------------------------------------------------------------------------------
 def detect_boost(env):
     """Detect boost on clusters"""
     if 'BOOST_ROOT' in env['ENV']:
         env['BOOST_INCLUDE_DIR'] = join(env['ENV']['BOOST_ROOT'], 'include')
     env.Tool(boost)
 
 
 # ------------------------------------------------------------------------------
 def detect_thrust(env):
     """Detect cuda on clusters"""
     if 'CUDA_ROOT' in env['ENV']:
         env['THRUST_INCLUDE_DIR'] = join(env['ENV']['CUDA_ROOT'], 'include')
     else:
         env['THRUST_INCLUDE_DIR'] = '/opt/cuda/include'
     if basename(env['CXX']) == "clang++":
       env.AppendUnique(CXXFLAGS=["-Wno-unused-local-typedef"])
     env.Tool(thrust)
 
 
 # ------------------------------------------------------------------------------
 def detect_dependencies(env):
     """Detect all dependencies"""
     detect_fftw(env)
     detect_boost(env)
     detect_thrust(env)
 
     # Activate cuda if needed
     if env['backend'] == 'cuda':
 	detect_cuda(env)
 
 
 # ------------------------------------------------------------------------------
 def gen_print(action_string, color_string, env):
     """Generic function for creating pretty compile output"""
     if env['verbose']:
         return None
 
     def print_fun(command, target, source, env):
         colors = env['COLOR_DICT']
         print("{}[{}] {}{}".format(colors[color_string],
                                    action_string,
                                    colors['end'],
                                    target[0]))
     return print_fun
 
 
 # ------------------------------------------------------------------------------
 # Main compilation
 # ------------------------------------------------------------------------------
 
 # Compilation colors
 colors = {
   'cyan': '\033[96m',
   'purple': '\033[95m',
   'blue': '\033[94m',
   'green': '\033[92m',
   'yellow': '\033[93m',
   'red': '\033[91m',
   'end': '\033[0m'
 }
 
 # Inherit all environment variables (for CXX detection, etc.)
 main_env = Environment(ENV=os.environ)
 main_env['COLOR_DICT'] = colors
 
 # Compiler detection
 compiler_default = 'g++'
 if 'CXX' in os.environ:
     compiler_default = os.environ['CXX']
     print("Detected default compiler {} from environment".format(compiler_default))
 
 # Build variables
 vars = Variables('build-setup.conf')
 vars.Add(EnumVariable('build_type', 'Build type', 'release',
                       allowed_values=('release', 'profiling', 'debug'),
                       ignorecase=2))
 vars.Add(EnumVariable('backend', 'Thrust backend', 'omp',
                       allowed_values=('omp', 'cuda'),
                       ignorecase=2))
 vars.Add(EnumVariable('sanitizer', 'Sanitizer type', 'none',
                       allowed_values=('none', 'memory', 'leaks', 'address'),
                       ignorecase=2))
 vars.Add('prefix', 'Prefix where to install', '/usr/local')
 vars.Add('CXX', 'Compiler', compiler_default)
 vars.Add('py_exec', 'Python executable', 'python')
 
 # Cosmetic
 vars.Add(BoolVariable('verbose', 'Activate verbosity', False))
 vars.Add(BoolVariable('color', 'Color the non-verbose compilation output', False))
 
 # Tamaas components
 vars.Add(BoolVariable('build_doc', 'Build documentation', False))
 vars.Add(BoolVariable('build_tests', 'Build test suite', False))
 vars.Add(BoolVariable('build_python', 'Build python wrapper', False))
 
 # Dependencies
 vars.Add(BoolVariable('use_googletest', 'Build tests using GTest', False))
 
 vars.Update(main_env)
 
 Help(vars.GenerateHelpText(main_env))
 
 # Save all options, not just those that differ from default
 with open('build-setup.conf', 'w') as setup:
     for key in vars.keys():
         setup.write("{} = '{}'\n".format(key, main_env[key]))
 
 main_env['should_configure'] = not main_env.GetOption('clean') and not main_env.GetOption('help')
 
 build_type = main_env['build_type']
 build_dir = 'build-' + main_env['build_type']
 
 if main_env['should_configure']:
     print("Building in " + build_dir)
 
 verbose = main_env['verbose']
 
 # Remove colors if not set
 if not main_env['color']:
     for key in colors:
         colors[key] = ''
 
 # Setting object suffix
 main_env['SHOBJSUFFIX'] = '.o'
 
 if not verbose:
     main_env['SHCXXCOMSTR'] = u'{0}[Compiling ($SHCXX)] {1}$SOURCE'.format(colors['green'],
                                                                            colors['end'])
     main_env['SHLINKCOMSTR'] = u'{0}[Linking] {1}$TARGET'.format(colors['purple'],
                                                                  colors['end'])
 
 # Include paths
 main_env.AppendUnique(CPPPATH=['#/src',
                                '#/src/core',
                                '#/src/bem',
                                '#/src/surface',
                                '#/src/python',
                                '#/src/percolation',
                                '#/src/model',
                                '#/src/solvers',
                                '#/src/gpu',
                                '#/python'])
 
 # Changing the shared object extension
 main_env['SHOBJSUFFIX'] = '.o'
 
 # Back to gcc if cuda is activated
 if main_env['backend'] == "cuda":
     main_env['CXX'] = "g++"
 
 # OpenMP flags - compiler dependent
 omp_libs = {
     "g++": ["gomp"],
     "c++": ["gomp"],
     "clang++": [],
     "icpc": []
 }
 
 omp_flags = {
     "g++": ["-fopenmp"],
     "c++": ["-fopenmp"],
     "clang++": ["-fopenmp=libomp"],
     "icpc": ["-qopenmp"]
 }
 
 omp_lib = omp_libs[main_env['CXX']]
 omp_flag = omp_flags[main_env['CXX']]
 
 main_env.AppendUnique(LIBS=omp_lib)
 main_env.AppendUnique(LINKFLAGS=omp_flag)
 
 # Flags and options
 main_env.AppendUnique(CXXFLAGS=['-std=c++11',
                                 '-Wall',
                                 omp_flag])
 
 if main_env['backend'] == 'omp':
     main_env.AppendUnique(CPPDEFINES=["THRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_OMP"])
 elif main_env['backend'] == 'cuda':
     main_env.AppendUnique(CPPDEFINES=["THRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_CUDA"])
     main_env.AppendUnique(CPPDEFINES=['USE_CUDA'])
 
+if main_env['build_python']:
+    main_env.AppendUnique(CPPDEFINES=['USE_PYTHON'])
+
 # Adding compile flags defined in evironment
 if 'CXXFLAGS' in os.environ:
     main_env.AppendUnique(CXXFLAGS=Split(os.environ['CXXFLAGS']))
 
 if build_type == 'debug':
     main_env.AppendUnique(CPPDEFINES=['TAMAAS_DEBUG'])
 
 # Compilation flags
 cxxflags_dict = {
     "debug": Split("-g -O0"),
     "profiling": Split("-g -pg -O2 -fno-omit-frame-pointer"),
     "release": Split("-O3")
 }
 
 # Link flags for shared libs
 shlinkflags_dict = {
     "debug": Split(""),
     "profiling": ['-pg'],
     "release": []
 }
 
 if main_env['sanitizer'] != 'none':
     if main_env['backend'] == 'cuda':
         print("Sanitizers with cuda are not yet supported!")
         Exit(1)
     cxxflags_dict[build_type].append('-fsanitize={}'.format(main_env['sanitizer']))
     shlinkflags_dict[build_type].append('-fsanitize={}'.format(main_env['sanitizer']))
 
 main_env.AppendUnique(CXXFLAGS=cxxflags_dict[build_type])
 main_env.AppendUnique(SHLINKFLAGS=shlinkflags_dict[build_type])
 main_env.AppendUnique(LINKFLAGS=shlinkflags_dict[build_type])
 
 main_env['LIBPATH'] = [abspath(join(build_dir, 'src'))]
 main_env['RPATH'] = "$LIBPATH"
 
 if main_env['should_configure']:
     detect_dependencies(main_env)
 
 # Writing information file
 write_info_file("src/tamaas_info.cpp")
 
 # Saving the env file
 env_content = """export PYTHONPATH=$PYTHONPATH:{0}/python
 export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:{0}/src
 """
 
 
 def write_env_file(target, source, env):
     """Builder to write content to file"""
     with open(str(target[0]), 'w') as env_file:
         env_file.write(env_content.format(abspath(build_dir)))
 
 
 main_env['gen_print'] = gen_print
 env_file_env = main_env.Clone(PRINT_CMD_LINE_FUNC=gen_print("Writing", "cyan", main_env))
 # Need to have a command and manage tamaas_environement.sh as target because
 # the build directory does not always exist
 env_file_env.Command(join(build_dir, 'tamaas_environement.sh'), None, write_env_file)
 
 Export('main_env')
 
 
 # Building subdirs
 def subdir(dir):
     return SConscript(join(dir, 'SConscript'),
                       variant_dir=join(build_dir, dir),
                       duplicate=True)
 
 # Building Tamaas library
 subdir('src')
 
 # Building Tamaas extra components
 for dir in ['python', 'tests', 'doc']:
     if main_env['build_{}'.format(dir)] and not main_env.GetOption('help'):
         subdir(dir)
diff --git a/python/SConscript b/python/SConscript
index 558466b..a3ddd05 100644
--- a/python/SConscript
+++ b/python/SConscript
@@ -1,55 +1,34 @@
 from __future__ import print_function
 
 import os
 from os.path import abspath, basename
-from subprocess import check_output, STDOUT
 import re
 
 Import('main_env')
 
 # Pybind11 wrapper
 env_pybind = main_env.Clone(SHLIBPREFIX='')
-
-# Run code below to get versions from user preference and not scons
-versions_script = """
-from __future__ import print_function
-import distutils.sysconfig as dsys
-from numpy.distutils.misc_util import get_numpy_include_dirs as numpy_dirs
-
-print(dsys.get_python_inc())
-for d in numpy_dirs():
-    print(d)"""
-
-includes = check_output([main_env['py_exec'], "-c", versions_script],
-                        universal_newlines=True).split('\n')
-
-includes += [Dir('#third-party/pybind11/include')]
-
-# Extension of shared library for python
-extension = check_output([main_env['py_exec'] + '-config', "--extension-suffix"],
-                         universal_newlines=True).split('\n')[0]
-env_pybind.AppendUnique(CPPPATH=includes)
 env_pybind.Tool(pybind11)
 
 pybind_sources = Split("""
 tamaas_module.cpp
 wrap/core.cpp
 wrap/percolation.cpp
 wrap/surface.cpp
 wrap/model.cpp
 wrap/solvers.cpp
 wrap/bem.cpp
 """)
 
 # Building the pybind library
 env_pybind.SharedLibrary(
-    target='tamaas/_tamaas' + extension,
+    target='tamaas/_tamaas',
     source=pybind_sources,
     LIBS=['Tamaas'],
     RPATH=[abspath('../src')]
 )
 
 # Copying the __init__.py file with extra python classes
 copy_env = env_pybind.Clone(
     PRINT_CMD_LINE_FUNC=main_env['gen_print']("Copying", "red", main_env))
 copy_env.Command("tamaas/__init__.py", "#/python/tamaas.py", Copy("$TARGET", "$SOURCE"))
diff --git a/site_scons/site_init.py b/site_scons/site_init.py
index 7a33330..3536579 100644
--- a/site_scons/site_init.py
+++ b/site_scons/site_init.py
@@ -1,101 +1,123 @@
 from SCons.Script import *
 from os.path import *
+from subprocess import check_output, STDOUT
 
 import os
 
 
 def fftw(env):
     """A Tool to search for fftw headers and libraries"""
 
     if env.GetOption('clean'):
         return
     env.SetDefault(FFTW_VERSION='3')
     env.SetDefault(FFTW_LIBRARY_WISH=[])
 
     if 'FFTW_INCLUDE_DIR' in env:
         env.AppendUnique(CPPPATH=env['FFTW_INCLUDE_DIR'])
 
     if 'FFTW_LIBRARY_DIR' in env:
         env.AppendUnique(LIBPATH=env['FFTW_LIBRARY_DIR'])
 
     version = env['FFTW_VERSION']
     if version == "2":
         lib_names = {'main': 'fftw'}
         inc_names = ['fftw.h']
     else:
         lib_names = {'main': 'fftw3',
                      'thread': 'fftw3_threads',
                      'omp': 'fftw3_omp'}
         inc_names = ['fftw3.h']
 
     try:
         libs = [lib_names[i] for i in env['FFTW_LIBRARY_WISH']]
     except:
         raise SCons.Errors.StopError(
            'Incompatible wishlist {0} from version {1}'.format(
              env['FFTW_LIBRARY_WISH'], env['FFTW_VERSION']))
 
     env.AppendUnique(LIBS=libs)
     if version == "2":
         env.Append(LIBS='m')
 
     conf = Configure(env)
     if not conf.CheckLibWithHeader(libs, inc_names, 'c++'):
         raise SCons.Errors.StopError(
            'Failed to find libraries {0} or '
            'headers {1}.'.format(str(lib_names), str(inc_names)))
 
     env = conf.Finish()
 
 
 def boost(env):
     """A tool to confugure for boost"""
     if env.GetOption('clean'):
         return
 
     if 'BOOST_INCLUDE_DIR' in env:
         env.AppendUnique(CPPPATH=env['BOOST_INCLUDE_DIR'])
     conf = Configure(env)
     if not conf.CheckCXXHeader('boost/preprocessor/seq.hpp'):
         raise SCons.Errors.StopError(
             'Failed to find Boost header')
     env = conf.Finish()
 
 def pybind11(env):
     """A tool to configure pybind11"""
     if env.GetOption('clean'):
         return
 
+    # Run code below to get versions from user preference and not scons
+    versions_script = """
+from __future__ import print_function
+import distutils.sysconfig as dsys
+from numpy.distutils.misc_util import get_numpy_include_dirs as numpy_dirs
+
+print(dsys.get_python_inc())
+for d in numpy_dirs():
+    print(d)"""
+
+    includes = check_output([env['py_exec'], "-c", versions_script],
+                            universal_newlines=True).split('\n')
+
+    includes += [Dir('#third-party/pybind11/include')]
+
+    # Extension of shared library for python
+    extension = check_output([env['py_exec'] + '-config', "--extension-suffix"],
+                            universal_newlines=True).split('\n')[0]
+    env['SHLIBSUFFIX'] = extension + '.so'
+    env.AppendUnique(CPPPATH=includes)
+
     conf = Configure(env)
     if not conf.CheckCXXHeader('pybind11/pybind11.h'):
         raise SCons.Errors.StopError(
             'Failed to find pybind11 header\n' +
             "Run 'git submodule update --init --recursive third-party/pybind11'")
     env = conf.Finish()
 
 
 def gtest(env):
     """A tool to configure GoogleTest"""
     if env.GetOption('clean'):
         return
 
     conf = Configure(env)
     if not conf.CheckCXXHeader('gtest/gtest.h'):
         raise SCons.Errors.StopError(
             'Failed to find GoogleTest header\n' +
             "Run 'git submodule update --init --recursive third-party/googletest'")
     env = conf.Finish()
 
 
 def thrust(env):
     """A tool to confugure for thrust"""
     if env.GetOption('clean'):
         return
 
     if 'THRUST_INCLUDE_DIR' in env:
         env.AppendUnique(CPPPATH=env['THRUST_INCLUDE_DIR'])
     conf = Configure(env)
     if not conf.CheckCXXHeader('thrust/version.h'):
         raise SCons.Errors.StopError(
             'Failed to find Thrust library')
     env = conf.Finish()
diff --git a/src/bem/bem_fft_base.cpp b/src/bem/bem_fft_base.cpp
index 8de76d4..838ec4f 100644
--- a/src/bem/bem_fft_base.cpp
+++ b/src/bem/bem_fft_base.cpp
@@ -1,199 +1,200 @@
 /**
  * @file
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #include "bem_fft_base.hh"
 /* -------------------------------------------------------------------------- */
 #include "bem_meta_functional.hh"
 #include "elastic_energy_functional.hh"
 #include "fft_plan_manager.hh"
 #include "surface_statistics.hh"
+#include "types.hh"
 /* -------------------------------------------------------------------------- */
 #include <omp.h>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /* -------------------------------------------------------------------------- */
 
 BemFFTBase::BemFFTBase(Surface<Real>& p)
     : BemInterface(p), functional(new MetaFunctional(*this)),
       surface_tractions(p.size(), p.getL()),
       surface_displacements(p.size(), p.getL()),
       surface_spectral_input(p.size(), p.getL()),
       surface_spectral_output(p.size(), p.getL()),
       surface_spectral_influence_disp(p.size(), p.getL()),
       true_displacements(p.size(), p.getL()),
       old_displacements(p.size(), p.getL()), gap(p.size(), p.getL()),
       e_star(std::nan("")), tmp_coeff(1. / M_PI) {
   //    this->setNumberOfThreads(this->nthreads);
   max_iterations = 5000;
   dump_freq = 100;
   this->functional->addFunctionalTerm(new ElasticEnergyFunctional(*this));
   std::cerr << this << " is setting up the surfaces";
   std::cerr << &p << " has size " << p.size() << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 
 BemFFTBase::~BemFFTBase() { FFTPlanManager::get().clean(); }
 
 /* -------------------------------------------------------------------------- */
 
 const Surface<Real>& BemFFTBase::getTractions() const {
   return surface_tractions;
 }
 /* -------------------------------------------------------------------------- */
 const Surface<Real>& BemFFTBase::getDisplacements() const {
   return surface_displacements;
 }
 /* -------------------------------------------------------------------------- */
 const Surface<Real>& BemFFTBase::getSpectralInfluenceOverDisplacement() {
   return surface_spectral_influence_disp;
 }
 /* -------------------------------------------------------------------------- */
 void BemFFTBase::setEffectiveModulus(Real e_star) {
   this->e_star = e_star;
   this->computeSpectralInfluenceOverDisplacement();
 }
 /* -------------------------------------------------------------------------- */
 Real BemFFTBase::getEffectiveModulus() { return this->e_star; }
 /* -------------------------------------------------------------------------- */
 void BemFFTBase::computeSpectralInfluenceOverDisplacement() {
   // if (std::isnan(tmp_coeff)) SURFACE_FATAL("constant tmp_coeff is nan");
   // if (std::isnan(e_star))    SURFACE_FATAL("constant e_star is nan");
 
   const Real L = surface.getL();
 
   auto wavevectors = FFTransform<Real, 2>::computeFrequencies<true>(
       surface_spectral_influence_disp.sizes());
 
   auto inf_it = surface_spectral_influence_disp.begin(),
        inf_begin = surface_spectral_influence_disp.begin();
   auto end = wavevectors.end(2);
   //#pragma omp parallel for firstprivate(inf_it)
   for (auto it = wavevectors.begin(2); it < end; ++it) {
     inf_it = inf_begin;
     inf_it += it - wavevectors.begin(2);
     VectorProxy<Real> q_vec(&(*it), 2);
     q_vec *= 2 * M_PI / L;
     Real q = q_vec.l2norm();
 
     // Influence function
     *inf_it = 2. / (q * e_star);
   }
 
   surface_spectral_influence_disp(0) = 0.;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::computeDisplacementsFromTractions() {
   this->applyInfluenceFunctions(surface_tractions, surface_displacements);
 }
 
 /* -------------------------------------------------------------------------- */
 void BemFFTBase::computeTractionsFromDisplacements() {
   this->applyInverseInfluenceFunctions(surface_displacements,
                                        surface_tractions);
 }
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::applyInfluenceFunctions(Surface<Real>& input,
                                          Surface<Real>& output) {
   FFTransform<Real, 2>& plan1 =
       FFTPlanManager::get().createPlan(input, surface_spectral_output);
 
   plan1.forwardTransform();
 
   surface_spectral_output *= surface_spectral_influence_disp;
 
   FFTransform<Real, 2>& plan2 =
       FFTPlanManager::get().createPlan(output, surface_spectral_output);
 
   plan2.backwardTransform();
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::applyInverseInfluenceFunctions(Surface<Real>& input,
                                                 Surface<Real>& output) {
 
   FFTransform<Real, 2>& plan1 =
       FFTPlanManager::get().createPlan(input, surface_spectral_output);
 
   plan1.forwardTransform();
 
   surface_spectral_output /= surface_spectral_influence_disp;
   surface_spectral_output(0) = 0;
 
   FFTransform<Real, 2>& plan2 =
       FFTPlanManager::get().createPlan(output, surface_spectral_output);
 
   plan2.backwardTransform();
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::computeGaps() {
   UInt size = surface.size();
 
 #pragma omp parallel for
   for (UInt i = 0; i < size * size; i++) {
     gap(i) = true_displacements(i) - surface(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::computeTrueDisplacements() {
   this->true_displacements = this->surface_displacements;
 
   UInt n = surface.size();
   UInt size = n * n;
 
   Real shift = 1e300;
 
   for (UInt i = 0; i < size; ++i) {
     if (surface_displacements(i) - shift - surface(i) < 0.) {
       shift = surface_displacements(i) - surface(i);
     }
   }
 
   for (UInt i = 0; i < size; ++i) {
     true_displacements(i) = surface_displacements(i) - shift;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::addFunctional(Functional* new_funct) {
   this->functional->addFunctionalTerm(new_funct);
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
diff --git a/src/core/fftransform.hh b/src/core/fftransform.hh
index dfd2ba4..a22890c 100644
--- a/src/core/fftransform.hh
+++ b/src/core/fftransform.hh
@@ -1,121 +1,115 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #ifndef FFTRANSFORM_HH
 #define FFTRANSFORM_HH
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "grid_hermitian.hh"
 #include "tamaas.hh"
-#include "types.hh"
+#include "static_types.hh"
 #include <vector>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 /**
  * @brief Generic class for FFT interface wrapping
  */
 template <typename T, UInt dim>
 class FFTransform {
 public:
   /// Constructor
   FFTransform(Grid<T, dim>& real, GridHermitian<T, dim>& spectral);
 
   /// Destructor
   virtual ~FFTransform();
 
 public:
   /// Perform FFT
   virtual void forwardTransform() = 0;
 
   /// Perform IFFT
   virtual void backwardTransform() = 0;
 
   /// Normalize the real surface after backward transform
   virtual void normalize();
 
 public:
   /// Fill a grid with modevectors: boolean if grid has hermitian dimensions
   template <bool hermitian>
   static Grid<T, dim> computeFrequencies(const std::array<UInt, dim>& sizes);
 
 protected:
   Grid<T, dim>& real;
   GridHermitian<T, dim>& spectral;
 };
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <bool hermitian>
 Grid<T, dim>
 FFTransform<T, dim>::computeFrequencies(const std::array<UInt, dim>& sizes) {
   // If hermitian is true, we suppose the dimensions of freq are
   // reduced based on hermitian symetry and that it has dim components
 
   auto& n = sizes;
   Grid<T, dim> freq(n, dim);
-  std::array<UInt, dim> tuple{{0}};
-
-#pragma omp parallel for firstprivate(tuple)
-  for (auto it = freq.begin(dim); it < freq.end(dim); ++it) {
-    VectorProxy<T> wavevector(&(*it), dim);
-    UInt index = freq.getNbPoints() - (freq.end(dim) - it);
-
-    /// Computing tuple from index
-    for (UInt d = dim - 1; d > 0; d--) {
-      tuple[d] = index % n[d];
-      index -= tuple[d];
-      index /= n[d];
-    }
-
-    tuple[0] = index;
-
-    UInt dmax = dim;
-
-    if (hermitian) {
-      dmax = dim - 1;
-      wavevector(dim - 1) = tuple[dim - 1];
-    }
-
-    for (UInt d = 0; d < dmax; d++) {
-      T td = tuple[d];
-      T nd = n[d];
-      T q = (tuple[d] <= n[d] / 2) ? td : td - nd;
-      wavevector(d) = q;
-    }
-  }
+  constexpr UInt dmax = dim - static_cast<UInt>(hermitian);
+
+  Loop::stridedLoop(
+      [dmax, &n](StaticVector<T, dim>&& wavevector, UInt index) {
+        std::array<UInt, dim> tuple{{0}};
+        /// Computing tuple from index
+        for (Int d = dim - 1; d >= 0; d--) {
+          tuple[d] = index % n[d];
+          index -= tuple[d];
+          index /= n[d];
+        }
+
+        if (hermitian)
+          wavevector(dim - 1) = tuple[dim - 1];
+
+        for (UInt d = 0; d < dmax; d++) {
+	  // Type conversion
+          T td = tuple[d];
+          T nd = n[d];
+          T q = (tuple[d] < n[d] / 2) ? td : td - nd;
+          wavevector(d) = q;
+        }
+      },
+      freq, Loop::range(freq.getNbPoints()));
 
   return freq;
 }
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 #endif  // FFTRANSFORM_HH
diff --git a/src/core/grid_base.hh b/src/core/grid_base.hh
index 6b91019..0d17874 100644
--- a/src/core/grid_base.hh
+++ b/src/core/grid_base.hh
@@ -1,317 +1,298 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __GRID_BASE_HH__
 #define __GRID_BASE_HH__
 /* -------------------------------------------------------------------------- */
 #include "array.hh"
 #include "iterator.hh"
 #include "loop.hh"
 #include "tamaas.hh"
 #include <cstddef>
 #include <limits>
 #include <utility>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
-/// Namespace for view index classes
-namespace view {
-/// Index class
-struct index {
-  index(Int i = -1) : i(i) {}
-  Int i;
-};
-
-/// Blocked index class
-struct blocked : public index {
-  blocked(Int i) : index(i) {}
-};
-
-/// Free index class
-struct free : public index {
-  free() : index(-1) {}
-};
-}  // namespace view
-
 /**
  * @brief Dimension agnostic grid with components stored per points
  */
 template <typename T>
 class GridBase {
 public:
   /// Constructor by default
   GridBase() = default;
   /// Copy constructor
   GridBase(const GridBase& o)
       : data(o.data), nb_components(o.nb_components), offset(o.offset) {}
   /// Move constructor (transfers data ownership)
   GridBase(GridBase&& o)
       : data(std::move(o.data)),
         nb_components(std::exchange(o.nb_components, 1)),
         offset(std::exchange(o.offset, 0)) {}
 
   /// Destructor
   virtual ~GridBase() = default;
 
   /* ------------------------------------------------------------------------ */
   /* Iterator class */
   /* ------------------------------------------------------------------------ */
 public:
   using iterator = iterator_::iterator<T>;
   using const_iterator = iterator_::iterator<const T>;
 
 public:
   /// Get grid dimension
   virtual UInt getDimension() const { return 1; }
   /// Get internal data pointer (const)
   const T* getInternalData() const { return this->data.data(); }
   /// Get internal data pointer (non-const)
   T* getInternalData() { return this->data.data(); }
   /// Get number of components
   UInt getNbComponents() const { return nb_components; }
   /// Set number of components
   void setNbComponents(UInt n) { nb_components = n; }
   /// Get offset
   UInt getOffset() const { return offset; }
   /// Get total size
   virtual UInt dataSize() const { return this->data.size(); }
   /// Get number of points
   UInt getNbPoints() const {
     return this->dataSize() / this->getNbComponents();
   }
   /// Set components
   void uniformSetComponents(const GridBase<T>& vec);
 
   /* ------------------------------------------------------------------------ */
   /* Iterator methods */
   /* ------------------------------------------------------------------------ */
   virtual iterator begin(UInt n = 1) {
     return iterator(this->getInternalData() + this->offset, 0, n);
   }
 
   virtual iterator end(UInt n = 1) {
     return iterator(this->getInternalData() + this->offset, this->dataSize(),
                     n);
   }
 
   virtual const_iterator begin(UInt n = 1) const {
     return const_iterator(this->getInternalData() + this->offset, 0, n);
   }
 
   virtual const_iterator end(UInt n = 1) const {
     return const_iterator(this->getInternalData() + this->offset,
                           this->dataSize(), n);
   }
 
   /* ------------------------------------------------------------------------ */
   /* Operators */
   /* ------------------------------------------------------------------------ */
   inline T& operator()(UInt i) { return this->data[i]; }
 
   inline const T& operator()(UInt i) const { return this->data[i]; }
 
 #define VEC_OPERATOR(op)                                                       \
   template <typename T1>                                                       \
   GridBase& operator op(const GridBase<T1>& other)
 
   VEC_OPERATOR(+=);
   VEC_OPERATOR(*=);
   VEC_OPERATOR(-=);
   VEC_OPERATOR(/=);
 #undef VEC_OPERATOR
 
   template <typename T1>
   T dot(const GridBase<T1>& other) const;
 
 #define SCALAR_OPERATOR(op) GridBase& operator op(const T& e)
 
   SCALAR_OPERATOR(+=);
   SCALAR_OPERATOR(*=);
   SCALAR_OPERATOR(-=);
   SCALAR_OPERATOR(/=);
   SCALAR_OPERATOR(=);
 #undef SCALAR_OPERATOR
 
   GridBase& operator=(const GridBase& o);
 
   /* ------------------------------------------------------------------------ */
   /* Min/max */
   /* ------------------------------------------------------------------------ */
   T min() const;
   T max() const;
   T sum() const;
   T mean() const { return sum() / static_cast<T>(dataSize()); }
   T var() const;
 
   /* ------------------------------------------------------------------------ */
   /* Move/Copy */
   /* ------------------------------------------------------------------------ */
 
 public:
   template <typename T1>
   void copy(const GridBase<T1>& other) {
     data = other.data;
     nb_components = other.nb_components;
     offset = other.offset;
   }
 
   template <typename T1>
   void move(GridBase<T1>&& other) {
     data = std::move(other.data);
     nb_components = std::exchange(other.nb_components, 1);
     offset = std::exchange(other.offset, 0);
   }
 
   GridBase& wrap(GridBase& o) {
     data.wrap(o.data);
     nb_components = o.nb_components;
     offset = o.offset;
     return *this;
   }
 
   GridBase& wrap(const GridBase& o) {
     data.wrap(o.data);
     nb_components = o.nb_components;
     offset = o.offset;
     return *this;
   }
 
 protected:
   Array<T> data;
   UInt nb_components = 1;
   UInt offset = 0;
 };
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void GridBase<T>::uniformSetComponents(const GridBase<T>& vec) {
   if (!(vec.dataSize() == this->nb_components))
     TAMAAS_EXCEPTION("Cannot set grid field with values of vector");
 
   auto begin_it(begin());
   auto end_it(end());  // silencing nvcc warnings
 
 #pragma omp parallel for
   for (auto it = begin(); it < end_it; ++it) {
     UInt i = it - begin_it;
     *it = vec.data[i % this->nb_components];
   }
 }
 
   /* --------------------------------------------------------------------------
    */
 
 #define SCALAR_OPERATOR_IMPL(op)                                               \
   template <typename T>                                                        \
   inline GridBase<T>& GridBase<T>::operator op(const T& e) {                   \
     Loop::loop([e] CUDA_LAMBDA(T& val) { val op e; }, *this);                  \
     return *this;                                                              \
   }
 
 SCALAR_OPERATOR_IMPL(+=);
 SCALAR_OPERATOR_IMPL(*=);
 SCALAR_OPERATOR_IMPL(-=);
 SCALAR_OPERATOR_IMPL(/=);
 SCALAR_OPERATOR_IMPL(=);
 
 #undef SCALAR_OPERATOR_IMPL
 
 /* -------------------------------------------------------------------------- */
 
 #define VEC_OPERATOR_IMPL(op)                                                  \
   template <typename T>                                                        \
   template <typename T1>                                                       \
   inline GridBase<T>& GridBase<T>::operator op(const GridBase<T1>& other) {    \
     TAMAAS_ASSERT(other.dataSize() == this->dataSize(),                        \
                   "surface size does not match");                              \
     Loop::loop([] CUDA_LAMBDA(T& x, const T1& y) { x op y; }, *this, other);   \
     return *this;                                                              \
   }
 
 VEC_OPERATOR_IMPL(+=);
 VEC_OPERATOR_IMPL(-=);
 VEC_OPERATOR_IMPL(*=);
 VEC_OPERATOR_IMPL(/=);
 #undef VEC_OPERATOR
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 GridBase<T>& GridBase<T>::operator=(const GridBase& o) {
   this->copy(o);
   return *this;
 }
 
 /* -------------------------------------------------------------------------- */
 
 using op = operation;
 
 template <typename T>
 inline T GridBase<T>::min() const {
   // const auto id = [] CUDA_LAMBDA (const T&x){return x;};
   return Loop::reduce<op::min>([] CUDA_LAMBDA(const T& x) { return x; }, *this);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::max() const {
   return Loop::reduce<op::max>([] CUDA_LAMBDA(const T& x) { return x; }, *this);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::sum() const {
   return Loop::reduce<op::plus>([] CUDA_LAMBDA(const T& x) { return x; },
                                 *this);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::var() const {
   const T mu = mean();
   const T var = Loop::reduce<op::plus>(
       [mu] CUDA_LAMBDA(const T& x) { return (x - mu) * (x - mu); }, *this);
   return var / (dataSize() - 1);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 template <typename T1>
 T GridBase<T>::dot(const GridBase<T1>& other) const {
   return Loop::reduce<op::plus>(
       [] CUDA_LAMBDA(const T& x, const T1& y) { return x * y; }, *this, other);
 }
 
 __END_TAMAAS__
 
 #endif  // __GRID_BASE_HH__
diff --git a/src/core/grid_tmpl.hh b/src/core/grid_tmpl.hh
index abaf75e..2f17992 100644
--- a/src/core/grid_tmpl.hh
+++ b/src/core/grid_tmpl.hh
@@ -1,173 +1,170 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __GRID_TMPL_HH__
 #define __GRID_TMPL_HH__
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "tamaas.hh"
 
 __BEGIN_TAMAAS__
 
 template <typename T, UInt dim>
 template <typename T1>
 void Grid<T, dim>::init(const T1& n, UInt nb_components) {
   std::copy(n.begin(), n.end(), this->n.begin());
   this->nb_components = nb_components;
   this->resize(this->n);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 inline UInt Grid<T, dim>::computeSize() const {
   UInt size = 1;
   for (UInt i = 0; i < dim; i++)
     size *= n[i];
   size *= this->nb_components;
   return size;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index_pos, UInt index,
                                        T1... rest) const {
   offset += index * strides[index_pos];
   return unpackOffset(offset, index_pos + 1, rest...);  // tail-rec bb
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index_pos,
                                        UInt index) const {
   return offset + index * strides[index_pos];
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline T& Grid<T, dim>::operator()(T1... args) {
   /// Checking access integrity
   constexpr UInt nargs = sizeof...(T1);
   static_assert(nargs == dim + 1 || nargs == 1 || nargs == dim,
                 "number of arguments in operator() does not match dimension");
 
   constexpr UInt start = (nargs == 1) ? dim : 0;
   UInt offset = unpackOffset(this->offset, start, args...);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline const T& Grid<T, dim>::operator()(T1... args) const {
   /// Checking access integrity
   constexpr UInt nargs = sizeof...(T1);
   static_assert(nargs == dim + 1 || nargs == 1 || nargs == dim,
                 "number of arguments in operator() does not match dimension");
 
   constexpr UInt start = (nargs == 1) ? dim : 0;
   UInt offset = unpackOffset(this->offset, start, args...);
   return this->data[offset];
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 inline UInt Grid<T, dim>::computeOffset(std::array<UInt, dim + 1> tuple) const {
   UInt offset = this->offset;
-  for (UInt i = 0; i < dim + 1; i++) {
-    offset += tuple[i] * strides[i];
-  }
-
-  return offset;
+  return std::inner_product(tuple.begin(), tuple.end(), strides.begin(),
+                            offset);
 }
 
 template <typename T, UInt dim>
 inline T& Grid<T, dim>::operator()(std::array<UInt, dim + 1> tuple) {
   UInt offset = computeOffset(tuple);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 inline const T& Grid<T, dim>::
 operator()(std::array<UInt, dim + 1> tuple) const {
   UInt offset = computeOffset(tuple);
   return this->data[offset];
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim>& Grid<T, dim>::operator=(const Grid<T, dim>& other) {
   this->copy(other);
   return *this;
 }
 
 template <typename T, UInt dim>
 Grid<T, dim>& Grid<T, dim>::operator=(Grid<T, dim>& other) {
   this->copy(other);
   return *this;
 }
 
 template <typename T, UInt dim>
 Grid<T, dim>& Grid<T, dim>::operator=(Grid<T, dim>&& other) {
   this->move(std::forward<Grid<T, dim>>(other));
   return *this;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename T1>
 void Grid<T, dim>::copy(const Grid<T1, dim>& other) {
   GridBase<T>::copy(other);
   this->n = other.n;
   this->strides = other.strides;
 }
 
 template <typename T, UInt dim>
 template <typename T1>
 void Grid<T, dim>::move(Grid<T1, dim>&& other) {
   GridBase<T>::move(std::forward<Grid<T1, dim>>(other));
   this->n = std::move(other.n);
   this->strides = std::move(other.strides);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Stream operator                                                            */
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 inline std::ostream& operator<<(std::ostream& stream,
                                 const Grid<T, dim>& _this) {
   _this.printself(stream);
   return stream;
 }
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 #endif  // __GRID_TMPL_HH__
diff --git a/src/core/grid_view.hh b/src/core/grid_view.hh
index 9bbe3d7..966c300 100644
--- a/src/core/grid_view.hh
+++ b/src/core/grid_view.hh
@@ -1,244 +1,113 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __GRID_VIEW_HH__
 #define __GRID_VIEW_HH__
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "tamaas.hh"
 #include <vector>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 /**
  * @brief View type on grid
- * /!\ iterators do not work on this type of grid
- *
+ * This is a view on a *contiguous* chunk of data defined by a grid
  */
 template <template <typename, UInt> class Base, typename T, UInt base_dim,
           UInt dim>
 class GridView : public Base<T, dim> {
 public:
   /// Constructor
-  GridView(GridBase<T>& grid_base, const std::vector<view::index>& multi_index);
+  GridView(GridBase<T>& grid_base, const std::vector<UInt>& multi_index);
   /// Move constructor
   GridView(GridView&& o)
       : Base<T, dim>(std::forward<Base<T, dim>>(o)),
         grid(std::exchange(o.grid, nullptr)), indexes(std::move(o.indexes)) {}
 
   /// Destructor
   virtual ~GridView() = default;
 
-public:
-  /// Compute strides and offset
-  void computeStrides() override;
-  /// Return view size
-  UInt dataSize() const override { return this->computeSize(); }
-  /// Delete resize
-  void resize(const std::array<UInt, dim>&) = delete;
-  /// Get a view only for the last component
-  GridView lastComponentView() const;
-
-public:
-  using iterator = typename GridBase<T>::iterator;
-  using const_iterator = typename GridBase<T>::const_iterator;
-
-  iterator begin(UInt n[[gnu::unused]] = 1) override {
-    TAMAAS_EXCEPTION("iterators on views do not work as expected");
-  }
-
-  iterator end(UInt n[[gnu::unused]] = 1) override {
-    TAMAAS_EXCEPTION("iterators on views do not work as expected");
-  }
-
-  const_iterator begin(UInt n[[gnu::unused]] = 1) const override {
-    TAMAAS_EXCEPTION("const_iterators on views do not work as expected");
-  }
-
-  const_iterator end(UInt n[[gnu::unused]] = 1) const override {
-    TAMAAS_EXCEPTION("const_iterators on views do not work as expected");
-  }
-
-public:
-  /// Access operator (1D access is extremely slow)
-  template <typename... T1>
-  T& operator()(T1... args);
-  /// Const access operator (1D access is extremely slow)
-  template <typename... T1>
-  const T& operator()(T1... args) const;
-
 protected:
   Base<T, base_dim>* grid;
-  std::vector<view::index> indexes;
+  std::vector<UInt> indexes;
 };
 
 /* -------------------------------------------------------------------------- */
-/* Template implementation                                                    */
-/* -------------------------------------------------------------------------- */
+template <template <typename, UInt> class Base, typename T, UInt base_dim,
+          typename... Args>
+GridView<Base, T, base_dim, base_dim - sizeof...(Args)>
+make_view(Base<T, base_dim>& base, Args... indices) {
+  std::vector<UInt> multi_index = {static_cast<UInt>(indices)...};
+  return GridView<Base, T, base_dim, base_dim - sizeof...(Args)>(base,
+                                                                 multi_index);
+}
+
+/* --------------------------------------------------------------------------
+   */
+/* Template implementation */
+/* --------------------------------------------------------------------------
+   */
 
 template <template <typename, UInt> class Base, typename T, UInt base_dim,
           UInt dim>
-GridView<Base, T, base_dim, dim>::GridView(
-    GridBase<T>& grid_base, const std::vector<view::index>& multi_index)
+GridView<Base, T, base_dim, dim>::GridView(GridBase<T>& grid_base,
+                                           const std::vector<UInt>& multi_index)
     : Base<T, dim>(), grid(nullptr), indexes(multi_index) {
   static_assert(base_dim >= dim,
                 "view dimension must be >= than the base class");
 
   // Checking view type
   grid = dynamic_cast<Base<T, base_dim>*>(&grid_base);
 
   if (!grid)
     TAMAAS_EXCEPTION("given base type is incompatible with view");
 
-  // Checking view integrity
-  if (multi_index.size() == base_dim) {
-    this->nb_components = grid->getNbComponents();
-  } else if (multi_index.size() == base_dim + 1) {
-    this->nb_components = 1;
-  } else {
-    TAMAAS_EXCEPTION("Multi-index has wrong number of components");
-  }
-
-  UInt view_dim = 0;
-  for (auto& it : multi_index) {
-    if (it.i == -1)
-      view_dim += 1;
-  }
-
-  if (view_dim != dim)
-    TAMAAS_EXCEPTION("Multi-index has the wrong number of free components");
-
-  if (indexes.size() == base_dim && this->nb_components == 1)
-    indexes.push_back(view::blocked(0));
-
-  this->data.wrapMemory(grid->getInternalData(), grid->dataSize());
-  computeStrides();
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <template <typename, UInt> class Base, typename T, UInt base_dim,
-          UInt dim>
-void GridView<Base, T, base_dim, dim>::computeStrides() {
-  /// Copying sizes with number of components in last position
-  auto stride_it = this->strides.begin();
-  auto size_it = this->n.begin();
-  auto grid_stride_it = grid->getStrides().begin();
-  auto grid_size_it = grid->sizes().begin();
-
-  for (auto& index : indexes) {
-    if (index.i == -1) {
-      *stride_it = *grid_stride_it;
-      *size_it = *grid_size_it;
-      ++stride_it, ++size_it;
-    } else
-      this->offset += index.i * *grid_stride_it;
-    ++grid_stride_it, ++grid_size_it;
-  }
-}
+  constexpr UInt dim_offset = base_dim - dim;
 
-/* -------------------------------------------------------------------------- */
+  if (dim_offset != multi_index.size())
+    TAMAAS_EXCEPTION("Number of blocked indices does not match view dimension");
 
-template <template <typename, UInt> class Base, typename T, UInt base_dim,
-          UInt dim>
-template <typename... T1>
-T& GridView<Base, T, base_dim, dim>::operator()(T1... args) {
-  constexpr UInt nargs = sizeof...(args);
-
-  static_assert(nargs == 1 || nargs == dim || nargs == dim + 1,
-                "number of arguments in operator() does not match dimension");
-
-  if (nargs == 1) {
-    UInt index = UInt(args...);
-    std::array<UInt, dim + 1> tuple = {0};
-
-    tuple[dim] = index % this->nb_components;
-    index -= tuple[dim];
-    index /= this->nb_components;
-
-    /// Computing tuple from index
-    for (UInt d = dim - 1; d > 0; d--) {
-      tuple[d] = index % this->n[d];
-      index -= tuple[d];
-      index /= this->n[d];
-    }
-
-    tuple[0] = index;
-    return Base<T, dim>::operator()(tuple);
-  } else {
-    return Base<T, dim>::operator()(args...);
-  }
-}
+  std::copy(grid->sizes().begin() + dim_offset, grid->sizes().end(),
+            this->n.begin());
+  std::copy(grid->getStrides().begin() + dim_offset, grid->getStrides().end(),
+            this->strides.begin());
 
-/* -------------------------------------------------------------------------- */
+  this->offset =
+      std::inner_product(multi_index.begin(), multi_index.end(),
+                         grid->getStrides().begin(), grid->getOffset());
 
-template <template <typename, UInt> class Base, typename T, UInt base_dim,
-          UInt dim>
-template <typename... T1>
-const T& GridView<Base, T, base_dim, dim>::operator()(T1... args) const {
-  constexpr UInt nargs = sizeof...(args);
-
-  static_assert(nargs == 1 || nargs == dim || nargs == dim + 1,
-                "number of arguments in operator() does not match dimension");
-
-  if (nargs == 1) {
-    UInt index = UInt(args...);
-    std::array<UInt, dim + 1> tuple = {0};
-
-    tuple[dim] = index % this->nb_components;
-    index -= tuple[dim];
-    index /= this->nb_components;
-
-    /// Computing tuple from index
-    for (UInt d = dim - 1; d > 0; d--) {
-      tuple[d] = index % this->n[d];
-      index -= tuple[d];
-      index /= this->n[d];
-    }
-
-    tuple[0] = index;
-    return Base<T, dim>::operator()(tuple);
-  } else {
-    return Base<T, dim>::operator()(args...);
-  }
+  this->nb_components = grid->getNbComponents();
+  this->data.wrapMemory(grid->getInternalData(), this->computeSize());
 }
 
 /* -------------------------------------------------------------------------- */
 
-template <template <typename, UInt> class Base, typename T, UInt base_dim,
-          UInt dim>
-GridView<Base, T, base_dim, dim>
-GridView<Base, T, base_dim, dim>::lastComponentView() const {
-  GridView<Base, T, base_dim, dim> tmp(*grid, indexes);
-  tmp.indexes.back() = view::blocked(this->nb_components - 1);
-  tmp.computeStrides();
-  return tmp;
-}
 
 __END_TAMAAS__
 
 #endif  // __GRID_VIEW_HH__
diff --git a/src/core/loop.hh b/src/core/loop.hh
index 28f2037..68309da 100644
--- a/src/core/loop.hh
+++ b/src/core/loop.hh
@@ -1,177 +1,183 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __LOOP_HH__
 #define __LOOP_HH__
 /* -------------------------------------------------------------------------- */
 #include "loops/apply.hh"
 #include "loops/loop_utils.hh"
 #include "tamaas.hh"
 #include <thrust/for_each.h>
 #include <thrust/iterator/counting_iterator.h>
 #include <thrust/iterator/zip_iterator.h>
 #include <thrust/transform_reduce.h>
 #include <thrust/tuple.h>
 #include <type_traits>
 
 __BEGIN_TAMAAS__
 
 /**
  * @brief Singleton class for automated loops using lambdas
 
  * This class is sweet candy :) It provides abstraction of the paralelism
  * paradigm used in loops and allows simple and less erro-prone loop syntax,
  * with minimum boiler plate. I love it <3
  */
 class Loop {
 public:
   /// Backends enumeration
   enum backend {
     omp,   ///< [OpenMP](http://www.openmp.org/specifications/) backend
     cuda,  ///< [Cuda](http://docs.nvidia.com/cuda/index.html) backend
   };
 
   /// Helper class to count iterations within lambda-loop
   template <typename T>
   class arange {
   public:
     using it_type = thrust::counting_iterator<T>;
     arange(T size) : range_size(size) {}
-    it_type begin() const { return it_type(T(0)); }
-    it_type end() const { return it_type(range_size); }
+    it_type begin(UInt i = 0) const { return it_type(T(0)); }
+    it_type end(UInt i = 0) const { return it_type(range_size); }
+    UInt getNbComponents() const { return 1; }
 
   private:
     T range_size;
   };
 
+  template <typename T>
+  static arange<T> range(T size) {
+    return arange<T>(size);
+  }
+
   /// Loop functor over any number of grids
   template <typename Functor, typename... Grids>
   static void loop(Functor&& func, Grids&&... containers);
 
   /// Strided loop over any number of grids
   template <typename Functor, typename... Grids>
   static void stridedLoop(Functor&& func, Grids&&... containers);
 
   /// Reduce over any number of grids
   template <operation op, typename Functor, typename... Grids>
   static auto reduce(Functor&& func, Grids&&... containers)
       -> decltype(func(containers(0)...));
 
   /// Strided reduce over any number of grids
   template <operation op, typename Functor, typename... Grids>
   static auto stridedReduce(Functor&& func, Grids&&... containers)
       -> decltype(func(containers(0)...));
 
   /// Constructor
   Loop() = delete;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Template implementation                                                    */
 /* -------------------------------------------------------------------------- */
 
 template <typename Functor, typename... Grids>
 void Loop::loop(Functor&& func, Grids&&... containers) {
   auto begin =
       thrust::make_zip_iterator(thrust::make_tuple(containers.begin()...));
   auto end = thrust::make_zip_iterator(thrust::make_tuple(containers.end()...));
 
   thrust::for_each(begin, end, detail::ApplyFunctor<Functor>(func));
 
 #ifdef USE_CUDA
   cudaDeviceSynchronize();
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename Functor, typename... Grids>
 void Loop::stridedLoop(Functor&& func, Grids&&... containers) {
   auto begin = thrust::make_zip_iterator(
       thrust::make_tuple(containers.begin(containers.getNbComponents())...));
   auto end = thrust::make_zip_iterator(
       thrust::make_tuple(containers.end(containers.getNbComponents())...));
 
   thrust::for_each(begin, end, detail::ApplyFunctor<Functor>(func));
 #ifdef USE_CUDA
   cudaDeviceSynchronize();
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <operation op, typename Functor, typename... Grids>
 auto Loop::reduce(Functor&& func, Grids&&... containers)
     -> decltype(func(containers(0)...)) {
   auto begin =
       thrust::make_zip_iterator(thrust::make_tuple(containers.begin()...));
   auto end = thrust::make_zip_iterator(thrust::make_tuple(containers.end()...));
 
   using reduce_type = decltype(func(containers(0)...));
   using apply_type = detail::ApplyFunctor<Functor, reduce_type>;
 
   auto red_helper = detail::reduction_helper<op, apply_type>(apply_type(func));
   auto result = thrust::reduce(
       begin, end, red_helper.template init<reduce_type>(), red_helper);
 
 #ifdef USE_CUDA
   cudaDeviceSynchronize();
 #endif
   return result;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <operation op, typename Functor, typename... Grids>
 auto Loop::stridedReduce(Functor&& func, Grids&&... containers)
     -> decltype(func(containers(0)...)) {
   auto begin = thrust::make_zip_iterator(
       thrust::make_tuple(containers.begin(containers.getNbComponents())...));
   auto end = thrust::make_zip_iterator(
       thrust::make_tuple(containers.end(containers.getNbComponents())...));
 
   using reduce_type = decltype(func(containers(0)...));
   using apply_type = detail::ApplyFunctor<Functor, reduce_type>;
 
   auto red_helper = detail::reduction_helper<op, apply_type>(apply_type(func));
 
   auto result = thrust::reduce(
       begin, end, red_helper.template init<reduce_type>(), red_helper);
 
 #ifdef USE_CUDA
   cudaDeviceSynchronize();
 #endif
   return result;
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
 
 #undef EXEC_CASE_MACRO
 #undef REDUCE_CASE_MACRO
 
 #endif  // __LOOP_HH__
diff --git a/src/model/model_type.cpp b/src/model/model_type.cpp
index 5fc2163..65861e4 100644
--- a/src/model/model_type.cpp
+++ b/src/model/model_type.cpp
@@ -1,47 +1,39 @@
 /**
  * @file
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016-2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "model_type.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
-std::vector<view::index> model_type_traits<model_type::basic_1d>::indices = {
-    view::free()};
-std::vector<view::index> model_type_traits<model_type::basic_2d>::indices = {
-    view::free(), view::free()};
-
-std::vector<view::index> model_type_traits<model_type::surface_1d>::indices = {
-    view::free()};
-std::vector<view::index> model_type_traits<model_type::surface_2d>::indices = {
-    view::free(), view::free()};
-
-std::vector<view::index> model_type_traits<model_type::volume_1d>::indices = {
-    view::blocked(0), view::free()};
-std::vector<view::index> model_type_traits<model_type::volume_2d>::indices = {
-    view::blocked(0), view::free(), view::free()};
+std::vector<UInt> model_type_traits<model_type::basic_1d>::indices = {};
+std::vector<UInt> model_type_traits<model_type::basic_2d>::indices = {};
+std::vector<UInt> model_type_traits<model_type::surface_1d>::indices = {};
+std::vector<UInt> model_type_traits<model_type::surface_2d>::indices = {};
+std::vector<UInt> model_type_traits<model_type::volume_1d>::indices = {0};
+std::vector<UInt> model_type_traits<model_type::volume_2d>::indices = {0};
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
diff --git a/src/model/model_type.hh b/src/model/model_type.hh
index 42b7bfc..bb4adcc 100644
--- a/src/model/model_type.hh
+++ b/src/model/model_type.hh
@@ -1,100 +1,100 @@
 /**
  * @file
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016-2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __MODEL_TYPE_HH__
 #define __MODEL_TYPE_HH__
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "grid_base.hh"
 #include "tamaas.hh"
 #include <boost/preprocessor/seq.hpp>
 #include <memory>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 
 /// Types for grid dimensions and number of components
 enum class model_type {
   basic_1d,    ///< one component line
   basic_2d,    ///< one component surface
   surface_1d,  ///< two components line
   surface_2d,  ///< three components surface
   volume_1d,   ///< two components volume
   volume_2d    ///< three components volume
 };
 
 /// Trait class to store physical dimension of domain, of boundary and
 /// number of components
 template <model_type type>
 struct model_type_traits {};
 
 #define MODEL_TYPE_TRAITS_MACRO(type, dim, comp, bdim)                         \
   template <>                                                                  \
   struct model_type_traits<model_type::type> {                                 \
     static constexpr UInt dimension = dim;                                     \
     static constexpr UInt components = comp;                                   \
     static constexpr UInt boundary_dimension = bdim;                           \
-    static std::vector<view::index> indices;                                   \
+    static std::vector<UInt> indices;                                          \
   }
 
 MODEL_TYPE_TRAITS_MACRO(basic_1d, 1, 1, 1);
 MODEL_TYPE_TRAITS_MACRO(basic_2d, 2, 1, 2);
 
 MODEL_TYPE_TRAITS_MACRO(surface_1d, 1, 2, 1);
 MODEL_TYPE_TRAITS_MACRO(surface_2d, 2, 3, 2);
 
 MODEL_TYPE_TRAITS_MACRO(volume_1d, 2, 2, 1);
 MODEL_TYPE_TRAITS_MACRO(volume_2d, 3, 3, 2);
 
 #undef MODEL_TYPE_TRAITS_MACRO
 
 #define TAMAAS_MODEL_TYPES                                                     \
   (model_type::basic_1d)(model_type::basic_2d)(model_type::surface_1d)(        \
       model_type::surface_2d)(model_type::volume_1d)(model_type::volume_2d)
 
 #define ALLOC_GRID_CASE_MACRO(a, b, type)                                      \
   case type: {                                                                 \
     constexpr UInt dim = (boundary)                                            \
                              ? model_type_traits<type>::boundary_dimension     \
                              : model_type_traits<type>::dimension;             \
     ptr = new Grid<T, dim>(n, nc);                                             \
     break;                                                                     \
   }
 
 /// Helper function for grid allocation
 template <bool boundary, typename T>
 GridBase<T>* allocateGrid(model_type type, const std::vector<UInt> n, UInt nc) {
   GridBase<T>* ptr = nullptr;
   switch (type) {
     BOOST_PP_SEQ_FOR_EACH(ALLOC_GRID_CASE_MACRO, ~, TAMAAS_MODEL_TYPES);
   }
 
   return ptr;
 }
 
 #undef ALLOC_GRID_CASE_MACRO
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 #endif  // __MODEL_TYPE_HH__
diff --git a/tests/SConscript b/tests/SConscript
index 014e436..42e6ca7 100644
--- a/tests/SConscript
+++ b/tests/SConscript
@@ -1,106 +1,112 @@
 from __future__ import print_function
 import os
 from os.path import join, abspath, basename
 
 
 # ------------------------------------------------------------------------------
 def copyComStr(env, main):
   if 'SHCXXCOMSTR' in main:
     env['CXXCOMSTR'] = main['SHCXXCOMSTR']
   if 'SHLINKCOMSTR' in main:
     env['LINKCOMSTR'] = main['SHLINKCOMSTR']
 
 
 # ------------------------------------------------------------------------------
 def make_python_tests(env):
     """Copy python tests to build directory"""
     test_env = env.Clone(
         PRINT_CMD_LINE_FUNC=main_env['gen_print']("Copying", "red", main_env))
 
     test_files = Split("""
     run_tests.sh
     test_hertz_pressure.py
     test_westergaard.py
     test_patch_westergaard.py
     test_surface.py
     test_autocorrelation.py
     test_hertz_disp.py
     test_hertz_kato.py
     test_saturated_pressure.py
     test_bem_grid.py
     test_flood_fill.py
+    fftfreq.py
     """)
 
     src_dir = "#/tests"
     build_dir = 'build-' + main_env['build_type'] + '/tests'
 
     for file in test_files:
         source = join(src_dir, file)
         test_env.Command(file, source, Copy("$TARGET", "$SOURCE"))
 
 
 # ------------------------------------------------------------------------------
 def compile_google_test(env, gtest_path):
     gtest_obj = env.Object('gtest.o', [join(gtest_path, "src/gtest-all.cc")])
     env.StaticLibrary('gtest', gtest_obj)
 
 # ------------------------------------------------------------------------------
 def make_google_tests(env):
     gtest_dir = Dir('#/third-party/googletest/googletest')
 
     gtest_env = env.Clone(CPPPATH=[gtest_dir],
                           CXXFLAGS=['-pthread', '-isystem',
                                     join(str(gtest_dir), "include")])
 
     colors = env['COLOR_DICT']
 
     if not env['verbose']:
         gtest_env['ARCOMSTR'] = u'{}[Ar]{} $TARGET'.format(colors['purple'],
                                                            colors['end'])
         gtest_env['RANLIBCOMSTR'] = u'{}[Randlib]{} $TARGET'.format(colors['purple'],
                                                                     colors['end'])
     gtest_env.Tool(gtest)
 
     gtest_path = str(gtest_dir)
     compile_google_test(gtest_env, gtest_path)
 
     env.AppendUnique(LIBS=['Tamaas'],
                      CXXFLAGS=gtest_env['CXXFLAGS'])
 
     google_test_files = Split("""
                               test_fft.cpp
                               test_grid.cpp
                               test_loop.cpp
                               test_model.cpp
                               """)
 
     google_test_files += ['libgtest.a']
 
-    gtest_main = env.Object("gtest_main.o", join(gtest_path, "src/gtest_main.cc"))
+    gtest_main = env.Object("gtest_main.o", 'tamaas_gtest_main.cc')
     env.Program('test_gtest_all', google_test_files + [gtest_main])
   
 
 # ------------------------------------------------------------------------------
 def make_bare_tests(env):
     env.Program("test_rough_surface.cpp")
 
 # ------------------------------------------------------------------------------
 Import('main_env')
 
 # Setup of test environment
 test_env = main_env.Clone()
 test_env.AppendUnique(LIBS=['Tamaas'],
                       LIBPATH=[abspath('build-' + main_env['build_type'] + '/src')])
 copyComStr(test_env, main_env)
 
 
 # Building tests that do not require any third party
 make_bare_tests(test_env)
 
 if test_env['build_python']:
+  test_env.Tool(pybind11)
+  #test_env.AppendUnique(LIBS=['python3.6m'])
+  test_env.ParseConfig("{}-config --ldflags".format(test_env['py_exec']))
+  test_env['CCFLAGS'] = []
+
   make_python_tests(test_env)
 
 # Building google tests
 if test_env['use_googletest']:
   make_google_tests(test_env)
 
diff --git a/tests/fftfreq.py b/tests/fftfreq.py
new file mode 100644
index 0000000..a58f780
--- /dev/null
+++ b/tests/fftfreq.py
@@ -0,0 +1,32 @@
+from __future__ import division
+
+__author__ = "Lucas Frérot"
+
+import numpy as np
+from numpy.fft import fftfreq
+
+
+def frequencies1D(frequencies):
+    n = len(frequencies)
+    frequencies[:] = n * fftfreq(n)
+
+
+def assign2DFrequencies(qx, qy, frequencies):
+    qx, qy = np.meshgrid(qx, qy, indexing='ij')
+
+    frequencies[:, :, 0] = qx
+    frequencies[:, :, 1] = qy
+
+
+def frequencies2D(frequencies):
+    n = frequencies.shape
+    qx = fftfreq(n[0]) * n[0]
+    qy = fftfreq(n[1]) * n[1]
+    assign2DFrequencies(qx, qy, frequencies)
+
+
+def hfrequencies2D(frequencies):
+    n = frequencies.shape
+    qx = fftfreq(n[0]) * n[0]
+    qy = np.arange(n[1])
+    assign2DFrequencies(qx, qy, frequencies)
diff --git a/src/model/model_type.cpp b/tests/tamaas_gtest_main.cc
similarity index 55%
copy from src/model/model_type.cpp
copy to tests/tamaas_gtest_main.cc
index 5fc2163..641ab76 100644
--- a/src/model/model_type.cpp
+++ b/tests/tamaas_gtest_main.cc
@@ -1,47 +1,57 @@
 /**
- * @file
+ *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
- * Copyright (©)  2016-2017 EPFL  (Ecole Polytechnique  Fédérale de
+ * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
-#include "model_type.hh"
-/* -------------------------------------------------------------------------- */
-__BEGIN_TAMAAS__
+#include "test.hh"
+#ifdef USE_PYTHON
+#include <pybind11/pybind11.h>
+#include <pybind11/embed.h>
+#include <pybind11/numpy.h>
+#endif
 /* -------------------------------------------------------------------------- */
 
-std::vector<view::index> model_type_traits<model_type::basic_1d>::indices = {
-    view::free()};
-std::vector<view::index> model_type_traits<model_type::basic_2d>::indices = {
-    view::free(), view::free()};
-
-std::vector<view::index> model_type_traits<model_type::surface_1d>::indices = {
-    view::free()};
-std::vector<view::index> model_type_traits<model_type::surface_2d>::indices = {
-    view::free(), view::free()};
+class TamaasEnvironment : public ::testing::Environment {
+public:
+  void SetUp() override {
+    tamaas::initialize();
+  }
 
-std::vector<view::index> model_type_traits<model_type::volume_1d>::indices = {
-    view::blocked(0), view::free()};
-std::vector<view::index> model_type_traits<model_type::volume_2d>::indices = {
-    view::blocked(0), view::free(), view::free()};
+  void TearDown() override {
+    tamaas::finalize();
+  }
+};
 
 /* -------------------------------------------------------------------------- */
-__END_TAMAAS__
+int main(int argc, char ** argv) {
+#ifdef USE_PYTHON
+  std::cerr << "Initializing python interpreter" << std::endl;
+  pybind11::scoped_interpreter guard{};
+#endif
+  ::testing::InitGoogleTest(&argc, argv);
+  ::testing::AddGlobalTestEnvironment(new TamaasEnvironment());
+
+  return RUN_ALL_TESTS();
+}
+
+
diff --git a/tests/test_fft.cpp b/tests/test_fft.cpp
index eec9114..0a37e93 100644
--- a/tests/test_fft.cpp
+++ b/tests/test_fft.cpp
@@ -1,81 +1,140 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "test.hh"
 #include "grid.hh"
 #include "grid_hermitian.hh"
 #include "fft_plan_manager.hh"
 
+#ifdef USE_PYTHON
+#include <pybind11/pybind11.h>
+#include <pybind11/embed.h>
+#include <pybind11/numpy.h>
+#endif
+
 using namespace tamaas;
 
 /* -------------------------------------------------------------------------- */
 
 TEST(TestFFTInterface, FFT1D) {
   constexpr UInt size = 20;
 
   double data[size] = {0};
   fftw_complex solution[size/2 + 1] = {{0}};
 
   std::iota(std::begin(data), std::end(data), 0);
   fftw_plan solution_plan = fftw_plan_dft_r2c_1d(size, data, solution, FFTW_ESTIMATE);
   fftw_execute(solution_plan);
 
   Grid<Real, 1> grid({size}, 1);
   std::iota(grid.begin(), grid.end(), 0);
   GridHermitian<Real, 1> result({size/2 + 1}, 1);
 
   FFTPlanManager::get().createPlan(grid, result).forwardTransform();
 #ifdef USE_CUDA
   cudaDeviceSynchronize();
 #endif
 
   ASSERT_TRUE(compare(result, solution, AreComplexEqual())) << "1D FFTW transform failed";
   fftw_destroy_plan(solution_plan);
   FFTPlanManager::get().destroyPlan(grid, result);
 }
 
+/* -------------------------------------------------------------------------- */
 TEST(TestFFTInterface, FFT2D) {
   constexpr UInt size = 10;
 
   double data[size * size] = {0};
   fftw_complex solution[size * (size/2 + 1)] = {{0}};
 
   std::iota(std::begin(data), std::end(data), 0);
   fftw_plan solution_plan = fftw_plan_dft_r2c_2d(size, size, data, solution, FFTW_ESTIMATE);
   fftw_execute(solution_plan);
 
   Grid<Real, 2> grid({size, size}, 1);
   std::iota(grid.begin(), grid.end(), 0);
   GridHermitian<Real, 2> result({size, size/2 + 1}, 1);
 
   FFTPlanManager::get().createPlan(grid, result).forwardTransform();
 #ifdef USE_CUDA
   cudaDeviceSynchronize();
 #endif
 
   ASSERT_TRUE(compare(result, solution, AreComplexEqual())) << "2D FFTW transform failed";
   fftw_destroy_plan(solution_plan);
   FFTPlanManager::get().destroyPlan(grid, result);
 }
+
+
+/* -------------------------------------------------------------------------- */
+#ifdef USE_PYTHON
+
+namespace py = pybind11;
+
+TEST(TestFFTInterface, Frequencies1D) {
+  std::array<UInt, 1> sizes = {{10}};
+  auto freq = FFTransform<Real, 1>::computeFrequencies<false>(sizes);
+
+  py::module fftfreq = py::module::import("fftfreq");
+
+  std::vector<Real> reference(freq.dataSize());
+  py::array py_arg(reference.size(), reference.data(), py::none());
+
+  fftfreq.attr("frequencies1D")(py_arg);
+
+  ASSERT_TRUE(compare(reference, freq, AreFloatEqual()))
+      << "Non hermitian frequencies are wrong";
+
+  auto hfreq = FFTransform<Real, 1>::computeFrequencies<true>(sizes);
+  std::iota(reference.begin(), reference.end(), 0);
+
+  for (Real r : reference)
+    std::cout << r << std::endl;
+  std::cout << hfreq << std::endl;
+  ASSERT_TRUE(compare(reference, hfreq, AreFloatEqual()))
+      << "Hermitian frequencies are wrong";
+}
+
+TEST(TestFFTInterface, Frequencies2D) {
+  std::array<UInt, 2> sizes = {{10, 10}};
+  auto freq = FFTransform<Real, 2>::computeFrequencies<false>(sizes);
+
+  py::module fftfreq = py::module::import("fftfreq");
+
+  std::vector<Real> reference(freq.dataSize());
+  py::array py_arg({10, 10, 2}, reference.data(), py::none());
+
+  fftfreq.attr("frequencies2D")(py_arg);
+  ASSERT_TRUE(compare(reference, freq, AreFloatEqual()))
+      << "Non hermitian frequencies are wrong";
+  
+
+  auto hfreq = FFTransform<Real, 2>::computeFrequencies<true>(sizes);
+  fftfreq.attr("hfrequencies2D")(py_arg);
+  ASSERT_TRUE(compare(reference, hfreq, AreFloatEqual()))
+      << "Hermitian frequencies are wrong";
+}
+
+#endif
diff --git a/tests/test_grid.cpp b/tests/test_grid.cpp
index ec558cf..57bac89 100644
--- a/tests/test_grid.cpp
+++ b/tests/test_grid.cpp
@@ -1,127 +1,143 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "test.hh"
 #include "grid.hh"
+#include "grid_view.hh"
 
 using namespace tamaas;
 
 /* -------------------------------------------------------------------------- */
 
 // Testing 1D creation
 TEST(TestGridCreation, Creation1d) {
   Grid<Real, 1> grid{{5}, 2};
   std::array<UInt, 1> correct_size{{5}};
   ASSERT_TRUE(grid.sizes() == correct_size) << "Wrong sizes";
 
   std::array<UInt, 2> correct_strides{{2, 1}};
   ASSERT_TRUE(grid.getStrides() == correct_strides) << "Wrong strides";
 }
 
 // Testing 2D creation
 TEST(TestGridCreation, Creation2d) {
   Grid<Real, 2> grid({5, 2}, 3);
   std::array<UInt, 2> correct_size{{5, 2}};
   ASSERT_TRUE(grid.sizes() == correct_size) << "Wrong sizes";
 
   std::array<UInt, 3> correct_strides{{6, 3, 1}};
   ASSERT_TRUE(grid.getStrides() == correct_strides) << "Wrong strides";
 }
 
 // Testing 3D creation
 TEST(TestGridCreation, Creation3d) {
   Grid<Real, 3> grid({8, 5, 2}, 3);
   std::array<UInt, 3> correct_size{{8, 5, 2}};
   ASSERT_TRUE(grid.sizes() == correct_size) << "Wrong sizes";
 
   std::array<UInt, 4> correct_strides{{30, 15, 3, 1}};
   ASSERT_TRUE(grid.getStrides() == correct_strides) << "Wrong strides";
 }
 
 /* -------------------------------------------------------------------------- */
 
 // Testing iterators in STL function iota
 TEST(TestGridIterators, Iota) {
   constexpr UInt N = 20;
   Grid<UInt, 1> grid({N}, 1);
   std::iota(grid.begin(), grid.end(), 0);
 
   const UInt * p = grid.getInternalData();
   for (UInt i = 0 ; i < N ; ++i)
     ASSERT_TRUE(p[i] == i) << "Iota fill failed";
 }
 
 // Testing filling grid with OpenMP loop
 TEST(TestGridIterators, OpenMPLoops) {
   Grid<UInt, 1> grid({20}, 1);
 
 #ifndef USE_CUDA
 #pragma omp parallel for
 #endif
   for (auto it = grid.begin() ; it < grid.end() ; ++it) {
     UInt i = it - grid.begin();
     *it = i;
   }
 
   std::vector<UInt> correct(20);
   std::iota(correct.begin(), correct.end(), 0);
 
   ASSERT_TRUE(compare(grid, correct)) << "Fill using OpenMP loop failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 // Testing scalar simple loop-based operators
 TEST(TestGridOperators, LoopOperators) {
   Grid<UInt, 1> grid({20}, 1);
   grid = 1;
   EXPECT_TRUE(std::all_of(grid.begin(), grid.end(), [](UInt x) {
     return x == 1;
   })) << "Assignement operator failed";
   grid += 1;
   EXPECT_TRUE(std::all_of(grid.begin(), grid.end(), [](UInt x) {
     return x == 2;
   })) << "Plus-equal operator failed";
 }
 
 // Testing loop-based functions with reductions
 TEST(TestGridOperators, Stats) {
   constexpr UInt N = 20;
   Grid<Real, 1> grid({N}, 1);
   std::iota(grid.begin(), grid.end(), 0);
 
   auto b = grid.begin(), e = grid.end();
   Real min = *std::min_element(b, e);
   Real max = *std::max_element(b, e);
   Real acc = std::accumulate(b, e, 0);
   Real mea = acc / N;
   Real var = 35;
 
-  EXPECT_TRUE(grid.min() == min) << "Minimum function failed";
-  EXPECT_TRUE(grid.max() == max) << "Maximum function failed";
-  EXPECT_TRUE(grid.sum() == acc) <<"Sum function failed";
-  EXPECT_TRUE(grid.mean() == mea) << "Mean function failed";
-  EXPECT_TRUE(grid.var() == var) << "Var function failed";
+  EXPECT_DOUBLE_EQ(grid.min(), min) << "Minimum function failed";
+  EXPECT_DOUBLE_EQ(grid.max(), max) << "Maximum function failed";
+  EXPECT_DOUBLE_EQ(grid.sum(), acc) <<"Sum function failed";
+  EXPECT_DOUBLE_EQ(grid.mean(), mea) << "Mean function failed";
+  EXPECT_DOUBLE_EQ(grid.var(), var) << "Var function failed";
+}
+
+/* -------------------------------------------------------------------------- */
+
+TEST(TestGridView, View2D) {
+  Grid<UInt, 3> grid3d({10, 10, 10}, 3);
+  std::iota(grid3d.begin(), grid3d.end(), 0);
+
+  for (UInt i : Loop::range(10)) {
+    auto grid_view = make_view(grid3d, i);
+    std::vector<UInt> reference(300);
+    std::iota(reference.begin(), reference.end(), 300 * i);
+    EXPECT_TRUE(
+        std::equal(grid_view.begin(), grid_view.end(), reference.begin()));
+  }
 }