#ifndef DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
#define DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH


#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/microstructure/matrix_operations.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>

#include <dune/common/parametertree.hh>

#include <dune/fufem/dunepython.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;


// constexpr unsigned int str2int(const std::string* str, int h = 0)
// {
//     return !str[h] ? 5381 : (str2int(str, h+1) * 33) ^ str[h];
// }

  template<class Domain>
  int indicatorFunction_material_1(const Domain& x)
  {
    double theta=0.25;
    if (x[0] < (-1/2+theta) && x[2]<(-1/2+theta))
        return 1;    //#Phase1
    else if (x[1]< (-1/2+theta) && x[2]>(1/2-theta))
        return 2;    //#Phase2
    else 
        return 0;    //#Phase3
  }
 MatrixRT material_1(const MatrixRT& G, const int& phase)
  {
      const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
      const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
      if (phase == 1)
        return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
      if (phase == 2)
        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 CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
  // using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
  using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
  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&) >;


protected:

  const GridView& gridView_;
  const ParameterTree& parameterSet_;

  int Phases_;

  MatrixFunc L1_;
  MatrixFunc L2_;
  MatrixFunc L3_;



  // const FieldVector<double , ...number of mu-Values/Phases> .. schwierig zur compile-time

  // const FuncScalar mu_; 
  // const FuncScalar lambda_; 
  // double gamma_;

  std::string materialFunctionName_;
  // std::string materialFunctionName_ = parameterSet_.get<std::string>("materialFunction", "material");

  // --- Number of material phases?
  // const int phases_;

  // Func2Tensor materialFunction_;   //actually not needed?? 

  // Func2Tensor elasticityTensor_;
  // Func2TensorParam elasticityTensor_;
  Func2TensorPhase elasticityTensor_;

  // FuncScalar indicatorFunction_;

  // GridViewFunction<double(const Domain&), GridView> indicatorFunction_;
  GridViewFunction<int(const Domain&), GridView> indicatorFunctionGVF_;
  Func2int indicatorFunction_;
  // static const auto indicatorFunction_;

//   VectorCT x_1_, x_2_, x_3_;            // (all) Corrector coefficient vectors 
//   VectorCT phi_1_, phi_2_, phi_3_;      // Corrector phi_i coefficient vectors 
//   FieldVector<double,3> m_1_, m_2_, m_3_;  // Corrector m_i coefficient vectors 

