diff --git a/src/cell/CMakeLists.txt b/src/cell/CMakeLists.txt
index 3dba371..86afb67 100644
--- a/src/cell/CMakeLists.txt
+++ b/src/cell/CMakeLists.txt
@@ -1,35 +1,36 @@
 # =============================================================================
 # file   CMakeLists.txt
 #
 # @author Till Junge <till.junge@epfl.ch>
 #
 # @date   08 Jan 2018
 #
 # @brief  configuration for cell implementations
 #
 # @section LICENSE
 #
 # Copyright © 2018 Till Junge
 #
 # µSpectre is free software; you can redistribute it and/or
 # modify it under the terms of the GNU General Public License as
 # published by the Free Software Foundation, either version 3, or (at
 # your option) any later version.
 #
 # µSpectre 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
 # General Public License for more details.
 #
 # You should have received a copy of the GNU General Public License
 # along with GNU Emacs; see the file COPYING. If not, write to the
 # Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 # Boston, MA 02111-1307, USA.
 # =============================================================================
 
 set (cells_SRC
   ${CMAKE_CURRENT_SOURCE_DIR}/cell_base.cc
   ${CMAKE_CURRENT_SOURCE_DIR}/cell_split.cc
+  ${CMAKE_CURRENT_SOURCE_DIR}/cell_laminate.cc
   )
 
 target_sources(muSpectre PRIVATE ${cells_SRC})
