#ifndef DUNE_MICROSTRUCTURE_VTK_FILLER_HH
#define DUNE_MICROSTRUCTURE_VTK_FILLER_HH

#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
//#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/microstructure/prestrain_material_geometry.hh>


	using std::string;
	using namespace Dune;
	using namespace Functions;
    using namespace Functions::BasisFactory;

	//static const int dimworld = 3;

	//using GT = UGGrid<dim>; //typename MicrostructuredRodGrid::GT<dim>;
	//using GV = typename GT::LeafGridView;
	/*
	using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
	using VectorRT = FieldVector< double, dimworld>;
	using ScalarRT = double;


	using ScalarCT = std::vector<double>;
	using VectorCT = BlockVector<FieldVector<double,1> >;

	using HierarchicVectorView = Dune::Functions::HierarchicVectorWrapper< VectorCT, double>;
*/
	/*
	template<class Basis>
	static auto scalarDiscGlobalBasisFunc(
		const Basis& feBasis,
		std::function< double(const typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate&) >& scalarFunc)
	{
		ScalarCT values;
    	Functions::interpolate(feBasis, values, scalarFunc);
    	auto discGlobalBasisFunc = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(feBasis, values);
    	return discGlobalBasisFunc;
	}

	template<class Basis>
	static auto vectorDiscGlobalBasisFunc(
		const Basis& feBasis,
		std::function< FieldVector< double, dimworld>(const typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate&) >& vectorFunc)
	{
		VectorCT values;
    	Functions::interpolate(feBasis, values, vectorFunc);
    	auto discGlobalBasisFunc = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
    		(Functions::subspaceBasis(feBasis), HierarchicVectorView(values));
    	return discGlobalBasisFunc;
	}
*/

template <int dim>
class FTKfillerContainer {

public:

	static const int dimworld = 3;

// 	using GT = UGGrid<dim>; //typename MicrostructuredRodGrid::GT<dim>;
    using GT =YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
	using GV = typename GT::LeafGridView;
	using Domain = typename GV::template Codim<0>::Geometry::GlobalCoordinate;
	using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
	using VectorRT = FieldVector< double, dimworld>;
	using ScalarRT = double;
	using FuncMatrix = std::function< MatrixRT(const Domain&) >;
	using FuncVector = std::function< VectorRT(const Domain&) >;
	using FuncScalar = std::function< ScalarRT(const Domain&) >;

	using ScalarCT = std::vector<double>;
	using VectorCT = BlockVector<FieldVector<double,1> >;
	using HierarchicVectorView = Dune::Functions::HierarchicVectorWrapper< VectorCT, double>;


    void vtkPrestrainNorm(const GV& gridView, const FuncMatrix& B, string filename)          // Updated
    {
	    VTKWriter<typename GT::LeafGridView>  vtkWriter(gridView);
	    Functions::LagrangeBasis<typename GT::LeafGridView, 1> feBasis(gridView);
	    FuncScalar normB = PrestrainImp<dim>().getNormPrestrain(B);

// 	    auto discGlobalBasisFunc = scalarDiscGlobalBasisFunc(feBasis, normB);
        auto discGlobalBasisFunc  = Dune::Functions::makeGridViewFunction(normB, gridView);
//     auto muLocal = localFunction(muGridF);


	    ScalarCT values;
    	Functions::interpolate(feBasis, values, normB);
    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(feBasis, values);
	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
    	vtkWriter.write(filename);

    }

    void vtkPrestrainDirector(const typename GT::LeafGridView& gridView, const FuncVector& b, string filename)
    {
    	VTKWriter<typename GT::LeafGridView>  vtkWriter(gridView);

        using namespace Functions::BasisFactory;

    	auto basis_CE = makeBasis(
	    gridView,
	    power<dimworld>(
      	lagrange<1>(),
      	flatLexicographic()));//flatInterleaved()

    	VectorCT values;
    	Functions::interpolate(basis_CE, values, b);
    	auto b_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_CE, values);

