 
     err = nc_put_var_double( ncid, coordsID[0], periodify(X, g2d_periodic).data());
     err = nc_put_var_double( ncid, coordsID[1], periodify(Y, g2d_periodic).data());
+
+    dg::SparseTensor<dg::HVec> metric = g2d.metric();
+    dg::HVec g_xx = metric.value(0,0), g_xy = metric.value(0,1), g_yy=metric.value(1,1);
+    dg::SparseElement<dg::HVec> vol_ = dg::tensor::volume(metric);
+    dg::HVec vol = vol_.value();
     //err = nc_put_var_double( ncid, coordsID[2], g.z().data());
     //compute and write deformation into netcdf
-    dg::blas1::pointwiseDivide( g2d.g_xy(), g2d.g_xx(), temp0);
+    dg::blas1::pointwiseDivide( g_xy, g_xx, temp0);
     const dg::HVec ones = dg::evaluate( dg::one, g2d);
-    X=g2d.g_yy();
+    X=g_yy;
     err = nc_put_var_double( ncid, defID, periodify(X, g2d_periodic).data());
     //compute and write ribeiroratio into netcdf
-    dg::blas1::pointwiseDivide( g2d.g_yy(), g2d.g_xx(), temp0);
+    dg::blas1::pointwiseDivide( g_yy, g_xx, temp0);
     X=temp0;
 
     err = nc_put_var_double( ncid, confID, periodify(X, g2d_periodic).data());
     std::cout << "Construction successful!\n";
 
     //compute error in volume element (in ribeiro grid g^xx is the volume element)
-    dg::blas1::pointwiseDot( g2d.g_xx(), g2d.g_yy(), temp0);
-    dg::blas1::pointwiseDot( g2d.g_xy(), g2d.g_xy(), temp1);
+    dg::blas1::pointwiseDot( g_xx, g_yy, temp0);
+    dg::blas1::pointwiseDot( g_xy, g_xy, temp1);
     dg::blas1::axpby( 1., temp0, -1., temp1, temp0);
-    dg::blas1::transfer( g2d.g_xx(),  temp1);
+    dg::blas1::transfer( g_xx,  temp1);
     dg::blas1::pointwiseDot( temp1, temp1, temp1);
     dg::blas1::axpby( 1., temp1, -1., temp0, temp0);
     double error = sqrt( dg::blas2::dot( temp0, w2d, temp0)/dg::blas2::dot( temp1, w2d, temp1));
     std::cout<< "Rel Error in Determinant is "<<error<<"\n";
 
     //compute error in determinant vs volume form
-    dg::blas1::pointwiseDot( g2d.g_xx(), g2d.g_yy(), temp0);
-    dg::blas1::pointwiseDot( g2d.g_xy(), g2d.g_xy(), temp1);
+    dg::blas1::pointwiseDot( g_xx, g_yy, temp0);
+    dg::blas1::pointwiseDot( g_xy, g_xy, temp1);
     dg::blas1::axpby( 1., temp0, -1., temp1, temp0);
     dg::blas1::transform( temp0, temp0, dg::SQRT<double>());
     dg::blas1::pointwiseDivide( ones, temp0, temp0);
     dg::blas1::transfer( temp0, X);
     err = nc_put_var_double( ncid, volID, periodify(X, g2d_periodic).data());
-    dg::blas1::axpby( 1., temp0, -1., g2d.vol(), temp0);
-    error = sqrt(dg::blas2::dot( temp0, w2d, temp0)/dg::blas2::dot( g2d.vol(), w2d, g2d.vol()));
+    dg::blas1::axpby( 1., temp0, -1., vol, temp0);
+    error = sqrt(dg::blas2::dot( temp0, w2d, temp0)/dg::blas2::dot( vol, w2d, vol));
     std::cout << "Rel Consistency  of volume is "<<error<<"\n";
 
     //compare g^xx to volume form
-    dg::blas1::transfer( g2d.g_xx(), temp0);
+    dg::blas1::transfer( g_xx, temp0);
     dg::blas1::pointwiseDivide( ones, temp0, temp0);