diff --git a/src/cell/cell_base.hh b/src/cell/cell_base.hh
index 5f7046e..0580bf7 100644
--- a/src/cell/cell_base.hh
+++ b/src/cell/cell_base.hh
@@ -1,523 +1,524 @@
 /**
  * @file   cell_base.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   01 Nov 2017
  *
  * @brief Base class representing a unit cell cell with single
  *        projection operator
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre 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
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #ifndef CELL_BASE_H
 #define CELL_BASE_H
 
 #include "common/common.hh"
 #include "common/ccoord_operations.hh"
 #include "common/field.hh"
 #include "common/utilities.hh"
 #include "materials/material_base.hh"
 #include "fft/projection_base.hh"
 #include "cell/cell_traits.hh"
 
 #include <vector>
 #include <memory>
 #include <tuple>
 #include <functional>
 #include <array>
 
 namespace muSpectre {
   /**
    * Cell adaptors implement the matrix-vector multiplication and
    * allow the system to be used like a sparse matrix in
    * conjugate-gradient-type solvers
    */
   template <class Cell>
   class CellAdaptor;
 
   /**
    * Base class for cells that is not templated and therefore can be
    * in solvers that see cells as runtime-polymorphic objects. This
    * allows the use of standard
    * (i.e. spectral-method-implementation-agnostic) solvers, as for
    * instance the scipy solvers
    */
 
   class Cell
   {
   public:
     //! sparse matrix emulation
     using Adaptor = CellAdaptor<Cell>;
 
     //! dynamic vector type for interactions with numpy/scipy/solvers etc.
     using Vector_t = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
 
     //! dynamic matrix type for setting strains
     using Matrix_t = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
 
     //! ref to constant vector
     using ConstVector_ref = Eigen::Map<const Vector_t>;
 
     //! output vector reference for solvers
     using Vector_ref = Eigen::Map<Vector_t>;
 
     //! Default constructor
     Cell() = default;
 
     //! Copy constructor
     Cell(const Cell &other) = default;
 
     //! Move constructor
     Cell(Cell &&other) = default;
 
     //! Destructor
     virtual ~Cell()  = default;
 
     //! Copy assignment operator
     Cell& operator=(const Cell &other) = default;
 
     //! Move assignment operator
     Cell& operator=(Cell &&other) = default;
 
     //! for handling double initialisations right
     bool is_initialised() const {return this->initialised;}
 
     //! returns the number of degrees of freedom in the cell
     virtual Dim_t get_nb_dof() const = 0;
 
     //! return the communicator object
     virtual const Communicator & get_communicator() const = 0;
 
     /**
      * formulation is hard set by the choice of the projection class
      */
     virtual const Formulation & get_formulation() const = 0;
 
     /**
      * returns the material dimension of the problem
      */
     virtual Dim_t get_material_dim() const = 0;
 
     /**
      * returns the number of rows and cols for the strain matrix type
      * (for full storage, the strain is stored in material_dim ×
      * material_dim matrices, but in symmetriy storage, it is a column
      * vector)
      */
     virtual std::array<Dim_t, 2> get_strain_shape() const = 0;
 
     /**
      * returns a writable map onto the strain field of this cell. This
      * corresponds to the unknowns in a typical solve cycle.
      */
     virtual Vector_ref get_strain_vector() = 0;
 
     /**
      * returns a read-only map onto the stress field of this
      * cell. This corresponds to the intermediate (and finally, total)
      * solution in a typical solve cycle
      */
     virtual ConstVector_ref get_stress_vector() const = 0;
 
 
     /**
      * evaluates and returns the stress for the currently set strain
      */
     virtual ConstVector_ref evaluate_stress() = 0;
 
     /**
      * evaluates and returns the stress and stiffness for the currently set strain
      */
     virtual std::array<ConstVector_ref, 2> evaluate_stress_tangent() = 0;
 
 
     /**
      * applies the projection operator in-place on the input vector
      */
     virtual void apply_projection(Eigen::Ref<Vector_t> vec) = 0;
 
     /**
      * freezes all the history variables of the materials
      */
     virtual void save_history_variables() = 0;
 
     /**
      * evaluates the directional and projected stiffness (this
      * corresponds to G:K:δF in de Geus 2017,
      * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that
      * this operation needs to be implemented with a copy in oder to
      * be compatible with scipy and EigenCG etc (At the very least,
      * the copy is only made once)
      */
     virtual Vector_ref evaluate_projected_directional_stiffness
       (Eigen::Ref<const Vector_t> delF) = 0;
 
 
     /**
      * set uniform strain (typically used to initialise problems
      */
     virtual void set_uniform_strain(const Eigen::Ref<const Matrix_t> &) = 0;
 
     //! get a sparse matrix view on the cell
     virtual Adaptor get_adaptor() = 0;
   protected:
     bool initialised{false}; //!< to handle double initialisation right
 
   private:
   };
   //! DimS spatial dimension (dimension of problem
   //! DimM material_dimension (dimension of constitutive law)
   template <Dim_t DimS, Dim_t DimM=DimS>
   class CellBase: public Cell
   {
   public:
     using Parent = Cell;
     using Ccoord = Ccoord_t<DimS>; //!< cell coordinates type
     using Rcoord = Rcoord_t<DimS>; //!< physical coordinates type
     //! global field collection
     using FieldCollection_t = GlobalFieldCollection<DimS>;
     //! the collection is handled in a `std::unique_ptr`
     using Collection_ptr = std::unique_ptr<FieldCollection_t>;
     //! polymorphic base material type
     using Material_t = MaterialBase<DimS, DimM>;
     //! materials handled through `std::unique_ptr`s
     using Material_ptr = std::unique_ptr<Material_t>;
     //! polymorphic base projection type
     using Projection_t = ProjectionBase<DimS, DimM>;
     //! projections handled through `std::unique_ptr`s
     using Projection_ptr = std::unique_ptr<Projection_t>;
     //! expected type for strain fields
     using StrainField_t =
       TensorField<FieldCollection_t, Real, secondOrder, DimM>;
     //! expected type for stress fields
     using StressField_t =
       TensorField<FieldCollection_t, Real, secondOrder, DimM>;
     //! expected type for tangent stiffness fields
     using TangentField_t =
       TensorField<FieldCollection_t, Real, fourthOrder, DimM>;
     //! combined stress and tangent field
     using FullResponse_t =
       std::tuple<const StressField_t&, const TangentField_t&>;
     //! iterator type over all cell pixel's
     using iterator = typename CcoordOps::Pixels<DimS>::iterator;
 
     //! dynamic vector type for interactions with numpy/scipy/solvers etc.
     using Vector_t = typename Parent::Vector_t;
 
     //! ref to constant vector
     using ConstVector_ref = typename Parent::ConstVector_ref;
 
     //! output vector reference for solvers
     using Vector_ref = typename Parent::Vector_ref;
 
 
     //! sparse matrix emulation
     using Adaptor = Parent::Adaptor;
 
     //! Default constructor
     CellBase() = delete;
 
     //! constructor using sizes and resolution
     CellBase(Projection_ptr projection);
 
     //! Copy constructor
     CellBase(const CellBase &other) = delete;
 
     //! Move constructor
     CellBase(CellBase &&other);
 
     //! Destructor
     virtual ~CellBase() = default;
 
     //! Copy assignment operator
     CellBase& operator=(const CellBase &other) = delete;
 
     //! Move assignment operator
     CellBase& operator=(CellBase &&other) = default;
 
     /**
      * Materials can only be moved. This is to assure exclusive
      * ownership of any material by this cell
      */
     Material_t & add_material(Material_ptr mat);
 
 
     /**
      * returns a writable map onto the strain field of this cell. This
      * corresponds to the unknowns in a typical solve cycle.
      */
     virtual Vector_ref get_strain_vector() override;
 
     /**
      * returns a read-only map onto the stress field of this
      * cell. This corresponds to the intermediate (and finally, total)
      * solution in a typical solve cycle
      */
     virtual ConstVector_ref get_stress_vector() const override;
 
     /**
      * evaluates and returns the stress for the currently set strain
      */
     virtual ConstVector_ref evaluate_stress() override;
 
     /**
      * evaluates and returns the stress and stiffness for the currently set strain
      */
     virtual std::array<ConstVector_ref, 2> evaluate_stress_tangent() override;
 
     /**
      * evaluates the directional and projected stiffness (this
      * corresponds to G:K:δF in de Geus 2017,
      * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that
      * this operation needs to be implemented with a copy in oder to
      * be compatible with scipy and EigenCG etc. (At the very least,
      * the copy is only made once)
      */
     virtual Vector_ref evaluate_projected_directional_stiffness
       (Eigen::Ref<const Vector_t> delF) override;
 
     //! return the template param DimM (required for polymorphic use of `Cell`
     Dim_t get_material_dim() const override final {return DimM;}
 
     /**
      * returns the number of rows and cols for the strain matrix type
      * (for full storage, the strain is stored in material_dim ×
      * material_dim matrices, but in symmetriy storage, it is a column
      * vector)
      */
     std::array<Dim_t, 2> get_strain_shape() const override final;
 
     /**
      * applies the projection operator in-place on the input vector
      */
     void apply_projection(Eigen::Ref<Vector_t> vec) override final;
 
 
 
     /**
      * set uniform strain (typically used to initialise problems
      */
     void set_uniform_strain(const Eigen::Ref<const Matrix_t> &) override;
 
 
     /**
      * evaluate all materials
      */
     virtual FullResponse_t evaluate_stress_tangent(StrainField_t & F);
 
     /**
      * evaluate directional stiffness (i.e. G:K:δF or G:K:δε)
      */
 
     StressField_t & directional_stiffness(const TangentField_t & K,
                                           const StrainField_t & delF,
                                           StressField_t & delP);
     /**
      * vectorized version for eigen solvers, no copy, but only works
      * when fields have ArrayStore=false
      */
     Vector_ref directional_stiffness_vec(const Eigen::Ref<const Vector_t> & delF);
     /**
      * Evaluate directional stiffness into a temporary array and
      * return a copy. This is a costly and wasteful interface to
      * directional_stiffness and should only be used for debugging or
      * in the python interface
      */
     Eigen::ArrayXXd directional_stiffness_with_copy
       (Eigen::Ref<Eigen::ArrayXXd> delF);
 
     /**
      * Convenience function circumventing the neeed to use the
      * underlying projection
      */
     StressField_t & project(StressField_t & field);
 
     //! returns a ref to the cell's strain field
     StrainField_t & get_strain();
 
     //! returns a ref to the cell's stress field
     const StressField_t & get_stress() const;
 
     //! returns a ref to the cell's tangent stiffness field
     const TangentField_t & get_tangent(bool create = false);
 
 
 
     //! returns a ref to a temporary field managed by the cell
     StrainField_t & get_managed_field(std::string unique_name);
 
     /**
      * general initialisation; initialises the projection and
      * fft_engine (i.e. infrastructure) but not the materials. These
      * need to be initialised separately
      */
     void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate);
 
     /**
      * for materials with state variables, these typically need to be
      * saved/updated an the end of each load increment, this function
      * calls this update for each material in the cell
      */
     void save_history_variables() override final;
 
     iterator begin(); //!< iterator to the first pixel
     iterator end();  //!< iterator past the last pixel
     //! number of pixels in the cell
     size_t size() const {return pixels.size();}
 
     //! return the subdomain resolutions of the cell
     const Ccoord & get_subdomain_resolutions() const {
       return this->subdomain_resolutions;}
     //! return the subdomain locations of the cell
     const Ccoord & get_subdomain_locations() const {
       return this->subdomain_locations;}
     //! return the domain resolutions of the cell
     const Ccoord & get_domain_resolutions() const {
       return this->domain_resolutions;}
     //! return the sizes of the cell
     const Rcoord & get_domain_lengths() const {return this->domain_lengths;}
     //! return the real space coordinates of the center of pixel
     Rcoord get_pixel_coordinate(Ccoord & pixel) const;
     //! return the physical dimension of a pixel
     Rcoord get_pixel_lengths() const;
     //! check if the pixel is inside of the cell
     bool is_inside(Rcoord point);
     //! check if the point is inside of the cell
     bool is_inside(Ccoord pixel);
     //calculate the volume of each pixel:
     Real get_pixel_volume();
 
     /**
      * formulation is hard set by the choice of the projection class
      */
     const Formulation & get_formulation() const override final {
       return this->projection->get_formulation();}
 
     /**
      * get a reference to the projection object. should only be
      * required for debugging
      */
     Eigen::Map<Eigen::ArrayXXd> get_projection() {
       return this->projection->get_operator();}
 
     //! returns the spatial size
     constexpr static Dim_t get_sdim() {return DimS;};
 
     //! return a sparse matrix adaptor to the cell
     Adaptor get_adaptor() override;
     //! returns the number of degrees of freedom in the cell
     Dim_t get_nb_dof() const override {return this->size()*ipow(DimS, 2);};
 
     //! return the communicator object
     virtual const Communicator & get_communicator() const override {
       return this->projection->get_communicator();
     }
 
   protected:
     //! make sure that every pixel is assigned to one and only one material
     virtual void check_material_coverage();
 
     const Ccoord & subdomain_resolutions; //!< the cell's subdomain resolutions
     const Ccoord & subdomain_locations; //!< the cell's subdomain resolutions
     const Ccoord & domain_resolutions; //!< the cell's domain resolutions
     CcoordOps::Pixels<DimS> pixels; //!< helper to iterate over the pixels
     const Rcoord & domain_lengths; //!< the cell's lengths
     Collection_ptr fields; //!< handle for the global fields of the cell
     StrainField_t & F; //!< ref to strain field
     StressField_t & P; //!< ref to stress field
     //! Tangent field might not even be required; so this is an
     //! optional ref_wrapper instead of a ref
     optional<std::reference_wrapper<TangentField_t>> K{};
     //! container of the materials present in the cell
     std::vector<Material_ptr> materials{};
     Projection_ptr projection; //!< handle for the projection operator
     bool initialised{false}; //!< to handle double initialisation right
     SplittedCell is_cell_splitted {SplittedCell::no};
