#ifndef DUNE_MICROSTRUCTURE_PRESTRAINIMP_HH
#define DUNE_MICROSTRUCTURE_PRESTRAINIMP_HH

// #include <dune/microstructure/microstructuredrodgrid.hh>

#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>

//template<class MicrostructuredRodGrid> reicht include

using namespace Dune;


class PrestrainImp {

public:

    static const int dim = 3;
	static const int nCompo = 3;

// 	using GridType_Ce = typename MicrostructuredRodGrid::CellGridType;
    using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
    using GridView = CellGridType::LeafGridView;
    using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;

    using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
    using Func2Tensor = std::function< MatrixRT(const Domain&) >;
// 	using Domain = typename GridType_Ce::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
// 	using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
// 	using Func2Tensor = std::function< MatrixRT(const Domain&) >;

protected:

	double p1, p2, theta;
	double width; //cell geometry

public:

    PrestrainImp(double p1, double p2, double theta, double width) : p1(p1), p2(p2), theta(theta), width(width){}

    Func2Tensor getPrestrain(unsigned int imp)
    {
    	using std::pow;
	    using std::abs;
	    using std::sqrt;
	    using std::sin;
	    using std::cos;

    	if (imp==1)
    	{

            Func2Tensor B1_ = [this] (const Domain& x) {              //  ISOTROPIC PRESSURE
                if (abs(x[0]) > (theta/2)  && x[2] > 0)
                    return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
                if (abs(x[0]) < (theta/2)  && x[2] < 0)
                    return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
                else
                    return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};

            };
            std::cout <<" Prestrain Type: 1 "<< std::endl;
            return B1_;
        }
        else if (imp==2)
        {
            Func2Tensor B2_ = [this] (const Domain& x) {              // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE    

                if (abs(x[0]) < (theta/2) && x[2] < 0 && x[2] > -(1.0/2.0) )
                    return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
                else
                    return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};


            };
            std::cout <<" Prestrain Type: 2 "<< std::endl;
            return B2_;


        }




            // TODO ANISOTROPIC PRESSURE

//         else if (imp==2)
//     	{
//
// 	    	Func2Tensor B2_ = [this] (const Domain& z)
// 	    	{
// 	    		auto pi = std::acos(-1.0);
// 	      		double beta = pi/4.0 + pi/4.0*z[1]; //z[1]=x3
// 	      		MatrixRT Id(0);
// 	      		for (int i=0;i<nCompo;i++)
// 	        		Id[i][i]=1.0;
// 	      		MatrixRT pressure = Id;
// 	      		pressure *= 1.0/6.0;
// 	      		MatrixRT n_ot_n = {
// 	      			{cos(beta)*cos(beta), 0.0, cos(beta)*sin(beta)},
// 	      			{0.0,                           0.0, 0.0						  },
// 	      			{cos(beta)*sin(beta), 0.0, sin(beta)*sin(beta)}
// 	      		};
// 	      		n_ot_n*=0.5;
// 	      		pressure += n_ot_n;
// 	      		pressure *= this->a;
// 	      		return pressure;
// 	    	};
//     		return B2_;
//     	}
    }

};


#endif