-    dg::blas1::axpby( 1., temp0, -1., g2d.vol(), temp0);
-    error=sqrt(dg::blas2::dot( temp0, w2d, temp0))/sqrt( dg::blas2::dot(g2d.vol(), w2d, g2d.vol()));
+    dg::blas1::axpby( 1., temp0, -1., vol, temp0);
+    error=sqrt(dg::blas2::dot( temp0, w2d, temp0))/sqrt( dg::blas2::dot(vol, w2d, vol));
     std::cout << "Rel Error of volume form is "<<error<<"\n";
 
-    const dg::HVec vol = dg::create::volume( g3d);
+    vol = dg::create::volume( g3d);
     dg::HVec ones3d = dg::evaluate( dg::one, g3d);
     double volume = dg::blas1::dot( vol, ones3d);
 
     std::cout << "TEST VOLUME IS:\n";
     if( psi_0 < psi_1) gp.psipmax = psi_1, gp.psipmin = psi_0;
     else               gp.psipmax = psi_0, gp.psipmin = psi_1;
-    dg::geo::Iris<Psip> iris(c.psip, gp.psipmin, gp.psipmax);
+    dg::geo::Iris iris(psip.f(), gp.psipmin, gp.psipmax);
     //dg::CylindricalGrid3d<dg::HVec> g3d( gp.R_0 -2.*gp.a, gp.R_0 + 2*gp.a, -2*gp.a, 2*gp.a, 0, 2*M_PI, 3, 2200, 2200, 1, dg::PER, dg::PER, dg::PER);
 //     dg::CartesianGrid2d g2dC( gp.R_0 -1.2*gp.a, gp.R_0 + 1.2*gp.a, -1.2*gp.a, 1.2*gp.a, 1, 1e3, 1e3, dg::PER, dg::PER);
     dg::CartesianGrid2d g2dC( gp.R_0 -2.0*gp.a, gp.R_0 + 2.0*gp.a, -2.0*gp.a, 2.0*gp.a, 1, 2e3, 2e3, dg::PER, dg::PER);

    err = nc_def_var( ncid, "volume", NC_DOUBLE, 2, dim3d, &volID);
    err = nc_def_var( ncid, "divB", NC_DOUBLE, 2, dim3d, &divBID);

    thrust::host_vector<double> psi_p = dg::pullback( psip, g2d);
    //g.display();
    err = nc_put_var_double( ncid, onesID, periodify(psi_p, g2d_periodic).data());
    dg::HVec X( g2d.size()), Y(X); //P = dg::pullback( dg::coo3, g);
    for( unsigned i=0; i<g2d.size(); i++)
    {
        X[i] = g2d.r()[i];
        Y[i] = g2d.z()[i];
    }

    dg::HVec temp0( g2d.size()), temp1(temp0);

    }
    //write parameters from file into variables
    dg::geo::solovev::GeomParameters gp(js);
    dg::geo::BinaryFunctorsLvl2 psip = dg::geo::solovev::createPsip( gp);
    std::cout << "Psi min "<<psip.f()(gp.R_0, 0)<<"\n";
    std::cout << "Type psi_0 and psi_1\n";
    double psi_0, psi_1;
    std::cin >> psi_0>> psi_1;
    gp.display( std::cout);
    dg::Timer t;
    //solovev::detail::Fpsi fpsi( gp, -10);
    std::cout << "Constructing ribeiro grid ... \n";
    t.tic();
    dg::geo::Ribeiro ribeiro( psip, psi_0, psi_1, gp.R_0, 0., 1);
    dg::CurvilinearProductGrid3d g3d(ribeiro, n, Nx, Ny,Nz, dg::DIR);
    dg::CurvilinearGrid2d g2d = g3d.perp_grid();
    dg::Grid2d g2d_periodic(g2d.x0(), g2d.x1(), g2d.y0(), g2d.y1(), g2d.n(), g2d.Nx(), g2d.Ny()+1); 
    t.toc();
    std::cout << "Construction took "<<t.diff()<<"s"<<std::endl;