+    LaminateCell is_cell_laminate {LaminateCell::no}; 
   private:
   };
 
 
   /**
    * lightweight resource handle wrapping a `muSpectre::Cell` or
    * a subclass thereof into `Eigen::EigenBase`, so it can be
    * interpreted as a sparse matrix by Eigen solvers
    */
   template <class Cell>
   class CellAdaptor: public Eigen::EigenBase<CellAdaptor<Cell>> {
 
   public:
     using Scalar = double;     //!< sparse matrix traits
     using RealScalar = double; //!< sparse matrix traits
     using StorageIndex = int;  //!< sparse matrix traits
     enum {
       ColsAtCompileTime = Eigen::Dynamic,
       MaxColsAtCompileTime = Eigen::Dynamic,
       RowsAtCompileTime = Eigen::Dynamic,
       MaxRowsAtCompileTime = Eigen::Dynamic,
       IsRowMajor = false
     };
 
     //! constructor
     CellAdaptor(Cell & cell):cell{cell}{}
     //!returns the number of logical rows
     Eigen::Index rows() const {return this->cell.get_nb_dof();}
     //!returns the number of logical columns
     Eigen::Index cols() const {return this->rows();}
 
     //! implementation of the evaluation
     template<typename Rhs>
     Eigen::Product<CellAdaptor,Rhs,Eigen::AliasFreeProduct>
     operator*(const Eigen::MatrixBase<Rhs>& x) const {
       return Eigen::Product<CellAdaptor,Rhs,Eigen::AliasFreeProduct>
         (*this, x.derived());
     }
     Cell & cell; //!< ref to the cell
   };
 
 }  // muSpectre
 
 
 namespace Eigen {
   namespace internal {
     //! Implementation of `muSpectre::CellAdaptor` * `Eigen::DenseVector` through a
     //! specialization of `Eigen::internal::generic_product_impl`:
     template<typename Rhs, class CellAdaptor>
     struct generic_product_impl<CellAdaptor, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
       : generic_product_impl_base<CellAdaptor,Rhs,generic_product_impl<CellAdaptor,Rhs> >
     {
       //! undocumented
       typedef typename Product<CellAdaptor,Rhs>::Scalar Scalar;
 
       //! undocumented
       template<typename Dest>
       static void scaleAndAddTo(Dest& dst, const CellAdaptor& lhs, const Rhs& rhs, const Scalar& /*alpha*/)
       {
         // This method should implement "dst += alpha * lhs * rhs" inplace,
         // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it.
         // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs,
         dst.noalias() += const_cast<CellAdaptor&>(lhs).cell.
           evaluate_projected_directional_stiffness(rhs);
       }
     };
   }
 }
 
 
 #endif /* CELL_BASE_H */
