diff --git a/src/materials/laminate_homogenisation.cc b/src/materials/laminate_homogenisation.cc
index 3192598..8a08039 100644
--- a/src/materials/laminate_homogenisation.cc
+++ b/src/materials/laminate_homogenisation.cc
@@ -1,621 +1,697 @@
 /**
  * @file   lamminate_hemogenization.hh
  *
  * @author Ali Falsafi <ali.falsafi@epfl.ch>
  *
  * @date   28 Sep 2018
  *
  * @brief
  *
  * 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/geometry.hh"
 #include "common/common.hh"
 #include "laminate_homogenisation.hh"
 #include "material_anisotropic.hh"
 #include "materials_toolbox.hh"
 
 #include <tuple>
 
 namespace muSpectre {
 
   /* ---------------------------------------------------------------------- */
   template<>
+  constexpr auto LamHomogen<threeD>::get_equation_indexes()
+    ->std::array<Dim_t, threeD>{
+    constexpr Dim_t Dim {3};
+    return std::array<Dim_t, Dim>{0, 1, 2};
+  }
+  template<>
+  constexpr auto LamHomogen<twoD>::get_equation_indexes()
+    ->std::array<Dim_t, twoD>{
+    constexpr Dim_t Dim {2};
+    return  std::array<Dim_t, Dim>{0, 1};
+  }
+
+  /* ---------------------------------------------------------------------- */
+  template<>
+  constexpr auto LamHomogen<threeD>::get_serial_non_equation_indexes()
+    ->std::array<Dim_t, threeD-1>{
+    constexpr Dim_t Dim {3};
+    return std::array<Dim_t, Dim-1>{3, 4};
+  }
+
+  template<>
+  constexpr auto LamHomogen<twoD>::get_serial_non_equation_indexes()
+    ->std::array<Dim_t, twoD-1>{
+    constexpr Dim_t Dim {2};
+    return  std::array<Dim_t, Dim-1>{2};
+  }
+
+  /*---------------------------------------------------------------------- */
+  template<>
   constexpr auto LamHomogen<threeD>::get_serial_indexes()
     ->std::array<Dim_t, 2*threeD-1>{
     constexpr Dim_t Dim {3};
     return std::array<Dim_t, 2*Dim-1>{0, 1, 2, 3, 6};
   }
   template<>
   constexpr auto LamHomogen<twoD>::get_serial_indexes()
     ->std::array<Dim_t, 2*twoD-1>{
     constexpr Dim_t Dim {2};
     return  std::array<Dim_t, 2*Dim-1>{0, 1, 2};
   }
 
   /* ---------------------------------------------------------------------- */
   template<>
   constexpr auto LamHomogen<threeD>::get_parallel_indexes()
     ->std::array<Dim_t, (threeD-1)*(threeD-1)>{
     constexpr Dim_t Dim {3};
     return std::array<Dim_t, (Dim-1)*(Dim-1)>{4, 5, 7, 8};;
   }
   template<>
   constexpr auto LamHomogen<twoD>::get_parallel_indexes()
     ->std::array<Dim_t, (twoD-1)*(twoD-1)>{
     constexpr Dim_t Dim {2};
     return  std::array<Dim_t, (Dim-1)*(Dim-1)>{3};
   }
 
   /* ---------------------------------------------------------------------- */
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   get_serial_strain(const Eigen::Ref<Strain_t>& Fs)->Serial_strain_t{
     Serial_strain_t ret_F;
     auto serial_indexes{get_serial_indexes()};
 
     for (int i{0}; i < ret_F.size(); ++i){
       ret_F(i) = Fs(serial_indexes[i]);
     }
     return ret_F;
   }
 
   /* ---------------------------------------------------------------------- */
 
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   get_parallel_strain(const Eigen::Ref<Strain_t>& Fs)->Parallel_strain_t{
     Parallel_strain_t ret_F;
     auto parallel_indexes(get_parallel_indexes());
 
     for (int i{0}; i < ret_F.size(); ++i){
       ret_F(i) = Fs(parallel_indexes[i]);
     }
     return ret_F;
   }
   /* ---------------------------------------------------------------------- */
-  
+
+  template<Dim_t Dim>
+  auto LamHomogen<Dim>::
+  get_equation_strain_from_serial(const Eigen::Ref<Serial_strain_t>& Fs)
+    -> Equation_strain_t{
+    Equation_strain_t ret_F;
+    auto equation_indexes(get_equation_indexes());
+
+    for (int i{0}; i < ret_F.size(); ++i){
+      ret_F(i) = Fs(equation_indexes[i]);
+    }
+    return ret_F;
+  }
+
+  template<Dim_t Dim>
+  auto LamHomogen<Dim>::
+  get_serial_strain_from_equation(const Eigen::Ref<Equation_strain_t>& Fs)
+    -> Serial_strain_t{
+    Serial_strain_t ret_F;
+    auto equation_indexes(get_equation_indexes());
+    auto non_equation_indexes(get_serial_non_equation_indexes());
+    for (unsigned int i{0}; i < equation_indexes.size(); ++i){
+      ret_F(i) = Fs(equation_indexes[i]);
+    }
+    for (unsigned int i{0}; i < non_equation_indexes.size(); ++i){
+      ret_F(i) = Fs(non_equation_indexes[i] - Dim);
+    }
+    return ret_F;
+  }
   /* ---------------------------------------------------------------------- */
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   get_total_strain(const Eigen::Ref<Serial_strain_t>& F_seri,
                    const Eigen::Ref<Parallel_strain_t>& F_para)
     ->Strain_t{
     Strain_t ret_F;
     auto serial_indexes{get_serial_indexes()};
     auto parallel_indexes{get_parallel_indexes()};
 
     for (unsigned int i{0}; i < serial_indexes.size(); ++i){
       ret_F(serial_indexes[i]) = F_seri(i);
     }
     for (unsigned int i{0}; i < parallel_indexes.size(); ++i){
       ret_F(parallel_indexes[i]) = F_para(i);
     }
     return ret_F;
   }
 
   /* ---------------------------------------------------------------------- */
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   get_serial_stiffness(const Eigen::Ref<Stiffness_t>& K)
     ->Serial_stiffness_t{
     Serial_stiffness_t ret_K{};
     auto serial_indexes{get_serial_indexes()};
     for (int i{0}; i < (2*Dim-1); ++i){
       for (int j{0}; j < (2*Dim-1); ++j){
         ret_K(i,j) = K(serial_indexes[i], serial_indexes[j]);
       }
     }
     return ret_K;
   }
 
   /* ---------------------------------------------------------------------- */
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   get_parallel_stiffness(const Eigen::Ref<Stiffness_t>& K)
     ->Parallel_stiffness_t{
     Parallel_stiffness_t ret_K{};
     auto parallel_indexes{get_parallel_indexes()};
 
     for (int i{0}; i < (Dim-1)*(Dim-1); ++i){
       for (int j{0}; j < (Dim-1)*(Dim-1); ++j){
         ret_K(i,j) = K(parallel_indexes[i], parallel_indexes[j]);
       }
     }
     return ret_K;
   }
 
   /* ---------------------------------------------------------------------- */
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   get_total_stiffness(const Eigen::Ref<Serial_stiffness_t>& K_seri,
                       const Eigen::Ref<Parallel_stiffness_t>& K_para)
     ->Stiffness_t{
     Stiffness_t ret_K{};
     auto serial_indexes{get_serial_indexes()};
     auto parallel_indexes{get_parallel_indexes()};
 
     for (int i{0}; i < (Dim-1)*(Dim-1); ++i){
       for (int j{0}; j < (Dim-1)*(Dim-1); ++j){
         ret_K(parallel_indexes[i],parallel_indexes[j]) = K_para(i, j);
       }
     }
 
     for (int i{0}; i < (Dim-1)*(Dim-1); ++i){
       for (int j{0}; j < (Dim-1)*(Dim-1); ++j){
         ret_K(serial_indexes[i],serial_indexes[j]) = K_seri(i, j);
       }
     }
     return ret_K;
   }
