From 7e709f8cfe642fda197c3733f6e0048007adab54 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Sat, 28 Aug 2021 18:07:56 +0200
Subject: [PATCH] Update MatrixMaterial and Plot MaterialFunction
---
.../prestrain_material_geometry.hh | 364 ++++++++++++++++--
inputs/cellsolver.parset | 22 +-
src/dune-microstructure.cc | 59 +++
3 files changed, 399 insertions(+), 46 deletions(-)
diff --git a/dune/microstructure/prestrain_material_geometry.hh b/dune/microstructure/prestrain_material_geometry.hh
index c52107ca..ad03b2ae 100644
--- a/dune/microstructure/prestrain_material_geometry.hh
+++ b/dune/microstructure/prestrain_material_geometry.hh
@@ -53,8 +53,6 @@ public:
return muTerm;
}
-
-
else if (imp == "analytical_Example"){
double mu1 = parameters.get<double>("mu1",10.0);
double beta = parameters.get<double>("beta",2.0);
@@ -71,34 +69,69 @@ public:
return muTerm;
}
- else if (imp == "matrix_material") // Matrix material with prestrained Fiber inclusions
- {
- double mu1 = parameters.get<double>("mu1",10.0);
+ 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 == "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", 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;
+// 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;
// determine if point is in upper/lower Layer
- if ((height/2) <= s[2]<= height) // upper 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((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ 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 = { (width/(2*nF))+((width/nF)*i) , (3/4)*height };
+ 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 )
return mu2; //richtig so? Fiber hat tendenziell größeres mu?
@@ -111,11 +144,14 @@ public:
{
for(size_t i=0; i<nF ;i++) // running through all the Fibers..
{
- if((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ 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 = { (width/(2*nF))+((width/nF)*i) , height/4.0 };
-
+ 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 )
return mu2; //richtig so? Fiber hat tendenziell größeres mu?
else
@@ -128,6 +164,73 @@ public:
};
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))
+
+
+
+ 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;
+
+ // 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 )
+ 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(-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 )
+ 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);
@@ -295,37 +398,73 @@ public:
return lambdaTerm;
}
- else if (imp == "matrix_material") // Matrix material with prestrained Fiber inclusions
+ 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;
+ };
+
+ 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;
+ };
+
+ 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", 3); //number of Fibers in each Layer
- double rF = parameters.get<double>("rF", 0.1);
+ 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
- auto lambdaTerm = [lambda1,lambda2, theta,nF,rF,height,width] (const Domain& x)
+ 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;
+// 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-LAMBDA is used" << std::endl;
+// std::cout << "matrixMaterial-MU is used" << std::endl;
// determine if point is in upper/lower Layer
- if ((height/2) <= s[2]<= height) // upper 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((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ 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 = { (width/(2*nF))+((width/nF)*i) , (3/4)*height };
+ 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 )
- return lambda2;
+ return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
else
return lambda1;
}
@@ -335,13 +474,16 @@ public:
{
for(size_t i=0; i<nF ;i++) // running through all the Fibers..
{
- if((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ 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 = { (width/(2*nF))+((width/nF)*i) , height/4.0 };
-
+ 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 )
- return lambda2;
+ return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
else
return lambda1;
}
@@ -352,8 +494,73 @@ public:
};
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))
+
+
+
+ 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;
+
+ // 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 )
+ return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
+ else
+ 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 )
+ return lambda2; //richtig so? Fiber hat tendenziell größeres mu?
+ 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);
@@ -569,31 +776,54 @@ public:
std::cout <<" Prestrain Type: analytical_Example "<< std::endl;
return B;
}
- else if (imp == "matrix_material") // Matrix material with prestrained Fiber inclusions
+ 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}};
+ 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"){
+
+ 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}};
+ 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.1);
+ 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
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 = 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};
-
-
-
+
// determine if point is in upper/lower Layer
- if ((height/2) <= s[2]<= height) // upper 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((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ 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 = { (width/(2*nF))+((width/nF)*i) , (3/4)*height };
+ 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 )
return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
@@ -607,10 +837,10 @@ public:
{
for(size_t i=0; i<nF ;i++) // running through all the Fibers..
{
- if((width/nF)*i<= s[0] <= (width/nF)*(i+1))
+ 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 = { (width/(2*nF))+((width/nF)*i) , height/4.0 };
+ 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 )
return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
@@ -626,6 +856,62 @@ public:
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))
+
+ 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};
+
+ // 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 )
+ 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(-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 )
+ 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...
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index e0b00e86..af1a8964 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -11,6 +11,7 @@ outputPath = "../../outputs/output.txt"
#print_debug = true #default = false
+
#############################################
# Cell Domain
#############################################
@@ -28,8 +29,8 @@ cellDomain = 1
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
########################################################################
-numLevels = 1 3 # computes all levels from first to second entry
-#numLevels = 3 3 # 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
@@ -58,17 +59,23 @@ gamma = 1.0
# Material parameters
#############################################
+write_materialFunctions = true # VTK mu-functions , lambda-functions
+
beta = 2.0 # ratio between material parameters mu1 & mu2 .... beta = 1.0 corresponds to homogeneous case
mu1 = 10.0
-lambda1 = 0.0
+lambda1 = 20.0
#lambda1 = 20.0
#lambda1 = 5.0
#material_implementation = "analytical_Example"
-material_implementation = "matrix_material"
#material_implementation = "isotropic_bilayer"
+#material_implementation = "matrix_material_circles"
+material_implementation = "matrix_material_squares"
+
+#material_implementation = "circle_fiber" #TEST
+#material_implementation = "square_fiber" #TEST
#############################################
# Prestrain parameters
@@ -86,10 +93,11 @@ alpha = 5.0 # ratio between prestrain parameters rho1 & rho2
------------Matrix Material --------------
-prestrainType = "matrix_material"
+#prestrainType = "matrix_material_circles"
+prestrainType = "matrix_material_squares"
-nF = 5 #number of Fibers (in each Layer)
-rF = 0.1 #Fiber radius
+nF = 8 #number of Fibers (in each Layer)
+#rF = 0.05 #Fiber radius max-fiber-radius = (width/(2.0*nF)
------------------------------------------
width = 1.0
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 510c1b94..8b4e0a04 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -47,6 +47,9 @@
#include <dune/microstructure/prestrain_material_geometry.hh>
#include <dune/microstructure/matrix_operations.hh>
+#include <dune/microstructure/vtk_filler.hh> //TEST
+
+
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
// #include <dune/fufem/discretizationerror.hh>
@@ -986,6 +989,8 @@ int main(int argc, char *argv[])
// Print Options
bool print_debug = parameterSet.get<bool>("print_debug", false);
+
+ bool write_materialFunctions = parameterSet.get<bool>("write_materialFunctions", false);
bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
bool write_corrector_phi2 = parameterSet.get<bool>("write_corrector_phi2", false);
bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
@@ -1543,9 +1548,63 @@ int main(int argc, char *argv[])
vtkWriter.write( VTKOutputName + "-level"+ std::to_string(level));
std::cout << "wrote data to file: " + VTKOutputName + "-level" + std::to_string(level) << std::endl;
+
+
+ if (write_materialFunctions )
+ {
+
+ using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+ VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
+ using GridViewVTK = VTKGridType::LeafGridView;
+ const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
+
+ auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
+ auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
+
+ std::vector<double> mu_CoeffP0;
+ Functions::interpolate(scalarP0FeBasis, mu_CoeffP0, muTerm);
+ auto mu_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, mu_CoeffP0);
+
+ std::vector<double> mu_CoeffP1;
+ Functions::interpolate(scalarP1FeBasis, mu_CoeffP1, muTerm);
+ auto mu_DGBF_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, mu_CoeffP1);
+
+
+ std::vector<double> lambda_CoeffP0;
+ Functions::interpolate(scalarP0FeBasis, lambda_CoeffP0, lambdaTerm);
+ auto lambda_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, lambda_CoeffP0);
+
+ std::vector<double> lambda_CoeffP1;
+ Functions::interpolate(scalarP1FeBasis, lambda_CoeffP1, lambdaTerm);
+ auto lambda_DGBF_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, lambda_CoeffP1);
+
+ VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
+
+ MaterialVtkWriter.addVertexData(
+ mu_DGBF_P1,
+ VTK::FieldInfo("mu_P1", VTK::FieldInfo::Type::scalar, 1));
+ MaterialVtkWriter.addCellData(
+ mu_DGBF_P0,
+ VTK::FieldInfo("mu_P0", VTK::FieldInfo::Type::scalar, 1));
+ MaterialVtkWriter.addVertexData(
+ lambda_DGBF_P1,
+ VTK::FieldInfo("lambda_P1", VTK::FieldInfo::Type::scalar, 1));
+ MaterialVtkWriter.addCellData(
+ lambda_DGBF_P0,
+ VTK::FieldInfo("lambda_P0", VTK::FieldInfo::Type::scalar, 1));
+
+ MaterialVtkWriter.write("MaterialFunctions-level"+ std::to_string(level) );
+ std::cout << "wrote data to file: MaterialFunctions-level" + std::to_string(level) << std::endl;
+ }
+
+
levelCounter++;
} // Level-Loop End
+
+
+
+
/////////////////////////////////////////
// Print Storage
--
GitLab