diff --git a/src/cell/cell_laminate.cc b/src/cell/cell_laminate.cc
new file mode 100644
index 0000000..d78576a
--- /dev/null
+++ b/src/cell/cell_laminate.cc
@@ -0,0 +1,69 @@
+/**
+ * @file   cell_split.cc
+ *
+ * @author Ali Falsafi <ali.faslafi@epfl.ch>
+ *
+ * @date   10 Sep 2018
+ *
+ * @brief  Implementation for cell base class
+ *
+ * Copyright © 2017 Till Junge
+ *
+ * µSpectre is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation, either version 3, or (at
+ * your option) any later version.
+ *
+ * µSpectre 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
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with GNU Emacs; see the file COPYING. If not, write to the
+ * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+ * Boston, MA 02111-1307, USA.
+ */
+
+#include "cell/cell_laminate.hh"
+#include "common/ccoord_operations.hh"
+#include "common/iterators.hh"
+#include "common/tensor_algebra.hh"
+#include "common/common.hh"
+
+
+
+namespace muSpectre {
+
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  CellLaminate<DimS, DimM>::CellLaminate(Projection_ptr projection, LaminateCell is_cell_laminate)
+    :Parent(std::move(projection), SplittedCell::yes)
+  {
+    this->is_cell_laminate = is_cell_laminate;
+  }
+  /* ---------------------------------------------------------------------- */
+  template <Dim_t DimS, Dim_t DimM>
+  auto CellLaminate<DimS, DimM>::
+  evaluate_stress_tangent(StrainField_t & F)-> FullResponse_t
+  {
+    if (this->initialised == false) {
+      this->initialise();
+    }
+    //! High level compatibility checks
+    if (F.size() != this->F.size()) {
+      throw std::runtime_error("Size mismatch");
+    }
+    constexpr bool create_tangent{true};
+    this->get_tangent(create_tangent);
+    for (auto &mat : this->materials) {
+      auto mat_stress_tangent = mat->compute_stresses_tangent(F);
+    }
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template class CellLaminate<twoD, twoD>;
+  template class CellLaminate<threeD, threeD>;
+  /* ---------------------------------------------------------------------- */
+
+}
diff --git a/src/cell/cell_laminate.hh b/src/cell/cell_laminate.hh
new file mode 100644
index 0000000..713aa9b
--- /dev/null
+++ b/src/cell/cell_laminate.hh
@@ -0,0 +1,105 @@
+/**
+ * @file   cell_split.hh
+ *
+ * @author Ali Falsafi <ali.falsafi@epfl.ch>
+ *
+ * @date   10 Sep 2018
+ *
+ * @brief Base class representing a unit cell able to handle
+ *        split material assignments
+ *
+ * Copyright © 2017 Till Junge
+ *
+ * µSpectre is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation, either version 3, or (at
+ * your option) any later version.
+ *
+ * µSpectre 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
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with GNU Emacs; see the file COPYING. If not, write to the
+ * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+ * Boston, MA 02111-1307, USA.
+ */
+
+
+
+#ifndef CELL_LAMINATE_H
+#define CELL_LAMINATE_H
+
+#include "common/common.hh"
+#include "common/ccoord_operations.hh"
+#include "common/field.hh"
+#include "common/utilities.hh"
+#include "materials/material_base.hh"
+#include "fft/projection_base.hh"
+#include "cell/cell_traits.hh"
+#include "cell/cell_base.hh"
+#include "cell/cell_split.hh"
+
+
+namespace muSpectre {
+  //! DimS spatial dimension (dimension of problem
+  //! DimM material_dimension (dimension of constitutive law)
+  /* This class handles the cells that has splitly assigned material to their pixels */
+  template <Dim_t DimS, Dim_t DimM=DimS>
+  class CellLaminate:public CellSplit<DimS, DimM> {
+  public:
+    using Parent = CellSplit<DimS, DimM>;
+    //! global field collection
+    using FieldCollection_t = GlobalFieldCollection<DimS>;
+    using Projection_t = ProjectionBase<DimS, DimM>;
+    //! projections handled through `std::unique_ptr`s
+    using Projection_ptr = std::unique_ptr<Projection_t>;
+    using StrainField_t =
+      TensorField<FieldCollection_t, Real, secondOrder, DimM>;
+    //! expected type for stress fields
+    using StressField_t =
+      TensorField<FieldCollection_t, Real, secondOrder, DimM>;
+    //! expected type for tangent stiffness fields
+    using TangentField_t =
+      TensorField<FieldCollection_t, Real, fourthOrder, DimM>;
+    //! combined stress and tangent field
+    using FullResponse_t =
+      std::tuple<const StressField_t&, const TangentField_t&>;
+
+    //! ref to constant vector
+    using ConstVector_ref = typename Parent::ConstVector_ref;
+    //! Default constructor
+    CellLaminate() = delete;
+
+    //! constructor using sizes and resolution
+    CellLaminate(Projection_ptr projection,
+                 LaminateCell is_cell_laminate= LaminateCell::yes);
+
+    //! Copy constructor
+    CellLaminate(const CellLaminate &other) = delete;
+
+    //! Move constructor
+    CellLaminate(CellLaminate &&other) = default;
+
+    //! Destructor
+    virtual ~CellLaminate() = default;
+
+    //! Copy assignment operator
+    CellLaminate& operator=(const CellLaminate &other) = delete;
+
+    //! Move assignment operator
+    CellLaminate& operator=(CellLaminate &&other) = default;
+    /* ---------------------------------------------------------------------- */
+  protected:
+    //full resppnse is composed of the stresses and tangent matrix
+    FullResponse_t evaluate_stress_tangent(StrainField_t & F) override;
+    std::array<ConstVector_ref, 2> evaluate_stress_tangent() override;
+    ConstVector_ref evaluate_stress() override;
+
+  };
+
+
+}
+
+#endif /* CELL_LAMINATE_H */
diff --git a/src/cell/cell_split.cc b/src/cell/cell_split.cc
index 8deaee3..ff8c305 100644
--- a/src/cell/cell_split.cc
+++ b/src/cell/cell_split.cc
@@ -1,315 +1,315 @@
 /**
  * @file   cell_split.cc
  *
  * @author Ali Falsafi <ali.faslafi@epfl.ch>
  *
  * @date   19 Apr 2018
  *
  * @brief  Implementation for cell base class
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre 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
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 #include "common/common.hh"
 #include "common/ccoord_operations.hh"
 #include "common/field.hh"
 #include "common/utilities.hh"
 #include "materials/material_base.hh"
 #include "fft/projection_base.hh"
 #include "cell/cell_traits.hh"
 #include "cell/cell_base.hh"
 #include "cell/cell_split.hh"
 
 #include <vector>
 #include <memory>
 #include <tuple>
 #include <sstream>
 #include <algorithm>
 
 namespace muSpectre {
 
-    /* ---------------------------------------------------------------------- */
+  /* ---------------------------------------------------------------------- */
 
   template <Dim_t DimS, Dim_t DimM>
   CellSplit<DimS, DimM>::CellSplit(Projection_ptr projection, SplittedCell is_cell_splitted )
     :Parent(std::move(projection))
   {
     this->is_cell_splitted = is_cell_splitted;
   }
 
   /* ---------------------------------------------------------------------- */
 
   template <Dim_t DimS, Dim_t DimM>
   std::vector<Real> CellSplit<DimS, DimM>::get_assigned_ratios(){
     auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions);
     std::vector<Real> pixel_assigned_ratios(nb_pixels, 0.0);
     for (auto & mat: this->materials) {
       for (auto & pixel: *mat) {
         auto index = CcoordOps::get_index(this->subdomain_resolutions,
                                           this->subdomain_locations,
                                           pixel);
         pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel);
       }
     }
     return pixel_assigned_ratios;
   }
 
    /* ---------------------------------------------------------------------- */
 
   template <Dim_t DimS, Dim_t DimM>
   std::vector<Real> CellSplit<DimS, DimM>::get_unassigned_ratios_incomplete_pixels(){
     auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions);
     std::vector<Real> pixel_assigned_ratios(nb_pixels, 0.0);
     std::vector<Real> pixel_assigned_ratios_incomplete_pixels;
     for (auto & mat: this->materials) {
       for (auto & pixel: *mat) {
         auto index = CcoordOps::get_index(this->subdomain_resolutions,
                                           this->subdomain_locations,
                                           pixel);
         pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel);
       }
     }
     for (auto && ratio:pixel_assigned_ratios){
       if (ratio<1){
         pixel_assigned_ratios_incomplete_pixels.push_back(1.0-ratio);
       }
     }
     return pixel_assigned_ratios_incomplete_pixels;
   }
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   auto CellSplit<DimS, DimM>::make_incomplete_pixels() -> IncompletePixels{
     return IncompletePixels(*this);
   }
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   std::vector<int> CellSplit<DimS, DimM>::get_index_incomplete_pixels(){
     auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions);
     std::vector<Real> pixel_assigned_ratios(nb_pixels, 0.0);
     std::vector<int> index_unassigned_pixels;
     for (auto & mat: this->materials) {
       for (auto & pixel: *mat) {
         auto index = CcoordOps::get_index(this->subdomain_resolutions,
                                           this->subdomain_locations,
                                           pixel);
         pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel);
       }
     }
     for (auto && tup: akantu::enumerate(this->pixels)){
       auto && i{std::get<0> (tup)};
       if (pixel_assigned_ratios[i] < 1){
         index_unassigned_pixels.push_back(i);
       }
     }
     return index_unassigned_pixels;
   }
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   std::vector<Ccoord_t<DimS>> CellSplit<DimS, DimM>::get_unassigned_pixels(){
     auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions);
   std::vector<Real> pixel_assigned_ratios(nb_pixels, 0.0);
   std::vector<Ccoord_t<DimS>> unassigned_pixels;
   for (auto & mat: this->materials) {
     for (auto & pixel: *mat) {
       auto index = CcoordOps::get_index(this->subdomain_resolutions,
                                         this->subdomain_locations,
                                         pixel);
         pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel);
       }
     }
   for (auto && tup: akantu::enumerate(this->pixels)){
     auto && index{std::get<0> (tup)};
     auto && pix{std::get<1> (tup)};
     if (pixel_assigned_ratios[index] < 1){
       unassigned_pixels.push_back(pix);
     }
   }
   return unassigned_pixels;
 }
 
   /* ---------------------------------------------------------------------- */
   /*template <Dim_t DimS, Dim_t DimM>
   template<class MaterialType>
   void CellSplit<DimS, DimM>::complete_material_assignment(MaterialType& material){
     std::vector<Real> pixel_assigned_ratio(this->get_assigned_ratios());
     for (auto && tup: akantu::enumerate(*this)) {
       auto && pixel = std::get<1>(tup);
       auto iterator = std::get<0>(tup);
       if (pixel_assigned_ratio[iterator] < 1.0){
         material.add_pixel_split(pixel, 1.0-pixel_assigned_ratio[iterator]);
       }
     }
   }*/
 
 
   template <Dim_t DimS, Dim_t DimM>
   void CellSplit<DimS, DimM>::check_material_coverage(){
     auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions);
     std::vector<Real> pixel_assigned_ratio(nb_pixels, 0.0);
     for (auto & mat: this->materials) {
       for (auto & pixel: *mat) {
         auto index = CcoordOps::get_index(this->subdomain_resolutions,
                                           this->subdomain_locations,
                                           pixel);
         pixel_assigned_ratio.at(index) += mat->get_assigned_ratio(pixel);
       }
     }
     std::vector<Ccoord_t<DimM>> over_assigned_pixels;
     std::vector<Ccoord_t<DimM>> under_assigned_pixels;
     for (size_t i = 0; i < nb_pixels; ++i) {
       if (pixel_assigned_ratio.at(i) > 1.0) {
         over_assigned_pixels.push_back(CcoordOps::get_ccoord(this->subdomain_resolutions,
                                                              this->subdomain_locations, i));
       }else if (pixel_assigned_ratio.at(i) < 1.0) {
         under_assigned_pixels.push_back(CcoordOps::get_ccoord(this->subdomain_resolutions,
                                                               this->subdomain_locations, i));
       }
     }
      if (over_assigned_pixels.size() != 0) {
        std::stringstream err {};
        err << "Execesive material is assigned to the following pixels: ";
        for (auto & pixel: over_assigned_pixels) {
          err << pixel << ", ";
        }
        err << "and that cannot be handled";
        throw std::runtime_error(err.str());
      }
      if (under_assigned_pixels.size() != 0) {
        std::stringstream err {};
        err << "Insufficient material is assigned to the following pixels: ";
        for (auto & pixel: under_assigned_pixels) {
          err << pixel << ", ";
        }
        err << "and that cannot be handled";
        throw std::runtime_error(err.str());
      }
   }