-
+  /* ---------------------------------------------------------------------- */
+  // TODO: complete this dunction and its inverse
+  template<Dim_t Dim>
+  auto LamHomogen<Dim>::
+  get_equation_stiffness_from_serial(const Eigen::Ref<Serial_stiffness_t>& Ks)
+    -> Equation_stiffness_t{
+    Equation_stiffness_t ret_K{};
+    for (unsigned int i = 0 ; i < Dim; ++i) {
+      for (unsigned int j = 0 ; j < Dim; ++j) {
+        ret_K(i,j) = Ks(i,j);
+      }
+      for (unsigned int j = 1 ; j < Dim; ++j) {
+        ret_K(i,j) += Ks(i,j+(Dim-1));
+      }
+    }
+    return ret_K;
+  }
   /* ---------------------------------------------------------------------- */
   // It does not to be especialised because all functions used init are
   // especialised. this function computes the strain in the second
   // layer of laminate material from strain in the first layer and total strain
 
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   linear_eqs_seri (Real ratio,
                    const Eigen::Ref<Strain_t> & F_lam,
                    const Eigen::Ref<Serial_strain_t> & F_1)
     ->Serial_strain_t{
     Serial_strain_t F_2_seri;
     Serial_strain_t F_lam_seri = get_serial_strain(F_lam);
     F_2_seri = (F_lam_seri - ratio * F_1) / (1-ratio);
     return F_2_seri;
   }
 
   /* ---------------------------------------------------------------------- */
   // this function recieves the strain and returns ths delta_stress and the
   // jaaconbian in general coordinates
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   delta_stress_eval (const Function_t &mat_1_stress_eval,
                      const Function_t &mat_2_stress_eval,
                      const Eigen::Ref<StrainMat_t>& F_1,
                      const Eigen::Ref<StrainMat_t>& F_2,
                      RotatorNormal<Dim> rotator,
                      Real ratio)
     ->std::tuple<Serial_stress_t, Serial_stiffness_t>{
     auto P_K_1 = mat_1_stress_eval(F_1);
     auto P_K_2 = mat_2_stress_eval(F_2);
 
     StressMat_t del_P_mat; del_P_mat = std::get<0>(P_K_1) - std::get<0>(P_K_2);
     del_P_mat = rotator.rotate(del_P_mat);
     auto del_P_seri = get_serial_strain(Eigen::Map<Stress_t>(del_P_mat.data()));
 
     Stiffness_t del_K; del_K =  (std::get<1>(P_K_1) -
                                  (-ratio/(1-ratio))* std::get<1>(P_K_2));
     del_K = rotator.rotate(del_K);
     auto del_K_seri = get_serial_stiffness(del_K);
 
     return std::make_tuple(del_P_seri, del_K_seri);
   }
 
 
   /* ---------------------------------------------------------------------- */
   // linear equation are solved in rot coordiantes to obtain F_2
   // form known F_1 and then these two are passed to calculate
   // stress in the general corrdinates
 
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   delta_stress_eval_F_1 (const Function_t &mat_1_stress_eval,
                          const Function_t &mat_2_stress_eval,
                          Eigen::Ref<Strain_t> F_0,
                          Eigen::Ref<Strain_t> F_1_rot,
                          RotatorNormal<Dim> rotator,
                          Real ratio)
     ->std::tuple<Serial_stress_t, Serial_stiffness_t>{
     StrainMat_t F_0_rot =
       rotator.rotate(Eigen::Map<StrainMat_t> (F_0.data()));
 
     Parallel_strain_t F_0_par_rot =
       get_parallel_strain(Eigen::Map<Strain_t>(F_0_rot.data()));
 
     StrainMat_t F_1 =
       rotator.rotate_back(Eigen::Map<StrainMat_t>(F_1_rot.data()));
 
     Serial_strain_t F_1_seri_rot =
       get_serial_strain(Eigen::Map<Strain_t> (F_1_rot.data()));
 
     Serial_strain_t F_2_seri_rot =
       linear_eqs_seri(ratio,
                       Eigen::Map<Strain_t> (F_0_rot.data()),
                       F_1_seri_rot);
     Strain_t F_2_rot =
       get_total_strain(F_2_seri_rot, F_0_par_rot);
 
     StrainMat_t F_2  =
       rotator.rotate_back(Eigen::Map<StrainMat_t> (F_2_rot.data()));
 
     return delta_stress_eval(mat_1_stress_eval,
                              mat_2_stress_eval,
                              Eigen::Map<StrainMat_t> (F_1.data()),
                              Eigen::Map<StrainMat_t> (F_2.data()),
                              rotator,
                              ratio);
   }
 
   /* ---------------------------------------------------------------------- */
   template<Dim_t Dim>
   inline auto LamHomogen<Dim>::
   del_energy_eval (const Real del_F_norm,
                    const Real delta_P_norm)
     ->Real{
     return del_F_norm * delta_P_norm;
   }
   /* ---------------------------------------------------------------------- */
   template<>
   auto LamHomogen<twoD>::
   lam_K_combine(const Eigen::Ref<Stiffness_t>& K_1,
                 const Eigen::Ref<Stiffness_t>& K_2,
                 const Real ratio)
     ->Stiffness_t{
     using Mat_A_t = Eigen::Matrix<Real, twoD, twoD>;
     using Vec_A_t = Eigen::Matrix<Real, twoD, 1>;
     using Vec_AT_t = Eigen::Matrix<Real,   1, twoD>;
 
     std::array<double, 3> cf = {1.0, sqrt(2.0), 2.0};
     auto get_A11 = [cf](const Eigen::Ref<Stiffness_t>& K){
       Mat_A_t A11 = Mat_A_t::Zero() ;
       A11<<
       cf[0]*get(K,0,0,0,0), cf[1]*get(K,0,0,0,1),
       cf[1]*get(K,0,0,0,1), cf[2]*get(K,0,1,0,1);
       return A11;
     };
     Mat_A_t A11c;
     A11c<<
       cf[0], cf[1],
       cf[1], cf[2];
     auto get_A12 = [cf](const Eigen::Ref<Stiffness_t>& K){
       Vec_A_t A12 = Vec_A_t::Zero() ;
       A12<<
       cf[0]*get(K,0,0,1,1),
       cf[1]*get(K,1,1,0,1);
       return A12;
     };
     Vec_A_t A12c = Vec_A_t::Zero();
     A12c <<
       cf[0],
       cf[1];
 
     Mat_A_t A11_1 {get_A11(K_1)};
     Mat_A_t A11_2 {get_A11(K_2)};
 
     Vec_A_t A12_1 = {get_A12(K_1)}; Vec_AT_t A21_1 = A12_1.transpose();
     Vec_A_t A12_2 = {get_A12(K_2)}; Vec_AT_t A21_2 = A12_2.transpose();
 
     Real A22_1 = get(K_1,1,1,1,1);
     Real A22_2 = get(K_2,1,1,1,1);
 
 
     auto get_inverse_average =
       [&ratio](const Eigen::Ref<Mat_A_t> &matrix_1,
                const Eigen::Ref<Mat_A_t> &matrix_2){
       return ((      ratio * matrix_1.inverse() +
                      (1-ratio)* matrix_2.inverse()).inverse());
     };
     auto get_average =
       [&ratio](Real A_1,
                Real A_2){
       return ratio*A_1 + (1-ratio)*A_2;
     };
 
     auto get_average_vec =
       [&ratio](Vec_A_t A_1,
                Vec_A_t A_2){
       return ratio*A_1 + (1-ratio)*A_2;};
 
     auto get_average_vecT =
       [&ratio](Vec_AT_t A_1,
                Vec_AT_t A_2){
       return ratio*A_1 + (1-ratio)*A_2;};
 
     // calculating average of A matrices  of the materials
     Mat_A_t A11 = get_inverse_average(A11_1, A11_2);
     Vec_A_t A12 = A11 * get_average_vec(A11_1.inverse()*A12_1,
                                         A11_2.inverse()*A12_2);
     Real A22 = get_average(A22_1 - A21_1 * A11_1.inverse() * A12_1,
                            A22_2 - A21_2 * A11_2.inverse() * A12_2) +
       get_average_vecT(A21_1 * A11_1.inverse(),
                        A21_2 * A11_2.inverse()) *
       A11 * get_average_vec(A11_1.inverse() * A12_1,
                             A11_2.inverse() * A12_2);
 
     std::vector<Real> c_maker_inp =
       {A11(0,0)/A11c(0,0),A12(0,0)/A12c(0,0),A11(0,1)/A11c(0,1),
        A22               ,A12(1,0)/A12c(1,0),
        A11(1,1)/A11c(1,1)};
     // now the resultant stiffness is calculated fro 6 elements obtained from
     // A matrices averaging routine:
     Stiffness_t ret_K = Stiffness_t::Zero();
     ret_K = MaterialAnisotropic<twoD>::c_maker(c_maker_inp);
     return ret_K ;
 
   }
 
   /*------------------------------------------------------------------*/
   template<>
   auto LamHomogen<threeD>::
   lam_K_combine(const Eigen::Ref<Stiffness_t>& K_1,
                 const Eigen::Ref<Stiffness_t>& K_2,
                 const Real ratio)
     ->Stiffness_t{
 
     // the combination method is obtained form P. 163 of
     // "Theory of Composites"
     //  Author : Milton_G_W
 
     //constructing "A" matrices( A11, A12, A21,  A22) according to the
     //procedure from the book:
 
     // this type of matrix will be used in calculating the combinatio of the
     // Stiffness matrixes
     using Mat_A_t = Eigen::Matrix<Real, threeD, threeD>;
 
     // these coeffs are used in constructing A matrices from k and vice versa
     std::array<double, 3> cf = {1.0, sqrt(2.0), 2.0};
     auto get_A11 = [cf](const Eigen::Ref<Stiffness_t>& K){
       Mat_A_t A11 = Mat_A_t::Zero() ;
       A11<<
       cf[0]*get(K,0,0,0,0), cf[1]*get(K,0,0,0,2), cf[1]*get(K,0,0,0,1),
       cf[1]*get(K,0,0,0,2), cf[2]*get(K,0,2,0,2), cf[2]*get(K,0,2,0,1),
       cf[1]*get(K,0,0,0,1), cf[2]*get(K,0,2,0,1), cf[2]*get(K,0,1,0,1);
       return A11;
     };
     Mat_A_t A11c;
     A11c<<
       cf[0], cf[1], cf[1],
       cf[1], cf[2], cf[2],
       cf[1], cf[2], cf[2];
 
     auto get_A12 = [cf](const Eigen::Ref<Stiffness_t>& K){
       Mat_A_t A12 = Mat_A_t::Zero();
       A12 <<
       cf[0]*get(K,0,0,1,1), cf[0]*get(K,0,0,2,2), cf[1]*get(K,0,0,1,2),
       cf[1]*get(K,1,1,0,2), cf[1]*get(K,2,2,0,2), cf[2]*get(K,1,2,0,2),
       cf[1]*get(K,1,1,0,1), cf[1]*get(K,2,2,0,1), cf[2]*get(K,1,2,0,1);
       return A12;
     };
     Mat_A_t A12c;
     A12c<<
       cf[0], cf[0], cf[1],
       cf[1], cf[1], cf[2],
       cf[1], cf[1], cf[2];
 
 
     auto get_A22 = [cf](const Eigen::Ref<Stiffness_t>& K){
       Mat_A_t A22 = Mat_A_t::Zero();
       A22 <<
       cf[0]*get(K,2,2,2,2), cf[0]*get(K,2,2,1,1), cf[1]*get(K,2,2,1,2),
       cf[0]*get(K,2,2,1,1), cf[0]*get(K,1,1,1,1), cf[1]*get(K,1,1,1,2),
       cf[1]*get(K,2,2,1,2), cf[1]*get(K,1,1,1,2), cf[2]*get(K,2,1,2,1);
       return A22;
     };
     Mat_A_t A22c;
     A22c <<
       cf[0], cf[0], cf[1],
       cf[0], cf[0], cf[1],
       cf[1], cf[1], cf[2];
 
 
     Mat_A_t A11_1 {get_A11(K_1)};
     Mat_A_t A11_2 {get_A11(K_2)};
 
     Mat_A_t A12_1 {get_A12(K_1)}; Mat_A_t A21_1 {A12_1.transpose()};
     Mat_A_t A12_2 {get_A12(K_2)}; Mat_A_t A21_2 {A12_2.transpose()};
 
     Mat_A_t A22_1 {get_A22(K_1)};
     Mat_A_t A22_2 {get_A22(K_2)};
 
     // these two functions are routines to compute average of "A" matrices
     auto get_inverse_average =
       [&ratio](const Eigen::Ref<Mat_A_t> &matrix_1,
                const Eigen::Ref<Mat_A_t> &matrix_2){
       return ((      ratio * matrix_1.inverse() +
                      (1-ratio)* matrix_2.inverse()).inverse());
     };
 
     auto get_average =
       [&ratio](const Mat_A_t &matrix_1,
                const Mat_A_t &matrix_2){
       return (ratio * matrix_1 + (1-ratio) * matrix_2);
     };
 
     // calculating average of A matrices  of the materials
     Mat_A_t A11 = get_inverse_average(A11_1, A11_2);
     Mat_A_t A12 = A11 * get_average(A11_1.inverse()*A12_1,
                                     A11_2.inverse()*A12_2);
     Mat_A_t A22 = get_average(A22_1 - A21_1 * A11_1.inverse() * A12_1,
                               A22_2 - A21_2 * A11_2.inverse() * A12_2) +
       get_average(A21_1 * A11_1.inverse(),
                   A21_2 * A11_2.inverse()) *
       A11 * get_average(A11_1.inverse() * A12_1,
                         A11_2.inverse() * A12_2);
 
     std::vector<Real> c_maker_inp =
       {A11(0,0)/A11c(0,0), A12(0,0)/A12c(0,0), A12(0,1)/A12c(0,1),
        A12(0,2)/A12c(0,2), A11(1,0)/A11c(1,0), A11(0,2)/A11c(0,2),
 
        A22(1,1)/A22c(1,1), A22(1,0)/A22c(1,0), A22(2,1)/A22c(2,1),
        A12(1,0)/A12c(1,0), A12(2,0)/A12c(2,0),
 
        A22(0,0)/A22c(0,0), A22(2,0)/A22c(2,0),
        A12(1,1)/A12c(1,1), A12(2,1)/A12c(2,1),
 
        A22(2,2)/A22c(2,2), A12(1,2)/A12c(1,2), A12(2,2)/A12c(2,2),
        A11(1,1)/A11c(1,1), A11(2,1)/A11c(2,1),
        A11(2,2)/A11c(2,2)};
     // now the resultant stiffness is calculated fro 21 elements obtained from
     // A matrices averaging routine :
     Stiffness_t ret_K = Stiffness_t::Zero();
     ret_K = MaterialAnisotropic<threeD>::c_maker(c_maker_inp);
     return ret_K ;
   }
 
   /* ---------------------------------------------------------------------- */
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   lam_P_combine(Eigen::Ref<StressMat_t> P_1,
                 Eigen::Ref<StressMat_t> P_2,
                 const Real ratio)
     ->StressMat_t{
     StressMat_t P_tmp{};
     //return P_tmp ;
     Serial_stress_t P_1_seri =
       get_serial_strain(Eigen::Map<Stress_t>(P_1.data()));
 
     Parallel_stress_t P_1_para =
       get_parallel_strain(Eigen::Map<Stress_t>(P_1.data()));
 
     Parallel_stress_t P_2_para =
       get_parallel_strain(Eigen::Map<Stress_t>(P_2.data()));
 
     Parallel_stress_t P_para = ratio * P_1_para + (1-ratio)* P_2_para ;
     auto ret_P = get_total_strain(P_1_seri, P_para);
     return Eigen::Map<StressMat_t> (ret_P.data());
   }
   /* ---------------------------------------------------------------------- */
   template<Dim_t Dim>
   auto LamHomogen<Dim>::
   laminate_solver(Eigen::Ref<Strain_t> F_coord,
                   Function_t mat_1_stress_eval,
                   Function_t mat_2_stress_eval,
                   Real ratio,
                   const Eigen::Ref<Vec_t>& normal_vec,
                   Real tol,
                   Dim_t max_iter)
     ->std::tuple<StressMat_t, Stiffness_t>{
     /*
      * here we rotate the strain such that the laminate intersection normal
      * would align with the x-axis F_lam is the total strain in the new
      * coordinates.
      */
     Real del_energy;
     Dim_t iter{0};
 
     RotatorNormal<Dim> rotator(normal_vec);
     ///////////temporary :
     StrainMat_t aaaa = StrainMat_t::Identity();
     rotator.set_rot_mat(aaaa);
     ///////////
     Strain_t F_1{F_coord},F_2{F_coord};
     StrainMat_t F_1_rot_mat {rotator.rotate(Eigen::Map<StrainMat_t>(F_1.data()))};
     Strain_t F_1_rot = Eigen::Map<Strain_t>(F_1_rot_mat.data());
     Parallel_strain_t F_para{get_parallel_strain(F_coord)};
     Serial_strain_t F_1_new_seri{get_serial_strain(F_1_rot)};
     Serial_strain_t F_1_seri {F_1_new_seri};
+
+    Equation_strain_t F_1_eq     =   get_equation_strain_from_serial(F_1_seri);
+    Equation_strain_t F_1_new_eq =   get_equation_strain_from_serial(F_1_seri);
     std::cout << std::endl <<  "F_1 Initial : " << F_1.transpose() << std::endl;
 
     std::tuple<Serial_stress_t, Serial_stiffness_t> delta_P_K_seri;
     Serial_stress_t delta_P_seri{};
     Serial_stiffness_t delta_K_seri{};
+    Equation_stress_t delta_P_eq{};
+    Equation_stiffness_t delta_K_eq{};
 
     std::tuple<StressMat_t, Stiffness_t> P_K_1, P_K_2;
     StressMat_t P_1{}, P_2{}, ret_P{};
     Stiffness_t K_1{}, K_2{}, ret_K{};
 
 
     /*
      * (P_1 - P_2) (as functionsback to life  of F_1)
      * and its Jacobian are calculated here in the rotated coordiantes
      */
     delta_P_K_seri = delta_stress_eval_F_1(mat_1_stress_eval,
                                            mat_2_stress_eval,
                                            F_coord,
                                            F_1,
                                            rotator,
                                            ratio);
     delta_P_seri = std::get<0> (delta_P_K_seri);
     delta_K_seri = std::get<1> (delta_P_K_seri);
-
+    delta_P_eq = get_equation_strain_from_serial(delta_P_seri);
+    delta_K_eq = get_equation_stiffness_from_serial(delta_K_seri);
 
     // solution loop:
     do {
       // updating variables:
-      F_1_seri =  F_1_new_seri;
-
+      F_1_eq =  F_1_new_eq;
+      F_1_seri = F_1_new_seri;
       //solving for ∇(delP) * [x_n+1 - x_n] = -delP
-      F_1_new_seri = F_1_seri - delta_K_seri.ldlt().solve(delta_P_seri);
-
+      F_1_new_eq = F_1_eq - delta_K_eq.ldlt().solve(delta_P_eq);
+      F_1_new_seri = get_serial_strain_from_equation(F_1_new_eq);
       // updating F_1 before stress and stiffness evaluation
       F_1_rot = get_total_strain(F_1_new_seri, F_para);
 
       // stress and stiffenss evaluation in the rotated coordinates
       delta_P_K_seri = delta_stress_eval_F_1(mat_1_stress_eval,
                                              mat_2_stress_eval,
                                              F_coord,
                                              Eigen::Map<Strain_t>(F_1_rot.data()),
                                              rotator,
                                              ratio);
 
       delta_P_seri = std::get<0> (delta_P_K_seri);
       delta_K_seri = std::get<1> (delta_P_K_seri);
+      delta_P_eq = get_equation_strain_from_serial(delta_P_seri);
+      delta_K_eq = get_equation_stiffness_from_serial(delta_K_seri);
 
-      std::cout << std::endl << "iter : " <<iter << std::endl
-                << "res : "<< delta_P_seri.transpose() << std::endl
-                << "delta K : "  << std::endl << delta_K_seri << std::endl
-                << "F_1 : "<< F_1_new_seri.transpose() << std::endl;
 
       // computing variable used for loop termination criterion
       // energy norm of the residual
-      del_energy = std::abs(delta_P_seri.dot(F_1_new_seri - F_1_seri));
+      del_energy = std::abs(delta_P_eq.dot(F_1_new_eq - F_1_eq));
       iter ++;
     } while(del_energy > tol && iter < max_iter);
 
     // computing stress and stiffness of each layer according to it's
     // correspondent strin calculated in the loop
     P_K_1 = mat_1_stress_eval(Eigen::Map<StrainMat_t> (F_1.data()));
     P_K_2 = mat_2_stress_eval(Eigen::Map<StrainMat_t> (F_2.data()));
     P_1 = std::get<0>(P_K_1); P_2 = std::get<0>(P_K_2);
     K_1 = std::get<1>(P_K_1); K_2 = std::get<1>(P_K_2);
-    std::cout << "K_1" << std::endl << K_1;
 
     // check if the loop has lead in convergence or not:
     if (iter == max_iter) {
       throw std::runtime_error("Error: the solver has not converged!!!!");
     }
     // storing the resultant strain in each layer in F_1 and F_2
     auto F_2_new_seri = linear_eqs_seri(ratio, F_coord, F_1_new_seri);
     F_1 = get_total_strain(F_1_new_seri, F_para);
     F_2 = get_total_strain(F_2_new_seri, F_para);
     std::cout << "steps : " << iter << std::endl;
     std::cout <<  "F_1 Final : " << F_1.transpose() << std::endl;
     std::cout <<  "F_2 Final : " << F_2.transpose() << std::endl;
 
 
     // combine the computed strains and tangents to have the resultant
     // stress and tangent of the pixel
     ret_P = lam_P_combine(P_1, P_2, ratio);
     ret_K = lam_K_combine(K_1, K_2, ratio);
 
     return std::make_tuple(rotator.rotate_back(ret_P),
                            rotator.rotate_back(ret_K));
   }
 
   /* ---------------------------------------------------------------------- */
   template struct LamHomogen<twoD>;
   template struct LamHomogen<threeD>;
 
 } //muSpectre
 