//   MatrixRT M1_, M2_, M3_;  // (assembled) corrector matrices M_i
//   const std::array<MatrixRT*, 3 > mContainer = {&M1_ , &M2_, &M3_};
//   const std::array<VectorCT, 3> phiContainer = {phi_1_,phi_2_,phi_3_};

  // ---- Basis for R_sym^{2x2}
  MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
  MatrixRT G2_ {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
  MatrixRT G3_ {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
  std::array<MatrixRT,3 > MatrixBasisContainer_ = {G1_, G2_, G3_};

  Func2Tensor x3G_1_ = [] (const Domain& x) {
                            return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};    //TODO könnte hier sign übergeben?
                        };

  Func2Tensor x3G_2_ = [] (const Domain& x) {
                            return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*x[2], 0.0}, {0.0, 0.0, 0.0}};
                        };

  Func2Tensor x3G_3_ = [] (const Domain& x) {                                                                               
                            return MatrixRT{{0.0, (1.0/sqrt(2.0))*x[2], 0.0}, {(1.0/sqrt(2.0))*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
                        };

  const std::array<Func2Tensor, 3> x3MatrixBasisContainer_ = {x3G_1_, x3G_2_, x3G_3_};

  
public: 
    ///////////////////////////////
    // constructor
    ///////////////////////////////
    prestrainedMaterial(const GridView gridView,
                        const ParameterTree& parameterSet)    // string: "name of material"? // mu_(mu), muValues? müsste Anzahl Phasen bereits kennen..
                       : gridView_(gridView), 
                         parameterSet_(parameterSet)
    {
	  std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material");


    std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;

    // Python::Module module = Python::import(materialFunctionName_);

    // elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));

    // elasticityTensor_ = setupElasticityTensor(materialFunctionName_);
    setupElasticityTensor(materialFunctionName_);
    // elasticityTensor_ = material_1;
    // elasticityTensor_ = setupElasticityTensor();

    // module.get("Phases").toC<int>(Phases_);

    // auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction"));
    // indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
    // indicatorFunction_ =  Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);


    // L1_ = Python::make_function<MatrixRT>(module.get("L1"));
    // L2_ = Python::make_function<MatrixRT>(module.get("L2"));
    // L3_ = Python::make_function<MatrixRT>(module.get("L3"));


    // Func2TensorParam elasticityTensor_ = Python::make_function<double>(module.get("L"));
    // Func2Tensor materialFunction_ = Python::make_function<double>(module.get("f"));
    // bool isotropic_ = true; // read from module File TODO 
    } 
/////////////////////////////////////////////////////////////////////////////////////////////////////////



//   static int indicatorFunction_material_1(const Domain& x)
//   {
//     double theta=0.25;
//     if (x[0] < (-1/2+theta) && x[2]<(-1/2+theta))
//         return 1;    //#Phase1
//     else if (x[1]< (-1/2+theta) && x[2]>(1/2-theta))
//         return 2;    //#Phase2
//     else 
//         return 0;    //#Phase3
//   }
//  static MatrixRT material_1(const MatrixRT& G, const int& phase)
//   {
//       FieldVector<double,3> mu = {80.0, 80.0, 60.0};
//       FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
//       if (phase == 1)
//         return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
//       if (phase == 2)
//         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();
//   }


//---function that determines elasticity Tensor 
    // auto setupElasticityTensor(const std::string name)  // cant use materialFunctionName_ here!?
    // {
    //     if(name == "material")
    //     {
    //       return material_1;
    //     }
    //     else
    //      DUNE_THROW(Exception, "There exists no material in materialDefinitions.hh with this name "); 
    // }



  //---function that determines elasticity Tensor 
    void setupElasticityTensor(const std::string name)  // cant use materialFunctionName_ here!?
    {
        if(name == "material")
        {
          // indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
          // indicatorFunction_ = indicatorFunction_material_1;
          indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1<Domain>, gridView_);
          indicatorFunction_ = indicatorFunction_material_1<Domain>;
          elasticityTensor_ = material_1;
        }
        else
         DUNE_THROW(Exception, "There exists no material in materialDefinitions.hh with this name "); 
    }


  //--- apply elasticityTensor_ to input Matrix G at position x
  MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
  {
    // 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,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_;}



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


    // shared_ptr<Func2TensorParam> getElasticityTensor(){return make_shared<Func2TensorParam>(elasticityTensor_);}


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




    // shared_ptr<MatrixCT> getStiffnessMatrix(){return make_shared<MatrixCT>(stiffnessMatrix_);}
    // shared_ptr<VectorCT> getLoad_alpha1(){return make_shared<VectorCT>(load_alpha1_);}
    // shared_ptr<VectorCT> getLoad_alpha2(){return make_shared<VectorCT>(load_alpha2_);}
    // shared_ptr<VectorCT> getLoad_alpha3(){return make_shared<VectorCT>(load_alpha3_);}

    // shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
    // shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}


    // --- Get Correctors
    // shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
    // auto getMcontainer(){return make_shared<std::array<MatrixRT*, 3 > >(mContainer);}
    // auto getMcontainer(){return mContainer;}
    // shared_ptr<std::array<VectorCT, 3>> getPhicontainer(){return make_shared<std::array<VectorCT, 3>>(phiContainer);}

    
    // // shared_ptr<std::array<VectorRT, 3>> getBasiscontainer(){return make_shared<std::array<VectorRT, 3>>(basisContainer_);}
    // auto getMatrixBasiscontainer(){return make_shared<std::array<MatrixRT,3 >>(MatrixBasisContainer_);}
    // // auto getx3MatrixBasiscontainer(){return make_shared<std::array<Func2Tensor, 3>>(x3MatrixBasisContainer_);}
    // auto getx3MatrixBasiscontainer(){return x3MatrixBasisContainer_;}


  

    
    // shared_ptr<VectorCT> getCorr_a(){return make_shared<VectorCT>(x_a_);}
    // shared_ptr<VectorCT> getCorr_phi1(){return make_shared<VectorCT>(phi_1_);}
    // shared_ptr<VectorCT> getCorr_phi2(){return make_shared<VectorCT>(phi_2_);}
    // shared_ptr<VectorCT> getCorr_phi3(){return make_shared<VectorCT>(phi_3_);}






}; // end class




#endif