      	vtkWriter.addVertexData(b_DGBF, VTK::FieldInfo("prestrain_director", VTK::FieldInfo::Type::vector, dimworld));
    	vtkWriter.write(filename);
    }


    void vtkMaterialCell(const GV& gridView_CE,
    	const FuncScalar& mu_Func,
    	std::string filename)
    {

    	std::cout << "vtkMaterialCell ...\n";

	    VTKWriter<GV>  vtkWriter(gridView_CE);

	    // Scalar data
	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);

	    ScalarCT material_Coeff;
    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, material_Coeff);
	    vtkWriter.addCellData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));

	    vtkWriter.write(filename);

	    std::cout << "vtkMaterialCell done.\n";
    }


 	  void vtkMaterial3DRod(const GV& gridView_R3D, double eps,
    	const FuncVector& domain_Func, const FuncScalar& mu_Func,
    	std::string filename)
    {

    	std::cout << "vtkProblem ...\n";

	    VTKWriter<GV>  vtkWriter(gridView_R3D);

	    // Vector data
	    const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()

	    VectorCT domain_Coeff;
	    Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
	    auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
	    vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));

	    // Scalar data
	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);

	    ScalarCT material_Coeff;
    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, material_Coeff);
	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));

	    vtkWriter.write(filename);

	    std::cout << "vtkProblem done.\n";
    }


    void vtkProblemCell(const GV& gridView_CE,
    	const FuncMatrix& B_Func, const FuncScalar& mu_Func,
    	std::string filename)
    {

    	std::cout << "vtkProblemCell ...\n";

	    VTKWriter<GV>  vtkWriter(gridView_CE);

	    // Vector data
	    const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()


    	FuncVector Be1_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[0];
    	};
    	VectorCT Be1_Coeff;
	    Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
	    auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
	    vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));

	    FuncVector Be2_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[1];
    	};
    	VectorCT Be2_Coeff;
	    Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
	    auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
	    vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be3_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[2];
    	};
    	VectorCT Be3_Coeff;
	    Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
	    auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
	    vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));

	    // Scalar data
	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
	    FuncScalar normB_Func = PrestrainImp<dim>().getNormPrestrain(B_Func);
	    ScalarCT prestrain_Coeff;
    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, prestrain_Coeff);
	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));

	    ScalarCT material_Coeff;
    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, material_Coeff);
	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));

	    vtkWriter.write(filename);

	    std::cout << "vtkProblem done.\n";
    }

    void vtkProblemCellBVec(const GV& gridView_CE,
    	const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
    	std::string filename)
    {

    	std::cout << "vtkProblemCellBVec ...\n";

	    VTKWriter<GV>  vtkWriter(gridView_CE);

	    // Vector data
	    const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()


    	FuncVector Be1_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[0];
    	};
    	VectorCT Be1_Coeff;
	    Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
	    auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
	    vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be2_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[1];
    	};
    	VectorCT Be2_Coeff;
	    Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
	    auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
	    vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be3_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[2];
    	};
    	VectorCT Be3_Coeff;
	    Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
	    auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
	    vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));


    	VectorCT B_vec_Coeff;
		Functions::interpolate(basis_CE, B_vec_Coeff, B_vec_Func);
		auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_CE, B_vec_Coeff);
  		vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));


	    // Scalar data
	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);

	    FuncScalar normB_Func = PrestrainImp<dim>().getNormPrestrain(B_Func);
	    ScalarCT prestrain_Coeff;
    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, prestrain_Coeff);
	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));

	    ScalarCT material_Coeff;
    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, material_Coeff);
	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));

	    vtkWriter.write(filename);

	    std::cout << "vtkProblem done.\n";
    }


    void vtkProblem3DRod(const GV& gridView_R3D, double eps,
    	const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncScalar& mu_Func,
    	std::string filename)
    {

    	std::cout << "vtkProblem ...\n";

	    VTKWriter<GV>  vtkWriter(gridView_R3D);

	    // Vector data
	    const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()

	    VectorCT domain_Coeff;
	    Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
	    auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
	    vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));


    	FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[0];
    	};
    	VectorCT Be1_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
	    auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
	    vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[1];
    	};
    	VectorCT Be2_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
	    auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
	    vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[2];
    	};
    	VectorCT Be3_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
	    auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
	    vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));



	    // Scalar data
	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);

	    FuncScalar normB_Func = PrestrainImp<dim>().getNormPrestrain(B_Func);
	    FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
	    {
	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
	    	return normB_Func(x_Ce);
	    };
	    ScalarCT prestrain_Coeff;
    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, prestrain_Coeff);
	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));

	    ScalarCT material_Coeff;
    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, material_Coeff);
	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));

	    vtkWriter.write(filename);

	    std::cout << "vtkProblem done.\n";
    }

     void vtkProblem3DRodBVec(const GV& gridView_R3D, double eps,
    	const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
    	std::string filename)
    {

    	std::cout << "vtkProblem ...\n";

	    VTKWriter<GV>  vtkWriter(gridView_R3D);

	    // Vector data
	    const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()

	    VectorCT domain_Coeff;
	    Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
	    auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
	    vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));


    	FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[0];
    	};
    	VectorCT Be1_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
	    auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
	    vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[1];
    	};
    	VectorCT Be2_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
	    auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
	    vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[2];
    	};
    	VectorCT Be3_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
	    auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
	    vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));

	    FuncVector B_vec_R3D_Func = [B_vec_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_vec_Func(x_Ce);
    	};
	    VectorCT B_vec_Coeff;
		Functions::interpolate(basis_R3D, B_vec_Coeff, B_vec_R3D_Func);
		auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_R3D, B_vec_Coeff);
  		vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));



	    // Scalar data
	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);

	    FuncScalar normB_Func = PrestrainImp<dim>().getNormPrestrain(B_Func);
	    FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
	    {
	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
	    	return normB_Func(x_Ce);
	    };
	    ScalarCT prestrain_Coeff;
    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, prestrain_Coeff);
	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));

	    FuncScalar mu_R3D_Func = [mu_Func, eps](const auto& x)
	    {
	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
	    	return mu_Func(x_Ce);
	    };
	    ScalarCT material_Coeff;
    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_R3D_Func);
    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, material_Coeff);
	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));

	    vtkWriter.write(filename);

	    std::cout << "vtkProblem done.\n";
    }

};