diff --git a/src/materials/laminate_homogenisation.hh b/src/materials/laminate_homogenisation.hh
index 114fb50..478a004 100644
--- a/src/materials/laminate_homogenisation.hh
+++ b/src/materials/laminate_homogenisation.hh
@@ -1,206 +1,217 @@
 /**
  * @file   lamminate_hemogenization.hh
  *
  * @author Ali Falsafi <ali.falsafi@epfl.ch>
  *
  * @date   28 Sep 2018
  *
  * @brief Laminatehomogenisation enables one to obtain the resulting stress
  *        and stiffness tensors of a laminate pixel that is consisted of two
  *        materialswith a certain normal vector of their interface plane.
  *        note that it is supposed to be used in static way. so it does note
  *        any data member. It is merely a collection of functions used to
  *        calculate effective stress and stiffness.
  *
  *
  * 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 LAMINATE_HOMOGENISATION_H
 #define LAMINATE_HOMOGENISATION_H
 
 #include "common/geometry.hh"
 #include "common/common.hh"
 #include "common/field_map.hh"
 
 
 namespace muSpectre{
 
   // TODO: to template with the formulation as well
 
   template<Dim_t Dim>
   struct LamHomogen{
     //! typedefs for data handled by this interface
     using Vec_t = Eigen::Matrix<Real, Dim, 1>;
+
     using Stiffness_t = T4Mat<Real, Dim>;
     using Serial_stiffness_t = Eigen::Matrix<Real, 2*Dim-1, 2*Dim-1>;
     using Parallel_stiffness_t = T4Mat<Real, (Dim-1)>;
-
+    using Equation_stiffness_t = Eigen::Matrix<Real, Dim, Dim>;
 
     using Strain_t = Eigen::Matrix<Real, Dim*Dim, 1>;
     using Serial_strain_t = Eigen::Matrix<Real, 2*Dim-1, 1>;
     using Parallel_strain_t = Eigen::Matrix<Real, (Dim-1)*(Dim-1), 1>;
+    using Equation_strain_t = Eigen::Matrix<Real, Dim, 1>;
 
     using Stress_t = Strain_t;
     using Serial_stress_t = Serial_strain_t;
     using Parallel_stress_t = Parallel_strain_t;
+    using Equation_stress_t = Equation_strain_t ;
 
     using StrainMat_t = Eigen::Matrix<Real, Dim, Dim>;
     using StressMat_t = StrainMat_t;
 
     using Function_t = std::function<std::tuple<StressMat_t,
                                                 Stiffness_t>
                                      (const Eigen::Ref<StrainMat_t> &)>;
 
     // these two functions return the indexes in column major strain and stress
     // tensors that the behavior is either serial or parallel concerning
     // combining the laminate layers
+    static constexpr auto get_equation_indexes()
+      -> std::array<Dim_t, Dim>;
+    static constexpr auto get_serial_non_equation_indexes()
+      -> std::array<Dim_t, Dim-1>;
     static constexpr auto get_serial_indexes()
       -> std::array<Dim_t, 2*Dim-1>;
     static constexpr auto get_parallel_indexes()
       -> std::array<Dim_t, (Dim-1)*(Dim-1)>;
-    static constexpr auto get_serial_elimination_indexes()
-      -> std::array<Dim_t, (Dim-1)*(Dim-1)>;
-    static constexpr auto get_parallel_elimination_indexes()
-      -> std::array<Dim_t, 2*Dim-1>;
     static constexpr auto get_if_serial()
       -> std::array<bool, Dim*Dim>;
 
 
     // these functions return the parts of a stress or strain tensor that
     // behave either serial or parallel in a laminate structure
     // they are used to obtain a certian part of tensor which is needed by the
     // solver used in this struct to calculate the effective stress and stiffness
     // or compose the complete tensor from its sub sets.
 
     // obtain the serial part of stress or strain tensor
     static auto  get_serial_strain(const Eigen::Ref<Strain_t>& F) -> Serial_strain_t;
+    // obtain the serial part of stress or strain tensor
+    static auto  get_equation_strain_from_serial(const Eigen::Ref<Serial_strain_t>& F)
+      -> Equation_strain_t;
+    static auto get_serial_strain_from_equation(const Eigen::Ref<Equation_strain_t>& Fs)
+      -> Serial_strain_t;
     // obtain the parallel part of stress or strain tensor
     static auto  get_parallel_strain(const Eigen::Ref<Strain_t>& F)
       -> Parallel_strain_t;
 
+
     // compose the complete strain or stress tensor from
     //its serial and parallel parts
     static auto  get_total_strain(const Eigen::Ref<Serial_strain_t>& F_seri,
                                   const Eigen::Ref<Parallel_strain_t>& F_para)->Strain_t;
 
     // obtain the serial part of stifness tensor
     static auto get_serial_stiffness(const Eigen::Ref<Stiffness_t>& K)
-      ->Serial_stiffness_t;
+      -> Serial_stiffness_t;
     // obtain the parallel part of stifness  tensor
     static auto get_parallel_stiffness(const Eigen::Ref<Stiffness_t>& K)
-      ->Parallel_stiffness_t;
+      -> Parallel_stiffness_t;
 
+    static auto get_equation_stiffness_from_serial(const Eigen::Ref<Serial_stiffness_t>& Ks)
+      -> Equation_stiffness_t;
     // compose the complete stiffness tensor from
-    //its serial and parallel parts
+    // its serial and parallel parts
     static auto get_total_stiffness(const Eigen::Ref<Serial_stiffness_t>& K_seri,
                                     const Eigen::Ref<Parallel_stiffness_t>& K_para)
       ->Stiffness_t;
 
     /**
      * this functions calculate the strain in the second materials considering
      * the fact that:
      * 1. In the directions that the laminate layers are in series the weighted
      * average of the strains should be equal to the total strain of the laminate
      * 2. In the directions that the laminate layers are in parallel the strain
      * is constant and equal in all of the layers.
      */
     static auto linear_eqs(Real ratio,
                            const Eigen::Ref<Strain_t>& F_lam,
                            const Eigen::Ref<Strain_t> & F_1)->Strain_t;
 
     static auto linear_eqs_seri(Real ratio,
                                 const Eigen::Ref<Strain_t> & F_lam,
                                 const Eigen::Ref<Serial_strain_t> & F_1)
       ->Serial_strain_t;
 
     /**
      * the objective in hemgenisation of a single laminate pixel is equating the
      * stress in the serial directions so the difference of stress between their
      * layers should tend to zero. this function return the stress difference and
      * the difference of Stiffness matrices which is used as the Jacobian in the
      * solution process
      */
     static auto delta_stress_eval(const Function_t &mat_1_stress_eval,
                                   const Function_t &mat_2_stress_eval,
                                   const Eigen::Ref<StrainMat_t>& F_1,
                                   const Eigen::Ref<StrainMat_t>& F_2,
                                   RotatorNormal<Dim> rotator,
                                   Real ratio)
       ->std::tuple<Serial_stress_t, Serial_stiffness_t>;
 
 
     static auto delta_stress_eval_F_1(const Function_t &mat_1_stress_eval,
                                       const Function_t &mat_2_stress_eval,
                                       Eigen::Ref<Strain_t> F_0,
                                       Eigen::Ref<Strain_t> F_1_rot,
                                       RotatorNormal<Dim> rotator,
                                       Real ratio)
       ->std::tuple<Serial_stress_t, Serial_stiffness_t>;
     /**
      * the following function s claculates the energy dimension error of the
      * solution. it will beused in each step of the solution to detremine the
      * relevant difference that implementation of that step has had on
      * convergence to the solution.
      */
     inline static auto del_energy_eval(const Real del_F_norm,
                                        const Real delta_P_norm)
       ->Real;
 
     /**
      * This functions calculate the resultant stress and tangent matrices
      * according to the computed F_1 and F_2 from the solver.
      *
      */
     static auto lam_P_combine(Eigen::Ref<StressMat_t> P_1,
                               Eigen::Ref<StressMat_t> P_2,
                               const Real ratio)
       ->StressMat_t;
 
     static auto lam_K_combine(const Eigen::Ref<Stiffness_t>& K_1,
                               const Eigen::Ref<Stiffness_t>& K_2,
                               const Real ratio)
       ->Stiffness_t;
     /**
      * This is the main solver function that might be called staically from
      * an external file. this will return the resultant stress and stiffness
      * tensor according to interanl "equilibrium" of the lamiante.
      * The inputs are :
      * 1- global Strain
      * 2- stress calculation function of the layer 1
      * 3- stress calculation function of the layer 2
      * 4- the ratio of the first material in the laminate sturucture of the pixel
      * 5- the normal vector of the interface of two layers
      * 6- the tolerance error for the internal solution of the laminate pixel
      * 7- the maximum iterations for the internal solution of the laminate pixel
      */
     static auto laminate_solver(Eigen::Ref<Strain_t> F_coord,
                                 Function_t mat_1_stress_eval,
                                 Function_t mat_2_stress_eval,
                                 Real ratio,
                                 const Eigen::Ref<Vec_t>& normal_vec,
                                 Real tol = 1e-6,
                                 Dim_t max_iter = 1000)
       ->std::tuple<StressMat_t, Stiffness_t>;
 
   }; // LamHomogen
 
 } //muSpectre
 
 
 #endif /* LAMINATE_HOMOGENISATION_H */