-    /* ---------------------------------------------------------------------- */
+  /* ---------------------------------------------------------------------- */
   //this piece of code handles the evaluation of stress an dtangent matrix
   //in case the cells have materials in which pixels are partially composed of
   //diffferent materials.
 
   template <Dim_t DimS, Dim_t DimM>
   typename CellSplit<DimS, DimM>::FullResponse_t
   CellSplit<DimS, DimM>::evaluate_stress_tangent(StrainField_t & F){
     if (this->initialised == false) {
       this->initialise();
     }
     //! High level compatibility checks
     if (F.size() != this->F.size()) {
       throw std::runtime_error("Size mismatch");
     }
     constexpr bool create_tangent{true};
     this->get_tangent(create_tangent);
 
     // Here we should first set P and K matrixes equal to zeros first because we want to add up contribution
     //of the partial influence of different materials assigend to each pixel. Therefore, this values should be
     // initiialised as zero filled tensors
     this->set_p_k_zero();
 
       for (auto & mat: this->materials) {
         mat->compute_stresses_tangent(F, this->P, this->K.value(),
                                       this->get_formulation(), this->is_cell_splitted);
       }
       return std::tie(this->P, this->K.value());
   }
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   auto CellSplit<DimS, DimM>::
   evaluate_stress() -> ConstVector_ref {
     if (not this->initialised) {
       this->initialise();
     }
     this->P.set_zero();
     for (auto & mat: this->materials) {
       mat->compute_stresses(this->F, this->P, this->get_formulation());
     }
 
     return this->P.const_eigenvec();
   }
 /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   auto CellSplit<DimS, DimM>::
   evaluate_stress_tangent() -> std::array<ConstVector_ref, 2> {
     if (not this->initialised) {
       this->initialise();
     }
 
     constexpr bool create_tangent{true};
     this->get_tangent(create_tangent);
 
     // Here we should first set P and K matrixes equal to zeros first because we want to add up contribution
     //of the partial influence of different materials assigend to each pixel. Therefore, this values should be
     // initiialised as zero filled tensors
     this->set_p_k_zero();
     for (auto & mat: this->materials) {
       mat->compute_stresses_tangent(this->F, this->P, this->K.value(),
                                     this->get_formulation(), this->is_cell_splitted);
     }
     const TangentField_t & k = this->K.value();
     return std::array<ConstVector_ref, 2>{
       this->P.const_eigenvec(), k.const_eigenvec()};
 
   }
     /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void CellSplit<DimS, DimM>::set_p_k_zero(){
     this->P.set_zero();
     this->K.value().get().set_zero();
   }
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   CellSplit<DimS, DimM>::IncompletePixels::IncompletePixels(CellSplit<DimS, DimM>& cell)
     :cell(cell),incomplete_assigned_ratios(cell.get_unassigned_ratios_incomplete_pixels()),
     index_incomplete_pixels(cell.get_index_incomplete_pixels())
   {}
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   CellSplit<DimS, DimM>::IncompletePixels::iterator::iterator
     (const IncompletePixels & pixels, bool begin)
     : incomplete_pixels(pixels),
       index{begin? 0:this->incomplete_pixels.index_incomplete_pixels.size()}
   {}
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   auto CellSplit<DimS, DimM>::IncompletePixels::iterator::operator++()->iterator&{
     this->index++;
     return *this;
   }
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   auto CellSplit<DimS, DimM>::IncompletePixels::iterator::operator!=(const iterator & other)->bool{
     //        return (this->incomplete_pixels.index_incomplete_pixels[index] !=
     //                other.incomplete_pixels.index_incomplete_pixels[index]);
     return (this->index != other.index);
   }
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   auto CellSplit<DimS, DimM>::IncompletePixels::iterator::operator*()->value_type const{
     auto ccoord = CcoordOps::get_ccoord(this->incomplete_pixels.cell.get_domain_resolutions(),
                                         this->incomplete_pixels.cell.get_subdomain_locations(),
                                         this->incomplete_pixels.index_incomplete_pixels[index]);
     auto ratio = this->incomplete_pixels.incomplete_assigned_ratios[index];
     return std::make_tuple(ccoord, ratio);
   }
   /* ---------------------------------------------------------------------- */
   template class CellSplit<twoD, twoD>;
   template class CellSplit<threeD, threeD>;
   /* ---------------------------------------------------------------------- */
 
 } //muspectre
diff --git a/src/cell/cell_split.hh b/src/cell/cell_split.hh
index 808ac8d..ac354b2 100644
--- a/src/cell/cell_split.hh
+++ b/src/cell/cell_split.hh
@@ -1,197 +1,215 @@
 /**
  * @file   cell_split.hh
  *
  * @author Ali Falsafi <ali.falsafi@epfl.ch>
  *
  * @date   19 Apr 2018
  *
  * @brief Base class representing a unit cell able to handle
  *        split material assignments
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre 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
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #ifndef CELL_SPLIT_H
 #define CELL_SPLIT_H
 
 #include "common/common.hh"
 #include "common/ccoord_operations.hh"
 #include "common/field.hh"
 #include "common/utilities.hh"
 #include "materials/material_base.hh"
 #include "fft/projection_base.hh"
 #include "cell/cell_traits.hh"
 #include "cell/cell_base.hh"
 #include "common/intersection_octree.hh"
 
 #include <vector>
 #include <memory>
 #include <tuple>
 #include <functional>
 namespace muSpectre {
 
   //! DimS spatial dimension (dimension of problem
   //! DimM material_dimension (dimension of constitutive law)
-  /* This class handles the cells that has splitly assigned material to their pixels */
+  /* This class handles the cells that has
+     splitly assigned material to their pixels */
   template <Dim_t DimS, Dim_t DimM=DimS>
   class CellSplit: public CellBase<DimS, DimM> {
     friend class CellBase<DimS, DimM>;
   public:
     using Parent = CellBase<DimS, DimM>; //!< base class
     //! global field collection
     using FieldCollection_t = GlobalFieldCollection<DimS>;
     using Projection_t = ProjectionBase<DimS, DimM>;
     //! projections handled through `std::unique_ptr`s
     using Projection_ptr = std::unique_ptr<Projection_t>;
     using StrainField_t =
       TensorField<FieldCollection_t, Real, secondOrder, DimM>;
     //! expected type for stress fields
     using StressField_t =
       TensorField<FieldCollection_t, Real, secondOrder, DimM>;
     //! expected type for tangent stiffness fields
     using TangentField_t =
       TensorField<FieldCollection_t, Real, fourthOrder, DimM>;
     //! combined stress and tangent field
     using FullResponse_t =
       std::tuple<const StressField_t&, const TangentField_t&>;
 
     //! ref to constant vector
     using ConstVector_ref = typename Parent::ConstVector_ref;
 
 
     //! Default constructor
     CellSplit() = delete;
 
     //! constructor using sizes and resolution
-    CellSplit(Projection_ptr projection, SplittedCell is_cell_splitted = SplittedCell::yes);
+    CellSplit(Projection_ptr projection,
+              SplittedCell is_cell_splitted = SplittedCell::yes);
 
     //! Copy constructor
     CellSplit(const CellSplit &other) = delete;
 
     //! Move constructor
     CellSplit(CellSplit &&other) = default;
 
     //! Destructor
     virtual ~CellSplit() = default;
 
     //! Copy assignment operator
     CellSplit& operator=(const CellSplit &other) = delete;
 
     //! Move assignment operator
     CellSplit& operator=(CellSplit &&other) = default;
 
-
-    //completes the assignmnet of material with a specific material so all the under-assigned pixels would be assigned to a material
+    /**
+     *completes the assignmnet of material with a specific material so
+     *all the under-assigned pixels would be assigned to a material.
+     */
+    //
     void complete_material_assignment(MaterialBase<DimS, DimM>& material);
 
     //returns the assigend ratios to each pixel
     std::vector<Real> get_assigned_ratios();
 
     //
     //template<class MaterialType>
-    void make_automatic_precipitate(std::vector<Rcoord_t<DimS>> preciptiate_vertices, MaterialBase<DimS, DimM>& material);
+    void make_automatic_precipitate(std::vector<Rcoord_t<DimS>>
+                                    preciptiate_vertices,
+                                    MaterialBase<DimS, DimM>& material);
     //
     std::vector<Real> get_unassigned_ratios_incomplete_pixels();
     std::vector<int> get_index_incomplete_pixels();
     std::vector<Ccoord_t<DimS>> get_unassigned_pixels();
     class IncompletePixels{
     public:
       //! constructor
       IncompletePixels(CellSplit<DimS, DimM>& cell);
       //! copy constructor
       IncompletePixels(const IncompletePixels & other) = default;
       //! move constructor
       IncompletePixels(IncompletePixels & other) = default;
       //Deconstructor
       virtual ~IncompletePixels() = default;
 
       //! iterator type over all incompletetedly assigned pixel's
       class iterator{
       public:
         using value_type = std::tuple<Ccoord_t<DimS>, Real>; //!< stl conformance
         using const_value_type = const value_type; //!< stl conformance
         using pointer = value_type*; //!< stl conformance
         using difference_type = std::ptrdiff_t; //!< stl conformance
         using iterator_category = std::forward_iterator_tag;//!<stl conformance
         using reference = value_type; //!< stl conformance
         //! constructor
         iterator(const IncompletePixels & pixels, bool begin = true);
         //deconstructor
         virtual ~iterator() = default;
         //! dereferencing
         auto operator*()->value_type const;
         //! pre-increment
         auto operator++()->iterator&;
         //! inequality
         auto operator!=(const iterator & other)->bool;
         //! equality
         inline bool operator==(const iterator & other) const;
       private:
         const IncompletePixels& incomplete_pixels;
         size_t index;
       };
       inline iterator begin() const {return iterator(*this);}
       //! stl conformance
       inline iterator end() const {return iterator(*this, false);}
       //! stl conformance
-      inline size_t size() const {return CcoordOps::get_size(this->cell.get_domain_resolutions());}
+      inline size_t size() const {
+        return CcoordOps::get_size(this->cell.get_domain_resolutions());
+      }
     private:
       CellSplit<DimS, DimM>& cell;
       std::vector<Real> incomplete_assigned_ratios;
       std::vector<int> index_incomplete_pixels;
     };
 
     auto make_incomplete_pixels() -> IncompletePixels;
   protected:
     void check_material_coverage() override final;
     void set_p_k_zero();
-    //full resppnse is consisted of the stresses and tangent matrix
-    FullResponse_t evaluate_stress_tangent(StrainField_t & F) override final;
-    std::array<ConstVector_ref, 2> evaluate_stress_tangent() override final;
-    ConstVector_ref evaluate_stress() override final;
+    //full resppnse is composed of the stresses and tangent matrix
+    FullResponse_t evaluate_stress_tangent(StrainField_t & F) override;
+    std::array<ConstVector_ref, 2> evaluate_stress_tangent() override;
+    ConstVector_ref evaluate_stress() override;
   private:
   };
  /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
