diff --git a/src/core/grid.hh b/src/core/grid.hh
index 8361b3f..74a37b9 100644
--- a/src/core/grid.hh
+++ b/src/core/grid.hh
@@ -1,188 +1,184 @@
 /**
  * @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_HH__
 #define __GRID_HH__
 /* -------------------------------------------------------------------------- */
 #include "grid_base.hh"
 #include "tamaas.hh"
 #include <array>
 #include <numeric>
 #include <utility>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /**
  * @brief Multi-dimensional & multi-component array class
  *
  * This class is a container for multi-component data stored on a multi-
  * dimensional grid.
  *
  * The access function is the parenthesis operator. For a grid of dimension d,
  * the operator takes d+1 arguments: the first d arguments are the position on
  * the grid and the last one is the component of the stored data.
  *
  * It is also possible to use the access operator with only one argument, it is
  * then considering the grid as a flat array, accessing the given cell of the
  * array.
  */
 template <typename T, UInt dim>
 class Grid : public GridBase<T> {
 public:
   /* ------------------------------------------------------------------------ */
   /* Types */
   /* ------------------------------------------------------------------------ */
   using value_type = T;
   static constexpr UInt dimension = dim;
 
   /* ------------------------------------------------------------------------ */
   /* Constructors */
   /* ------------------------------------------------------------------------ */
 public:
   /// Constructor by default (empty array)
   Grid();
   /// Constructor
   Grid(const std::array<UInt, dim>& n, UInt nb_components);
   /// Wrap on given data
   Grid(const std::array<UInt, dim>& n, UInt nb_components, value_type * data);
   /// Constructor with vectors
   Grid(const std::vector<UInt>& n, UInt nb_components);
   /// Constructor with initializer list
   Grid(const std::initializer_list<UInt>& n, UInt nb_components);
   /// Copy constructor
   Grid(const Grid& o);
   /// Move constructor (transfers data ownership)
   Grid(Grid&& o);
   /// Destructor
   ~Grid() override = default;
 
 private:
   /// Init from standard container
   template <typename Container>
   void init(const Container& n, UInt nb_components);
 
 public:
   /* ------------------------------------------------------------------------ */
   /* Common operations */
   /* ------------------------------------------------------------------------ */
   /// Resize array
   void resize(const std::array<UInt, dim>& n);
   /// Compute size
   inline UInt computeSize() const;
   /// Get grid dimension
   inline UInt getDimension() const override { return dim; }
   /// Compute strides
   virtual void computeStrides();
   /// Print
   virtual void printself(std::ostream& str) const;
   /// Get sizes
   const std::array<UInt, dim>& sizes() const { return n; }
   /// Get strides
   const std::array<UInt, dim + 1>& getStrides() const { return this->strides; }
 
   /* ------------------------------------------------------------------------ */
   /* Access operators */
   /* ------------------------------------------------------------------------ */
   /// Variadic access operator (non-const)
   template <typename... T1>
   inline T& operator()(T1... args);
   /// Variadic access operator
   template <typename... T1>
   inline const T& operator()(T1... args) const;
 
   /// Tuple index access operator
   inline T& operator()(std::array<UInt, dim + 1> tuple);
   inline const T& operator()(std::array<UInt, dim + 1> tuple) const;
 
 private:
   /// Unpacking the arguments of variadic ()
   template <typename... T1>
   inline UInt unpackOffset(UInt offset, UInt index_pos, UInt index,
                            T1... rest) const;
   /// End case for recursion
   template <typename... T1>
   inline UInt unpackOffset(UInt offset, UInt index_pos, UInt index) const;
   /// Computing offset for a tuple index
   inline UInt computeOffset(std::array<UInt, dim + 1> tuple) const;
 
   /* ------------------------------------------------------------------------ */
   /* Move/Copy operators */
   /* ------------------------------------------------------------------------ */
 public:
+  using GridBase<T>::operator=;
   // = operator
   Grid& operator=(const Grid& other);
   // = operator (not const input, otherwise gcc is confused)
   Grid& operator=(Grid& other);
   // = operator (move)
   Grid& operator=(Grid&& other);
-  // = operator for scalar
-  Grid& operator=(T val) {
-    GridBase<T>::operator=(val);
-    return *this;
-  }
 
   // Copy data from another grid
   template <typename T1>
   void copy(const Grid<T1, dim>& other);
   // Move data from another grid
   template <typename T1>
   void move(Grid<T1, dim>&& other);
 
   // Wrap memory (non-const)
   void wrap(Grid & other) {
     GridBase<T>::wrap(other);
     std::copy(other.n.begin(), other.n.end(), n.begin());
     this->computeStrides();
   }
 
   // Wrap memory
   void wrap(const Grid & other) {
     GridBase<T>::wrap(other);
     std::copy(other.n.begin(), other.n.end(), n.begin());
     this->computeStrides();
   }
 
   /* ------------------------------------------------------------------------ */
   /* Member variables */
   /* ------------------------------------------------------------------------ */
 protected:
   std::array<UInt, dim> n;            ///< shape of grid: size per dimension
   std::array<UInt, dim + 1> strides;  ///< strides for access
 };
 
 __END_TAMAAS__
 
 /* -------------------------------------------------------------------------- */
 /* Inline/template function definitions                                       */
 /* -------------------------------------------------------------------------- */
 
 #include "grid_tmpl.hh"
 
 /* -------------------------------------------------------------------------- */
 
 #endif  // __GRID_HH__
