#ifndef DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
#define DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH


#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/common/parametertree.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/microstructure/matrix_operations.hh>
#include <dune/microstructure/materialDefinitions.hh>


using namespace Dune;
using namespace MatrixOperations;
using std::pow;
using std::abs;
using std::sqrt;
using std::sin;
using std::cos;

using std::shared_ptr;
using std::make_shared;





//  template<class Domain>
//  MatrixRT MAT(const MatrixRT& G, const Domain& x)      //unfortunately not possible as a GridViewFunction? 
//   {
//       const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
//       const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
//       double theta=0.25;
//       if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
//          return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
//       else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
//          return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
//       else 
//          return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
//   }

template <class GridView>     // needed for GridViewFunctions
class prestrainedMaterial
{

public:
  static const int dimworld = 3; //GridView::dimensionworld;
  static const int dim = 3; //const int dim = Domain::dimension;

  using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
  using LocalDomain = typename GridView::template Codim<0>::Geometry::LocalCoordinate;
  using Element = typename GridViewEntitySet<GridView, 0>::Element;
  using ScalarRT = FieldVector< double, 1>;
  using VectorRT = FieldVector< double, dimworld>;
  using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
  using FuncScalar = std::function< double(const Domain&) >;
  using Func2int = std::function< int(const Domain&) >;
  using Func2Tensor = std::function< MatrixRT(const Domain&) >;
  using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
  using Func2TensorPhase = std::function< MatrixRT(const MatrixRT& ,const int&) >;
  using MatrixFunc  = std::function< MatrixRT(const MatrixRT&) >;
  using MatrixPhase = std::function< MatrixRT(const int&) >;

protected:

  const GridView& gridView_;
  const ParameterTree& parameterSet_;
  std::string materialFunctionName_;

  // MatrixFunc L1_;
  // MatrixFunc L2_;
  // MatrixFunc L3_;
  // const FuncScalar mu_; 
  // const FuncScalar lambda_; 
  // const int phases_;
  // Func2Tensor elasticityTensor_;
  // Func2TensorParam elasticityTensor_;
  Func2TensorPhase elasticityTensor_;
  MatrixPhase prestrain_;


  Func2int indicatorFunction_;
  GridViewFunction<int(const Domain&), GridView> indicatorFunctionGVF_;
  LocalFunction<int(const Domain&), typename GridViewEntitySet<GridView, 0>::Element >  localIndicatorFunction_;

public: 
    ///////////////////////////////
    // constructor
    ///////////////////////////////
    prestrainedMaterial(const GridView gridView,
                        const ParameterTree& parameterSet)  
                       : gridView_(gridView), 
                         parameterSet_(parameterSet)
    {
	  std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material_1");
    std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;

    setupMaterial(materialFunctionName_);

    indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
    localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);