/*
	template<int dim>
	void vtkProblemCell(const typename UGGrid<dim>::LeafGridView& gridView_CE,
    	const std::function< MatrixRT(const typename UGGrid<dim>::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate&)& B_Func,
    	const std::function< ScalarRT(const typename UGGrid<dim>::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate&)& mu_Func,
    	std::string filename)
    {
    	using GT = UGGrid<dim>;
    	using GV = typename GT::LeafGridView;
    	using Domain = typename GV::template Codim<0>::Geometry::GlobalCoordinate;
    	using FuncMatrix = std::function< MatrixRT(const Domain&) >;
		using FuncVector = std::function< VectorRT(const Domain&) >;
		using FuncScalar = std::function< ScalarRT(const Domain&) >;
    	std::cout << "vtkProblemCell ...\n";

	    VTKWriter<GV>  vtkWriter(gridView_CE);

	    // Vector data
	    const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()


    	FuncVector Be1_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[0];
    	};
    	VectorCT Be1_Coeff;
	    Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
	    auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
	    vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be2_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[1];
    	};
    	VectorCT Be2_Coeff;
	    Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
	    auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
	    vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be3_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[2];
    	};
    	VectorCT Be3_Coeff;
	    Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
	    auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
	    vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));


	    // Scalar data
	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);

	    FuncScalar normB_Func = getNormPrestrain(B_Func);
	    ScalarCT prestrain_Coeff;
    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, prestrain_Coeff);
	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));

	    ScalarCT material_Coeff;
    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, material_Coeff);
	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));

	    vtkWriter.write(filename);

	    std::cout << "vtkProblem done.\n";
    }

    void vtkProblemCellBVec(const GV& gridView_CE,
    	const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
    	std::string filename)
    {

    	std::cout << "vtkProblemCell ...\n";

	    VTKWriter<GV>  vtkWriter(gridView_CE);

	    // Vector data
	    const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()


    	FuncVector Be1_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[0];
    	};
    	VectorCT Be1_Coeff;
	    Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
	    auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
	    vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be2_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[1];
    	};
    	VectorCT Be2_Coeff;
	    Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
	    auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
	    vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be3_Func = [B_Func](const auto& x)
    	{
    		return B_Func(x)[2];
    	};
    	VectorCT Be3_Coeff;
	    Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
	    auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
	    vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));


    	VectorCT B_vec_Coeff;
		Functions::interpolate(basis_CE, B_vec_Coeff, B_vec_Func);
		auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_CE, B_vec_Coeff);
  		vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));


	    // Scalar data
	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);

	    FuncScalar normB_Func = getNormPrestrain(B_Func);
	    ScalarCT prestrain_Coeff;
    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, prestrain_Coeff);
	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));

	    ScalarCT material_Coeff;
    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, material_Coeff);
	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));

	    vtkWriter.write(filename);

	    std::cout << "vtkProblem done.\n";
    }


    void vtkProblem3DRod(const GV& gridView_R3D, double eps,
    	const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncScalar& mu_Func,
    	std::string filename)
    {

    	std::cout << "vtkProblem ...\n";

	    VTKWriter<GV>  vtkWriter(gridView_R3D);

	    // Vector data
	    const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()

	    VectorCT domain_Coeff;
	    Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
	    auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
	    vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));


    	FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[0];
    	};
    	VectorCT Be1_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
	    auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
	    vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[1];
    	};
    	VectorCT Be2_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
	    auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
	    vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[2];
    	};
    	VectorCT Be3_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
	    auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
	    vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));



	    // Scalar data
	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);

	    FuncScalar normB_Func = getNormPrestrain(B_Func);
	    FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
	    {
	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
	    	return normB_Func(x_Ce);
	    };
	    ScalarCT prestrain_Coeff;
    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, prestrain_Coeff);
	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));

	    ScalarCT material_Coeff;
    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, material_Coeff);
	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));

	    vtkWriter.write(filename);

	    std::cout << "vtkProblem done.\n";
    }

     void vtkProblem3DRodBVec(const GV& gridView_R3D, double eps,
    	const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
    	std::string filename)
    {

    	std::cout << "vtkProblem ...\n";

	    VTKWriter<GV>  vtkWriter(gridView_R3D);

	    // Vector data
	    const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()

	    VectorCT domain_Coeff;
	    Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
	    auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
	    vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));


    	FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[0];
    	};
    	VectorCT Be1_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
	    auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
	    vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[1];
    	};
    	VectorCT Be2_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
	    auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
	    vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));


	    FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_Func(x_Ce)[2];
    	};
    	VectorCT Be3_R3D_Coeff;
	    Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
	    auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
	        (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
	    vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));

	    FuncVector B_vec_R3D_Func = [B_vec_Func, eps](const auto& x)
    	{
    		auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
    		return B_vec_Func(x_Ce);
    	};
	    VectorCT B_vec_Coeff;
		Functions::interpolate(basis_R3D, B_vec_Coeff, B_vec_R3D_Func);
		auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_R3D, B_vec_Coeff);
  		vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));



	    // Scalar data
	    Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);

	    FuncScalar normB_Func = getNormPrestrain(B_Func);
	    FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
	    {
	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
	    	return normB_Func(x_Ce);
	    };
	    ScalarCT prestrain_Coeff;
    	Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
    	auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, prestrain_Coeff);
	    vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));

	    FuncScalar mu_R3D_Func = [mu_Func, eps](const auto& x)
	    {
	    	auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
	    	return mu_Func(x_Ce);
	    };
	    ScalarCT material_Coeff;
    	Functions::interpolate(scalarFeBasis, material_Coeff, mu_R3D_Func);
    	auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
    		(scalarFeBasis, material_Coeff);
	    vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));

	    vtkWriter.write(filename);

	    std::cout << "vtkProblem done.\n";
    }
*/

#endif