diff --git a/src/core/grid_base.hh b/src/core/grid_base.hh
index 0d17874..bcd296b 100644
--- a/src/core/grid_base.hh
+++ b/src/core/grid_base.hh
@@ -1,298 +1,338 @@
 /**
  * @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 "static_types.hh"
 #include <cstddef>
 #include <limits>
 #include <utility>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /**
  * @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);
+/// Broadcast operators
+#define BROADCAST_OPERATOR(op)                                                 \
+  template <typename DT, typename ST, UInt size>                               \
+  GridBase& operator op(const StaticArray<DT, ST, size>& b)
+
+  BROADCAST_OPERATOR(+=);
+  BROADCAST_OPERATOR(-=);
+  BROADCAST_OPERATOR(*=);
+  BROADCAST_OPERATOR(/=);
+  BROADCAST_OPERATOR(=);
+#undef BROADCAST_OPERATOR
 
   /* ------------------------------------------------------------------------ */
   /* 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:
+  GridBase& operator=(const GridBase& o) {
+    this->copy(o);
+    return *this;
+  }
+
+  GridBase& operator=(GridBase& o) {
+    this->copy(o);
+    return *this;
+  }
+
+  GridBase& operator=(GridBase&& o) {
+    this->move(std::move(o));
+    return *this;
+  }
+
   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;
-}
+#define BROADCAST_OPERATOR_IMPL(op)                                            \
+  template <typename T>                                                        \
+  template <typename DT, typename ST, UInt size>                               \
+  GridBase<T>& GridBase<T>::operator op(const StaticArray<DT, ST, size>& b) {  \
+    TAMAAS_ASSERT(this->dataSize() % size == 0,                                \
+                  "Broadcast operator cannot broadcast"                        \
+                      << size << " on array of size " << this->dataSize());    \
+    Tensor<StaticArray, DT, size> a(b);                                        \
+    Loop::stridedLoop([a] CUDA_LAMBDA(VectorProxy<T, size>&& g) { g op a; },   \
+                      *this);                                                  \
+    return *this;                                                              \
+  }
 
-/* -------------------------------------------------------------------------- */
+BROADCAST_OPERATOR_IMPL(+=);
+BROADCAST_OPERATOR_IMPL(-=);
+BROADCAST_OPERATOR_IMPL(*=);
+BROADCAST_OPERATOR_IMPL(/=);
+BROADCAST_OPERATOR_IMPL(=);
+#undef BROADCAST_OPERATOR_IMPL
 
-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);
+  return Loop::reduce<operation::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);
+  return Loop::reduce<operation::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);
+  return Loop::reduce<operation::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>(
+  const T var = Loop::reduce<operation::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>(
+  return Loop::reduce<operation::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 2f17992..922d217 100644
--- a/src/core/grid_tmpl.hh
+++ b/src/core/grid_tmpl.hh
@@ -1,170 +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;
   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));