diff --git a/tests/split_test_laminate_solver.cc b/tests/split_test_laminate_solver.cc
index 827c6a5..2810e9f 100644
--- a/tests/split_test_laminate_solver.cc
+++ b/tests/split_test_laminate_solver.cc
@@ -1,416 +1,416 @@
 /**
  * @file   test_laminate_solver.cc
  *
  * @author Till Junge <ali.falsafi@epfl.ch>
  *
  * @date   19 Oct 2018
  *
  * @brief  Tests for the large-strain, laminate homogenisation algorithm
  *
  * 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 <type_traits>
 
 #include <boost/mpl/list.hpp>
 #include <boost/range/combine.hpp>
 
 #include "materials/material_linear_elastic1.hh"
 #include "materials/material_linear_elastic2.hh"
 #include "materials/material_orthotropic.hh"
 #include "materials/laminate_homogenisation.hh"
 #include "tests.hh"
 #include "tests/test_goodies.hh"
 #include "common/field_collection.hh"
 #include "common/iterators.hh"
 #include "materials/materials_toolbox.hh"
 
 namespace muSpectre{
   BOOST_AUTO_TEST_SUITE(laminate_homogenisation);
 
   template <class Mat_t, Dim_t Dim, Dim_t c>
   struct MaterialFixture{
     //constructor
     MaterialFixture();
   };
 
   template<Dim_t Dim, Dim_t c>
   struct MaterialFixture<MaterialAnisotropic<Dim, Dim>, Dim, c>{
     MaterialFixture():
       mat("Name_aniso", aniso_inp_maker()){}
 
     std::vector<Real> aniso_inp_maker(){
       std::vector<Real> aniso_inp{};
       switch (Dim){
       case (twoD):{
         aniso_inp =
           {c*1.1,c*2.1,c*3.1,c*4.1,c*5.1,c*6.1};;
         break;
       }
       case (threeD):{
         aniso_inp =
           {c*1.1,c*2.1,c*3.1,c*4.1,c*5.1,c*6.1,c*7.1,c*8.1,c*9.1,c*10.1,
            c*11.1,c*12.1,c*13.1,c*14.1,c*15.1,c*16.1,c*17.1,c*18.1,c*19.1,
            c*20.1,c*21.1};
         break;
       }
       default:
         throw std::runtime_error("Unknown dimension");
       }
       return aniso_inp;
     }
 
     MaterialAnisotropic<Dim, Dim> mat;
   };
 
 
   template<Dim_t Dim, Dim_t c>
   struct MaterialFixture<MaterialOrthotropic<Dim, Dim>, Dim, c>{
     MaterialFixture():
       mat("Name_orthro", ortho_inp_maker()){}
     std::vector<Real> ortho_inp_maker(){
       std::vector<Real> ortho_inp{};
       switch (Dim){
       case (twoD):{
         ortho_inp = {c*1.1,c*2.1,c*3.1,c*4.1};
         break;
       }
       case (threeD):{
         ortho_inp = {c*1.1,c*2.1,c*3.1,c*4.1,
                      c*5.1,c*6.1,c*7.1,c*8.1,c*9.1};
         break;
       }
       default:
         throw std::runtime_error("Unknown dimension");
       }
       return ortho_inp;
     }
 
     MaterialOrthotropic<Dim, Dim> mat;
   };
 
 
 
   template<Dim_t Dim, Dim_t c>
   struct MaterialFixture<MaterialLinearElastic1<Dim, Dim>, Dim, c>{
     MaterialFixture():
     mat("Name_LinElastic1", young_maker(), poisson_maker()){}
 
     Real young_maker() {
       Real lambda{c*2}, mu{c*1.5};
       Real young{mu*(3*lambda + 2*mu)/(lambda + mu)};
       return young;
     }
 
     Real poisson_maker() {
       Real lambda{c*2}, mu{c*1.5};
       Real poisson{lambda/(2*(lambda + mu))};
       return poisson;
     }
 
     MaterialLinearElastic1<Dim, Dim> mat;
   };
 
 
 
   template <class Mat1_t, class Mat2_t, Dim_t Dim, Dim_t c1, Dim_t c2>
   struct MaterialPairFixture
   {
     using Vec_t = Eigen::Matrix<Real, Dim, 1>;
     using Stiffness_t = T4Mat<Real, Dim>;
     using Serial_stiffness_t = Eigen::Matrix<Real, 2*Dim-1, 2*Dim-1>;
     using Parallel_stiffness_t = T4Mat<Real, (Dim-1)>;
 
     using Strain_t = Eigen::Matrix<Real, Dim*Dim, 1>;
     using Serial_strain_t = Eigen::Matrix<Real, 2*Dim-1, 1>;
     using Parallel_strain_t = Eigen::Matrix<Real, (Dim-1)*(Dim-1), 1>;
 
     using Stress_t = Strain_t;
     using Serial_stress_t = Serial_strain_t;
     using Paralle_stress_t = Parallel_strain_t;
 
     using StrainMat_t = Eigen::Matrix<Real, Dim, Dim>;
     using StressMat_t = StrainMat_t;
     using T4MatMat_t = Eigen::Map<T4Mat<Real, Dim>>;
 
     using Output_t = std::tuple<StressMat_t,
                                 Stiffness_t>;
     using Function_t = std::function<Output_t(const Eigen::Ref<StrainMat_t> &)>;
 
     // constructor :
     MaterialPairFixture(){
       auto mat_fix_1 = std::unique_ptr<MaterialFixture<Mat1_t, Dim, c1>>();
       auto mat_fix_2 = std::unique_ptr<MaterialFixture<Mat2_t, Dim, c2>>();
       this->normal_vec = this->normal_vec / this->normal_vec.norm();
       //std::cout << this-> normal_vec;
     }
 
     constexpr Dim_t get_dim  ()  {return Dim;}
 
 
 
     MaterialFixture<Mat1_t, Dim, c1> mat_fix_1{};
     MaterialFixture<Mat2_t, Dim, c2> mat_fix_2{};
     // Vec_t normal_vec{Vec_t::Random()};
     Vec_t normal_vec{ Vec_t::UnitX()};
-    //Real ratio {(double) std::rand() / (RAND_MAX)};
-    Real ratio  = 0.5;
+    Real ratio {(double) std::rand() / (RAND_MAX)};
+    //Real ratio  = 0.5;
     static constexpr Dim_t fix_dim{Dim};
     };
 
   /*
   using fix_list = boost::mpl::list
     <MaterialPairFixture<MaterialLinearElastic1<threeD, threeD>,
                          MaterialLinearElastic1<threeD, threeD>,
                          threeD, 1, 1>,
      MaterialPairFixture<MaterialOrthotropic<threeD, threeD>,
                          MaterialOrthotropic<threeD, threeD>,
                          threeD, 2, 2>,
      MaterialPairFixture<MaterialAnisotropic<threeD, threeD>,
                          MaterialAnisotropic<threeD, threeD>,
                          threeD, 2, 2>,
      MaterialPairFixture<MaterialLinearElastic1<twoD, twoD>,
                          MaterialLinearElastic1<twoD, twoD>,
                          twoD, 2, 2>,
      MaterialPairFixture<MaterialOrthotropic<twoD, twoD>,
                          MaterialOrthotropic<twoD, twoD>,
                          twoD, 1, 1>,
      MaterialPairFixture<MaterialAnisotropic<twoD, twoD>,
                          MaterialAnisotropic<twoD, twoD>,
                          twoD, 2, 2>>;
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(identical_material, Fix, fix_list, Fix) {
     auto & mat1{Fix::mat_fix_1.mat};
     auto & mat2{Fix::mat_fix_2.mat};
 
     using Strain_t = typename Fix::Strain_t;
     using StrainMat_t = typename Fix::StrainMat_t;
     StrainMat_t F;
     F = StrainMat_t::Random();
     using Fucntion_t = typename Fix::Function_t;
 
     Fucntion_t mat1_evaluate_stress_func =
       [&mat1](const Eigen::Ref<StrainMat_t> & strain){
       return mat1.evaluate_stress_tangent(std::move(strain));
     };
 
     Fucntion_t mat2_evaluate_stress_func =
       [&mat2](const Eigen::Ref<StrainMat_t> & strain){
       return mat2.evaluate_stress_tangent(std::move(strain));
     };
 
     auto P_K_lam =
       LamHomogen<Fix::fix_dim>::laminate_solver(Eigen::Map<Strain_t>(F.data()),
                                                 mat1_evaluate_stress_func,
                                                 mat2_evaluate_stress_func,
                                                 Fix::ratio,
                                                 Fix::normal_vec,
                                                 1e-8,
                                                 1e5);
     auto P_lam = std::get<0> (P_K_lam);
     auto K_lam = std::get<1> (P_K_lam);
 
     auto P_K_ref = mat1_evaluate_stress_func(F);
     auto P_ref = std::get<0> (P_K_ref);
     auto K_ref = std::get<1> (P_K_ref);
 
     auto err_P{(P_lam - P_ref).norm()};
     auto err_K{(K_lam - K_ref).norm()};
 
 
     BOOST_CHECK_LT(err_P, tol);
     BOOST_CHECK_LT(err_K, tol);
   };
   */
