From 09782f0da1c88f2a25cecabe85ada07eb11e57a7 Mon Sep 17 00:00:00 2001
From: Stefan Neukamm <stefan.neukamm@tu-dresden.de>
Date: Thu, 7 Jul 2022 09:40:41 +0200
Subject: [PATCH] =?UTF-8?q?material=5Fneukamm=20hinzugef=C3=BCgt?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
.../prestrain_material_geometry.hh | 2136 +++++++++--------
geometries/material_neukamm.py | 16 +
inputs/cellsolver.parset | 9 +-
outputs/Results/1/BMatrix.txt | 3 +
outputs/Results/1/QMatrix.txt | 9 +
outputs/Results/1/output.txt | 61 +
outputs/Results/1/parameter | 14 +
7 files changed, 1201 insertions(+), 1047 deletions(-)
create mode 100644 geometries/material_neukamm.py
create mode 100644 outputs/Results/1/BMatrix.txt
create mode 100644 outputs/Results/1/QMatrix.txt
create mode 100644 outputs/Results/1/output.txt
create mode 100644 outputs/Results/1/parameter
diff --git a/dune/microstructure/prestrain_material_geometry.hh b/dune/microstructure/prestrain_material_geometry.hh
index d39bdbe8..e26cf1d4 100644
--- a/dune/microstructure/prestrain_material_geometry.hh
+++ b/dune/microstructure/prestrain_material_geometry.hh
@@ -23,860 +23,893 @@ class IsotropicMaterialImp
{
public:
- static const int dimWorld = dim;
- using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
- using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
- using Func2Tensor = std::function< MatrixRT(const Domain&) >;
- using FuncScalar = std::function< double(const Domain&) >;
+ static const int dimWorld = dim;
+ using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+ using Domain = typename CellGridType::LeafGridView::template 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){
+ FuncScalar getMu(ParameterTree parameters){
- std::string imp = parameters.get<std::string>("material_prestrain_imp", "analytical_Example");
- //std::array<std::string,5> imps({"homogen_poisson" "bilayer_poisson" "chess_poisson" "3Dchess_poisson" "vertical_bilayer_poisson"});
+ std::string imp = parameters.get<std::string>("material_prestrain_imp", "analytical_Example");
+ //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 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);
+ double width = parameters.get<double>("width", 1.0);
+ double height = parameters.get<double>("height", 1.0);
- if (imp == "homogeneous"){
+ if (imp == "homogeneous"){
- double mu1 = parameters.get<double>("mu1", 10);
+ double mu1 = parameters.get<double>("mu1", 10);
- auto muTerm = [mu1] (const Domain& z) {return mu1;};
+ auto muTerm = [mu1] (const Domain& z) {return mu1;};
- return muTerm;
-
- }
- else if (imp == "two_phase_material_1"){
- std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
-
-// Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
- Python::Module module = Python::import("two_phase_material_1");
- auto indicatorFunction = Python::make_function<double>(module.get("f"));
-
- auto muTerm = [mu,indicatorFunction] (const Domain& z)
- {
-// std::cout << "Test:" << indicatorFunction(z) << std::endl;
- if (indicatorFunction(z) == 1)
- return mu[0];
- else
- return mu[1];
- };
- return muTerm;
- }
- else if (imp == "two_phase_material_2"){
- std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
-
-// Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
- Python::Module module = Python::import("two_phase_material_2");
- auto indicatorFunction = Python::make_function<double>(module.get("f"));
-
- auto muTerm = [mu,indicatorFunction] (const Domain& z)
- {
-// std::cout << "Test:" << indicatorFunction(z) << std::endl;
- if (indicatorFunction(z) == 1)
- return mu[0];
- else
- return mu[1];
- };
- return muTerm;
- }
- else if (imp == "two_phase_material_3"){
- std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
-
-// Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
- Python::Module module = Python::import("two_phase_material_3");
- auto indicatorFunction = Python::make_function<double>(module.get("f"));
-
- auto muTerm = [mu,indicatorFunction] (const Domain& z)
- {
-// std::cout << "Test:" << indicatorFunction(z) << std::endl;
- if (indicatorFunction(z) == 1)
- return mu[0];
- else
- return mu[1];
- };
- return muTerm;
- }
- else if (imp == "isotropic_bilayer")
- {
- double mu1 = parameters.get<double>("mu1",10.0);
- double beta = parameters.get<double>("beta",2.0);
- double mu2 = beta*mu1;
+ return muTerm;
+
+ }
+ else if (imp == "material_neukamm"){
+ std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("material_neukamm"));
+ Python::Module module = Python::import("material_neukamm");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+ auto muTerm = [mu,indicatorFunction] (const Domain& z)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(z) == 1)
+ return mu[0];
+ else
+ return mu[1];
+ };
+ return muTerm;
+ }
+ else if (imp == "two_phase_material_1"){
+ std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+ Python::Module module = Python::import("two_phase_material_1");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+ auto muTerm = [mu,indicatorFunction] (const Domain& z)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(z) == 1)
+ return mu[0];
+ else
+ return mu[1];
+ };
+ return muTerm;
+ }
+ else if (imp == "two_phase_material_2"){
+ std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+ Python::Module module = Python::import("two_phase_material_2");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+ auto muTerm = [mu,indicatorFunction] (const Domain& z)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(z) == 1)
+ return mu[0];
+ else
+ return mu[1];
+ };
+ return muTerm;
+ }
+ else if (imp == "two_phase_material_3"){
+ std::array<double,2> mu = parameters.get<std::array<double,2>>("mu", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+ Python::Module module = Python::import("two_phase_material_3");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+ auto muTerm = [mu,indicatorFunction] (const Domain& z)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(z) == 1)
+ return mu[0];
+ else
+ return mu[1];
+ };
+ return muTerm;
+ }
+ else if (imp == "isotropic_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, theta] (const Domain& z) {
- if (z[2]>0)
- return mu1;
- else
- return mu2;
- };
+ auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+ if (z[2]>0)
+ return mu1;
+ else
+ return mu2;
+ };
- 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;
+ 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))
+ 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
+ else
return mu2;
};
- return muTerm;
- }
- else if (imp == "parametrized_Laminate"){
- double mu1 = parameters.get<double>("mu1",1.0);
- double beta = parameters.get<double>("beta",2.0);
- double mu2 = beta*mu1;
+ return muTerm;
+ }
+ else if (imp == "parametrized_Laminate"){
+ double mu1 = parameters.get<double>("mu1",1.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))
- if (abs(z[0]) > (theta/2.0))
+ auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+ // std::cout << "Analytical-MU is used" << std::endl;
+ // if (abs(z[0]) >= (theta/2.0))
+ if (abs(z[0]) > (theta/2.0))
return mu1;
- else
+ else
return mu2;
};
- return muTerm;
- }
- else if (imp == "circle_fiber"){
- double mu1 = parameters.get<double>("mu1",10.0);
- double beta = parameters.get<double>("beta",2.0);
- double mu2 = beta*mu1;
-
- auto muTerm = [mu1,mu2,width] (const Domain& z) {
- if (z[0] < 0 && sqrt(pow(z[1],2) + pow(z[2],2) ) < width/4.0 || 0 < z[0] && sqrt(pow(z[1],2) + pow(z[2],2) ) > width/4.0)
- return mu1;
- else
- return mu2;
+ return muTerm;
+ }
+ else if (imp == "circle_fiber"){
+ double mu1 = parameters.get<double>("mu1",10.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double mu2 = beta*mu1;
+
+ auto muTerm = [mu1,mu2,width] (const Domain& z) {
+ if (z[0] < 0 && sqrt(pow(z[1],2) + pow(z[2],2) ) < width/4.0 || 0 < z[0] && sqrt(pow(z[1],2) + pow(z[2],2) ) > width/4.0)
+ return mu1;
+ else
+ return mu2;
};
- return muTerm;
- }
- else if (imp == "square_fiber"){
- double mu1 = parameters.get<double>("mu1",10.0);
- double beta = parameters.get<double>("beta",2.0);
- double mu2 = beta*mu1;
-
- auto muTerm = [mu1,mu2,width] (const Domain& z) {
- if (z[0] < 0 && std::max(abs(z[1]), abs(z[2])) < width/4.0 || 0 < z[0] && std::max(abs(z[1]), abs(z[2])) > width/4.0)
- return mu1;
- else
- return mu2;
+ return muTerm;
+ }
+ else if (imp == "square_fiber"){
+ double mu1 = parameters.get<double>("mu1",10.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double mu2 = beta*mu1;
+
+ auto muTerm = [mu1,mu2,width] (const Domain& z) {
+ if (z[0] < 0 && std::max(abs(z[1]), abs(z[2])) < width/4.0 || 0 < z[0] && std::max(abs(z[1]), abs(z[2])) > width/4.0)
+ return mu1;
+ else
+ return mu2;
};
- return muTerm;
- }
- else if (imp == "matrix_material_circles") // Matrix material with prestrained Fiber inclusions (circular cross-section)
- {
- double mu1 = parameters.get<double>("mu1",10.0);
- double beta = parameters.get<double>("beta",2.0);
- double mu2 = beta*mu1;
+ return muTerm;
+ }
+ else if (imp == "matrix_material_circles") // Matrix material with prestrained Fiber inclusions (circular cross-section)
+ {
+ 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", 2); //number of Fibers in each Layer
- double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
+ int nF = parameters.get<int>("nF", 2); //number of Fibers in each Layer
+ double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
- assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+ 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;
- // TEST
- double domainShift = 0.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;
+ 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;
+ // TEST
+ double domainShift = 0.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;
- double output;
+ double output;
- // determine if point is in upper/lower Layer
- if ( 0 <= s[2] && s[2] <= 1.0/2.0) // upper Layer
- {
- for(size_t i=0; i<nF ;i++) // running through all the Fibers..
- {
- if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1)) // have to evenly space fibers...
- {
-// std::cout << " i : " << i << std::endl;
-// printvector(std::cout, s, "s" , "--");
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
-// printvector(std::cout, Fcenter, "Fcenter" , "--");
-//
- if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
- output = mu2;
-// return mu2; //richtig so? Fiber hat tendenziell größeres mu?
- else
- output = mu1;
-// return mu1;
- }
- else {}
- }
- }
- else // lower Layer
- {
- for(size_t i=0; i<nF ;i++) // running through all the Fibers..
- {
- if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1))
- {
-// std::cout << " i : " << i << std::endl;
-// printvector(std::cout, s, "s" , "--");
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
-// printvector(std::cout, Fcenter, "Fcenter" , "--");
+ // determine if point is in upper/lower Layer
+ if ( 0 <= s[2] && s[2] <= 1.0/2.0) // upper Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1)) // have to evenly space fibers...
+ {
+ // std::cout << " i : " << i << std::endl;
+ // printvector(std::cout, s, "s" , "--");
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+ // printvector(std::cout, Fcenter, "Fcenter" , "--");
+ //
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ output = mu2;
+ // return mu2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ output = mu1;
+ // return mu1;
+ }
+ else {}
+ }
+ }
+ else // lower Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1))
+ {
+ // std::cout << " i : " << i << std::endl;
+ // printvector(std::cout, s, "s" , "--");
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
+ // printvector(std::cout, Fcenter, "Fcenter" , "--");
- if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
- output = mu2;
-// return mu2; //richtig so? Fiber hat tendenziell größeres mu?
- else
- output = mu1;
-// return mu1;
- }
- else{}
- }
-
- }
- return output;
- };
- return muTerm;
- }
- else if (imp == "matrix_material_squares") // Matrix material with prestrained Fiber inclusions (square-shaped cross-section)
- {
- double mu1 = parameters.get<double>("mu1",10.0);
- double beta = parameters.get<double>("beta",2.0);
- double mu2 = beta*mu1;
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ output = mu2;
+ // return mu2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ output = mu1;
+ // return mu1;
+ }
+ else{}
+ }
+
+ }
+ return output;
+ };
+ return muTerm;
+ }
+ else if (imp == "matrix_material_squares") // Matrix material with prestrained Fiber inclusions (square-shaped cross-section)
+ {
+ 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", 2); //number of Fibers in each Layer
- double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
+ int nF = parameters.get<int>("nF", 2); //number of Fibers in each Layer
+ double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
- assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+ 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;
- // TEST
- double domainShift = 0.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;
- double output;
+ 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;
+ // TEST
+ double domainShift = 0.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;
+ double output;
- // determine if point is in upper/lower Layer
- if ( 0 <= s[2] && s[2] <= 1.0/2.0) // upper Layer
- {
- for(size_t i=0; i<nF ;i++) // running through all the Fibers..
- {
- if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1)) // have to evenly space fibers...
- {
-// std::cout << " i : " << i << std::endl;
-// printvector(std::cout, s, "s" , "--");
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
-// printvector(std::cout, Fcenter, "Fcenter" , "--");
-//
- if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
- output = mu2;
-// return mu2; //richtig so? Fiber hat tendenziell größeres mu?
- else
- output = mu1;
-// return mu1;
- }
- }
- }
- else // lower Layer
- {
- for(size_t i=0; i<nF ;i++) // running through all the Fibers..
- {
- if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1))
- {
-// std::cout << " i : " << i << std::endl;
-// printvector(std::cout, s, "s" , "--");
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
-// printvector(std::cout, Fcenter, "Fcenter" , "--");
+ // determine if point is in upper/lower Layer
+ if ( 0 <= s[2] && s[2] <= 1.0/2.0) // upper Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1)) // have to evenly space fibers...
+ {
+ // std::cout << " i : " << i << std::endl;
+ // printvector(std::cout, s, "s" , "--");
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+ // printvector(std::cout, Fcenter, "Fcenter" , "--");
+ //
+ if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
+ output = mu2;
+ // return mu2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ output = mu1;
+ // return mu1;
+ }
+ }
+ }
+ else // lower Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1))
+ {
+ // std::cout << " i : " << i << std::endl;
+ // printvector(std::cout, s, "s" , "--");
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
+ // printvector(std::cout, Fcenter, "Fcenter" , "--");
- if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
- output = mu2;
-// return mu2; //richtig so? Fiber hat tendenziell größeres mu?
- else
- output = mu1;
-// return mu1;
- }
- }
- }
- return output;
- };
- 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;
+ if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
+ output = mu2;
+ // return mu2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ output = mu1;
+ // return mu1;
+ }
+ }
+ }
+ return output;
+ };
+ 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) {
+ auto muTerm = [mu1, mu2, phi] (const Domain& z) {
if (isInRotatedPlane(phi, z[dim-2], z[dim-1]))
return mu1;
else
return mu2;
- };
+ };
- return muTerm;
- }
+ 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;
+ 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");
+ 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;
+ 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 == "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_prestrain_imp", "analytical_Example");
- //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?
+ 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_prestrain_imp", "analytical_Example");
+ //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);
+ double width = parameters.get<double>("width", 1.0);
+ double height = parameters.get<double>("height", 1.0);
- if (imp == "homogeneous"){
- double lambda1 = parameters.get<double>("lambda1",0.0);
+ if (imp == "homogeneous"){
+ double lambda1 = parameters.get<double>("lambda1",0.0);
- auto lambdaTerm = [lambda1] (const Domain& z) {return lambda1;};
+ auto lambdaTerm = [lambda1] (const Domain& z) {return lambda1;};
- return lambdaTerm;
- }
- if (imp == "two_phase_material_1"){
- std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
-
-// Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
- Python::Module module = Python::import("two_phase_material_1");
- auto indicatorFunction = Python::make_function<double>(module.get("f"));
-
- auto lambdaTerm = [lambda,indicatorFunction] (const Domain& z)
- {
-// std::cout << "Test:" << indicatorFunction(z) << std::endl;
- if (indicatorFunction(z) == 1)
- return lambda[0];
- else
- return lambda[1];
- };
- return lambdaTerm;
- }
- if (imp == "two_phase_material_2"){
- std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
-
-// Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
- Python::Module module = Python::import("two_phase_material_2");
- auto indicatorFunction = Python::make_function<double>(module.get("f"));
-
- auto lambdaTerm = [lambda,indicatorFunction] (const Domain& z)
- {
-// std::cout << "Test:" << indicatorFunction(z) << std::endl;
- if (indicatorFunction(z) == 1)
- return lambda[0];
- else
- return lambda[1];
- };
- return lambdaTerm;
- }
- if (imp == "two_phase_material_3"){
- std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
-
-// Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
- Python::Module module = Python::import("two_phase_material_3");
- auto indicatorFunction = Python::make_function<double>(module.get("f"));
-
- auto lambdaTerm = [lambda,indicatorFunction] (const Domain& z)
- {
-// std::cout << "Test:" << indicatorFunction(z) << std::endl;
- if (indicatorFunction(z) == 1)
- return lambda[0];
- else
- return lambda[1];
- };
- return lambdaTerm;
- }
- else if (imp == "isotropic_bilayer")
- {
- double lambda1 = parameters.get<double>("lambda1",0.0);
- double beta = parameters.get<double>("beta",2.0);
- double lambda2 = beta*lambda1;
+ return lambdaTerm;
+ }
+ if (imp == "material_neukamm"){
+ std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
+
+ Python::Module module = Python::import("material_neukamm");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+ auto lambdaTerm = [lambda,indicatorFunction] (const Domain& z)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(z) == 1)
+ return lambda[0];
+ else
+ return lambda[1];
+ };
+ return lambdaTerm;
+ }
+ if (imp == "two_phase_material_1"){
+ std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+ Python::Module module = Python::import("two_phase_material_1");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+ auto lambdaTerm = [lambda,indicatorFunction] (const Domain& z)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(z) == 1)
+ return lambda[0];
+ else
+ return lambda[1];
+ };
+ return lambdaTerm;
+ }
+ if (imp == "two_phase_material_2"){
+ std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+ Python::Module module = Python::import("two_phase_material_2");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+ auto lambdaTerm = [lambda,indicatorFunction] (const Domain& z)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(z) == 1)
+ return lambda[0];
+ else
+ return lambda[1];
+ };
+ return lambdaTerm;
+ }
+ if (imp == "two_phase_material_3"){
+ std::array<double,2> lambda = parameters.get<std::array<double,2>>("lambda", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+ Python::Module module = Python::import("two_phase_material_3");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+ auto lambdaTerm = [lambda,indicatorFunction] (const Domain& z)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(z) == 1)
+ return lambda[0];
+ else
+ return lambda[1];
+ };
+ return lambdaTerm;
+ }
+ else if (imp == "isotropic_bilayer")
+ {
+ 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) {
- if (z[2]>0)
- return lambda1;
- else
- return lambda2;
- };
- 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) {
+ if (z[2]>0)
+ return lambda1;
+ else
+ return lambda2;
+ };
+ 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) {
+ 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 == "parametrized_Laminate"){
- double lambda1 = parameters.get<double>("lambda1",0.0);
- double beta = parameters.get<double>("beta",2.0);
- double lambda2 = beta*lambda1;
+ // 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 == "parametrized_Laminate"){
+ 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) {
+ auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
-// std::cout << "Analytical-LAMBDA is used" << std::endl; //TEST
-// std::cout << "Lambda1:" << lambda1 << std::endl;
-// if (abs(z[0]) > (theta/2.0))
- if (abs(z[0]) > (theta/2.0))
- return lambda1;
- else
- return lambda2;
- };
- return lambdaTerm;
- }
- else if (imp == "circle_fiber"){
- double lambda1 = parameters.get<double>("lambda1",0.0);
- double beta = parameters.get<double>("beta",2.0);
- double lambda2 = beta*lambda1;
-
- auto lambdaTerm = [lambda1,lambda2,width] (const Domain& z) // Prestrain inducing twist (bisher nur extension)
+ // std::cout << "Analytical-LAMBDA is used" << std::endl; //TEST
+ // std::cout << "Lambda1:" << lambda1 << std::endl;
+ // if (abs(z[0]) > (theta/2.0))
+ if (abs(z[0]) > (theta/2.0))
+ return lambda1;
+ else
+ return lambda2;
+ };
+ return lambdaTerm;
+ }
+ else if (imp == "circle_fiber"){
+ double lambda1 = parameters.get<double>("lambda1",0.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double lambda2 = beta*lambda1;
+
+ auto lambdaTerm = [lambda1,lambda2,width] (const Domain& z) // Prestrain inducing twist (bisher nur extension)
{
- if (z[0] < 0 && sqrt(pow(z[1],2) + pow(z[2],2) ) < width/4.0 || 0 < z[0] && sqrt(pow(z[1],2) + pow(z[2],2) ) > width/4.0)
- return lambda1;
- else
- return lambda2;
- };
+ if (z[0] < 0 && sqrt(pow(z[1],2) + pow(z[2],2) ) < width/4.0 || 0 < z[0] && sqrt(pow(z[1],2) + pow(z[2],2) ) > width/4.0)
+ return lambda1;
+ else
+ return lambda2;
+ };
- return lambdaTerm;
- }
- else if (imp == "square_fiber"){
- double lambda1 = parameters.get<double>("lambda1",0.0);
- double beta = parameters.get<double>("beta",2.0);
- double lambda2 = beta*lambda1;
-
- auto lambdaTerm = [lambda1,lambda2,width] (const Domain& z) // Prestrain inducing twist (bisher nur extension)
+ return lambdaTerm;
+ }
+ else if (imp == "square_fiber"){
+ double lambda1 = parameters.get<double>("lambda1",0.0);
+ double beta = parameters.get<double>("beta",2.0);
+ double lambda2 = beta*lambda1;
+
+ auto lambdaTerm = [lambda1,lambda2,width] (const Domain& z) // Prestrain inducing twist (bisher nur extension)
{
- if (z[0] < 0 && std::max(abs(z[1]), abs(z[2])) < width/4.0 || 0 < z[0] && std::max(abs(z[1]), abs(z[2])) > width/4.0)
- return lambda1;
- else
- return lambda2;
- };
+ if (z[0] < 0 && std::max(abs(z[1]), abs(z[2])) < width/4.0 || 0 < z[0] && std::max(abs(z[1]), abs(z[2])) > width/4.0)
+ return lambda1;
+ else
+ return lambda2;
+ };
- return lambdaTerm;
- }
- else if (imp == "matrix_material_circles") // Matrix material with prestrained Fiber inclusions (circular cross-section)
- {
- double lambda1 = parameters.get<double>("lambda1",0.0);
- double beta = parameters.get<double>("beta",2.0);
- double lambda2 = beta*lambda1;
+ return lambdaTerm;
+ }
+ else if (imp == "matrix_material_circles") // Matrix material with prestrained Fiber inclusions (circular cross-section)
+ {
+ 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", 2); //number of Fibers in each Layer
- double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
+ int nF = parameters.get<int>("nF", 2); //number of Fibers in each Layer
+ double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
- assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+ 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;
- // TEST
- double domainShift = 0.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;
- double output;
+ 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;
+ // TEST
+ double domainShift = 0.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;
+ double output;
- // determine if point is in upper/lower Layer
- if ( 0 <= s[2] && s[2] <= 1.0/2.0) // upper Layer
- {
- for(size_t i=0; i<nF ;i++) // running through all the Fibers..
- {
- if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1)) // have to evenly space fibers...
- {
-// std::cout << " i : " << i << std::endl;
-// printvector(std::cout, s, "s" , "--");
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
-// printvector(std::cout, Fcenter, "Fcenter" , "--");
-//
- if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
- output = lambda2;
-// return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
- else
- output = lambda1;
-// return lambda1;
- }
- }
- }
- else // lower Layer
- {
- for(size_t i=0; i<nF ;i++) // running through all the Fibers..
- {
- if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1))
- {
-// std::cout << " i : " << i << std::endl;
-// printvector(std::cout, s, "s" , "--");
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
-// printvector(std::cout, Fcenter, "Fcenter" , "--");
+ // determine if point is in upper/lower Layer
+ if ( 0 <= s[2] && s[2] <= 1.0/2.0) // upper Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1)) // have to evenly space fibers...
+ {
+ // std::cout << " i : " << i << std::endl;
+ // printvector(std::cout, s, "s" , "--");
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+ // printvector(std::cout, Fcenter, "Fcenter" , "--");
+ //
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ output = lambda2;
+ // return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ output = lambda1;
+ // return lambda1;
+ }
+ }
+ }
+ else // lower Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1))
+ {
+ // std::cout << " i : " << i << std::endl;
+ // printvector(std::cout, s, "s" , "--");
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
+ // printvector(std::cout, Fcenter, "Fcenter" , "--");
- if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
- output = lambda2;
-// return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
- else
- output = lambda1;
-// return lambda1;
- }
- }
-
- }
- return output;
- };
- return lambdaTerm;
- }
- else if (imp == "matrix_material_squares") // Matrix material with prestrained Fiber inclusions (square-shaped cross-section)
- {
- double lambda1 = parameters.get<double>("lambda1",0.0);
- double beta = parameters.get<double>("beta",2.0);
- double lambda2 = beta*lambda1;
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ output = lambda2;
+ // return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ output = lambda1;
+ // return lambda1;
+ }
+ }
+
+ }
+ return output;
+ };
+ return lambdaTerm;
+ }
+ else if (imp == "matrix_material_squares") // Matrix material with prestrained Fiber inclusions (square-shaped cross-section)
+ {
+ 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", 2); //number of Fibers in each Layer
- double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
+ int nF = parameters.get<int>("nF", 2); //number of Fibers in each Layer
+ double rF = parameters.get<double>("rF", 0.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
- assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+ 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;
- // TEST
- double domainShift = 0.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;
- double output;
+ 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;
+ // TEST
+ double domainShift = 0.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;
+ double output;
- // determine if point is in upper/lower Layer
- if ( 0 <= s[2] && s[2] <= 1.0/2.0) // upper Layer
- {
- for(size_t i=0; i<nF ;i++) // running through all the Fibers..
- {
- if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1)) // have to evenly space fibers...
- {
-// std::cout << " i : " << i << std::endl;
-// printvector(std::cout, s, "s" , "--");
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
-// printvector(std::cout, Fcenter, "Fcenter" , "--");
-//
- if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
- output = lambda2;
-// return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
- else
- output = lambda1;
-// return lambda1;
- }
- }
- }
- else // lower Layer
- {
- for(size_t i=0; i<nF ;i++) // running through all the Fibers..
- {
- if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1))
- {
-// std::cout << " i : " << i << std::endl;
-// printvector(std::cout, s, "s" , "--");
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
-// printvector(std::cout, Fcenter, "Fcenter" , "--");
+ // determine if point is in upper/lower Layer
+ if ( 0 <= s[2] && s[2] <= 1.0/2.0) // upper Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1)) // have to evenly space fibers...
+ {
+ // std::cout << " i : " << i << std::endl;
+ // printvector(std::cout, s, "s" , "--");
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+ // printvector(std::cout, Fcenter, "Fcenter" , "--");
+ //
+ if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
+ output = lambda2;
+ // return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ output = lambda1;
+ // return lambda1;
+ }
+ }
+ }
+ else // lower Layer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if(-1.0/2.0 + (1.0/double(nF))*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/double(nF))*(i+1))
+ {
+ // std::cout << " i : " << i << std::endl;
+ // printvector(std::cout, s, "s" , "--");
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0 };
+ // printvector(std::cout, Fcenter, "Fcenter" , "--");
- if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
- output = lambda2;
-// return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
- else
- output = lambda1;
-// return lambda1;
- }
- }
- }
- return output;
- };
- return lambdaTerm;
- }
- else if (imp == "bilayer_lame"){
- double lambda1 = parameters.get<double>("mu1",384.615);
- double lambda2 = parameters.get<double>("mu2",384.615);
+ if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
+ output = lambda2;
+ // return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ output = lambda1;
+ // return lambda1;
+ }
+ }
+ }
+ return output;
+ };
+ 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; };
+ 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 == "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 == "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;
-// }
- }
+ 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;
+ // }
+ }
@@ -897,319 +930,338 @@ 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 = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
-// using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
- using ScalarRT = double;
- using Func2Tensor = std::function< MatrixRT(const Domain&) >;
- using FuncScalar = std::function< ScalarRT(const Domain&) >;
- using FuncMatrix = std::function< MatrixRT(const Domain&) >;
- using FuncVector = std::function< VectorRT(const Domain&) >;
+ // static const int dim = 3;
+ static const int dimWorld = 3;
+ using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+ // using GridView = CellGridType::LeafGridView;
+ using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
+ // using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+ using ScalarRT = double;
+ using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+ using FuncScalar = std::function< ScalarRT(const Domain&) >;
+ using FuncMatrix = std::function< MatrixRT(const Domain&) >;
+ using FuncVector = std::function< VectorRT(const Domain&) >;
protected:
-// double p1, p2, theta;
-// double width; //Cell geometry
+ // double p1, p2, theta;
+ // double width; //Cell geometry
public:
-// Func2Tensor getPrestrain(std::string imp)
- Func2Tensor getPrestrain(ParameterTree parameters)
- {
+ // Func2Tensor getPrestrain(std::string imp)
+ Func2Tensor getPrestrain(ParameterTree parameters)
+ {
- std::string imp = parameters.get<std::string>("material_prestrain_imp", "analytical_Example");
+ std::string imp = parameters.get<std::string>("material_prestrain_imp", "analytical_Example");
- double width = parameters.get<double>("width", 1.0);
- double height = parameters.get<double>("height", 1.0);
+ 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;
+ 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 == "homogeneous"){
- Func2Tensor B = [p1] (const Domain& x) // ISOTROPIC PRESSURE
- {
- return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
- };
- std::cout <<" Prestrain Type: homogeneous" << std::endl;
- return B;
- }
- else if (imp == "two_phase_material_1"){
- std::array<double,2> rho = parameters.get<std::array<double,2>>("rho", {1.0,3.0});
-
-// Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
- Python::Module module = Python::import("two_phase_material_1");
- auto indicatorFunction = Python::make_function<double>(module.get("f"));
+ if (imp == "homogeneous"){
+ Func2Tensor B = [p1] (const Domain& x) // ISOTROPIC PRESSURE
+ {
+ return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ };
+ std::cout <<" Prestrain Type: homogeneous" << std::endl;
+ return B;
+ }
+ else if (imp == "two_phase_material_1"){
+ std::array<double,2> rho = parameters.get<std::array<double,2>>("rho", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+ Python::Module module = Python::import("two_phase_material_1");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
- Func2Tensor B = [rho,indicatorFunction] (const Domain& x)
- {
-// std::cout << "Test:" << indicatorFunction(z) << std::endl;
- if (indicatorFunction(x) == 1)
- return MatrixRT{{rho[0], 0.0 , 0.0}, {0.0, rho[0], 0.0}, {0.0, 0.0, rho[0]}};
- else
- return MatrixRT{{rho[1], 0.0 , 0.0}, {0.0, rho[1], 0.0}, {0.0, 0.0, rho[1]}};
- };
- std::cout <<" Prestrain Type: two_phase_material_1" << std::endl;
- return B;
- }
- else if (imp == "two_phase_material_2"){
- std::array<double,2> rho = parameters.get<std::array<double,2>>("rho", {1.0,3.0});
-
-// Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
- Python::Module module = Python::import("two_phase_material_2");
- auto indicatorFunction = Python::make_function<double>(module.get("f"));
+ Func2Tensor B = [rho,indicatorFunction] (const Domain& x)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(x) == 1)
+ return MatrixRT{{rho[0], 0.0 , 0.0}, {0.0, rho[0], 0.0}, {0.0, 0.0, rho[0]}};
+ else
+ return MatrixRT{{rho[1], 0.0 , 0.0}, {0.0, rho[1], 0.0}, {0.0, 0.0, rho[1]}};
+ };
+ std::cout <<" Prestrain Type: two_phase_material_1" << std::endl;
+ return B;
+ }
+ else if (imp == "material_neukamm"){
+ std::array<double,2> rho = parameters.get<std::array<double,2>>("rho", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+ Python::Module module = Python::import("material_neukamm");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
- Func2Tensor B = [rho,indicatorFunction] (const Domain& x)
- {
-// std::cout << "Test:" << indicatorFunction(z) << std::endl;
- if (indicatorFunction(x) == 1)
- return MatrixRT{{rho[0], 0.0 , 0.0}, {0.0, rho[0], 0.0}, {0.0, 0.0, rho[0]}};
- else
- return MatrixRT{{rho[1], 0.0 , 0.0}, {0.0, rho[1], 0.0}, {0.0, 0.0, rho[1]}};
- };
- std::cout <<" Prestrain Type: two_phase_material_2" << std::endl;
- return B;
- }
- else if (imp == "two_phase_material_3"){
- std::array<double,2> rho = parameters.get<std::array<double,2>>("rho", {1.0,3.0});
-
-// Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
- Python::Module module = Python::import("two_phase_material_3");
- auto indicatorFunction = Python::make_function<double>(module.get("f"));
+ Func2Tensor B = [rho,indicatorFunction] (const Domain& x)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(x) == 1)
+ return MatrixRT{{rho[0], 0.0 , 0.0}, {0.0, rho[0], 0.0}, {0.0, 0.0, rho[0]}};
+ else
+ return MatrixRT{{rho[1], 0.0 , 0.0}, {0.0, rho[1], 0.0}, {0.0, 0.0, rho[1]}};
+ };
+ std::cout <<" Prestrain Type: two_phase_material_1" << std::endl;
+ return B;
+ }
+ else if (imp == "two_phase_material_2"){
+ std::array<double,2> rho = parameters.get<std::array<double,2>>("rho", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+ Python::Module module = Python::import("two_phase_material_2");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
- Func2Tensor B = [rho,indicatorFunction] (const Domain& x)
- {
-// std::cout << "Test:" << indicatorFunction(z) << std::endl;
- if (indicatorFunction(x) == 1)
- return MatrixRT{{rho[0], 0.0 , 0.0}, {0.0, rho[0], 0.0}, {0.0, 0.0, rho[0]}};
- else
- return MatrixRT{{rho[1], 0.0 , 0.0}, {0.0, rho[1], 0.0}, {0.0, 0.0, rho[1]}};
- };
- std::cout <<" Prestrain Type: two_phase_material_3" << std::endl;
- return B;
- }
- else if (imp == "isotropic_bilayer")
- {
- Func2Tensor B = [p1,p2,theta] (const Domain& x) // ISOTROPIC PRESSURE
- {
- if (x[2] > 0)
- return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
- else if (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;
- }
+ Func2Tensor B = [rho,indicatorFunction] (const Domain& x)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(x) == 1)
+ return MatrixRT{{rho[0], 0.0 , 0.0}, {0.0, rho[0], 0.0}, {0.0, 0.0, rho[0]}};
+ else
+ return MatrixRT{{rho[1], 0.0 , 0.0}, {0.0, rho[1], 0.0}, {0.0, 0.0, rho[1]}};
+ };
+ std::cout <<" Prestrain Type: two_phase_material_2" << std::endl;
+ return B;
+ }
+ else if (imp == "two_phase_material_3"){
+ std::array<double,2> rho = parameters.get<std::array<double,2>>("rho", {1.0,3.0});
+
+ // Python::Module module = Python::import(parameters.get<std::string>("two_phase_material_1"));
+ Python::Module module = Python::import("two_phase_material_3");
+ auto indicatorFunction = Python::make_function<double>(module.get("f"));
+
+
+ Func2Tensor B = [rho,indicatorFunction] (const Domain& x)
+ {
+ // std::cout << "Test:" << indicatorFunction(z) << std::endl;
+ if (indicatorFunction(x) == 1)
+ return MatrixRT{{rho[0], 0.0 , 0.0}, {0.0, rho[0], 0.0}, {0.0, 0.0, rho[0]}};
+ else
+ return MatrixRT{{rho[1], 0.0 , 0.0}, {0.0, rho[1], 0.0}, {0.0, 0.0, rho[1]}};
+ };
+ std::cout <<" Prestrain Type: two_phase_material_3" << std::endl;
+ return B;
+ }
+ else if (imp == "isotropic_bilayer")
+ {
+ Func2Tensor B = [p1,p2,theta] (const Domain& x) // ISOTROPIC PRESSURE
+ {
+ if (x[2] > 0)
+ return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ else if (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 == "parametrized_Laminate")
- {
- Func2Tensor B = [p1,p2,theta] (const Domain& x)
- {
- 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 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}};
- 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: parametrized_Laminate"<< std::endl;
- return B;
- }
- else if (imp == "circle_fiber"){
-
- Func2Tensor B = [p1,theta,width] (const Domain& x) // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE
- {
+ 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 == "parametrized_Laminate")
+ {
+ Func2Tensor B = [p1,p2,theta] (const Domain& x)
+ {
+ 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 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}};
+ 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: parametrized_Laminate"<< std::endl;
+ return B;
+ }
+ else if (imp == "circle_fiber"){
+
+ Func2Tensor B = [p1,theta,width] (const Domain& x) // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE
+ {
if (x[0] < 0 && sqrt(pow(x[1],2) + pow(x[2],2) ) < width/4.0 || 0 < x[0] && sqrt(pow(x[1],2) + pow(x[2],2) ) > width/4.0)
- return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ 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: circle_fiber"<< std::endl;
- return B;
- }
- else if (imp == "square_fiber"){
+ return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
+ std::cout <<" Prestrain Type: circle_fiber"<< std::endl;
+ return B;
+ }
+ else if (imp == "square_fiber"){
- Func2Tensor B = [p1,theta,width] (const Domain& x) // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE
- {
+ Func2Tensor B = [p1,theta,width] (const Domain& x) // Bilayer with one rectangular Fiber & ISOTROPIC PRESSURE
+ {
if (x[0] < 0 && std::max(abs(x[1]), abs(x[2])) < width/4.0 || 0 < x[0] && std::max(abs(x[1]), abs(x[2])) > width/4.0)
- return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ 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: square_fiber"<< std::endl;
- return B;
- }
- else if (imp == "matrix_material_circles") // 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.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
+ return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
+ std::cout <<" Prestrain Type: square_fiber"<< std::endl;
+ return B;
+ }
+ else if (imp == "matrix_material_circles") // 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.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
- assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+ 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;
- double domainShift = 0.0;
- // shift x to domain [0,1]^3
- FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+ 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;
+ double domainShift = 0.0;
+ // shift x to domain [0,1]^3
+ FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
- MatrixRT output;
-
- // determine if point is in upper/lower Layer
- if ((0.0 <= s[2] && s[2] <= 1.0/2.0)) // upperLayer
- {
- for(size_t i=0; i<nF ;i++) // running through all the Fibers..
- {
- if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/nF)*(i+1))
- {
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
-//
- if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
- output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
-// return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
- else
- output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-// 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(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/nF)*(i+1))
- {
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0};
+ MatrixRT output;
+
+ // determine if point is in upper/lower Layer
+ if ((0.0 <= s[2] && s[2] <= 1.0/2.0)) // upperLayer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/nF)*(i+1))
+ {
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+ //
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ // return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ else
+ output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ // 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(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/nF)*(i+1))
+ {
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0};
- if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
- output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
-// return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
- else
- output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-// 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}};
- }
- }
-
- }
- return output;
- };
- std::cout <<" Prestrain Type: matrix_material"<< std::endl;
- return B;
- }
+ if(sqrt(pow(s[0]-Fcenter[0],2)+pow(s[2]-Fcenter[1],2)) <= rF )
+ output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ // return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ else
+ output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ // 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}};
+ }
+ }
+
+ }
+ return output;
+ };
+ std::cout <<" Prestrain Type: matrix_material"<< std::endl;
+ return B;
+ }
- else if (imp == "matrix_material_squares") // 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.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
+ else if (imp == "matrix_material_squares") // 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.5*(width/(2.0*nF)) ); //default: half of the max-fiber-radius mrF = (width/(2.0*nF))
- assert( (2*rF)*nF <= width && (height/4)+rF <= height); //check that fibers are not bigger than Domain
+ 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;
- double domainShift = 0.0;
- // shift x to domain [0,1]^3
- FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
+ 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;
+ double domainShift = 0.0;
+ // shift x to domain [0,1]^3
+ FieldVector<double,3> s = {x[0]+domainShift, x[1]+domainShift, x[2]+domainShift};
- MatrixRT output;
-
- // determine if point is in upper/lower Layer
- if ((0.0 <= s[2] && s[2] <= 1.0/2.0)) // upperLayer
- {
- for(size_t i=0; i<nF ;i++) // running through all the Fibers..
- {
- if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/nF)*(i+1))
- {
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
-//
- if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
- output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
-// return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
- else
- output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-// 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(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/nF)*(i+1))
- {
- // compute Fiber center
- FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0};
+ MatrixRT output;
+
+ // determine if point is in upper/lower Layer
+ if ((0.0 <= s[2] && s[2] <= 1.0/2.0)) // upperLayer
+ {
+ for(size_t i=0; i<nF ;i++) // running through all the Fibers..
+ {
+ if(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/nF)*(i+1))
+ {
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , 1.0/4.0};
+ //
+ if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
+ output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ // return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ else
+ output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ // 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(-1.0/2.0 + (1.0/nF)*i<= s[0] && s[0] <= -1.0/2.0 + (1.0/nF)*(i+1))
+ {
+ // compute Fiber center
+ FieldVector<double,2> Fcenter = { (1.0/(2.0*nF))+((1.0/double(nF))*i)-(1.0/2.0) , -1.0/4.0};
- if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
- output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
-// return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
- else
- output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
-// 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}};
- }
- }
-
- }
- return output;
- };
- std::cout <<" Prestrain Type: matrix_material"<< std::endl;
- return B;
- }
- else
- {
- // TODO other geometries etc...
+ if(std::max( abs(s[0]-Fcenter[0]), abs(s[2]-Fcenter[1]) ) <= rF )
+ output = MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ // return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ else
+ output = MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ // 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}};
+ }
+ }
+
+ }
+ return output;
+ };
+ std::cout <<" Prestrain Type: matrix_material"<< std::endl;
+ return B;
+ }
+ else
+ {
+ // TODO other geometries etc...
- }
- // TODO ANISOTROPIC PRESSURE
+ }
+ // TODO ANISOTROPIC PRESSURE
- }
+ }
- FuncScalar getNormPrestrain(const FuncMatrix& B)
- {
- FuncScalar normB = [&](const Domain& z) {
- return norm(B(z));
- };
- return normB;
- }
+ FuncScalar getNormPrestrain(const FuncMatrix& B)
+ {
+ FuncScalar normB = [&](const Domain& z) {
+ return norm(B(z));
+ };
+ return normB;
+ }
};
diff --git a/geometries/material_neukamm.py b/geometries/material_neukamm.py
new file mode 100644
index 00000000..d513ae98
--- /dev/null
+++ b/geometries/material_neukamm.py
@@ -0,0 +1,16 @@
+import math
+
+
+#Indicator function that determines both phases
+# x[0] : y1-component
+# x[1] : y2-component
+# x[2] : x3-component
+def f(x):
+ theta=0.25
+ # --- replace with your definition of indicatorFunction:
+ if ((abs(x[0]) < theta/2) and x[2]<-1/2+theta):
+ return 1 #Phase1
+ elif ((abs(x[1]) < theta/2) and x[2]>1/2-theta):
+ return 1 #Phase1
+ else :
+ return 0 #Phase2
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 7a95dd04..e977338d 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -64,9 +64,9 @@ lambda1=1.0
# better: pass material parameters as a vector
-mu=1.0 3.0
-lambda=1.0 3.0
-rho=1.0 8.0
+mu=80 60
+lambda=80 25
+rho=1.0 0
# ---volume fraction (default value = 1.0/4.0)
@@ -80,8 +80,7 @@ theta=0.25
#--- choose composite-Implementation:
geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
-#geometryFunction = two_phase_material_1
-material_prestrain_imp= "two_phase_material_1" #(Python-indicator-function with same name determines material phases)
+material_prestrain_imp= "material_neukamm" #(Python-indicator-function with same name determines material phases)
#material_prestrain_imp= "two_phase_material_2" #(Python-indicator-function with same name determines material phases)
#material_prestrain_imp= "two_phase_material_3" #(Python-indicator-function with same name determines material phases)
#material_prestrain_imp= "parametrized_Laminate"
diff --git a/outputs/Results/1/BMatrix.txt b/outputs/Results/1/BMatrix.txt
new file mode 100644
index 00000000..0fb3c029
--- /dev/null
+++ b/outputs/Results/1/BMatrix.txt
@@ -0,0 +1,3 @@
+1 1 -0.0144287202371896055
+1 2 0.014428690223795131
+1 3 -2.0320334857961725e-11
diff --git a/outputs/Results/1/QMatrix.txt b/outputs/Results/1/QMatrix.txt
new file mode 100644
index 00000000..9d8813d5
--- /dev/null
+++ b/outputs/Results/1/QMatrix.txt
@@ -0,0 +1,9 @@
+1 1 12.4813265250416485
+1 2 2.0831795458888962
+1 3 3.83811740666402797e-10
+2 1 2.0831795513179836
+2 2 12.481326525238547
+2 3 1.78942192122814843e-10
+3 1 -6.40561063029264521e-10
+3 2 -1.00945173229121783e-09
+3 3 10.3922965257922275
diff --git a/outputs/Results/1/output.txt b/outputs/Results/1/output.txt
new file mode 100644
index 00000000..215cb327
--- /dev/null
+++ b/outputs/Results/1/output.txt
@@ -0,0 +1,61 @@
+material_prestrain used: material_neukamm
+----- Input Parameters -----:
+alpha: 8
+gamma: 1
+theta: 0.25
+beta: 3
+material parameters:
+mu1: 1
+mu2: 3
+lambda1: 1
+lambda2: 3
+----------------------------:
+Number of Elements in each direction: [4,4,4]
+size of FiniteElementBasis: 240
+Solver-type used: CG-Solver
+---------- OUTPUT ----------
+ --------------------
+Corrector-Matrix M_1:
+-0.00922212 -2.49142e-10 0
+-2.49142e-10 -0.00494017 0
+0 0 0
+
+ --------------------
+Corrector-Matrix M_2:
+0.00494016 -2.03726e-10 0
+-2.03726e-10 0.00922211 0
+0 0 0
+
+ --------------------
+Corrector-Matrix M_3:
+1.90697e-10 -4.74653e-11 0
+-4.74653e-11 1.93256e-10 0
+0 0 0
+
+ --------------------
+Effective Matrix Q:
+12.4813 2.08318 3.83812e-10
+2.08318 12.4813 1.78942e-10
+-6.40561e-10 -1.00945e-09 10.3923
+
+--- Prestrain Output ---
+B_hat: -0.150032 0.150032 -2.16498e-10
+B_eff: -0.0144287 0.0144287 -2.03203e-11 (Effective Prestrain)
+------------------------
+q1=12.4813
+q2=12.4813
+q3=10.3923
+q12=2.08318
+b1=-0.0144287
+b2=0.0144287
+b3=-2.03203e-11
+b1_hat: -0.150032
+b2_hat: 0.150032
+b3_hat: -2.16498e-10
+mu_gamma=10.3923
+q_onetwo=2.083180
+---------------------------------------------------------------------------------------------------------
+ Levels | q1 | q2 | q3 | b1 | b2 | b3 |
+---------------------------------------------------------------------------------------------------------
+ 2 & 1.24813e+01 & 1.24813e+01 & 1.03923e+01 & -1.44287e-02 & 1.44287e-02 & -2.03203e-11 &
+---------------------------------------------------------------------------------------------------------
diff --git a/outputs/Results/1/parameter b/outputs/Results/1/parameter
new file mode 100644
index 00000000..f19d4939
--- /dev/null
+++ b/outputs/Results/1/parameter
@@ -0,0 +1,14 @@
+mu=80 60
+lambda=80 25
+rho=1.0 0
+
+def f(x):
+ theta=0.25
+ # --- replace with your definition of indicatorFunction:
+ if ((abs(x[0]) < theta/2) and x[2]<-1/2+theta):
+ return 1 #Phase1
+ elif ((abs(x[1]) < theta/2) and x[2]>1/2-theta):
+ return 1 #Phase1
+ else :
+ return 0 #Phase2
+
--
GitLab