diff --git a/test/triqs/gfs/g_nu_nup.cpp b/test/triqs/gfs/g_nu_nup.cpp
index 1fb9217b..8b319d9a 100644
--- a/test/triqs/gfs/g_nu_nup.cpp
+++ b/test/triqs/gfs/g_nu_nup.cpp
@@ -12,11 +12,8 @@ int main() {
  try {
   double beta = 1;
 
-  //  auto G = gf<cartesian_product<imfreq, imfreq>, tensor_valued<4>, matrix_operations>{
-  //    {{beta, Fermion, 100}, {beta, Fermion, 100}}, {2, 2, 2, 2}};
-
-  /*
-  auto G = gf<imfreq_x2>{{beta, Fermion, 100}, {2, 2}};
+  auto G = gf<cartesian_product<imfreq, imfreq>, tensor_valued<4>>{{{beta, Fermion, 100}, {beta, Fermion, 100}}, {2, 2, 2, 2},
+                                                                   make_memory_layout(0, 2, 3, 1, 4, 5)};
 
   // try to assign to expression 
   placeholder<0> nu_;
@@ -24,16 +21,19 @@ int main() {
 
   G(nu_, nup_) << 1 / (nu_ + nup_ + 1);
 
-  for (auto w : G.mesh())
-   std::cerr  << G[w](0,0) <<std::endl;
+  //for (auto w : G.mesh())
+  // std::cerr  << G[w](0,0,0,0) <<std::endl;
+
+  // Check print the memory_layout
+  std::cerr << G.data().indexmap().get_memory_layout()<<std::endl;
 
   // inverse 
-  auto Ginv = inverse(G);
+  //auto Ginv = inverse(G);
 
   // hdf5 
   h5::file file("ess_g_nu_nup.h5", H5F_ACC_TRUNC );
   h5_write(file, "g", G);
- */
+ 
 
  }
  TRIQS_CATCH_AND_ABORT;
diff --git a/test/triqs/gfs/g_om_nu_nup.cpp b/test/triqs/gfs/g_om_nu_nup.cpp
index 53cb8a94..974591df 100644
--- a/test/triqs/gfs/g_om_nu_nup.cpp
+++ b/test/triqs/gfs/g_om_nu_nup.cpp
@@ -12,8 +12,8 @@ int main() {
  try {
   double beta = 1;
 
-  /*
-  auto G = gf<imfreq_bff, tensor_valued<4>>{{{beta, Boson, 3}, {beta, Fermion, 10}, {beta, Fermion, 10}}, {2, 2, 2, 2}};
+  auto G = gf<cartesian_product<imfreq, imfreq, imfreq>, tensor_valued<4>>{
+      {{beta, Boson, 3}, {beta, Fermion, 10}, {beta, Fermion, 10}}, {2, 2, 2, 2}};
 
   // try to assign to expression 
   placeholder<2> om_;
@@ -38,7 +38,6 @@ int main() {
   //h5_write(file, "g[0]c", gg);
   h5_write(file, "g[0]", g0);
 
- */
   }
  TRIQS_CATCH_AND_ABORT;
 }
diff --git a/triqs/gfs.hpp b/triqs/gfs.hpp
index e87be1a1..1d1ebdb6 100644
--- a/triqs/gfs.hpp
+++ b/triqs/gfs.hpp
@@ -21,6 +21,8 @@
 #ifndef TRIQS_GFS_ALL_H
 #define TRIQS_GFS_ALL_H
 
+#define TRIQS_GF_INCLUDED
+
 #include <triqs/gfs/block.hpp>
 #include <triqs/gfs/imtime.hpp>
 #include <triqs/gfs/imfreq.hpp>
diff --git a/triqs/gfs/block.hpp b/triqs/gfs/block.hpp
index 21869af2..9bbf5167 100644
--- a/triqs/gfs/block.hpp
+++ b/triqs/gfs/block.hpp
@@ -73,9 +73,11 @@ namespace gfs {
    using mesh_t = gf_mesh<block_index>;
    using gf_t = gf<block_index, Target>;
    using gf_view_t = gf_view<block_index, Target>;
+   using aux_t = nothing;
+
    struct target_shape_t {};
 
-   static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t) { return std::vector<Target>(m.size()); }
+   static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t, aux_t) { return std::vector<Target>(m.size()); }
    static typename gf_t::singularity_t make_singularity(mesh_t const &m, target_shape_t) {
     return {};
    }