    // Python::Module module = Python::import(materialFunctionName_);
    // elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
    // elasticityTensor_ = setupElasticityTensor(materialFunctionName_);
    // auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction"));
    // indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
    // indicatorFunction_ =  Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);
    // module.get("Phases").toC<int>(Phases_);
    // L1_ = Python::make_function<MatrixRT>(module.get("L1"));
    // L2_ = Python::make_function<MatrixRT>(module.get("L2"));
    // L3_ = Python::make_function<MatrixRT>(module.get("L3"));
    // bool isotropic_ = true; // read from module File 
    // auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working 
    } 

  //---function that determines elasticity Tensor 
    void setupMaterial(const std::string name)  // cant use materialFunctionName_ here!?
    {
        // std::cout << "Using material definition:" << name << std::endl;
        if(name == "material_neukamm")
        {   
          indicatorFunction_ = indicatorFunction_material_neukamm<Domain>;
          elasticityTensor_ = material_neukamm;
          prestrain_ = prestrain_material_neukamm;
        }
        else if(name == "two_phase_material_1")
        {
          indicatorFunction_ = indicatorFunction_two_phase_material_1<Domain>;
          elasticityTensor_ = two_phase_material_1;
          prestrain_ = prestrain_two_phase_material_1;
        }
        else if(name == "two_phase_material_2")
        {
          indicatorFunction_ = indicatorFunction_two_phase_material_2<Domain>;
          elasticityTensor_ = two_phase_material_2;
          prestrain_ = prestrain_two_phase_material_2;
        }
        else if(name == "homogeneous")
        {
          indicatorFunction_ = indicatorFunction_material_homogeneous<Domain>;
          elasticityTensor_ = material_homogeneous;
          prestrain_ = prestrain_homogeneous;
        }
        else
         DUNE_THROW(Exception, "There exists no material with this name "); 
    }


     //--- apply elasticityTensor_ to input Matrix G at position x
 // input: Matrix G, material-phase
  MatrixRT ElasticityTensor(const MatrixRT& G, const int& phase) const  
  {
    return elasticityTensor_(G,phase);
  
  }

  int IndicatorFunction(const Domain& x) const
  {
    return indicatorFunction_(x);
  }

  // static int indicatorFunction_material_1(const Domain& x)
  // {
  //   // std::cout << "x[0] , x[1] , x[2]" << x[0] << " , " << x[1] << ", "<< x[2] << std::endl;
  //   double theta=0.25;
  //   std::cout << "-1/2" << -1/2 << std::endl;
  //   if ((x[1]< (-0.5+theta)) and (x[2]>(0.5-theta)))
  //   {
  //       return 2;    //#Phase1
  //   }
  //   else if ((x[0] < (-0.5+theta)) and (x[2]<(-0.5+theta)))
  //   {
  //       return 1;    //#Phase2
  //   }
  //   else 
  //   {
  //       return 0;    //#Phase3
  //   }
  // }



  // MatrixRT localElasticityTensor(const MatrixRT& G, const Domain& x, Element& element) const  //local Version ... muss binding öfters durchführen
  // {
  //   localIndicatorFunction_.bind(*element);
  //   return elasticityTensor_(G,localIndicatorFunction_(x));
  // }




  // MatrixRT ElasticityTensorGlobal(const MatrixRT& G, const Domain& x) const  //global Version (takes global coordinates)
  // {
  //   return MAT(G,x);
  // }


  // //--- apply elasticityTensor_ to input Matrix G at position x
  // MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const  //global Version (takes global coordinates)
  // {
  //   // Python::Module module = Python::import("material");
  //   // auto indicatorFunctionTest = Python::make_function<double>(module.get("indicatorFunction"));
  //   // return elasticityTensor_(G,indicatorFunctionTest(x));
  //   // auto IFunction = localFunction(indicatorFunction_);
  //   // auto IFunction = this->getIndicatorFunction();
  //   // auto tmp = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
  //   // auto tmpLocal = LocalF
  //   // return elasticityTensor_(G,IFunction(x));
  //   // return elasticityTensor_(G,localIndicatorFunction_(x));
  //   return elasticityTensor_(G,indicatorFunction_(x));
  //   // return elasticityTensor_(G,indicatorFunctionGVF_(x));
  //   // int phase = indicatorFunction_(x);
  //   // auto tmp = elasticityTensor_(G,phase);
  //   // auto tmp = material_1(G,indicatorFunction_(x));
  //   // printmatrix(std::cout, material_1(G,indicatorFunction_(x)), "material_1(G,indicatorFunction_(x))", "--");
  //   // printmatrix(std::cout, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--");
  //   // return tmp;
  //   // return material_1(G,indicatorFunction_(x));
  // }



//////////////////////////////////////////////////////////////////////////////
  // MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
  // {
  //   //--- apply elasticityTensor_ to input Matrix G at position x
  //   return elasticityTensor_(G,x);

  // }

  // MatrixRT applyElasticityTensorLocal(const MatrixRT& G, const Domain& x) const
  // {
  //   //--- apply elasticityTensor_ to input Matrix G at position x (local coordinates)
  //     // MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};

  //   	if (indicatorFunction_(x) == 1) 
	// 		  return L1_(G);
	// 	  else if (indicatorFunction_(x) == 2) 
	// 		  return L2_(G);
  //     else
  //       return L3_(G);

  // }


  // -----------------------------------------------------------------
  // --- write material (grid)functions to VTK
  void write_materialFunctions()
  {


	return;

  };
    ///////////////////////////////
    // Getter
    ///////////////////////////////
    ParameterTree getParameterSet() const {return parameterSet_;}
    // Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
    // Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
    Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;}
    MatrixPhase getPrestrainFunction() const {return prestrain_;} 

    auto getLocalIndicatorFunction() const {return localIndicatorFunction_;}   //get as localFunction
    auto getIndicatorFunctionGVF() const {return indicatorFunctionGVF_;}   //get as GridViewFunction
    auto getIndicatorFunction() const {return indicatorFunction_;}


    // shared_ptr<Func2int> getIndicatorFunction(){return make_shared<Func2int>(indicatorFunction_);}
    // auto getIndicatorFunction(){return make_shared<Func2int>(indicatorFunction_);}
    // shared_ptr<Func2TensorParam> getElasticityTensor(){return make_shared<Func2TensorParam>(elasticityTensor_);}

    //getLameParameters() .. Throw Exception if isotropic_ = false! 




}; // end class


#endif