diff --git a/doc/reference/c++/gf/concepts.rst b/doc/reference/c++/gf/concepts.rst
index f0edf3cc..256081bb 100644
--- a/doc/reference/c++/gf/concepts.rst
+++ b/doc/reference/c++/gf/concepts.rst
@@ -168,7 +168,8 @@ MeshPoint concept
 +------------------------------------------------+-----------------------------------------------------------------------------+
 | void reset()                                   | Reset the mesh point to the first point                                     |
 +------------------------------------------------+-----------------------------------------------------------------------------+
-| operator mesh_t::domain_t::point_t() const     | cast to the corresponding domain point                                      |
+| cast_t                                         | type of the corresponding domain point                                      |
+| operator cast_t() const                        | cast to the corresponding domain point                                      |
 +------------------------------------------------+-----------------------------------------------------------------------------+
 | Implements the basic operations on the domain  | Only for non composite mesh                                                 |
 | by using the cast operation                    |                                                                             |
diff --git a/test/triqs/gf/gf_2times.cpp b/test/triqs/gf/gf_2times.cpp
index 005e64b1..c9b49bc0 100644
--- a/test/triqs/gf/gf_2times.cpp
+++ b/test/triqs/gf/gf_2times.cpp
@@ -34,7 +34,7 @@ int main() {
  //G2(t_,tp_) << t_ + 3*tp_;
 
  TEST( G(1,1) );
- TEST( G(G.mesh()(1,1)) );
+ TEST( G[G.mesh()(1,1) ]);
  TEST( G.on_mesh(1,1));
 
  //G2(t_,tp_) << G(tp_,t_);
diff --git a/test/triqs/gf/gf_2times_b.cpp b/test/triqs/gf/gf_2times_b.cpp
index f15f6f8a..3dbbf88c 100644
--- a/test/triqs/gf/gf_2times_b.cpp
+++ b/test/triqs/gf/gf_2times_b.cpp
@@ -20,7 +20,7 @@ int main(){
  int nt=tmax/dt;
  auto  R= make_gf<two_real_times> (tmax,nt,make_shape(1,1));//results
 
- for(auto point:R.mesh()) R(point)=0;
+ for(auto point:R.mesh()) R[point]=0;
 
  const auto rg = on_mesh(R);
  R.on_mesh(1,1) = 10;
diff --git a/test/triqs/gf/gf_re_im_freq_time.cpp b/test/triqs/gf/gf_re_im_freq_time.cpp
index 14aa9c39..579113a3 100644
--- a/test/triqs/gf/gf_re_im_freq_time.cpp
+++ b/test/triqs/gf/gf_re_im_freq_time.cpp
@@ -47,7 +47,7 @@ int main() {
  G_w.singularity()(2)=triqs::arrays::matrix<double>{{0}};
  auto G_w2 = slice_mesh_imtime(G_w_tau, index);
  for(auto& w:G_w.mesh())
-  if ( std::abs(G_w(w)-G_w2(w)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_slice error : w="<< w <<" ,G_w="<< G_w(w)<<" ,G_w2="<< G_w2(w) <<"\n";
+  if ( std::abs(G_w[w]-G_w2[w]) > precision) TRIQS_RUNTIME_ERROR<<" fourier_slice error : w="<< w <<" ,G_w="<< G_w[w]<<" ,G_w2="<< G_w2[w] <<"\n";
  
  //test of the interpolation
  std::cout << G_t_tau(0.789,0.123) << std::endl; 
@@ -70,6 +70,6 @@ int main() {
  auto G_t = inverse_fourier(G_w);
  auto G_t2 = inverse_fourier(slice_mesh_imtime(G_w_tau, index) );
  for(auto& t:G_t.mesh())
-  if ( std::abs(G_t(t)-G_t2(t)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_slice_re_time error : t="<< t <<" ,G_t="<< G_t(t) <<" ,G_t2="<< G_t2(t) <<"\n";  
+  if ( std::abs(G_t[t]-G_t2[t]) > precision) TRIQS_RUNTIME_ERROR<<" fourier_slice_re_time error : t="<< t <<" ,G_t="<< G_t[t] <<" ,G_t2="<< G_t2[t] <<"\n";  
   
 }
diff --git a/test/triqs/gf/gf_retw.cpp b/test/triqs/gf/gf_retw.cpp
index 0c733374..24717188 100644
--- a/test/triqs/gf/gf_retw.cpp
+++ b/test/triqs/gf/gf_retw.cpp
@@ -26,8 +26,8 @@ int main() {
  auto Gt2 = make_gf<retime,scalar_valued> (0, tmax, N);
 
  int i =0;
- for (auto & t : Gt.mesh()) Gt(t) = 1.0*t;
- for (auto & w : Gw.mesh()) Gw(w) = 1.0*w;
+ for (auto & t : Gt.mesh()) Gt[t] = 1.0*t;
+ for (auto & w : Gw.mesh()) Gw[w] = 1.0*w;
  
  triqs::clef::placeholder<0> t_;
  triqs::clef::placeholder<1> w_;
diff --git a/test/triqs/gf/pb_affichage.cpp b/test/triqs/gf/pb_affichage.cpp
index 8bc81a44..a993af49 100644
--- a/test/triqs/gf/pb_affichage.cpp
+++ b/test/triqs/gf/pb_affichage.cpp
@@ -20,7 +20,7 @@ int main(){
  int nt=tmax/dt;
  auto  R= make_gf<two_real_times> (tmax,nt,make_shape(1,1));//results
 
- for(auto point:R.mesh()) R(point)=0;
+ for(auto point:R.mesh()) R[point]=0;
 
  const auto rg = on_mesh(R);
  R.on_mesh(1,1) = 10;
diff --git a/test/triqs/gf/test_fourier_matsubara.cpp b/test/triqs/gf/test_fourier_matsubara.cpp
index fa915d5f..2d93a573 100644
--- a/test/triqs/gf/test_fourier_matsubara.cpp
+++ b/test/triqs/gf/test_fourier_matsubara.cpp
@@ -27,7 +27,7 @@ int main() {
  auto Gw1 = make_gf<imfreq> (beta, Fermion, make_shape(1,1), N);
  Gw1(om_) << 1/(om_-E);
 //  for(auto const& w:Gw1.mesh()){
-//   std::cout<<"w="<<std::complex<double>(w)<<", Gw1=" << Gw1(w)(0,0)<<std::endl;
+//   std::cout<<"w="<<std::complex<double>(w)<<", Gw1=" << Gw1[w](0,0)<<std::endl;
 //  }
  h5_write(file, "Gw1", Gw1);   // the original lorentzian
  
@@ -44,14 +44,14 @@ int main() {
  for(auto const& w:Gw1.mesh()){
 //   std::cout<<"w="<<std::complex<double>(w)<<",Gw1b=" << Gw1b(w)(0,0)<<std::endl;
 //   std::cout<<"w="<<std::complex<double>(w)<<",Delta Gw1b=" << Gw1b(w)(0,0)-Gw1(w)(0,0)<<std::endl;
-  if ( std::abs(Gw1b(w)(0,0)-Gw1(w)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_matsubara error : w="<<std::complex<double>(w)<<" ,Gw1b="<<std::abs(Gw1b(w)(0,0))<<"\n";
+  if ( std::abs(Gw1b[w](0,0)-Gw1[w](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_matsubara error : w="<<std::complex<double>(w)<<" ,Gw1b="<<std::abs(Gw1b[w](0,0))<<"\n";
  }
  h5_write(file,"Gw1b",Gw1b); // must be 0
  
  ///verification that the TF is OK
  for(auto const & t:Gt1.mesh()){
-  Gt1(t)-= exp(-E*t) * ( (t>0?-1:0)+1/(1+exp(E*beta)) );
-  if ( std::abs(Gt1(t)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_matsubara error : t="<<t<<" ,G1="<<std::abs(Gt1(t)(0,0))<<"\n";
+  Gt1[t]-= exp(-E*t) * ( (t>0?-1:0)+1/(1+exp(E*beta)) );
+  if ( std::abs(Gt1[t](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_matsubara error : t="<<t<<" ,G1="<<std::abs(Gt1[t](0,0))<<"\n";
  }
  h5_write(file,"Gt1b",Gt1); // must be 0
  
diff --git a/test/triqs/gf/test_fourier_real_time.cpp b/test/triqs/gf/test_fourier_real_time.cpp
index 55a3e8ce..f67a3508 100644
--- a/test/triqs/gf/test_fourier_real_time.cpp
+++ b/test/triqs/gf/test_fourier_real_time.cpp
@@ -35,7 +35,7 @@ int main() {
   
   auto Gw1 = make_gf<refreq> (-wmax, wmax, Nw, make_shape(1,1),triqs::gfs::full_bins);
   double a = Gw1.mesh().delta() * sqrt( Gw1.mesh().size() );
-  for(auto const & w:Gw1.mesh()) Gw1(w)=lorentzian(w,a);
+  for(auto const & w:Gw1.mesh()) Gw1[w]=lorentzian(w,a);
   Gw1.singularity()(2)=triqs::arrays::matrix<double>{{2.0*a}};
   h5_write(file,"Gw1",Gw1);   // the original lorentzian
   
@@ -45,15 +45,15 @@ int main() {
   // verification that TF(TF^-1)=Id
   auto Gw1b = fourier(Gt1);
   for(auto const & w:Gw1b.mesh()){
-    Gw1b(w)-=Gw1(w);
-    if ( std::abs(Gw1b(w)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : w="<<w<<" ,G1="<<std::abs(Gw1b(w)(0,0))<<"\n";
+    Gw1b[w]-=Gw1[w];
+    if ( std::abs(Gw1b[w](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : w="<<w<<" ,G1="<<std::abs(Gw1b[w](0,0))<<"\n";
   }
   h5_write(file,"Gw1b",Gw1b); // must be 0
   
   // verification that TF is the lorentzian_inverse function
   for(auto const & t:Gt1.mesh()){
-    Gt1(t)-=lorentzian_inverse(t,a);
-    if ( std::abs(Gt1(t)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : t="<<t<<" ,G1="<<std::abs(Gt1(t)(0,0))<<"\n";
+    Gt1[t]-=lorentzian_inverse(t,a);
+    if ( std::abs(Gt1[t](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : t="<<t<<" ,G1="<<std::abs(Gt1[t](0,0))<<"\n";
   }
   h5_write(file,"Gt1b",Gt1); // must be 0
   
@@ -66,8 +66,8 @@ int main() {
   
   auto Gt2 = make_gf<retime> (-tmax, tmax, Nt, make_shape(1,1));
   a = 2*acos(-1.) / ( Gt2.mesh().delta() *sqrt( Gt2.mesh().size() ) );
-  for(auto const & t:Gt2.mesh()) Gt2(t) = 0.5 *I * ( lorentzian_inverse(-t,a)*theta(-t)-lorentzian_inverse(t,a)*theta(t) );
-  //for(auto const & t:Gt2.mesh()) Gt2(t) = 0.5_j * ( lorentzian_inverse(-t,a)*theta(-t)-lorentzian_inverse(t,a)*theta(t) );
+  for(auto const & t:Gt2.mesh()) Gt2[t] = 0.5 *I * ( lorentzian_inverse(-t,a)*theta(-t)-lorentzian_inverse(t,a)*theta(t) );
+  //for(auto const & t:Gt2.mesh()) Gt2[t] = 0.5_j * ( lorentzian_inverse(-t,a)*theta(-t)-lorentzian_inverse(t,a)*theta(t) );
   Gt2.singularity()(1)=triqs::arrays::matrix<double>{{1.0}};
   h5_write(file,"Gt2",Gt2);
   
@@ -75,9 +75,9 @@ int main() {
   h5_write(file,"Gw2",Gw2);
   
   for(auto const & w:Gw2.mesh()){
-    Gw2(w)-= 0.5/(w+a*I)+0.5/(w-a*I);
-    //Gw2(w)-= 0.5/(w+a*1_j)+0.5/(w-a*1_j);
-    if ( std::abs(Gw2(w)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : w="<<w<<" ,G2="<<std::abs(Gw2(w)(0,0))<<"\n";
+    Gw2[w]-= 0.5/(w+a*I)+0.5/(w-a*I);
+    //Gw2[w]-= 0.5/(w+a*1_j)+0.5/(w-a*1_j);
+    if ( std::abs(Gw2[w](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : w="<<w<<" ,G2="<<std::abs(Gw2[w](0,0))<<"\n";
   }
   h5_write(file,"Gw2b",Gw2);
 
@@ -88,8 +88,8 @@ int main() {
   tmax=4*acos(-1.);
   
   auto Gt3 = make_gf<retime> (-tmax, tmax, Nt, make_shape(1,1));
-  for(auto const & t:Gt3.mesh()) Gt3(t) = 1.0 * std::cos(10*t) + 0.25*std::sin(4*t) + 0.5 * I*std::sin(8*t+0.3*acos(-1.)) ;
-  //for(auto const & t:Gt3.mesh()) Gt3(t) = 1.0 * std::cos(10*t) + 0.25*std::sin(4*t) + 0.5_j*std::sin(8*t+0.3*acos(-1.)) ;
+  for(auto const & t:Gt3.mesh()) Gt3[t] = 1.0 * std::cos(10*t) + 0.25*std::sin(4*t) + 0.5 * I*std::sin(8*t+0.3*acos(-1.)) ;
+  //for(auto const & t:Gt3.mesh()) Gt3[t] = 1.0 * std::cos(10*t) + 0.25*std::sin(4*t) + 0.5_j*std::sin(8*t+0.3*acos(-1.)) ;
   h5_write(file,"Gt3",Gt3);
   
   auto Gw3 = fourier(Gt3);
diff --git a/triqs/clef/clef.hpp b/triqs/clef/clef.hpp
index 82e153f9..23a31c89 100644
--- a/triqs/clef/clef.hpp
+++ b/triqs/clef/clef.hpp
@@ -382,6 +382,20 @@ namespace triqs { namespace clef {
   make_expr_call(Obj&& obj, Args &&... args) 
   { return typename result_of::make_expr_call<Obj,Args...>::type (tags::function(),std::forward<Obj>(obj), std::forward<Args>(args)...);}
 
+ /* --------------------------------------------------------------------------------------------------
+  * Create a [] call (subscript) node of an object
+  * The object can be kept as a : a ref, a copy, a view
+  * --------------------------------------------------------------------------------------------------- */
+
+ namespace result_of { 
+  template< typename Obj, typename Arg> struct make_expr_subscript : 
+   std::enable_if< is_any_lazy<Arg>::value, expr<tags::subscript,typename remove_cv_ref<Obj>::type, typename remove_cv_ref<Arg>::type> > {};
+ }
+ template< typename Obj, typename Arg>
+  typename result_of::make_expr_subscript<Obj,Arg>::type 
+  make_expr_subscript(Obj&& obj, Arg && arg) 
+  { return typename result_of::make_expr_subscript<Obj,Arg>::type (tags::subscript(),std::forward<Obj>(obj), std::forward<Arg>(arg));}
+
  /* --------------------------------------------------------------------------------------------------
   *  function class : stores any expression polymorphically
   *  f(x_,y_ ) = an expression associates this expression dynamically to f, which 
diff --git a/triqs/gf/block.hpp b/triqs/gf/block.hpp
index 137f3a9e..2d057616 100644
--- a/triqs/gf/block.hpp
+++ b/triqs/gf/block.hpp
@@ -142,7 +142,7 @@ namespace triqs { namespace gfs {
     big_gf_t big_gf;
     mesh_iterator_t mesh_it;
 
-    typename Target::view_type const & dereference() const { return big_gf(*mesh_it);}
+    typename Target::view_type const & dereference() const { return big_gf[*mesh_it];}
     bool equal(block_gf_iterator const & other) const { return ((mesh_it == other.mesh_it));}
     public:
     block_gf_iterator(gf_view<block_index,Target,Opt> bgf, bool at_end = false): big_gf(std::move(bgf)), mesh_it(&big_gf.mesh(),at_end) {}
diff --git a/triqs/gf/gf.hpp b/triqs/gf/gf.hpp
index 1cc42b05..d35adfe6 100644
--- a/triqs/gf/gf.hpp
+++ b/triqs/gf/gf.hpp
@@ -172,7 +172,6 @@ namespace triqs { namespace gfs {
      typename std::add_const<
      typename boost::lazy_disable_if<  // disable the template if one the following conditions it true
      boost::mpl::or_< // starting condition [OR]
-     std::is_base_of< tag::mesh_point, typename std::remove_reference<Arg0>::type>,  // Arg0 is (a & or a &&) to a mesh_point_t
      clef::is_any_lazy<Arg0, Args...>                          // One of Args is a lazy expression
       , boost::mpl::bool_<(sizeof...(Args)!= evaluator_t::arity -1 ) && (evaluator_t::arity !=-1)> // if -1 : no check
       >,                                                       // end of OR
@@ -188,18 +187,8 @@ namespace triqs { namespace gfs {
      operator()(Arg0 arg0, Args... args) const {
       return clef::make_expr_call(view_type(*this),arg0, args...);
      }
-
-    typedef typename std::result_of<data_proxy_t(data_t       &,size_t)>::type r_type;
-    typedef typename std::result_of<data_proxy_t(data_t const &,size_t)>::type cr_type;
-
-    r_type  operator() (mesh_point_t const & x)       { return _data_proxy(_data, x.linear_index());}
-    cr_type operator() (mesh_point_t const & x) const { return _data_proxy(_data, x.linear_index());}
-
-    template<typename ... U>
-     r_type  operator() (closest_pt_wrap<U...> const & p)       { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point<Variable,Target,Opt>::invoke(this,p)));}
-    template<typename ... U>
-     cr_type operator() (closest_pt_wrap<U...> const & p) const { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point<Variable,Target,Opt>::invoke(this,p)));}
-
+   
+    /*
     // on mesh component for composite meshes
     // enable iif the first arg is a mesh_point_t for the first component of the mesh_t
     template<typename Arg0, typename ... Args, bool MeshIsComposite = std::is_base_of<tag::composite, mesh_t>::value >
@@ -211,6 +200,30 @@ namespace triqs { namespace gfs {
      typename std::enable_if< MeshIsComposite && std::is_base_of< tag::mesh_point, Arg0>::value,  cr_type>::type
      operator() (Arg0 const & arg0, Args const & ... args) const
      { return _data_proxy(_data, _mesh.mesh_pt_components_to_linear(arg0, args...));}
+*/
+
+    //// [] and access to the grid point
+    typedef typename std::result_of<data_proxy_t(data_t       &,size_t)>::type r_type;
+    typedef typename std::result_of<data_proxy_t(data_t const &,size_t)>::type cr_type;
+
+    r_type  operator[] (mesh_index_t const & arg)       { return _data_proxy(_data,_mesh.index_to_linear(arg));}
+    cr_type operator[] (mesh_index_t const & arg) const { return _data_proxy(_data,_mesh.index_to_linear(arg));}
+
+    r_type  operator[] (mesh_point_t const & x)       { return _data_proxy(_data, x.linear_index());}
+    cr_type operator[] (mesh_point_t const & x) const { return _data_proxy(_data, x.linear_index());}
+
+    template<typename ... U>
+     r_type  operator[] (closest_pt_wrap<U...> const & p)       { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point<Variable,Target,Opt>::invoke(this,p)));}
+    template<typename ... U>
+     cr_type operator[] (closest_pt_wrap<U...> const & p) const { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point<Variable,Target,Opt>::invoke(this,p)));}
+
+    // Interaction with the CLEF library : calling the gf with any clef expression as argument build a new clef expression
+    template<typename Arg>
+     typename boost::lazy_enable_if<    // enable the template if
+     clef::is_any_lazy<Arg>,  // One of Args is a lazy expression
+     clef::result_of::make_expr_subscript<view_type,Arg>
+      >::type     // end of lazy_enable_if
+      operator[](Arg && arg) const { return clef::make_expr_subscript(view_type(*this),std::forward<Arg>(arg));}
 
     /// A direct access to the grid point
     template<typename... Args>
@@ -231,18 +244,7 @@ namespace triqs { namespace gfs {
     };
     _on_mesh_wrapper_const friend on_mesh(gf_impl const & f) { return f;}
     _on_mesh_wrapper friend on_mesh(gf_impl & f) { return f;}
-
    public:
-    r_type  operator[] (mesh_index_t const & arg)       { return _data_proxy(_data,_mesh.index_to_linear(arg));}
-    cr_type operator[] (mesh_index_t const & arg) const { return _data_proxy(_data,_mesh.index_to_linear(arg));}
-
-    // Interaction with the CLEF library : calling the gf with any clef expression as argument build a new clef expression
-    template<typename Arg>
-     typename boost::lazy_enable_if<    // enable the template if
-     clef::is_any_lazy<Arg>,  // One of Args is a lazy expression
-     clef::result_of::make_expr_call<view_type,Arg>
-      >::type     // end of lazy_enable_if
-      operator[](Arg && arg) const { return clef::make_expr_call(view_type(*this),std::forward<Arg>(arg));}
 
     //----------------------------- HDF5 -----------------------------
 
@@ -362,10 +364,11 @@ namespace triqs { namespace gfs {
 
   private:
   template<typename RHS> void triqs_clef_auto_assign_impl (RHS const & rhs, std::integral_constant<bool,false>) {
-   for (auto const & w: this->mesh()) (*this)(w) = rhs(w);
+   for (auto const & w: this->mesh()) (*this)[w] = rhs(w);
+   //for (auto const & w: this->mesh()) (*this)[w] = rhs(typename B::mesh_t::mesh_point_t::cast_t(w));
   }
   template<typename RHS> void triqs_clef_auto_assign_impl (RHS const & rhs, std::integral_constant<bool,true>) {
-   for (auto const & w: this->mesh()) (*this)(w) = triqs::tuple::apply(rhs,w.components_tuple());
+   for (auto const & w: this->mesh()) (*this)[w] = triqs::tuple::apply(rhs,w.components_tuple());
    //for (auto w: this->mesh()) triqs::tuple::apply(*this,w.components_tuple()) = triqs::tuple::apply(rhs,w.components_tuple());
   }
 
@@ -389,14 +392,14 @@ namespace triqs { namespace gfs {
  template<typename Tag, typename D, typename Target = matrix_valued> struct gf_keeper{ gf_view<D,Target> g; gf_keeper (gf_view<D,Target> const & g_) : g(g_) {} };
 
  // ---------------------------------- slicing ------------------------------------
- 
+
  //slice
  template<typename Variable, typename Target, typename Opt, bool V, typename... Args>
   gf_view<Variable,matrix_valued,Opt> slice_target (gf_impl<Variable,Target,Opt,V> const & g, Args... args) {
    static_assert(std::is_same<Target,matrix_valued>::value, "slice_target only for matrix_valued GF's");
    using arrays::range;
    //auto sg=slice_target (g.singularity(),range(args,args+1)...);
-   return gf_view<Variable,matrix_valued,Opt>(g.mesh(), g.data()(tqa::range(), args... ),  slice_target (g.singularity(),args...) , g.symmetry());
+   return gf_view<Variable,matrix_valued,Opt>(g.mesh(), g.data()(range(), args... ),  slice_target (g.singularity(),args...) , g.symmetry());
   }
 
  template<typename Variable, typename Target, typename Opt, bool V, typename... Args>
@@ -404,14 +407,14 @@ namespace triqs { namespace gfs {
    static_assert(std::is_same<Target,matrix_valued>::value, "slice_target only for matrix_valued GF's");
    using arrays::range;
    auto sg=slice_target (g.singularity(),range(args,args+1)...);
-   return gf_view<Variable,scalar_valued,Opt>(g.mesh(), g.data()(tqa::range(), args... ), sg , g.symmetry());
+   return gf_view<Variable,scalar_valued,Opt>(g.mesh(), g.data()(range(), args... ), sg , g.symmetry());
   }
 
-/*
-  template<typename Variable1,typename Variable2, typename Target, typename Opt, bool V, typename... Args>
-  gf_view<Variable2,Target,Opt> slice_mesh (gf_impl<Variable1,Target,Opt,V> const & g, Args... args) {
-   return gf_view<Variable2,Target,Opt>(g.mesh().slice(args...), g.data()(g.mesh().slice_get_range(args...),arrays::ellipsis()), g.singularity(), g.symmetry());
-  }*/
+ /*
+    template<typename Variable1,typename Variable2, typename Target, typename Opt, bool V, typename... Args>
+    gf_view<Variable2,Target,Opt> slice_mesh (gf_impl<Variable1,Target,Opt,V> const & g, Args... args) {
+    return gf_view<Variable2,Target,Opt>(g.mesh().slice(args...), g.data()(g.mesh().slice_get_range(args...),arrays::ellipsis()), g.singularity(), g.symmetry());
+    }*/
 
 }}
 
diff --git a/triqs/gf/imfreq.hpp b/triqs/gf/imfreq.hpp
index 14f7379f..c680f88e 100644
--- a/triqs/gf/imfreq.hpp
+++ b/triqs/gf/imfreq.hpp
@@ -57,6 +57,9 @@ namespace triqs { namespace gfs {
     typedef typename std::conditional < std::is_same<Target, matrix_valued>::value, arrays::matrix_view<std::complex<double>>, std::complex<double>>::type rtype; 
     template<typename G>
      rtype operator() (G const * g, long n)  const {return g->data()(n, arrays::ellipsis()); }
+    // crucial because the mesh_point is cast in a complex, not an int !
+    template<typename G>
+     rtype operator() (G const * g, linear_mesh<matsubara_domain<true>>::mesh_point_t const & p)  const { return (*this)(g,p.index());}
     template<typename G>
      local::tail_view operator()(G const * g, freq_infty const &) const {return g->singularity();}
    };
diff --git a/triqs/gf/local/fourier_matsubara.cpp b/triqs/gf/local/fourier_matsubara.cpp
index 70209f6e..92c0a043 100644
--- a/triqs/gf/local/fourier_matsubara.cpp
+++ b/triqs/gf/local/fourier_matsubara.cpp
@@ -68,15 +68,15 @@ namespace triqs { namespace gfs {
    }
    if (gw.domain().statistic == Fermion){
     for (auto & t : gt.mesh())  
-     g_in(t.index()) = fact * exp(iomega*t) * ( gt(t) - ( oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) ) );
+     g_in[t.index()] = fact * exp(iomega*t) * ( gt[t] - ( oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) ) );
    }
    else {
     for (auto & t : gt.mesh()) 
-     g_in(t.index()) =  fact * ( gt(t) -  ( oneBoson(a1,b1,t,beta) + oneBoson(a2,b2,t,beta) + oneBoson(a3,b3,t,beta) ) );  
+     g_in[t.index()] =  fact * ( gt[t] -  ( oneBoson(a1,b1,t,beta) + oneBoson(a2,b2,t,beta) + oneBoson(a3,b3,t,beta) ) );  
    }
    details::fourier_base(g_in, g_out, L, true);
    for (auto & w : gw.mesh()) {
-    gw(w) = g_out( w.index() ) * exp(iomega2*w.index() ) + a1/(w-b1) + a2/(w-b2) + a3/(w-b3); 
+    gw[w] = g_out( w.index() ) * exp(iomega2*w.index() ) + a1/(w-b1) + a2/(w-b2) + a3/(w-b3); 
    }
    gw.singularity() = gt.singularity();// set tail
   }
@@ -110,7 +110,7 @@ namespace triqs { namespace gfs {
    }
    g_in() = 0;
    for (auto & w: gw.mesh()) {
-    g_in( w.index() ) =  fact * exp(w.index()*iomega2) * ( gw(w) - (a1/(w-b1) + a2/(w-b2) + a3/(w-b3)) );
+    g_in[ w.index() ] =  fact * exp(w.index()*iomega2) * ( gw[w] - (a1/(w-b1) + a2/(w-b2) + a3/(w-b3)) );
    }
    // for bosons GF(w=0) is divided by 2 to avoid counting it twice
    if (gw.domain().statistic == Boson && !Green_Function_Are_Complex_in_time ) g_in(0) *= 0.5;
@@ -122,13 +122,13 @@ namespace triqs { namespace gfs {
    //typedef typename gf<imtime>::mesh_type::gf_result_type gt_result_type;
    if (gw.domain().statistic == Fermion){
     for (auto & t : gt.mesh()){
-     gt(t) = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() ) * exp(-iomega*t)
+     gt[t] = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() ) * exp(-iomega*t)
      + oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) );
     }
    }
    else {
     for (auto & t : gt.mesh()) 
-     gt(t) = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() )
+     gt[t] = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() )
      + oneBoson(a1,b1,t,beta) + oneBoson(a2,b2,t,beta) + oneBoson(a3,b3,t,beta) );
    }
    if (gt.mesh().kind() == full_bins)
diff --git a/triqs/gf/local/fourier_real.cpp b/triqs/gf/local/fourier_real.cpp
index 02804c5a..504c43dc 100644
--- a/triqs/gf/local/fourier_real.cpp
+++ b/triqs/gf/local/fourier_real.cpp
@@ -58,12 +58,12 @@ namespace triqs { namespace gfs {
    dcomplex a1 = (t1 + I * t2/a )/2., a2 = (t1 - I * t2/a )/2.;
    
    for (auto const & t : gt.mesh())
-    g_in(t.index()) = (gt(t) - (a1*th_expo(t,a) + a2*th_expo_neg(t,a))) * std::exp(I*t*wmin);
+    g_in[t.index()] = (gt[t] - (a1*th_expo(t,a) + a2*th_expo_neg(t,a))) * std::exp(I*t*wmin);
    
    details::fourier_base(g_in, g_out, L, true);
    
    for (auto const & w : gw.mesh())
-    gw(w) = gt.mesh().delta() * std::exp(I*(w-wmin)*tmin) * g_out(w.index())
+    gw[w] = gt.mesh().delta() * std::exp(I*(w-wmin)*tmin) * g_out(w.index())
     + a1*th_expo_inv(w,a) + a2*th_expo_neg_inv(w,a);
    
    gw.singularity() = gt.singularity();// set tail
@@ -90,14 +90,14 @@ namespace triqs { namespace gfs {
    g_in() = 0;
    
    for (auto const & w: gw.mesh())
-    g_in(w.index()) = (gw(w) - a1*th_expo_inv(w,a) - a2*th_expo_neg_inv(w,a)  ) * std::exp(-I*w*tmin);
+    g_in(w.index()) = (gw[w] - a1*th_expo_inv(w,a) - a2*th_expo_neg_inv(w,a)  ) * std::exp(-I*w*tmin);
    
    details::fourier_base(g_in, g_out, L, false);
    
    const double corr = 1.0/(gt.mesh().delta()*L);
    for (auto const & t : gt.mesh())
-    gt(t) = corr * std::exp(I*wmin*(tmin-t)) *
-    g_out( t.index() ) + a1 * th_expo(t,a) + a2 * th_expo_neg(t,a) ;
+    gt[t] = corr * std::exp(I*wmin*(tmin-t)) *
+    g_out[ t.index() ] + a1 * th_expo(t,a) + a2 * th_expo_neg(t,a) ;
  
    // set tail
    gt.singularity() = gw.singularity();
diff --git a/triqs/gf/local/functions.cpp b/triqs/gf/local/functions.cpp
index 894c854e..f72982b6 100644
--- a/triqs/gf/local/functions.cpp
+++ b/triqs/gf/local/functions.cpp
@@ -44,7 +44,7 @@ namespace triqs { namespace gfs {
     dcomplex d= t(1)(n1,n2) , A=t(2)(n1,n2),B = t(3)(n1,n2) ;
     double b1 = 0,b2 =1, b3 =-1;
     dcomplex a1 = d-B, a2 = (A+B)/2, a3 = (B-A)/2;
-    for (auto & w : G.mesh())  dens_part(n1,n2)+= G(w)(n1,n2)  -  (a1/(w - b1) + a2 / (w-b2) + a3/(w-b3)); 
+    for (auto & w : G.mesh())  dens_part(n1,n2)+= G[w](n1,n2)  -  (a1/(w - b1) + a2 / (w-b2) + a3/(w-b3)); 
     dens_part(n1,n2) = dens_part(n1,n2)/Beta;
     dens_tail(n1,n2) = d + F(a1,b1,Beta) + F(a2,b2,Beta)+ F(a3,b3,Beta);
     // If  the Green function are NOT complex, then one use the symmetry property
@@ -77,7 +77,7 @@ namespace triqs { namespace gfs {
    res() = 0.0;
 
    for (auto l : gl.mesh()) {
-     res -= sqrt(2*l.index()+1) * gl(l);
+     res -= sqrt(2*l.index()+1) * gl[l];
    }
    res /= gl.domain().beta;
 
@@ -95,7 +95,7 @@ namespace triqs { namespace gfs {
 
    for (int p=1; p<=t.order_max(); p++)
      for (auto l : gl.mesh())
-       t(p) += (triqs::utility::legendre_t(l.index(),p)/pow(gl.domain().beta,p)) * gl(l);
+       t(p) += (triqs::utility::legendre_t(l.index(),p)/pow(gl.domain().beta,p)) * gl[l];
 
    return t;
 
@@ -113,7 +113,7 @@ namespace triqs { namespace gfs {
 
    tqa::array<double,2> corr(disc.shape()); corr() = 0;
    for (auto l : gl.mesh()) {
-     corr += t(l.index()) * gl(l);
+     corr += t(l.index()) * gl[l];
    }
 
    tqa::range R;
diff --git a/triqs/gf/local/legendre_matsubara.cpp b/triqs/gf/local/legendre_matsubara.cpp
index db014b10..f62ec6ef 100644
--- a/triqs/gf/local/legendre_matsubara.cpp
+++ b/triqs/gf/local/legendre_matsubara.cpp
@@ -37,7 +37,7 @@ void legendre_matsubara_direct(gf_view<imfreq> & gw, gf_view<legendre> const & g
   // Use the transformation matrix
   for (auto om: gw.mesh()) {
     for (auto l: gl.mesh()) {
-      gw(om) += legendre_T(om.index(),l.index()) * gl(l);
+      gw[om] += legendre_T(om.index(),l.index()) * gl[l];
     }
   }
 
@@ -72,7 +72,7 @@ void legendre_matsubara_direct (gf_view<imtime> & gt, gf_view<legendre> const &
   for (auto t : gt.mesh()) {
     L.reset( 2*t/gt.domain().beta-1 );
     for (auto l : gl.mesh()) {
-      gt(t) += sqrt(2*l.index()+1) / gt.domain().beta * gl(l) * L.next();
+      gt[t] += sqrt(2*l.index()+1) / gt.domain().beta * gl[l] * L.next();
     }
   }
 
@@ -89,7 +89,7 @@ void legendre_matsubara_inverse (gf_view<legendre> & gl, gf_view<imtime> const &
   for (auto t : gt.mesh()) {
     L.reset( 2*t/gt.domain().beta-1 );
     for (auto l : gl.mesh()) {
-      gl(l) += sqrt(2*l.index()+1) * L.next() * gt(t);
+      gl[l] += sqrt(2*l.index()+1) * L.next() * gt[t];
     }
   }
   gl.data() *= gt.mesh().delta();
diff --git a/triqs/gf/local/pade.cpp b/triqs/gf/local/pade.cpp
index 368b44b7..518c8616 100644
--- a/triqs/gf/local/pade.cpp
+++ b/triqs/gf/local/pade.cpp
@@ -53,7 +53,7 @@ namespace triqs { namespace gfs {
       gr() = 0.0;
       for (auto om : gr.mesh()) {
         dcomplex e = om + dcomplex(0.0,1.0)*freq_offset;
-        gr(om)(n1,n2) = PA(e);
+        gr[om](n1,n2) = PA(e);
       }
 
     }
diff --git a/triqs/gf/meshes/discrete.hpp b/triqs/gf/meshes/discrete.hpp
index 5fb64255..96585a4f 100644
--- a/triqs/gf/meshes/discrete.hpp
+++ b/triqs/gf/meshes/discrete.hpp
@@ -49,7 +49,8 @@ namespace triqs { namespace gfs {
     public:
     mesh_point_t( discrete_mesh const & mesh, index_t const & index_): m(&mesh), _index(index_) {}
     void advance() { ++_index;}
-    operator size_t () const { return m->index_to_point(_index);} 
+    typedef size_t cast_t;
+    operator cast_t() const { return m->index_to_point(_index);} 
     size_t linear_index() const { return _index;}
     size_t index() const { return _index;}
     bool at_end() const { return (_index == m->size());}
diff --git a/triqs/gf/meshes/linear.hpp b/triqs/gf/meshes/linear.hpp
index bbad7127..5d25380c 100644
--- a/triqs/gf/meshes/linear.hpp
+++ b/triqs/gf/meshes/linear.hpp
@@ -78,11 +78,12 @@ namespace triqs { namespace gfs {
    /// The wrapper for the mesh point
    class mesh_point_t : tag::mesh_point, public arith_ops_by_cast<mesh_point_t, domain_pt_t  > {
     linear_mesh const * m;  
-    index_t _index; 
+    index_t _index;
     public:
     mesh_point_t( linear_mesh const & mesh, index_t const & index_): m(&mesh), _index(index_) {}
     void advance() { ++_index;}
-    operator domain_pt_t () const { return m->index_to_point(_index);} 
+    typedef domain_pt_t cast_t; 
+    operator cast_t () const { return m->index_to_point(_index);} 
     size_t linear_index() const { return _index;}
     size_t index() const { return _index;}
     bool at_end() const { return (_index == m->size());}
diff --git a/triqs/gfs.hpp b/triqs/gfs.hpp
index 18dafa2a..0df2e09c 100644
--- a/triqs/gfs.hpp
+++ b/triqs/gfs.hpp
@@ -22,11 +22,12 @@
 #define TRIQS_GFS_ALL_H
 
 // The basic classes
+#include <triqs/gf/block.hpp>
 #include <triqs/gf/imtime.hpp>
 #include <triqs/gf/imfreq.hpp>
 #include <triqs/gf/retime.hpp>
 #include <triqs/gf/refreq.hpp>
-#include <triqs/gf/block.hpp>
+#include <triqs/gf/legendre.hpp>
 //#include <triqs/gf/product.hpp>