diff --git a/dune/microstructure/prestrain_material_geometry.hh b/dune/microstructure/prestrain_material_geometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c52107ca5bf2ef401c22d91cf9c10833017dacba
--- /dev/null
+++ b/dune/microstructure/prestrain_material_geometry.hh
@@ -0,0 +1,645 @@
+#ifndef DUNE_MICROSTRUCTURE_PRESTRAI_MATERIAL_GEOMETRY_HH
+#define DUNE_MICROSTRUCTURE_PRESTRAI_MATERIAL_GEOMETRY_HH
+
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.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;
+
+
+
+
+
+
+class IsotropicMaterialImp
+{
+
+public:
+ static const int dim = 3;
+ static const int dimWorld = 3;
+ using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+ using GridView = CellGridType::LeafGridView;
+ using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+ using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+ using FuncScalar = std::function< double(const Domain&) >;
+
+
+ FuncScalar getMu(ParameterTree parameters){
+
+ std::string imp = parameters.get<std::string>("material_implementation");
+ //std::array<std::string,5> imps({"homogen_poisson" "bilayer_poisson" "chess_poisson" "3Dchess_poisson" "vertical_bilayer_poisson"});
+
+ double phi = M_PI*parameters.get<double>("material_angle", 0.0)/180; //TODO
+ double theta = parameters.get<double>("theta",1.0/4.0); //TODO alle parameter hier laden?
+
+ double width = parameters.get<double>("width", 1.0);
+ double height = parameters.get<double>("height", 1.0);
+
+
+ if (imp == "homogeneous"){
+ double mu = parameters.get<double>("mu", 10);
+
+ auto muTerm = [mu] (const Domain& z) {
+ return mu;};
+
+ return muTerm;
+ }
+
+
+ else if (imp == "analytical_Example"){
+ double mu1 = parameters.get<double>("mu1",10.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double mu2 = beta*mu1;
+
+ auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+
+// std::cout << "Analytical-MU is used" << std::endl;
+ if (abs(z[0]) >= (theta/2.0))
+ return mu1;
+ else
+ return mu2;
+ };
+
+ return muTerm;
+ }
+ else if (imp == "matrix_material") // Matrix material with prestrained Fiber inclusions
+ {
+ double mu1 = parameters.get<double>("mu1",10.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double mu2 = beta*mu1;
+
+ int nF = parameters.get<int>("nF", 3); //number of Fibers in each Layer
+ double rF = parameters.get<double>("rF", 0.1);
+
+ assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+
+ auto muTerm = [mu1,mu2,theta,nF,rF,height,width] (const Domain& x)
+ {
+ //TODO if Domain == 1... if Domain == 2 ...
+ double domainShift = 1.0/2.0;
+ // shift x to domain [0,1]^3
+ FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+// std::cout << "matrixMaterial-MU is used" << std::endl;
+
+ // determine if point is in upper/lower Layer
+ if ((height/2) <= s[2]<= height) // upper Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ {
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (width/(2*nF))+((width/nF)*i) , (3/4)*height };
+//
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ return mu2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ return mu1;
+ }
+ }
+ }
+ else // lower Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ {
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (width/(2*nF))+((width/nF)*i) , height/4.0 };
+
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ return mu2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ return mu1;
+ }
+ }
+
+ }
+
+ };
+ return muTerm;
+ }
+ else if (imp == "bilayer"){
+ double mu1 = parameters.get<double>("mu1",10.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double mu2 = beta*mu1;
+
+ auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+ if ( isInRotatedPlane(phi, z[dim-2], z[dim-1]) )
+ return mu1;
+ else
+ return mu2; };
+
+ return muTerm;
+ }
+
+
+
+ else if (imp == "helix_poisson"){ // interessant!
+ double mu1 = parameters.get<double>("mu1",10.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double mu2 = beta*mu1;
+
+ double r = parameters.get<double>("radius");
+ double r2 = parameters.get<double>("radius_helix");
+
+ auto muTerm = [mu1, mu2, r, r2] (const Domain& z) {
+ std::array<double,2> h = {r*cos(2*M_PI*z[0]), r*sin(2*M_PI*z[0])};
+ if ( sqrt(pow(z[1]-h[0],2)+pow(z[2]-h[1],2)) < r2 ) // distance to the helix?
+ return mu1;
+ else
+ return mu2;
+ };
+
+ return muTerm;
+ }
+
+ else if (imp == "chess_poisson"){
+ double mu1 = parameters.get<double>("mu1",10.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double mu2 = beta*mu1;
+
+ auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+ if ( (isInRotatedPlane(phi, z[dim-2], z[dim-1]) and isInRotatedPlane(phi + M_PI/2 , z[dim-2], z[dim-1])) or
+ (isInRotatedPlane(phi + M_PI, z[dim-2], z[dim-1]) and isInRotatedPlane(phi + 3*M_PI/2, z[dim-2], z[dim-1])) )
+ return mu1;
+ else
+ return mu2;
+ };
+
+ return muTerm;
+ }
+
+// else if (imp == "3Dchess_poisson"){
+// double E1 = parameters.get<double>("E1", 17e6); //Faser
+// double nu1 = parameters.get<double>("nu1", 0.3);
+// double mu1 = lameMu(E1, nu1);
+//
+// double E2 = parameters.get<double>("E2", 17e6); //Polymer, weicher, Querkontraktion groß
+// double nu2 = parameters.get<double>("nu2", 0.5);
+// double mu2 = lameMu(E2, nu2);
+//
+// auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+// if (z[0]<0.5)
+// if ( (isInRotatedPlane(phi, z[1], z[2]) and isInRotatedPlane(phi + M_PI/2 , z[1], z[2])) or
+// (isInRotatedPlane(phi + M_PI, z[1], z[2]) and isInRotatedPlane(phi + 3*M_PI/2, z[1], z[2])) )
+// return mu1;
+// else
+// return mu2;
+// else
+// if ( (isInRotatedPlane(phi, z[1], z[2]) and isInRotatedPlane(phi + M_PI/2 , z[1], z[2])) or
+// (isInRotatedPlane(phi + M_PI, z[1], z[2]) and isInRotatedPlane(phi + 3*M_PI/2, z[1], z[2])) )
+// return mu2;
+// else
+// return mu1;
+// };
+//
+// return muTerm;
+// }
+// else if (imp == "vertical_bilayer_poisson"){
+// double E1 = parameters.get<double>("E1", 17e6); //material1
+// double nu1 = parameters.get<double>("nu1", 0.3);
+// double mu1 = lameMu(E1, nu1);
+// double E2 = parameters.get<double>("E2", 17e6); //material2
+// double nu2 = parameters.get<double>("nu2", 0.5);
+// double mu2 = lameMu(E2, nu2);
+//
+// auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+// if ( isInRotatedPlane(phi+M_PI/2.0, z[dim-2], z[dim-1]) )//x3
+// return mu1;
+// else
+// return mu2; };
+//
+// return muTerm;
+// }
+//
+// else if (imp == "microstructured_bilayer_poisson"){
+// double E1 = parameters.get<double>("E1", 17e6); //material1
+// double nu1 = parameters.get<double>("nu1", 0.3);
+// double mu1 = lameMu(E1, nu1);
+// double E2 = parameters.get<double>("E2", 17e6); //material2
+// double nu2 = parameters.get<double>("nu2", 0.5);
+// double mu2 = lameMu(E2, nu2);
+//
+// auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+// if ( isInRotatedPlane(phi, z[0]-0.5, z[1]) ) //TODO understand this
+// return mu1;
+// else
+// return mu2; };
+//
+// return muTerm;
+// }
+//
+// else { //(imp == "bilayer_poisson")
+// double E1 = parameters.get<double>("E1", 17e6); //material1
+// double nu1 = parameters.get<double>("nu1", 0.3);
+// double mu1 = lameMu(E1, nu1);
+// double E2 = parameters.get<double>("E2", 17e6); //material2
+// double nu2 = parameters.get<double>("nu2", 0.5);
+// double mu2 = lameMu(E2, nu2);
+//
+// auto muTerm = [mu1, mu2, phi] (const Domain& z) {
+// if ( isInRotatedPlane(phi, z[dim-2], z[dim-1]) ) //x2
+// return mu1;
+// else
+// return mu2; };
+//
+// return muTerm;
+// }
+ }
+
+ FuncScalar getLambda(ParameterTree parameters)
+ {
+
+ std::string imp = parameters.get<std::string>("material_implementation");
+ //std::array<std::string,5> imps({"homogen_poisson" "bilayer_poisson" "chess_poisson" "3Dchess_poisson" "vertical_bilayer_poisson"});
+ //int i_imp = parameters.get<int>("impnumber");
+ //std::string imp = imps[i_imp];
+ double phi = M_PI*parameters.get<double>("material_angle", 0.0)/180;
+ double theta = parameters.get<double>("theta",1.0/4.0); //TODO alle parameter hier laden?
+
+ double width = parameters.get<double>("width", 1.0);
+ double height = parameters.get<double>("height", 1.0);
+
+ if (imp == "homogenenous"){
+ double lambda = parameters.get<double>("lambda",384.615);
+
+ auto lambdaTerm = [lambda] (const Domain& z) {return lambda;};
+
+ return lambdaTerm;
+ }
+
+
+ else if (imp == "analytical_Example"){
+ double lambda1 = parameters.get<double>("lambda1",0.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double lambda2 = beta*lambda1;
+
+ auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
+
+// std::cout << "Analytical-LAMBDA is used" << std::endl; //TEST
+ if (abs(z[0]) >= (theta/2.0))
+ return lambda1;
+ else
+ return lambda2;
+ };
+
+ return lambdaTerm;
+ }
+ else if (imp == "matrix_material") // Matrix material with prestrained Fiber inclusions
+ {
+ double lambda1 = parameters.get<double>("lambda1",0.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double lambda2 = beta*lambda1;
+
+ int nF = parameters.get<int>("nF", 3); //number of Fibers in each Layer
+ double rF = parameters.get<double>("rF", 0.1);
+
+ assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+
+ auto lambdaTerm = [lambda1,lambda2, theta,nF,rF,height,width] (const Domain& x)
+ {
+ //TODO if Domain == 1... if Domain == 2 ...
+ double domainShift = 1.0/2.0;
+ // shift x to domain [0,1]^3
+ FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+// std::cout << "matrixMaterial-LAMBDA is used" << std::endl;
+
+ // determine if point is in upper/lower Layer
+ if ((height/2) <= s[2]<= height) // upper Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ {
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (width/(2*nF))+((width/nF)*i) , (3/4)*height };
+//
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ return lambda2;
+ else
+ return lambda1;
+ }
+ }
+ }
+ else // lower Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ {
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (width/(2*nF))+((width/nF)*i) , height/4.0 };
+
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ return lambda2;
+ else
+ return lambda1;
+ }
+ }
+
+ }
+
+ };
+ return lambdaTerm;
+ }
+
+
+ else if (imp == "bilayer_lame"){
+ double lambda1 = parameters.get<double>("mu1",384.615);
+ double lambda2 = parameters.get<double>("mu2",384.615);
+
+ auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+ if ( isInRotatedPlane(phi, z[dim-2], z[dim-1]) )
+ return lambda1;
+ else
+ return lambda2; };
+
+ return lambdaTerm;
+ }
+
+ else if (imp == "helix_poisson"){
+ double E1 = parameters.get<double>("E1", 17e6); //Faser
+ double nu1 = parameters.get<double>("nu1", 0.3);
+ double lambda1 = lameLambda(E1,nu1);
+
+ double E2 = parameters.get<double>("E2", 17e6); //Polymer, weicher, Querkontraktion groß
+ double nu2 = parameters.get<double>("nu2", 0.5);
+ double lambda2 = lameLambda(E2, nu2);
+
+ double r = parameters.get<double>("radius");
+ double r2 = parameters.get<double>("radius_helix");
+
+ auto lambdaTerm = [lambda1, lambda2, r, r2] (const Domain& z) {
+ std::array<double,2> h = {r*cos(2*M_PI*z[0]), r*sin(2*M_PI*z[0])};
+ MatrixRT ret(0.0);
+ if ( sqrt(pow(z[1]-h[0],2)+pow(z[2]-h[1],2)) < r2 )
+ return lambda1;
+ else
+ return lambda2;
+ };
+
+ return lambdaTerm;
+ }
+
+ else if (imp == "chess_poisson"){
+ double E1 = parameters.get<double>("E1", 17e6); //Faser
+ double nu1 = parameters.get<double>("nu1", 0.3);
+ double lambda1 = lameLambda(E1,nu1);
+
+ double E2 = parameters.get<double>("E2", 17e6); //Polymer, weicher, Querkontraktion groß
+ double nu2 = parameters.get<double>("nu2", 0.5);
+ double lambda2 = lameLambda(E2, nu2);
+
+ auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+ if ( (isInRotatedPlane(phi, z[dim-2], z[dim-1]) and isInRotatedPlane(phi + M_PI/2 , z[dim-2], z[dim-1])) or
+ (isInRotatedPlane(phi + M_PI, z[dim-2], z[dim-1]) and isInRotatedPlane(phi + 3*M_PI/2, z[dim-2], z[dim-1])) )
+ return lambda1;
+ else
+ return lambda2;
+ };
+
+ return lambdaTerm;
+ }
+
+// else if (imp == "3Dchess_poisson"){
+// double E1 = parameters.get<double>("E1", 17e6); //Faser
+// double nu1 = parameters.get<double>("nu1", 0.3);
+// double lambda1 = lameLambda(E1,nu1);
+//
+// double E2 = parameters.get<double>("E2", 17e6); //Polymer, weicher, Querkontraktion groß
+// double nu2 = parameters.get<double>("nu2", 0.5);
+// double lambda2 = lameLambda(E2, nu2);
+//
+// auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+// if (z[0]<0.5)
+// if ( (isInRotatedPlane(phi, z[1], z[2]) and isInRotatedPlane(phi + M_PI/2 , z[1], z[2])) or
+// (isInRotatedPlane(phi + M_PI, z[1], z[2]) and isInRotatedPlane(phi + 3*M_PI/2, z[1], z[2])) )
+// return lambda1;
+// else
+// return lambda2;
+// else
+// if ( (isInRotatedPlane(phi, z[1], z[2]) and isInRotatedPlane(phi + M_PI/2 , z[1], z[2])) or
+// (isInRotatedPlane(phi + M_PI, z[1], z[2]) and isInRotatedPlane(phi + 3*M_PI/2, z[1], z[2])) )
+// return lambda2;
+// else
+// return lambda1;
+// };
+//
+// return lambdaTerm;
+// }
+//
+// else if (imp == "vertical_bilayer_poisson"){
+// double E1 = parameters.get<double>("E1", 17e6); //material1
+// double nu1 = parameters.get<double>("nu1", 0.3);
+// double lambda1 = lameLambda(E1,nu1);
+// double E2 = parameters.get<double>("E2", 17e6); //material2
+// double nu2 = parameters.get<double>("nu2", 0.5);
+// double lambda2 = lameLambda(E2, nu2);
+//
+// auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+// if ( isInRotatedPlane(phi+M_PI/2.0, z[dim-2], z[dim-1]) )
+// return lambda1;
+// else
+// return lambda2; };
+//
+// return lambdaTerm;
+// }
+//
+// else if (imp == "microstructured_bilayer_poisson"){
+// double E1 = parameters.get<double>("E1", 17e6); //material1
+// double nu1 = parameters.get<double>("nu1", 0.3);
+// double lambda1 = lameLambda(E1,nu1);
+// double E2 = parameters.get<double>("E2", 17e6); //material2
+// double nu2 = parameters.get<double>("nu2", 0.5);
+// double lambda2 = lameLambda(E2, nu2);
+//
+// auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+// if ( isInRotatedPlane(phi, z[0]-0.5, z[1]) )
+// return lambda1;
+// else
+// return lambda2; };
+//
+// return lambdaTerm;
+// }
+//
+// else { //(imp == "bilayer_poisson")
+// double E1 = parameters.get<double>("E1", 17e6); //material1
+// double nu1 = parameters.get<double>("nu1", 0.3);
+// double lambda1 = lameLambda(E1, nu1);
+// double E2 = parameters.get<double>("E2", 17e6); //material2
+// double nu2 = parameters.get<double>("nu2", 0.5);
+// double lambda2 = lameLambda(E2, nu2);
+//
+// auto lambdaTerm = [lambda1, lambda2, phi] (const Domain& z) {
+// if ( isInRotatedPlane(phi, z[dim-2], z[dim-1]) )
+// return lambda1;
+// else
+// return lambda2; };
+//
+// return lambdaTerm;
+// }
+ }
+
+
+
+};
+
+
+
+
+
+
+
+
+
+
+
+
+class PrestrainImp
+{
+
+public:
+ static const int dim = 3;
+ static const int dimWorld = 3;
+ using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+ using GridView = CellGridType::LeafGridView;
+ using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+ 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(std::string imp)
+ Func2Tensor getPrestrain(ParameterTree parameters)
+ {
+
+ std::string imp = parameters.get<std::string>("prestrainType", "analytical_Example");
+
+ double width = parameters.get<double>("width", 1.0);
+ double height = parameters.get<double>("height", 1.0);
+
+ double theta = parameters.get<double>("theta",1.0/4.0);
+ double p1 = parameters.get<double>("rho1", 1.0);
+ double alpha = parameters.get<double>("alpha", 2.0);
+ double p2 = alpha*p1;
+
+ if (imp == "isotropic_bilayer")
+ {
+ Func2Tensor B = [p1,p2,theta] (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: isotropic_bilayer " << std::endl;
+ return B;
+ }
+ else if (imp == "analytical_Example")
+ {
+ Func2Tensor B = [p1,theta] (const Domain& x) // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE
+ {
+// if (abs(x[0]) < (theta/2) && x[2] < 0 ) //does not make a difference
+ 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: analytical_Example "<< std::endl;
+ return B;
+ }
+ else if (imp == "matrix_material") // Matrix material with prestrained Fiber inclusions
+ {
+ int nF = parameters.get<int>("nF", 3); //number of Fibers in each Layer
+ double rF = parameters.get<double>("rF", 0.1);
+
+ assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+
+ Func2Tensor B = [p1,p2,theta,nF,rF,height,width] (const Domain& x)
+ {
+ //TODO if Domain == 1... if Domain == 2 ...
+ double domainShift = 1.0/2.0;
+ // shift x to domain [0,1]^3
+ FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+
+
+
+ // determine if point is in upper/lower Layer
+ if ((height/2) <= s[2]<= height) // upper Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ {
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (width/(2*nF))+((width/nF)*i) , (3/4)*height };
+//
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ 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}};
+// return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+ }
+ }
+ }
+ else // lower Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ {
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (width/(2*nF))+((width/nF)*i) , height/4.0 };
+
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ 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}};
+// return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+ }
+ }
+
+ }
+
+ };
+ std::cout <<" Prestrain Type: matrix_material"<< std::endl;
+ return B;
+ }
+ else
+ {
+ // TODO other geometries etc...
+
+ }
+ // TODO ANISOTROPIC PRESSURE
+ }
+
+};
+
+
+
+
+
+
+
+#endif
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 19b4e6848d73f8d7a11083a2c3d23dc8165cf58e..e0b00e86dce0fdd86d60c92f922e671b59fa2c21 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -28,7 +28,7 @@ cellDomain = 1
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
########################################################################
-numLevels = 5 5 # computes all levels from first to second entry
+numLevels = 1 3 # computes all levels from first to second entry
#numLevels = 3 3 # computes all levels from first to second entry
#numLevels = 1 6
@@ -43,7 +43,6 @@ numLevels = 5 5 # computes all levels from first to second entry
#nElements_Cell = 100 100 2
#nElements_Cell = 100 100 100 // does not work
#nElements_Cell = 10 10 10
-
#nElements_Cell = 2 2 2
#nElements_Cell = 4 4 4
#nElements_Cell = 8 8 8
@@ -51,11 +50,9 @@ numLevels = 5 5 # computes all levels from first to second entry
#nElements_Cell = 32 32 32
#nElements_Cell = 64 64 64
-width = 1.0
-
-gamma = 50.0
-#gamma = 1.0
+#gamma = 50.0
+gamma = 1.0
#############################################
# Material parameters
@@ -65,21 +62,39 @@ beta = 2.0 # ratio between material parameters mu1 & mu2 .... beta = 1.0 corr
mu1 = 10.0
lambda1 = 0.0
+#lambda1 = 20.0
#lambda1 = 5.0
+#material_implementation = "analytical_Example"
+material_implementation = "matrix_material"
+#material_implementation = "isotropic_bilayer"
+
+
#############################################
# Prestrain parameters
#############################################
+rho1 = 1.0
+alpha = 5.0 # ratio between prestrain parameters rho1 & rho2
+
#theta = 0.3 # volume fraction #default = 1.0/4.0
#theta = 0.25 # volume fraction
#theta = 0.75 # volume fraction
-prestrainType = 2
+#prestrainType = "analytical_Example"
+#prestrainType = "isotropic_bilayer"
-rho1 = 1.0
+------------Matrix Material --------------
+
+prestrainType = "matrix_material"
+
+nF = 5 #number of Fibers (in each Layer)
+rF = 0.1 #Fiber radius
+------------------------------------------
+
+width = 1.0
+height = 1.0
-alpha = 2.0 # ratio between prestrain parameters rho1 & rho2
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 0b1a20b04230378933fa4b27236695e215b53912..510c1b949e6055338ba08c0a2f07abaf170da472 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -45,7 +45,7 @@
#include <dune/common/fmatrix.hh>
-#include <dune/microstructure/prestrainimp.hh>
+#include <dune/microstructure/prestrain_material_geometry.hh>
#include <dune/microstructure/matrix_operations.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
@@ -871,13 +871,19 @@ int main(int argc, char *argv[])
///////////////////////////////////
// Get Prestrain/Parameters
///////////////////////////////////
- unsigned int prestraintype = parameterSet.get<unsigned int>("prestrainType", 2);
- double rho1 = parameterSet.get<double>("rho1", 1.0);
- double rho2 = alpha*rho1;
- auto prestrainImp = PrestrainImp(rho1, rho2, theta, width);
- auto B_Term = prestrainImp.getPrestrain(prestraintype);
+// unsigned int prestraintype = parameterSet.get<unsigned int>("prestrainType", "analytical_Example"); //OLD
+// std::string prestraintype = parameterSet.get<std::string>("prestrainType", "analytical_Example");
+// double rho1 = parameterSet.get<double>("rho1", 1.0);
+// double rho2 = alpha*rho1;
+// auto prestrainImp = PrestrainImp(rho1, rho2, theta, width);
+// auto B_Term = prestrainImp.getPrestrain(prestraintype);
+
+
+
+ auto prestrainImp = PrestrainImp(); //NEW
+ auto B_Term = prestrainImp.getPrestrain(parameterSet);
- log << "prestrain imp: " << prestraintype << "\nrho1 = " << rho1 << "\nrho2 = " << rho2 << std::endl;
+// log << "prestrain imp: " << prestraintype << "\nrho1 = " << rho1 << "\nrho2 = " << rho2 << std::endl;
log << "alpha: " << alpha << std::endl;
log << "gamma: " << gamma << std::endl;
log << "theta: " << theta << std::endl;
@@ -945,6 +951,12 @@ int main(int argc, char *argv[])
///////////////////////////////////
// Create Lambda-Functions for material Parameters depending on microstructure
///////////////////////////////////
+
+ auto materialImp = IsotropicMaterialImp();
+ auto muTerm = materialImp.getMu(parameterSet);
+ auto lambdaTerm = materialImp.getLambda(parameterSet);
+
+/*
auto muTerm = [mu1, mu2, theta] (const Domain& z) {
if (abs(z[0]) >= (theta/2.0))
return mu1;
@@ -957,7 +969,7 @@ int main(int argc, char *argv[])
return lambda1;
else
return lambda2;
- };
+ };*/
auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
auto muLocal = localFunction(muGridF);
@@ -1098,7 +1110,7 @@ int main(int argc, char *argv[])
VectorCT x_2 = load_alpha2;
VectorCT x_3 = load_alpha3;
- auto load_alpha1BS = load_alpha1;
+// auto load_alpha1BS = load_alpha1;
// printvector(std::cout, load_alpha1, "load_alpha1 before SOLVER", "--" );
// printvector(std::cout, load_alpha2, "load_alpha2 before SOLVER", "--" );