+  this->move(std::move(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/static_types.hh b/src/core/static_types.hh
index 303093f..f15e0f3 100644
--- a/src/core/static_types.hh
+++ b/src/core/static_types.hh
@@ -1,328 +1,333 @@
 /**
  * @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 __STATIC_TYPES_HH__
 #define __STATIC_TYPES_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 
 __BEGIN_TAMAAS__
 
 /* -------------------------------------------------------------------------- */
 
 namespace detail {
 template <UInt acc, UInt n, UInt... ns>
 struct product_tail_rec : product_tail_rec<acc * n, ns...> {};
 
 template <UInt acc, UInt n>
 struct product_tail_rec<acc, n> : std::integral_constant<UInt, acc * n> {};
 
 template <UInt N, UInt n, UInt... ns>
 struct get_rec : get_rec<N - 1, ns...> {};
 
 template <UInt n, UInt... ns>
 struct get_rec<0, n, ns...> : std::integral_constant<UInt, n> {};
 }
 
 template <UInt... ns>
 struct product : detail::product_tail_rec<1, ns...> {};
 
 template <UInt N, UInt... ns>
 struct get : detail::get_rec<N, ns...> {};
 
 template <typename T>
 struct is_arithmetic : std::is_arithmetic<T> {};
 
 template <typename T>
 struct is_arithmetic<thrust::complex<T>> : std::true_type {};
 
 /**
  * @brief Static Array
  *
  * This class is meant to be a small and fast object for intermediate
  * calculations, possibly on wrapped memory belonging to a grid. Support type
  * show be either a pointer or a C array. It should not contain any virtual
  * method.
  */
 template <typename DataType, typename SupportType, UInt size>
 class StaticArray {
   static_assert(std::is_array<SupportType>::value ||
                     std::is_pointer<SupportType>::value,
                 "the support type of StaticArray should be either a pointer or "
                 "a C-array");
 
   using T = DataType;
 
 public:
   /// Access operator
   __device__ __host__ T& operator()(UInt i) {
     // TAMAAS_ASSERT(i < n, "Access out of bounds");
     return _mem[i];
   }
 
   /// Access operator
   __device__ __host__ const T& operator()(UInt i) const {
     // TAMAAS_ASSERT(i < n, "Access out of bounds");
     return _mem[i];
   }
 
   /// L2 norm squared
   __device__ __host__ T l2squared() const {
     T res = 0;
     for (UInt i = 0; i < size; ++i)
-      res += _mem[i] * _mem[i];
+      res += (*this)(i) * (*this)(i);
     return res;
   }
 
   /// L2 norm
   __device__ __host__ T l2norm() const { return std::sqrt(l2squared()); }
 
 #define VECTOR_OP(op)                                                          \
   template <typename DT, typename ST>                                          \
   __device__ __host__ void operator op(const StaticArray<DT, ST, size>& o) {   \
     for (UInt i = 0; i < size; ++i)                                            \
       (*this)(i) op o(i);                                                      \
   }
 
   VECTOR_OP(+=)
   VECTOR_OP(-=)
   VECTOR_OP(*=)
   VECTOR_OP(/=)
 #undef VECTOR_OP
 
 #define SCALAR_OP(op)                                                          \
   template <typename T1>                                                       \
   __device__ __host__                                                          \
       typename std::enable_if<is_arithmetic<T1>::value, StaticArray&>::type    \
       operator op(const T1& x) {                                               \
     for (UInt i = 0; i < size; ++i)                                            \
       (*this)(i) op x;                                                         \
     return *this;                                                              \
   }
 
   SCALAR_OP(+=)
   SCALAR_OP(-=)
   SCALAR_OP(*=)
   SCALAR_OP(/=)
   SCALAR_OP(=)
 #undef SCALAR_OP
 
-  /// Overriding the implit copy operator
+  /// Overriding the implicit copy operator
   StaticArray& operator=(const StaticArray& o) { return this->copy(o); }
 
+  template <typename DT, typename ST>
+  void operator=(const StaticArray<DT, ST, size>& o) {
+    this->copy(o);
+  }
+
   template <typename DT, typename ST>
   StaticArray& copy(const StaticArray<DT, ST, size>& o) {
     for (UInt i = 0; i < size; ++i)
       (*this)(i) = o(i);
     return *this;
   }
 
 protected:
   SupportType _mem;
 };
 
 /**
  * @brief Static Tensor
  *
  * This class implements a multi-dimensional tensor behavior.
  */
 template <typename DataType, typename SupportType = DataType*, UInt... dims>
 class StaticTensor
     : public StaticArray<DataType, SupportType, product<dims...>::value> {
   using parent = StaticArray<DataType, SupportType, product<dims...>::value>;
   using T = DataType;
 
 public:
   static constexpr UInt dim = sizeof...(dims);
   using parent::operator=;
 
 private:
   template <typename... Idx>
   static UInt unpackOffset(UInt offset, UInt index, Idx... rest) {
     constexpr UInt size = sizeof...(rest);
     offset += index;
     offset *= get<dim - size, dims...>::value;
     return unpackOffset(offset, rest...);
   }
 
   template <typename... Idx>
   static UInt unpackOffset(UInt offset, UInt index) {
     return offset + index;
   }
 
 public:
   template <typename... Idx>
   const T& operator()(Idx... idx) const {
-    return this->_mem[unpackOffset(0, idx...)];
+    return parent::operator()(unpackOffset(0, idx...));
   }
 
   template <typename... Idx>
   T& operator()(Idx... idx) {
-    return this->_mem[unpackOffset(0, idx...)];
+    return parent::operator()(unpackOffset(0, idx...));
   }
 };
 
 /* -------------------------------------------------------------------------- */
 /* Common Static Types */
 /* -------------------------------------------------------------------------- */
 template <typename DataType, typename SupportType, UInt n, UInt m>
 using StaticMatrix = StaticTensor<DataType, SupportType, n, m>;
 
 /// Vector class with size determined at compile-time
 template <typename DataType, typename SupportType, UInt n>
 class StaticVector : public StaticTensor<DataType, SupportType, n> {
   using T = DataType;
 
 public:
   using StaticTensor<DataType, SupportType, n>::operator=;
   /// Matrix-vector product
   template <typename DT1, typename ST1, typename DT2, typename ST2, UInt m>
   __device__ __host__ void mul(const StaticMatrix<DT1, ST1, n, m>& mat,
                                const StaticVector<DT2, ST2, m>& vec) {
     *this = T(0);
     for (UInt i = 0; i < n; ++i)
       for (UInt j = 0; j < m; ++j)
         (*this)(i) += mat(i, j) * vec(j);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 /* Proxy Static Types */
 /* -------------------------------------------------------------------------- */
 
 /// Proxy type for tensor
 template <template <typename, typename, UInt...> class StaticParent, typename T,
           UInt... dims>
 class TensorProxy : public StaticParent<T, T*, dims...> {
 public:
   /// Implicit wrap constructor
   __device__ __host__ TensorProxy(T& spot) { this->_mem = &spot; }
 
   using StaticTensor<T, T*, dims...>::operator=;
 };
 
 template <typename T, UInt n, UInt m>
 using MatrixProxy = TensorProxy<StaticMatrix, T, n, m>;
 
 template <typename T, UInt n>
 using VectorProxy = TensorProxy<StaticVector, T, n>;
 
 /* -------------------------------------------------------------------------- */
 /* On the stack static types */
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, typename, UInt...> class StaticParent, typename T,
           UInt... dims>
 class Tensor : public StaticParent<T, T[product<dims...>::value], dims...> {
   using parent = StaticParent<T, T[product<dims...>::value], dims...>;
 
 public:
   using parent::operator=;
   using parent::copy;
 
   /// Default constructor
   Tensor() = default;
   /// Construct with default value
   Tensor(T val) { *this = val; }
 
   /// Construct by copy from static tensor
   template <typename DT, typename ST>
   Tensor(const StaticParent<DT, ST, dims...>& o) {
     this->copy(o);
   }
 };
 
 template <typename T, UInt n, UInt m>
 using Matrix = Tensor<StaticMatrix, T, n, m>;
 
 template <typename T, UInt n>
 using Vector = Tensor<StaticVector, T, n>;
 
 /* -------------------------------------------------------------------------- */
 /* Arithmetic operators creating temporaries */
 /* -------------------------------------------------------------------------- */
 
 // template <template <typename, typename, UInt...> class StaticParent,
 //           typename DT1, typename ST1, typename DT2, typename ST2, UInt...
 //           dims>
 // Tensor<StaticParent, decltype(DT1(0) + DT2(0)), dims...>
 // operator+(const StaticParent<DT1, ST1, dims...>& t1,
 //           const StaticParent<DT2, ST2, dims...>& t2) {
 //   Tensor<StaticParent, decltype(DT1(0) + DT2(0)), dims...> a(t1);
 //   return a += t2;
 // }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt dim>
 Vector<decltype(DT1(0) + DT2(0)), dim>
 operator+(const StaticVector<DT1, ST1, dim>& a,
           const StaticVector<DT2, ST2, dim>& b) {
   Vector<decltype(DT1(0) + DT2(0)), dim> res(a);
   res += b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt dim>
 Vector<decltype(DT1(0) - DT2(0)), dim>
 operator-(const StaticVector<DT1, ST1, dim>& a,
           const StaticVector<DT2, ST2, dim>& b) {
   Vector<decltype(DT1(0) - DT2(0)), dim> res(a);
   res -= b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m>
 Matrix<decltype(DT1(0) + DT2(0)), n, m>
 operator+(const StaticMatrix<DT1, ST1, n, m>& a,
           const StaticMatrix<DT2, ST2, n, m>& b) {
   Matrix<decltype(DT1(0) + DT2(0)), n, m> res(a);
   res += b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m>
 Matrix<decltype(DT1(0) - DT2(0)), n, m>
 operator-(const StaticMatrix<DT1, ST1, n, m>& a,
           const StaticMatrix<DT2, ST2, n, m>& b) {
   Matrix<decltype(DT1(0) - DT2(0)), n, m> res(a);
   res -= b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m>
 Vector<decltype(DT1(0) * DT2(0)), n>
 operator*(const StaticMatrix<DT1, ST1, n, m>& a,
           const StaticVector<DT2, ST2, m>& b) {
   Vector<decltype(DT1(0) * DT2(0)), n> res;
   res.mul(a, b);
   return res;
 }
 
 __END_TAMAAS__
 
 #endif  // __STATIC_TYPES_HH__
diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp
index 7604742..6ac3a69 100644
--- a/src/model/westergaard.cpp
+++ b/src/model/westergaard.cpp
@@ -1,238 +1,236 @@
 /**
  * @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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "westergaard.hh"
 #include "fft_plan_manager.hh"
 #include "influence.hh"
 #include "loop.hh"
 #include "model.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 Westergaard<type>::Westergaard(Model* model) : BEEngine(model) {
   // Copy sizes
   std::array<UInt, trait::dimension> sizes;
   std::copy(model->getDiscretization().begin(),
             model->getDiscretization().end(), sizes.begin());
   auto hermitian_sizes =
       GridHermitian<Real, trait::dimension>::hermitianDimensions(sizes);
 
   // Copy boundary sizes
   auto sizes_it = sizes.begin();
   std::array<UInt, trait::boundary_dimension> boundary_hermitian_sizes;
   // Ignoring first dimension if model is volumetric
   if (trait::dimension > trait::boundary_dimension)
     ++sizes_it;
   std::copy(sizes_it, sizes.end(), boundary_hermitian_sizes.begin());
   boundary_hermitian_sizes =
       GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
           boundary_hermitian_sizes);
 
   constexpr UInt nb_components = trait::components;
 
   buffer.setNbComponents(nb_components);
   buffer.resize(boundary_hermitian_sizes);
   influence.setNbComponents(nb_components * nb_components);
   influence.resize(hermitian_sizes);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Westergaard<type>::initInfluence() {
   TAMAAS_EXCEPTION("Not implemented");
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define INFLUENCE_BASIC_MACRO                                                  \
   do {                                                                         \
     constexpr UInt bdim = trait::boundary_dimension;                           \
     auto wavevectors =                                                         \
         FFTransform<Real, bdim>::computeFrequencies<true>(influence.sizes());  \
     Real E_star = model->getHertzModulus();                                    \
-    VectorProxy<const Real, bdim> domain(model->getSystemSize()[0]);          \
+    VectorProxy<const Real, bdim> domain(model->getSystemSize()[0]);           \
     Loop::stridedLoop(influence::Basic<bdim>(E_star, domain), wavevectors,     \
                       influence);                                              \
     influence(0) = 0;                                                          \
   } while (0)
 
 template <>
 void Westergaard<model_type::basic_2d>::initInfluence() {
   INFLUENCE_BASIC_MACRO;
 }
 
 template <>
 void Westergaard<model_type::basic_1d>::initInfluence() {
   INFLUENCE_BASIC_MACRO;
 }
 
 #undef INFLUENCE_BASIC_MACRO
 
 #define INFLUENCE_SURFACE_MACRO                                                \
   do {                                                                         \
     constexpr UInt bdim = trait::boundary_dimension;                           \
     auto wavevectors =                                                         \
         FFTransform<Real, bdim>::computeFrequencies<true>(influence.sizes());  \
     Real E = model->getYoungModulus();                                         \
     Real nu = model->getPoissonRatio();                                        \
-    VectorProxy<const Real, bdim> domain(model->getSystemSize()[0]);          \
+    VectorProxy<const Real, bdim> domain(model->getSystemSize()[0]);           \
     Loop::stridedLoop(influence::Surface<bdim>(E, nu, domain), wavevectors,    \
                       influence);                                              \
+    MatrixProxy<Complex, bdim + 1, bdim + 1> inf_0(influence(0));              \
+    inf_0 = 0;                                                                 \
   } while (0)
 
 template <>
 void Westergaard<model_type::surface_1d>::initInfluence() {
   INFLUENCE_SURFACE_MACRO;
 }
 
 template <>
 void Westergaard<model_type::surface_2d>::initInfluence() {
   INFLUENCE_SURFACE_MACRO;
 }
 
 #undef INFLUENCE_SURFACE_MACRO
 
 /* ------------------------------------------------------------------------ */
 
 #define NEUMANN_BASIC(type)                                                    \
   template <>                                                                  \
   void Westergaard<type>::solveNeumann(GridBase<Real>& in,                     \
                                        GridBase<Real>& out) const {            \
     auto apply = [](decltype(buffer)& buffer,                                  \
                     const decltype(influence)& influence) {                    \
       buffer *= influence;                                                     \
       buffer(0) = 0;                                                           \
     };                                                                         \
                                                                                \
     fourierApply(apply, in, out);                                              \
   }
 
 #define DIRICHLET_BASIC(type)                                                  \
   template <>                                                                  \
   void Westergaard<type>::solveDirichlet(GridBase<Real>& in,                   \
                                          GridBase<Real>& out) const {          \
     auto apply = [](decltype(buffer)& buffer,                                  \
                     const decltype(influence)& influence) {                    \
       buffer /= influence;                                                     \
       buffer(0) = 0;                                                           \
     };                                                                         \
                                                                                \
     fourierApply(apply, in, out);                                              \
   }
 
 NEUMANN_BASIC(model_type::basic_1d);
 NEUMANN_BASIC(model_type::basic_2d);
 DIRICHLET_BASIC(model_type::basic_1d);
 DIRICHLET_BASIC(model_type::basic_2d);
 
 #undef NEUMANN_BASIC
 #undef DIRICHLET_BASIC
 
 /* -------------------------------------------------------------------------- */
 
 namespace detail {
 /// Class needed for cuda
 template <UInt d>
 struct SurfaceApply {
   void operator()(GridHermitian<Real, d>& buffer_,
                   const GridHermitian<Real, d>& influence) {
     Loop::stridedLoop(
         [] CUDA_LAMBDA(VectorProxy<Complex, d + 1> && buf,
                        MatrixProxy<const Complex, d + 1, d + 1> && inf) {
-          Vector<Complex, d + 1> copy(buf);
-          buf.mul(inf, copy);
+	  buf = inf * buf;
         },
         buffer_, influence);
-
-    VectorProxy<Complex, d + 1> buff_0(buffer_(0));
-    buff_0 = 0;
   }
 };
 }  // namespace detail
 
 #define NEUMANN_SURFACE(type)                                                  \
   template <>                                                                  \
   void Westergaard<type>::solveNeumann(GridBase<Real>& in,                     \
                                        GridBase<Real>& out) const {            \
     fourierApply(detail::SurfaceApply<trait::boundary_dimension>(), in, out);  \
   }
 
 /// \TODO Find a way to make it work with cuda
 NEUMANN_SURFACE(model_type::surface_1d);
 NEUMANN_SURFACE(model_type::surface_2d);
 
 #undef NEUMANN_SURFACE
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Westergaard<type>::solveNeumann(GridBase<Real>& /*neumann*/,
                                      GridBase<Real>& /*dirichlet*/) const {
   TAMAAS_EXCEPTION("Not yet implemented");
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Westergaard<type>::solveDirichlet(GridBase<Real>& /*neumann*/,
                                        GridBase<Real>& /*dirichlet*/) const {
   TAMAAS_EXCEPTION("Not yet implemented");
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 Real Westergaard<type>::getNeumannNorm() {
   constexpr UInt comp = model_type_traits<type>::components;
   Real _norm = 0.0;
 
   _norm = Loop::stridedReduce<op::plus>(
     [] CUDA_LAMBDA(MatrixProxy<Complex, comp, comp>&& m) {
       Real n = thrust::norm(m.l2squared());
       return std::isnan(n) ? 0 : n;
     },
     influence);
 
   return std::sqrt(_norm);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation                                                     */
 /* -------------------------------------------------------------------------- */
 template class Westergaard<model_type::basic_1d>;
 template class Westergaard<model_type::basic_2d>;
 template class Westergaard<model_type::surface_1d>;
 template class Westergaard<model_type::surface_2d>;
 template class Westergaard<model_type::volume_1d>;
 template class Westergaard<model_type::volume_2d>;
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
diff --git a/tests/test_grid.cpp b/tests/test_grid.cpp
index 57bac89..eb6b202 100644
--- a/tests/test_grid.cpp
+++ b/tests/test_grid.cpp
@@ -1,143 +1,199 @@
 /**
  *
  * @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_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()));
   }
 }
+
+/* -------------------------------------------------------------------------- */
+TEST(TestGridOperators, BroadcastOperators) {
+  Grid<UInt, 2> grid({10, 10}, 3), solution({10, 10}, 3);
+
+  // Filling up solution
+  Loop::stridedLoop([](VectorProxy<UInt, 3>&& v) {
+      v(0) = 1;
+      v(1) = 2;
+      v(2) = 3;
+    }, solution);
+
+  Vector<UInt, 3> v;
+  v(0) = 1;
+  v(1) = 2;
+  v(2) = 3;
+
+  grid = v;
+
+  EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
+      << "Broadcast assignement failure";
+
+  v(1) = 1;
+  v(2) = 1;
+
+  solution += 1;
+  grid += v;
+
+  EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
+      << "Broadcast += failure";
+
+  solution -= 1;
+  grid -= v;
+
+  EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
+      << "Broadcast -= failure";
+
+  v(1) = 2;
+  v(2) = 3;
+
+  solution *= solution;
+  grid *= v;
+
+  EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
+      << "Broadcast *= failure";
+
+  Loop::stridedLoop([](VectorProxy<UInt, 3>&& v) {
+      v(0) = 1;
+      v(1) = 2;
+      v(2) = 3;
+    }, solution);
+  grid /= v;
+
+  EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
+      << "Broadcast /= failure";
+}
diff --git a/tests/test_loop.cpp b/tests/test_loop.cpp
index df1a9d4..4124a30 100644
--- a/tests/test_loop.cpp
+++ b/tests/test_loop.cpp
@@ -1,255 +1,254 @@
 /**
  *
  * @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 "static_types.hh"
 #include "loop.hh"
 /* -------------------------------------------------------------------------- */
 /* WARNING: here we cannot use lambdas for tests because GoogleTest declares  */
 /* test functions as private members of classes which is incompatible with    */
 /* cuda's extended lambdas. I know... it's f*cking stupid                     */
 /* -------------------------------------------------------------------------- */
 
 using namespace tamaas;
 
 template <typename T>
 struct AddOneInplace {
   CUDA_LAMBDA void operator()(T& x) { x += 1; }
 };
 
 // Testing loops on one grid
 TEST(TestLoops, OneArgument) {
   Grid<Real, 1> grid({20}, 1);
   Grid<Real, 1> solution({20}, 1);
 
   auto add_one = [] (Real & x) { return x + 1; };
 
   std::iota(grid.begin(), grid.end(), 1);
 
   // Makeing solution
   std::transform(grid.begin(), grid.end(), solution.begin(), add_one);
 
   auto add_one_inplace = AddOneInplace<Real>();
 
   Loop::loop(add_one_inplace, grid);
 
   ASSERT_TRUE(compare(grid, solution, AreFloatEqual())) << "One argument loop failed";
 }
 
 struct PrimalTest {
   CUDA_LAMBDA void operator()(Int& primal, Int& val) {
     val = (primal > 0)? -1:1;
   }
 };
 
 // Testing loops on two grids
 TEST(TestLoops, TwoArguments) {
   // Why no ints?
   Grid<Int, 2> grid({20, 20}, 1);
   Grid<Int, 2> primal({20, 20}, 1);
   Grid<Int, 2> solution({20, 20}, 1);
 
   primal(0, 0) = 1; primal(0, 1) = 1;
   primal(1, 0) = 1; primal(1, 1) = 1;
 
   std::transform(primal.begin(), primal.end(), solution.begin(),
 		 [] (Int & primal) {
 		   return (primal > 0)? -1:1;
 		 });
 
   auto primal_test = PrimalTest();
 
   Loop::loop(primal_test, primal, grid);
 
   ASSERT_TRUE(compare(solution, grid)) << "Two argument loop failed";
 }
 
 struct AssignUInt {
   CUDA_LAMBDA void operator()(UInt& x, UInt i) {
     x = i;
   }
 };
 
 // Testing an enumeration
 TEST(TestLoops, Enumeration) {
   Grid<UInt, 1> grid({100}, 1);
   Grid<UInt, 1> solution({100}, 1);
   std::iota(solution.begin(), solution.end(), 0);
 
   auto assign_uint = AssignUInt();
 
   Loop::loop(assign_uint, grid, Loop::arange<UInt>(100));
 
   ASSERT_TRUE(compare(solution, grid)) << "Enumeration loop failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 struct Identity {
   CUDA_LAMBDA UInt operator()(UInt& x) const { return x; }
 };
 
 // Testing one grid reductions
 TEST(TestReductions, OneArgument) {
   Grid<UInt, 1> grid({6}, 1);
   std::iota(grid.begin(), grid.end(), 1);
 
   const auto id = Identity();
 
   // Sum reduction
   UInt sol = std::accumulate(grid.begin(), grid.end(), 0, std::plus<UInt>());
   UInt red = Loop::reduce<operation::plus>(id, grid);
   ASSERT_TRUE(sol == red) << "Addition reduction failed on one argument";
 
   // Product reduction
   sol = std::accumulate(grid.begin(), grid.end(), 1, std::multiplies<UInt>());
   red = Loop::reduce<operation::times>(id, grid);
   ASSERT_TRUE(sol == red) << "Multiplication reduction failed on one argument";
 
   // Min reduction
   sol = *std::min_element(grid.begin(), grid.end());
   red = Loop::reduce<operation::min>(id, grid);
   ASSERT_TRUE(sol == red) << "Min reduction failed on one argument";
 
   // Max reduction
   sol = *std::max_element(grid.begin(), grid.end());
   red = Loop::reduce<operation::max>(id, grid);
   ASSERT_TRUE(sol == red) << "Max reduction failed on one argument";
 }
 
 struct PrimalReduce {
   CUDA_LAMBDA UInt operator()(UInt& p, UInt& val) {
     return (p > 0)? val:0;
   }
 };
 
 TEST(TestReductions, TwoArguments) {
   Grid<UInt, 1> grid({20}, 1);
   Grid<UInt, 1> primal({20}, 1);
   grid = 1;
   primal(0) = 1;
   primal(1) = 1;
 
   auto primal_reduce = PrimalReduce();
 
   // Reduce on values where primal > 0
   UInt red = Loop::reduce<operation::plus>(primal_reduce, primal, grid);
   ASSERT_TRUE(red == 2) << "Two args reduction failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 TEST(TestStridedLoops, OneArgument) {
   Grid<UInt, 2> grid({10, 10}, 2);
   std::iota(grid.begin(), grid.end(), 1);
 
   Grid<UInt, 2> solution({10, 10}, 2);
   solution = grid;
 
   std::for_each(solution.begin(), solution.end(), [] (UInt & x) {
       if (x % 2 == 1)
 	x += 1;
     });
 
   auto add_one_inplace = AddOneInplace<UInt>();
 
   Loop::stridedLoop(add_one_inplace, grid);
 
   ASSERT_TRUE(compare(solution, grid)) << "One argument strided loop failed";
 }
 
 template <typename T> using WrapVector = VectorProxy<T, 2>;
 
 struct AddOneVector {
   CUDA_LAMBDA void operator()(WrapVector<UInt>&& x) {
     x(0) += 1;
   }
 };
 
 TEST(TestStridedLoops, VectorStride) {
   Grid<UInt, 2> grid({10, 10}, 2);
   std::iota(grid.begin(), grid.end(), 1);
 
   Grid<UInt, 2> solution({10, 10}, 2);
   solution = grid;
 
   std::for_each(solution.begin(), solution.end(), [] (UInt & x) {
       if (x % 2 == 1)
 	x += 1;
     });
 
   auto add_one_inplace = AddOneVector();
 
   Loop::stridedLoop(add_one_inplace, grid);
 
   ASSERT_TRUE(compare(solution, grid)) << "Static vector strided loop failed";
 }
 
 template <typename T> using WrapMatrix = MatrixProxy<T, 2, 2>;
 
 struct SetOneMatrix {
   CUDA_LAMBDA void operator()(WrapMatrix<UInt>&& x) {
     x(0, 0) = 1;
     x(1, 1) = 1;
   }
 };
 
 TEST(TestStridedLoops, MatrixStride) {
   Grid<UInt, 2> grid({10, 10}, 4);
   Grid<UInt, 2> solution({10, 10}, 4);
   std::iota(solution.begin(), solution.end(), 0);
 
   std::for_each(solution.begin(), solution.end(), [] (UInt & x) {
       if (x % 4 == 0 || x % 4 == 3)
 	x = 1;
       else
 	x = 0;
     });
 
   auto set_one = SetOneMatrix();
 
   Loop::stridedLoop(set_one, grid);
 
   ASSERT_TRUE(compare(solution, grid)) << "Static matrix strided loop failed";
 }
 
 TEST(TestStridedReduction, VectorReduce) {
   Grid<UInt, 2> grid({10, 10}, 3);
 
   Loop::stridedLoop([](VectorProxy<UInt, 3>&& v) {
       v(0) = 1;
       v(1) = 2;
       v(2) = 3;
     }, grid);
 
-  auto res = Loop::stridedReduce<operation::plus>([](VectorProxy<UInt, 3>&& v) -> Vector<UInt, 3> {
-      return v;
-    }, grid);
+  auto res = Loop::stridedReduce<operation::plus>(
+      [](VectorProxy<UInt, 3>&& v) -> Vector<UInt, 3> { return v; }, grid);
 
   ASSERT_EQ(res(0), 100);
   ASSERT_EQ(res(1), 200);
   ASSERT_EQ(res(2), 300);
 }