diff --git a/triqs/gfs/data_proxies.hpp b/triqs/gfs/data_proxies.hpp
index 35de0686..ebb32db9 100644
--- a/triqs/gfs/data_proxies.hpp
+++ b/triqs/gfs/data_proxies.hpp
@@ -104,19 +104,6 @@ namespace triqs { namespace gfs {
   AUTO_DECL operator()(S& data, Tu const& tu) const RETURN(_impl(data, tu, triqs::tuple::_get_seq<Tu>()));
  };
 
- //---------------------------- multi variable with index mixer----------------------------------
- template <typename T, int Nvar, int TargetDim, typename IndexMixer>
- struct data_proxy_array_index_mixer : data_proxy_array_common<T, Nvar + TargetDim> {
-
-  template <typename S1, typename S2> static utility::mini_vector<int, Nvar + TargetDim> join_shape(S1 const& s1, S2 const& s2) {
-   return tuple::apply_construct_parenthesis<utility::mini_vector<int, Nvar + TargetDim>>(IndexMixer::invoke(s1, s2));
-  }
-
-  template <typename S, typename Tu> auto operator()(S& data, Tu const& tu) const {
-   return tuple::apply(data, IndexMixer::invoke(tu, tuple::make_tuple_repeat<TargetDim>(arrays::range())));
-  }
- };
-
  //---------------------------- vector ----------------------------------
 
  template<typename V> struct view_proxy : public V { 
@@ -142,15 +129,6 @@ namespace triqs { namespace gfs {
   /// The data access
   template <typename S> AUTO_DECL operator()(S& data, size_t i) const  RETURN(data[i]);
 
-  /*
-  T        &  operator()(storage_t  &           data, size_t i) { return data[i];}
-  T  const &  operator()(storage_t      const & data, size_t i) const { return data[i];}
-  Tv       &  operator()(storage_view_t &       data, size_t i)       { return data[i];}
-  Tv const &  operator()(storage_view_t const & data, size_t i) const { return data[i];}
-  Tcv       &  operator()(storage_const_view_t &       data, size_t i)       { return data[i];}
-  Tcv const &  operator()(storage_const_view_t const & data, size_t i) const { return data[i];}
-*/
-
   template<typename S, typename RHS> static void assign_to_scalar   (S & data, RHS && rhs) {for (size_t i =0; i<data.size(); ++i) data[i] = rhs;}
   template <typename ST, typename RHS> static void rebind(ST& data, RHS&& rhs) { data.clear(); for (auto & x : rhs.data()) data.push_back(x);}
  };
diff --git a/triqs/gfs/gf.hpp b/triqs/gfs/gf.hpp
index 320a1245..e9004548 100644
--- a/triqs/gfs/gf.hpp
+++ b/triqs/gfs/gf.hpp
@@ -19,7 +19,6 @@
  *
  ******************************************************************************/
 #pragma once
-#define TRIQS_GF_INCLUDED
 #include <triqs/utility/first_include.hpp>
 #include <triqs/utility/std_vector_expr_template.hpp>
 #include <triqs/utility/factory.hpp>
@@ -73,6 +72,8 @@ namespace gfs {
            typename Singularity = gf_default_singularity_t<Variable, Target>, typename Opt = void>
  using gf_const_view = gf_view<Variable, Target, Singularity, Opt, true>;
 
+ template<typename Variable> constexpr int get_n_variables(Variable) { return 1;} // default.
+
  // the implementation class
  template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsView, bool IsConst> class gf_impl;
 
@@ -462,8 +463,20 @@ namespace gfs {
 
   using target_shape_t = typename factory::target_shape_t;
 
-  gf(typename B::mesh_t m, target_shape_t shape = target_shape_t{}, typename B::indices_t const &ind = typename B::indices_t{}, std::string name = "")
-     : B(std::move(m), factory::make_data(m, shape), factory::make_singularity(m, shape), typename B::symmetry_t{}, ind, name, // clean unncessary types
+  // with aux, and indices
+  gf(typename B::mesh_t m, target_shape_t shape, typename factory::aux_t aux,
+     typename B::indices_t const &ind = typename B::indices_t{}, std::string name = "")
+     : B(std::move(m), factory::make_data(m, shape, aux), factory::make_singularity(m, shape), typename B::symmetry_t{}, ind,
+         name, // clean unncessary types
+         typename B::evaluator_t{}) {
+   if (this->_indices.is_empty()) this->_indices = typename B::indices_t(shape);
+  }
+
+  // without the aux (defaulted)
+  gf(typename B::mesh_t m, target_shape_t shape = target_shape_t{}, typename B::indices_t const &ind = typename B::indices_t{},
+     std::string name = "")
+     : B(std::move(m), factory::make_data(m, shape, typename factory::aux_t{}), factory::make_singularity(m, shape),
+         typename B::symmetry_t{}, ind, name, // clean unncessary types
          typename B::evaluator_t{}) {
    if (this->_indices.is_empty()) this->_indices = typename B::indices_t(shape);
   }
@@ -805,9 +818,10 @@ namespace gfs {
    using gf_t = gf<Var, Target, Singularity, Opt>;
    using target_shape_t = arrays::mini_vector<int, VarDim>;
    using mesh_t = typename gf_t::mesh_t;
+   using aux_t = arrays::memory_layout< get_n_variables(Var{}) + VarDim>;
    //
-   static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape) {
-    typename gf_t::data_t A(gf_t::data_proxy_t::join_shape(m.size_of_components(), shape));
+   static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape, aux_t ml) {
+    typename gf_t::data_t A(gf_t::data_proxy_t::join_shape(m.size_of_components(), shape), ml);
     A() = 0;
     return A;
    }
diff --git a/triqs/gfs/imfreq_bff.hpp b/triqs/gfs/imfreq_bff.hpp
index 5803edbd..93e95fe6 100644
--- a/triqs/gfs/imfreq_bff.hpp
+++ b/triqs/gfs/imfreq_bff.hpp
@@ -23,7 +23,9 @@
 #include "./meshes/matsubara_freq.hpp"
 namespace triqs {
 namespace gfs {
+ 
 
+ WRONGE CODE 
  // short cut. Here only the change compare to default multi var implementation
  using imfreq_bff = cartesian_product<imfreq, imfreq, imfreq>;
  //struct imfreq_bff {}; // = cartesian_product<imfreq, imfreq, imfreq>;
diff --git a/triqs/gfs/imfreq_x2.hpp b/triqs/gfs/imfreq_x2.hpp
index 215419d1..b12f2a7c 100644
--- a/triqs/gfs/imfreq_x2.hpp
+++ b/triqs/gfs/imfreq_x2.hpp
@@ -24,6 +24,8 @@
 namespace triqs {
 namespace gfs {
 
+ DEAD CODE
+
  struct matrix_operations {};
 
  // in this version, we store all frequencies, to be able to write matrix operation simply.
diff --git a/triqs/gfs/m_tail.hpp b/triqs/gfs/m_tail.hpp
index 1d5ef4e3..b894ceb3 100644
--- a/triqs/gfs/m_tail.hpp
+++ b/triqs/gfs/m_tail.hpp
@@ -73,8 +73,9 @@ namespace gfs {
    using gf_t=gf<Variable, local::tail> ;
    //using gf_view_t=gf_view<block_index, Target> ;
    struct target_shape_t {};
+   using aux_t = nothing;
 
-   static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t) { return std::vector<local::tail>(m.size()); }
+   static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t, aux_t) { return std::vector<local::tail>(m.size()); }
    static nothing make_singularity(mesh_t const &m, target_shape_t) {
     return {};
    }
diff --git a/triqs/gfs/product.hpp b/triqs/gfs/product.hpp
index 69830510..b9b0d308 100644
--- a/triqs/gfs/product.hpp
+++ b/triqs/gfs/product.hpp
@@ -32,6 +32,8 @@ namespace gfs {
   static constexpr size_t size = sizeof...(Ms);
  };
 
+ template <typename... Ms>  constexpr int get_n_variables(cartesian_product<Ms...>) { return sizeof...(Ms);} 
+ 
  // use alias
  template <typename... Ms> struct cartesian_product<std::tuple<Ms...>> : cartesian_product<Ms...> {};