-  void CellSplit<DimS, DimM>::make_automatic_precipitate(std::vector<Rcoord_t<DimS>> precipitate_vertices, MaterialBase<DimS, DimM>& material){
-    RootNode<DimS> precipitate (*this, precipitate_vertices);
-    auto precipitate_intersects= precipitate.get_intersected_pixels();
-    auto precipitate_intersection_ratios = precipitate.get_intersection_ratios();
-    for (auto tup:akantu::zip(precipitate_intersects, precipitate_intersection_ratios)){
+  void CellSplit<DimS, DimM>::make_automatic_precipitate
+  (std::vector<Rcoord_t<DimS>> precipitate_vertices,
+   MaterialBase<DimS, DimM>& material){
+    RootNode<DimS> precipitate (*this,
+                                precipitate_vertices);
+    auto precipitate_intersects=
+      precipitate.get_intersected_pixels();
+
+    auto precipitate_intersection_ratios =
+      precipitate.get_intersection_ratios();
+
+    for (auto tup:akantu::zip(precipitate_intersects,
+                              precipitate_intersection_ratios)){
       auto pix = std::get<0> (tup);
       auto ratio = std::get<1> (tup);
       material.add_pixel_split(pix,ratio);
     }
   }
  /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
-  void CellSplit<DimS, DimM>::complete_material_assignment(MaterialBase<DimS, DimM>& material){
+  void CellSplit<DimS, DimM>::complete_material_assignment
+  (MaterialBase<DimS, DimM>& material){
     std::vector<Real> pixel_assigned_ratio(this->get_assigned_ratios());
     for (auto && tup: akantu::enumerate(*this)) {
       auto && pixel = std::get<1>(tup);
       auto iterator = std::get<0>(tup);
       if (pixel_assigned_ratio[iterator] < 1.0){
         material.add_pixel_split(pixel, 1.0-pixel_assigned_ratio[iterator]);
       }
     }
   }
 
 } //muspectre
 
 
 #endif /* CELL_SPLIT_H */
diff --git a/src/common/common.hh b/src/common/common.hh
index d30229f..78defc4 100644
--- a/src/common/common.hh
+++ b/src/common/common.hh
@@ -1,317 +1,323 @@
 /**
  * @file   common.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   01 May 2017
  *
  * @brief  Small definitions of commonly used types throughout µSpectre
  *
  * @section  LICENSE
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre 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
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include <array>
 #include <cmath>
 #include <complex>
 #include <iostream>
 #include <type_traits>
 #include <string>
 
 #ifndef COMMON_H
 #define COMMON_H
 
 
 namespace muSpectre {
 
   /**
    * Eigen uses signed integers for dimensions. For consistency,
    µSpectre uses them througout the code. needs to represent -1 for
    eigen
    */
   using Dim_t = int;
 
   constexpr Dim_t oneD{1}; //!< constant for a one-dimensional problem
   constexpr Dim_t twoD{2}; //!< constant for a two-dimensional problem
   constexpr Dim_t threeD{3}; //!< constant for a three-dimensional problem
   constexpr Dim_t firstOrder{1}; //!< constant for vectors
   constexpr Dim_t secondOrder{2}; //!< constant second-order tensors
   constexpr Dim_t fourthOrder{4}; //!< constant fourth-order tensors
 
   //@{
   //! @anchor scalars
   //! Scalar types used for mathematical calculations
   using Uint = unsigned int;
   using Int = int;
   using Real = double;
   using Complex = std::complex<Real>;
   //@}
 
   //! Ccoord_t are cell coordinates, i.e. integer coordinates
   template<Dim_t dim>
   using Ccoord_t = std::array<Dim_t, dim>;
   //! Real space coordinates
   template<Dim_t dim>
   using Rcoord_t = std::array<Real, dim>;
 
   /**
    * Allows inserting `muSpectre::Ccoord_t` and `muSpectre::Rcoord_t`
    * into `std::ostream`s
    */
   template<typename T, size_t dim>
   std::ostream & operator << (std::ostream & os,
                               const std::array<T, dim> & index) {
     os << "(";
     for (size_t i = 0; i < dim-1; ++i) {
       os << index[i] << ", ";
     }
     os << index.back() << ")";
     return os;
   }
 
   //! element-wise division
   template <size_t dim>
   Rcoord_t<dim> operator/(const Rcoord_t<dim> & a, const Rcoord_t<dim> & b) {
     Rcoord_t<dim> retval{a};
     for (size_t i = 0; i < dim; ++i) {
       retval[i]/=b[i];
     }
     return retval;
   }
 
   //! element-wise division
   template <size_t dim>
   Rcoord_t<dim> operator/(const Rcoord_t<dim> & a, const Ccoord_t<dim> & b) {
     Rcoord_t<dim> retval{a};
     for (size_t i = 0; i < dim; ++i) {
       retval[i]/=b[i];
     }
     return retval;
   }
 
   //! convenience definitions
   constexpr Real pi{3.1415926535897932384626433};
 
   //! compile-time potentiation required for field-size computations
   template <typename R, typename I>
   constexpr R ipow(R base, I exponent) {
     static_assert(std::is_integral<I>::value, "Type must be integer");
     R retval{1};
     for (I i = 0; i < exponent; ++i) {
       retval *= base;
     }
     return retval;
   }
 
   /**
    * Copyright banner to be printed to the terminal by executables
    * Arguments are the executable's name, year of writing and the name
    * + address of the copyright holder
    */
   void banner(std::string name, Uint year, std::string cpy_holder);
 
 
   /**
    * Planner flags for FFT (follows FFTW, hopefully this choice will
    * be compatible with alternative FFT implementations)
    * @enum muSpectre::FFT_PlanFlags
    */
   enum class FFT_PlanFlags {
     estimate, //!< cheapest plan for slowest execution
     measure,  //!< more expensive plan for fast execution
     patient   //!< very expensive plan for fastest execution
   };
 
   //! continuum mechanics flags
   enum class Formulation{
     finite_strain,   //!< causes evaluation in PK1(F)
     small_strain,    //!< causes evaluation in   σ(ε)
     small_strain_sym //!< symmetric storage as vector ε
   };
 
   //! split cell flags
   enum class SplittedCell{
     yes,
     no
   };
 
+  //! laminate cell flags
+  enum class LaminateCell{
+    yes,
+    no
+  };
+
   /**
    * compile time computation of voigt vector
    */
   template<bool sym=true>
   constexpr Dim_t vsize(Dim_t dim) {
     if (sym) {
       return (dim * (dim - 1) / 2 + dim);
     } else {
       return dim*dim;
     }
   }
 
   //! compute the number of degrees of freedom to store for the strain
   //! tenor given dimension dim
   constexpr Dim_t dof_for_formulation(const Formulation form,
                                       const Dim_t dim) {
     switch (form) {
     case Formulation::small_strain_sym: {
       return vsize(dim);
       break;
     }
     default:
       return ipow(dim, 2);
       break;
     }
   }
 
   //! inserts `muSpectre::Formulation`s into `std::ostream`s
   std::ostream & operator<<(std::ostream & os, Formulation f);
 
   /* ---------------------------------------------------------------------- */
   //! Material laws can declare which type of stress measure they provide,
   //! and µSpectre will handle conversions
   enum class StressMeasure {
     Cauchy,      //!< Cauchy stress σ
     PK1,         //!< First Piola-Kirchhoff stress
     PK2,         //!< Second Piola-Kirchhoff stress
     Kirchhoff,   //!< Kirchhoff stress τ
     Biot,        //!< Biot stress
     Mandel,      //!< Mandel stress
     no_stress_   //!< only for triggering static_asserts
   };
   //! inserts `muSpectre::StressMeasure`s into `std::ostream`s
   std::ostream & operator<<(std::ostream & os, StressMeasure s);
 
   /* ---------------------------------------------------------------------- */
   //! Material laws can declare which type of strain measure they require and
   //! µSpectre will provide it
   enum class StrainMeasure {
     Gradient,      //!< placement gradient (δy/δx)
     Infinitesimal, //!< small strain tensor .5(∇u + ∇uᵀ)
     GreenLagrange, //!< Green-Lagrange strain .5(Fᵀ·F - I)
     Biot,          //!< Biot strain
     Log,           //!< logarithmic strain
     Almansi,       //!< Almansi strain
     RCauchyGreen,  //!< Right Cauchy-Green tensor
     LCauchyGreen,  //!< Left Cauchy-Green tensor
     no_strain_     //!< only for triggering static_assert
   };
   //! inserts `muSpectre::StrainMeasure`s into `std::ostream`s
   std::ostream & operator<<(std::ostream & os, StrainMeasure s);
 
     /* ---------------------------------------------------------------------- */
     /**
      * all isotropic elastic moduli to identify conversions, such as E
      * = µ(3λ + 2µ)/(λ+µ). For the full description, see
      * https://en.wikipedia.org/wiki/Lam%C3%A9_parameters
      * Not all the conversions are implemented, so please add as needed
      */
     enum class ElasticModulus {
       Bulk,          //!< Bulk modulus K
       K = Bulk,      //!< alias for ``ElasticModulus::Bulk``
       Young,         //!< Young's modulus E
       E = Young,     //!< alias for ``ElasticModulus::Young``
       lambda,        //!< Lamé's first parameter λ
       Shear,         //!< Shear modulus G or µ
       G = Shear,     //!< alias for ``ElasticModulus::Shear``
       mu = Shear,    //!< alias for ``ElasticModulus::Shear``
       Poisson,       //!< Poisson's ratio ν
       nu = Poisson,  //!< alias for ``ElasticModulus::Poisson``
       Pwave,         //!< P-wave modulus M
       M=Pwave,       //!< alias for ``ElasticModulus::Pwave``
       no_modulus_};  //!< only for triggering static_asserts
 
     /**
      * define comparison in order to exploit that moduli can be
      * expressed in terms of any two other moduli in any order (e.g. K
      * = K(E, ν) = K(ν, E)
      */
     constexpr inline bool operator<(ElasticModulus A, ElasticModulus B) {
       return static_cast<int>(A) < static_cast<int>(B);
     }
   /* ---------------------------------------------------------------------- */
   /** Compile-time function to g strain measure stored by muSpectre
       depending on the formulation
    **/
   constexpr StrainMeasure get_stored_strain_type(Formulation form) {
     switch (form) {
     case Formulation::finite_strain: {
       return StrainMeasure::Gradient;
       break;
     }
     case Formulation::small_strain: {
       return StrainMeasure::Infinitesimal;
       break;
     }
     default:
       return StrainMeasure::no_strain_;
       break;
     }
   }
 
   /** Compile-time function to g stress measure stored by muSpectre
       depending on the formulation
    **/
   constexpr StressMeasure get_stored_stress_type(Formulation form) {
     switch (form) {
     case Formulation::finite_strain: {
       return StressMeasure::PK1;
       break;
     }
     case Formulation::small_strain: {
       return StressMeasure::Cauchy;
       break;
     }
     default:
       return StressMeasure::no_stress_;
       break;
     }
   }
 
   /* ---------------------------------------------------------------------- */
   /** Compile-time functions to get the stress and strain measures
       after they may have been modified by choosing a formulation.
 
       For instance, a law that expecs a Green-Lagrange strain as input
       will get the infinitesimal strain tensor instead in a small
       strain computation
    **/
   constexpr StrainMeasure get_formulation_strain_type(Formulation form,
                                                       StrainMeasure expected) {
     switch (form) {
     case Formulation::finite_strain: {
       return expected;
       break;
     }
     case Formulation::small_strain: {
       return get_stored_strain_type(form);
       break;
     }
     default:
       return StrainMeasure::no_strain_;
       break;
     }
   }
 
 }  // muSpectre
 
 
 #ifndef EXPLICITLY_TURNED_ON_CXX17
 #include "common/utilities.hh"
 #endif
 
 #endif /* COMMON_H */
diff --git a/src/materials/material_base.cc b/src/materials/material_base.cc
index f96bf03..2ffb1bb 100644
--- a/src/materials/material_base.cc
+++ b/src/materials/material_base.cc
@@ -1,94 +1,105 @@
 /**
  * @file   material_base.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   01 Nov 2017
  *
  * @brief  implementation of materi
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre 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
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "materials/material_base.hh"
 
 namespace muSpectre {
 
   //----------------------------------------------------------------------------//
   template <Dim_t DimS, Dim_t DimM>
   MaterialBase<DimS, DimM>::MaterialBase(std::string name)
     :name(name),assigned_ratio{make_field<ScalarField_t>("assigned ratio",
                                                          this->internal_fields)},
+     assigned_normal_vec{make_field<VectorField_t>("assigned normal vectors",
+                                                   this->internal_fields)},
      assigned_ratio_mapped{assigned_ratio}{
-    static_assert((DimM == oneD)||
-                  (DimM == twoD)||
-                  (DimM == threeD), "only 1, 2, or threeD supported");
-  }
+       static_assert((DimM == oneD)||
+                     (DimM == twoD)||
+                     (DimM == threeD), "only 1, 2, or threeD supported");
+     }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   const std::string & MaterialBase<DimS, DimM>::get_name() const {
     return this->name;
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void MaterialBase<DimS, DimM>::add_pixel(const Ccoord &ccoord) {
     this->internal_fields.add_pixel(ccoord);
   }
 
-    /* ---------------------------------------------------------------------- */
+  /* ---------------------------------------------------------------------- */
   template < Dim_t DimS, Dim_t DimM>
   void MaterialBase< DimS, DimM>::
   add_pixel_split(const Ccoord &local_ccoord, Real ratio){
     auto & this_mat = static_cast<MaterialBase&>(*this);
     this_mat.add_pixel(local_ccoord);
     this->assigned_ratio.push_back(ratio);
   }
   /* ---------------------------------------------------------------------- */
+  template < Dim_t DimS, Dim_t DimM>
+  void MaterialBase< DimS, DimM>::
+  add_pixel_laminate(const Ccoord &local_ccoord, Real ratio, Rcoord normal_vec){
+    auto & this_mat = static_cast<MaterialBase&>(*this);
+    this_mat.add_pixel(local_ccoord);
+    this->assigned_ratio.push_back(ratio);
+    this->assigned_normal_vec.push_back(normal_vec);
+  }
+  /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void MaterialBase<DimS, DimM>::compute_stresses(const Field_t & F,
                                                   Field_t & P,
                                                   Formulation form, SplittedCell is_cell_splitted) {
     this->compute_stresses(StrainField_t::check_ref(F),
                            StressField_t::check_ref(P),
                            form, is_cell_splitted);
   }
 
   template <Dim_t DimS, Dim_t DimM>
   Real MaterialBase<DimS, DimM>::get_assigned_ratio(Ccoord pixel){
     return this->assigned_ratio_mapped[pixel];
   }
 
   /* ---------------------------------------------------------------------- */
   template <Dim_t DimS, Dim_t DimM>
   void MaterialBase<DimS, DimM>::compute_stresses_tangent(const Field_t & F,
                                                           Field_t & P,
                                                           Field_t & K,
                                                           Formulation form, SplittedCell is_cell_splitted) {
     this->compute_stresses_tangent(StrainField_t::check_ref(F),
                                    StressField_t::check_ref(P),
                                    TangentField_t::check_ref(K),
                                    form, is_cell_splitted);
   }
 
   template class MaterialBase<2, 2>;
   template class MaterialBase<2, 3>;
   template class MaterialBase<3, 3>;
 
 }  // muSpectre
diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh
index 48c2eed..f2c1871 100644
--- a/src/materials/material_base.hh
+++ b/src/materials/material_base.hh
@@ -1,170 +1,177 @@
 /**
  * @file   material_base.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   25 Oct 2017
  *
  * @brief  Base class for materials (constitutive models)
  *
  * Copyright © 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre 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
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #ifndef MATERIAL_BASE_H
 #define MATERIAL_BASE_H
 
 #include "common/common.hh"
 #include "common/field.hh"
 #include "common/field_collection.hh"
 
 #include <string>
 
 namespace muSpectre {
 
   //! DimS spatial dimension (dimension of problem
   //! DimM material_dimension (dimension of constitutive law) and
   /**
    * @a DimM is the material dimension (i.e., the dimension of constitutive
    * law; even for e.g. two-dimensional problems the constitutive law could
    * live in three-dimensional space for e.g. plane strain or stress problems)
    */
   template<Dim_t DimS, Dim_t DimM>
   class MaterialBase
   {
   public:
     //! typedefs for data handled by this interface
     //! global field collection for cell-wide fields, like stress, strain, etc
     using GFieldCollection_t = GlobalFieldCollection<DimS>;
     //! field collection for internal variables, such as eigen-strains,
     //! plastic strains, damage variables, etc, but also for managing which
     //! pixels the material is responsible for
     using MFieldCollection_t = LocalFieldCollection<DimS>;
 
     using iterator = typename MFieldCollection_t::iterator; //!< pixel iterator
     //! polymorphic base class for fields only to be used for debugging
     using Field_t = internal::FieldBase<GFieldCollection_t>;
     //! Full type for stress fields
     using StressField_t = TensorField<GFieldCollection_t, Real, secondOrder, DimM>;
     //! Full type for strain fields
     using StrainField_t = StressField_t;
     //! Full type for tangent stiffness fields fields
     using TangentField_t = TensorField<GFieldCollection_t, Real, fourthOrder, DimM>;
     using Ccoord = Ccoord_t<DimS>; //!< cell coordinates type
+    using Rcoord = Rcoord_t<DimS>; //!< real coordinates type
     //! Default constructor
     MaterialBase() = delete;
 
     //! Construct by name
     MaterialBase(std::string name);
 
     //! Copy constructor
     MaterialBase(const MaterialBase &other) = delete;
 
     //! Move constructor
     MaterialBase(MaterialBase &&other) = delete;
 
     //! Destructor
     virtual ~MaterialBase() = default;
 
     //! Copy assignment operator
     MaterialBase& operator=(const MaterialBase &other) = delete;
 
     //! Move assignment operator
     MaterialBase& operator=(MaterialBase &&other) = delete;
 
 
     /**
      *  take responsibility for a pixel identified by its cell coordinates
      *  WARNING: this won't work for materials with additional info per pixel
      *  (as, e.g. for eigenstrain), we need to pass more parameters. Materials
      *  of this type need to overload add_pixel
      */
     virtual void add_pixel(const Ccoord & ccooord);
 
     virtual void add_pixel_split (const Ccoord &local_ccoord,
                                   Real ratio);
+
+    virtual void add_pixel_laminate (const Ccoord &local_ccoord,
+                                     Real ratio,
+                                     Rcoord normal_vec);
     //! allocate memory, etc, but also: wipe history variables!
     virtual void initialise() = 0;
 
     /**
      * for materials with state variables, these typically need to be
      * saved/updated an the end of each load increment, the virtual
      * base implementation does nothing, but materials with history
      * variables need to implement this
      */
     virtual void save_history_variables() {};
 
     //! return the material's name
     const std::string & get_name() const;
 
     //! spatial dimension for static inheritance
     constexpr static Dim_t sdim() {return DimS;}
     //! material dimension for static inheritance
     constexpr static Dim_t mdim() {return DimM;}
     //! computes stress
     virtual void compute_stresses(const StrainField_t & F,
                                   StressField_t & P,
                                   Formulation form, SplittedCell is_cell_splitted) = 0;
     /**
      * Convenience function to compute stresses, mostly for debugging and
      * testing. Has runtime-cost associated with compatibility-checking and
      * conversion of the Field_t arguments that can be avoided by using the
      * version with strongly typed field references
      */
     void compute_stresses(const Field_t & F,
                           Field_t & P,
                           Formulation form, SplittedCell is_cell_splitted = SplittedCell::no);
     //! computes stress and tangent moduli
     virtual void compute_stresses_tangent(const StrainField_t & F,
                                           StressField_t & P,
                                           TangentField_t & K,
                                           Formulation form, SplittedCell is_cell_splitted) = 0;
     /**
      * Convenience function to compute stresses and tangent moduli, mostly for
      * debugging and testing. Has runtime-cost associated with
      * compatibility-checking and conversion of the Field_t arguments that can
      * be avoided by using the version with strongly typed field references
      */
     void compute_stresses_tangent(const Field_t & F,
                                   Field_t & P,
                                   Field_t & K,
                                   Formulation form, SplittedCell is_cell_splitted = SplittedCell::no);
     // this function return the ratio of which the input pixel is consisted of this material
     Real get_assigned_ratio(Ccoord pixel);
     //! iterator to first pixel handled by this material
     inline iterator begin() {return this->internal_fields.begin();}
     //! iterator past the last pixel handled by this material
     inline iterator end()  {return this->internal_fields.end();}
     //! number of pixels assigned to this material
     inline size_t size() const {return this->internal_fields.size();}
 
     //! gives access to internal fields
     inline MFieldCollection_t& get_collection() {return this->internal_fields;}
   protected:
-
     const std::string name; //!< material's name (for output and debugging)
     MFieldCollection_t internal_fields{};//!< storage for internal variables
     using ScalarField_t = ScalarField<LocalFieldCollection<DimS>, Real> ;
     ScalarField_t & assigned_ratio; //!< field holding the assigning ratio of the material
+    using VectorField_t = TensorField<LocalFieldCollection<DimS>, Real, firstOrder, DimM>;
+    VectorField_t & assigned_normal_vec; //!< field holding the assigning ratio of the material
     using ScalarFieldMap_t = ScalarFieldMap<LocalFieldCollection<DimS>, Real, true>;
     ScalarFieldMap_t assigned_ratio_mapped;
 
   private:
   };
 }  // muSpectre
 
 #endif /* MATERIAL_BASE_H */
+