#ifndef DUNE_MICROSTRUCTURE_MATERIALDEFINITIONS.HH
#define DUNE_MICROSTRUCTURE_MATERIALDEFINITIONS.HH

#include <dune/common/parametertree.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/microstructure/matrix_operations.hh>



using namespace Dune;
using namespace MatrixOperations;
using std::pow;
using std::abs;
using std::sqrt;
using std::sin;
using std::cos;
using MatrixRT = FieldMatrix< double, 3, 3>;
using VectorRT = FieldVector< double, 3>;




// (30.8.22) DEPRECATED CLASS! 



  ////----------------------------------------------------------------------------------
//   template<class Domain>
//   int indicatorFunction_material_new(const Domain& x)
//   {
//     double theta=0.25;
//     if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
//         return 1;    //#Phase1
//     else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
//         return 2;    //#Phase2
//     else 
//         return 3;    //#Phase3
//   }
//  MatrixRT material_new(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 VoigtToMatrix(isotropic(mu[0],lambda[0])* matrixToVoigt(G));
//         // return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id(); //#Phase1 //Isotrop(mu,lambda)
//       else if (phase == 2)
//         return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id(); //#Phase2
//       else 
//         return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id(); //#Phase3
//   }
//   MatrixRT prestrain_material_new(const int& phase)
//   {
//       if (phase == 1)
//         return {{1.0, 0.0, 0.0},    //#Phase1
//                 {0.0, 1.0, 0.0},
//                 {0.0, 0.0, 1.0}
//                };
//       else if (phase == 2)
//         return {{1.0, 0.0, 0.0},    //#Phase2
//                 {0.0, 1.0, 0.0},
//                 {0.0, 0.0, 1.0}
//                };
//       else 
//         return {{0.0, 0.0, 0.0},    //#Phase3
//                 {0.0, 0.0, 0.0},
//                 {0.0, 0.0, 0.0}
//                };
//   }
  ////----------------------------------------------------------------------------------





  // ----------------------------------------------------------------------------------
  template<class Domain>
  int indicatorFunction_material_neukamm(const Domain& x)
  {
    double theta=0.25;
    if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
        return 1;    //#Phase1
    else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
        return 2;    //#Phase2
    else 
        return 3;    //#Phase3
  }
 MatrixRT material_neukamm(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(); //#Phase1 //Isotrop(mu,lambda)
      else if (phase == 2)
        return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id(); //#Phase2
      else 
        return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id(); //#Phase3
  }
  MatrixRT prestrain_material_neukamm(const int& phase)
  {
      if (phase == 1)
        return {{1.0, 0.0, 0.0},    //#Phase1
                {0.0, 1.0, 0.0},
                {0.0, 0.0, 1.0}
               };
      else if (phase == 2)
        return {{1.0, 0.0, 0.0},    //#Phase2
                {0.0, 1.0, 0.0},
                {0.0, 0.0, 1.0}
               };
      else 
        return {{0.0, 0.0, 0.0},    //#Phase3
                {0.0, 0.0, 0.0},
                {0.0, 0.0, 0.0}
               };
  }
  // ----------------------------------------------------------------------------------


  // ----------------------------------------------------------------------------------
  template<class Domain>
  int indicatorFunction_two_phase_material_1(const Domain& x)
  {
    double theta=0.25;
    if(abs(x[0]) > (theta/2.0))  
        return 1;    //#Phase1
    else 
        return 0;    //#Phase2
  }
 MatrixRT two_phase_material_1(const MatrixRT& G, const int& phase)
  {
      const FieldVector<double,2> mu     = {80.0, 160.0};
      const FieldVector<double,2> lambda = {80.0, 160.0};
      if (phase == 1)
        return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
      else 
        return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
  }
  MatrixRT prestrain_two_phase_material_1(const int& phase)
  {
      const FieldVector<double,2> rho = {3.0, 5.0};
      if (phase == 1)
        return {{rho[0], 0.0, 0.0},    //#Phase1
                {0.0, rho[0], 0.0},
                {0.0, 0.0, rho[0]}
               };
      else 
        return {{rho[1], 0.0, 0.0},    //#Phase2
                {0.0, rho[1], 0.0},
                {0.0, 0.0, rho[1]}
               };
  }
  // ----------------------------------------------------------------------------------


  // ----------------------------------------------------------------------------------
  template<class Domain>
  int indicatorFunction_two_phase_material_2(const Domain& x)
  {
    double theta=0.25;
    if(abs(x[0]) > (theta/2.0))  
        return 1;    //#Phase1
    else 
        return 0;    //#Phase2
  }
 MatrixRT two_phase_material_2(const MatrixRT& G, const int& phase)
  {
      const FieldVector<double,2> mu     = {5.0, 15.0};
      const FieldVector<double,2> lambda = {6.0, 8.0};
      if (phase == 1)
        return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
      else 
        return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
  }
  MatrixRT prestrain_two_phase_material_2(const int& phase)
  {
      const FieldVector<double,2> rho = {3.0, 5.0};
      if (phase == 1)
        return {{rho[0], 0.0, 0.0},    //#Phase1
                {0.0, rho[0], 0.0},
                {0.0, 0.0, rho[0]}
               };
      else 
        return {{rho[1], 0.0, 0.0},    //#Phase2
                {0.0, rho[1], 0.0},
                {0.0, 0.0, rho[1]}
               };
  }
  // ----------------------------------------------------------------------------------



  // ----------------------------------------------------------------------------------
  template<class Domain>
  int indicatorFunction_material_2(const Domain& x)
  {
    double theta=0.25;
    if(abs(x[0]) > (theta/2.0))  
        return 1;    //#Phase1
    else 
        return 0;    //#Phase2
  }
 MatrixRT material_2(const MatrixRT& G, const int& phase)
  {
      const FieldVector<double,2> mu     = {80.0, 160.0};
      const FieldVector<double,2> lambda = {80.0, 160.0};
      if (phase == 1)
        return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id(); //#Phase1
      else 
        return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id(); //#Phase2
  }
  MatrixRT prestrain_material_2(const int& phase)
  {
      const FieldVector<double,2> rho = {3.0, 5.0};
      if (phase == 1)
        return {{rho[0], 0.0, 0.0},    //#Phase1
                {0.0, rho[0], 0.0},
                {0.0, 0.0, rho[0]}
               };
      else 
        return {{rho[1], 0.0, 0.0},    //#Phase2
                {0.0, rho[1], 0.0},
                {0.0, 0.0, rho[1]}
               };
  }
  // ----------------------------------------------------------------------------------


  // ----------------------------------------------------------------------------------
  template<class Domain>
  int indicatorFunction_material_homogeneous(const Domain& x)
  {
      return 0;  
  }
 MatrixRT material_homogeneous(const MatrixRT& G, const int& phase)
 {
      const FieldVector<double,1> mu     = {80.0};
      const FieldVector<double,1> lambda = {80.0};
      return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
 }
   MatrixRT prestrain_homogeneous(const int& phase)
  {
      if (phase == 1)
        return {{1.0, 0.0, 0.0},    //#Phase1
                {0.0, 1.0, 0.0},
                {0.0, 0.0, 1.0}
               };
      else if (phase == 2)
        return {{1.0, 0.0, 0.0},    //#Phase2
                {0.0, 1.0, 0.0},
                {0.0, 0.0, 1.0}
               };
      else 
        return {{1.0, 0.0, 0.0},    //#Phase3
                {0.0, 1.0, 0.0},
                {0.0, 0.0, 1.0}
               };
  }
 // ----------------------------------------------------------------------------------
























#endif