-
+  /**
     using fix_list2 = boost::mpl::list
     <MaterialPairFixture<MaterialLinearElastic1<threeD, threeD>,
                          MaterialLinearElastic1<threeD, threeD>,
                          threeD, 7, 3>,
      MaterialPairFixture<MaterialLinearElastic1<threeD, threeD>,
                          MaterialLinearElastic1<threeD, threeD>,
                          threeD, 3, 7>,
      MaterialPairFixture<MaterialLinearElastic1<twoD, twoD>,
                          MaterialLinearElastic1<twoD, twoD>,
                          twoD, 3, 1>,
      MaterialPairFixture<MaterialLinearElastic1<twoD, twoD>,
                          MaterialLinearElastic1<twoD, twoD>,
                          twoD, 1, 3>>;
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(different_material_PK_1, Fix, fix_list2, Fix) {
     auto & mat1{Fix::mat_fix_1.mat};
     auto & mat2{Fix::mat_fix_2.mat};
     using Strain_t = typename Fix::Strain_t;
     using StrainMat_t = typename Fix::StrainMat_t;
     StrainMat_t F;
     F = StrainMat_t::Zero();
     F(0, 0) = 0.01;
     F(1, 0) = 0.005;
     F(0, 1) = 0.005;
     F(1, 1) = 0.01;
 
     using Function_t = typename Fix::Function_t;
     Dim_t Dim = Fix::get_dim();
     Function_t mat1_evaluate_stress_func =
       [Dim, &mat1](const Eigen::Ref<StrainMat_t> & strain){
       // here the material calculates the stress in its own traits conversion
       // from strain to stress:
       auto  mat_stress_tgt = mat1.evaluate_stress_tangent(std::move(strain));
 
       // here we extract the type of strain and stress that mat1 deals with
       using traits = typename std::remove_reference_t<decltype(mat1)>::traits;
       const auto strain_measure{traits::strain_measure};
       const auto stress_measure{traits::stress_measure};
 
       auto ret_fun = [Dim, mat_stress_tgt, stress_measure, strain_measure]
       (const Eigen::Ref<StrainMat_t> & strain){
 
         auto && ret_P_K = MatTB::PK1_stress<stress_measure, strain_measure>
         (strain, std::get<0>(mat_stress_tgt), std::get<1>(mat_stress_tgt));
 
         return ret_P_K;
       };
 
       return ret_fun(std::move(strain));
     };
 
     Function_t mat2_evaluate_stress_func =
       [Dim, &mat2](const Eigen::Ref<StrainMat_t> & strain){
       // here we extract the type of strain and stress that mat2 deals with
       using traits = typename std::remove_reference_t<decltype(mat2)>::traits;
       const auto strain_measure{traits::strain_measure};
       const auto stress_measure{traits::stress_measure};
 
       // here the material calculates the stress in its own traits conversion
       // from strain to stress:
       auto  mat_stress_tgt = mat2.evaluate_stress_tangent(std::move(strain));
 
 
 
       auto ret_fun = [Dim, mat_stress_tgt, stress_measure, strain_measure]
       (const Eigen::Ref<StrainMat_t> & strain){
 
         auto && ret_P_K = MatTB::PK1_stress<stress_measure, strain_measure>
         (strain, std::get<0>(mat_stress_tgt), std::get<1>(mat_stress_tgt));
 
         return ret_P_K;
       };
 
       return ret_fun(std::move(strain));
       };
     auto P_K_lam =
       LamHomogen<Fix::fix_dim>::laminate_solver(Eigen::Map<Strain_t>(F.data()),
                                                 mat1_evaluate_stress_func,
                                                 mat2_evaluate_stress_func,
                                                 Fix::ratio,
                                                 Fix::normal_vec,
                                                 1e-9,
                                                 100);
     //auto P_lam = std::get<0> (P_K_lam);
     auto K_lam = std::get<1> (P_K_lam);
 
     //auto P_K_1 = mat1_evaluate_stress_func(F);
     //auto P_1 = std::get<0> (P_K_1);
     //auto K_1 = std::get<1> (P_K_1);
 
     //auto P_K_2 = mat2_evaluate_stress_func(F);
     //auto P_2 = std::get<0> (P_K_2);
     //auto K_2 = std::get<1> (P_K_2);
 
     auto K_ref = K_lam;
     //auto P_ref = Fix::get_mean_P(P_1, P_2, Fix::ratio);
 
     auto err_K{(K_lam - K_ref).norm()};
 
     //BOOST_CHECK_LT(err_P, tol);
     BOOST_CHECK_LT(err_K, tol);
   }
 
-
+  **/
 
 
 
 
 
   using fix_list3 = boost::mpl::list
     <MaterialPairFixture<MaterialLinearElastic1<threeD, threeD>,
                          MaterialLinearElastic1<threeD, threeD>,
                          threeD, 7, 3>,
      MaterialPairFixture<MaterialLinearElastic1<threeD, threeD>,
                          MaterialLinearElastic1<threeD, threeD>,
                          threeD, 3, 7>,
      MaterialPairFixture<MaterialLinearElastic1<twoD, twoD>,
                          MaterialLinearElastic1<twoD, twoD>,
                          twoD, 3, 1>,
      MaterialPairFixture<MaterialLinearElastic1<twoD, twoD>,
                          MaterialLinearElastic1<twoD, twoD>,
                          twoD, 1, 3>>;
 
   BOOST_FIXTURE_TEST_CASE_TEMPLATE(different_material_PK_2, Fix, fix_list3, Fix) {
     auto & mat1{Fix::mat_fix_1.mat};
     auto & mat2{Fix::mat_fix_2.mat};
     using Strain_t = typename Fix::Strain_t;
     using StrainMat_t = typename Fix::StrainMat_t;
     StrainMat_t F;
     F = StrainMat_t::Zero();
     F(0, 0) = 0.01;
     F(1, 0) = 0.005;
     F(0, 1) = 0.005;
     F(1, 1) = 0.01;
 
     using Function_t = typename Fix::Function_t;
     Function_t mat1_evaluate_stress_func = [&mat1]
       (const Eigen::Ref<StrainMat_t>& strain){
       return mat1.evaluate_stress_tangent(std::move(strain));};
 
     Function_t mat2_evaluate_stress_func = [&mat2]
       (const Eigen::Ref<StrainMat_t>& strain){
       return mat2.evaluate_stress_tangent(std::move(strain));};
 
     auto P_K_lam =
       LamHomogen<Fix::fix_dim>::laminate_solver(Eigen::Map<Strain_t>(F.data()),
                                                 mat1_evaluate_stress_func,
                                                 mat2_evaluate_stress_func,
                                                 Fix::ratio,
                                                 Fix::normal_vec,
                                                 1e-9,
                                                 100);
     //auto P_lam = std::get<0> (P_K_lam);
     auto K_lam = std::get<1> (P_K_lam);
 
     //auto P_K_1 = mat1_evaluate_stress_func(F);
     //auto P_1 = std::get<0> (P_K_1);
     //auto K_1 = std::get<1> (P_K_1);
 
     //auto P_K_2 = mat2_evaluate_stress_func(F);
     //auto P_2 = std::get<0> (P_K_2);
     //auto K_2 = std::get<1> (P_K_2);
 
     auto K_ref = K_lam;
     //auto P_ref = Fix::get_mean_P(P_1, P_2, Fix::ratio);
 
     auto err_K{(K_lam - K_ref).norm()};
 
     //BOOST_CHECK_LT(err_P, tol);
     BOOST_CHECK_LT(err_K, tol);
   }
 
 
 
   BOOST_AUTO_TEST_SUITE_END();
 
 } //muSpectre