diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index cc567f41d32677272ac27ae675717443606540b5..06b3cff70acd4f22aaf532ebe3823f2393db1943 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -15,9 +15,15 @@ cellDomain = 1
#############################################
#levels = 2
-numLevels = {1,3} # computes all levels from first to second entry
-#
+#######################################################################
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+########################################################################
+
+numLevels = 1 3 # computes all levels from first to second entry
+
+
#Elements_Cell = 20 20 20 # number elements in each direction (y1 y2 x3)
#nElements_Cell = 30 30 30
#nElements_Cell = 30 30 30
@@ -53,7 +59,8 @@ lambda1 = 0.0
# Prestrain parameters
#############################################
-theta = 0.3 # volume fraction
+#theta = 0.3 # volume fraction
+theta = 0.25 # volume fraction
prestrainType = 2
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index b6710dcd7b301e11f6fae0e7506c637a8d94fc35..0ed65a791245516350dd113ba4d198f54c680b17 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -2237,7 +2237,7 @@ int main(int argc, char *argv[])
- std::array<int, dim> nElements = parameterSet.get<std::array<int,dim>>("nElements_Cell", {10, 10, 10});
+// std::array<int, dim> nElements = parameterSet.get<std::array<int,dim>>("nElements_Cell", {10, 10, 10});
///// LEVEL - APPROACH // TODO
@@ -2291,8 +2291,20 @@ int main(int argc, char *argv[])
std::cout << "output variant :" << std::get<std::string>(x) << std::endl;
- std::vector<std::variant<std::string, size_t , double>> Storage;
+
+
+
+ ///////////////////////////////////
+ // Create Data Storage
+ ///////////////////////////////////
+
+ // Storage:: #1 level #2 L2SymError #3 L2SymErrorOrder #4 L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
+ std::vector<std::variant<std::string, size_t , double>> Storage_Error;
+
+ // Storage:: #1 level #2 |q1_a-q1_c| #3 |q2_a-q2_c| #4 |q3_a-q3_c| #5 |b1_a-b1_c| #6 |b2_a-b2_c| #7 |b3_a-b3_c|
+ std::vector<std::variant<std::string, size_t , double>> Storage_Quantities; // TODO :better error ?
+
// Storage:: #1 level #2 L2SymError #3 L2SymErrorOrder #4 L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
@@ -2333,7 +2345,8 @@ int main(int argc, char *argv[])
std::cout << "Level: " << level << std::endl;
std::cout << " ----------------------------------" << std::endl;
- Storage.push_back(level);
+ Storage_Error.push_back(level);
+ Storage_Quantities.push_back(level);
std::array<int, dim> nElements = { (int)std::pow(2,level) , (int)std::pow(2,level) , (int)std::pow(2,level) };
@@ -2820,7 +2833,7 @@ int main(int argc, char *argv[])
double GGterm = 0.0;
GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
// TEST
- //TEST
+ //TEST
std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
std::setprecision(std::numeric_limits<float>::digits10);
@@ -2829,6 +2842,7 @@ int main(int argc, char *argv[])
// Q[a][b] = (coeffContainer[a]*loadContainer[b]) + GGterm; // TEST
std::cout << "GGTerm:" << GGterm << std::endl;
std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
+
// std::cout << "Test : " << coeffContainer[a]*loadContainer[b] << std::endl; //same...
}
printmatrix(std::cout, Q, "Matrix Q", "--");
@@ -2843,13 +2857,14 @@ int main(int argc, char *argv[])
auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B_Term); // <L i(x3G_alpha) , B>
B_hat[a] = (coeffContainer[a]*tmp2) + GGterm;
- std::cout << "check this Contribution: " << (coeffContainer[a]*tmp2) << std::endl; //see orthotropic.tex
- std::cout <<"B_hat[i]: " << B_hat[a] << std::endl;
+
+// std::cout << "check this Contribution: " << (coeffContainer[a]*tmp2) << std::endl; //see orthotropic.tex
+// std::cout <<"B_hat[i]: " << B_hat[a] << std::endl;
}
log << "Computed prestrain B_hat: " << std::endl;
log << B_hat << std::endl;
-
+ std::cout << "check this Contribution: " << (coeffContainer[2]*tmp2) << std::endl; //see orthotropic.tex
// Compute effective Prestrain
FieldVector<double, 3> Beff;
@@ -2858,11 +2873,11 @@ int main(int argc, char *argv[])
log << "Computed prestrain B_eff: " << std::endl;
log << Beff << std::endl;
- std::cout << "Beff : " << std::endl;
- for(size_t i=0; i<3; i++)
- {
- std::cout <<"Beff[i]: " << Beff[i] << std::endl;
- }
+// std::cout << "Beff : " << std::endl;
+// for(size_t i=0; i<3; i++)
+// {
+// std::cout <<"Beff[i]: " << Beff[i] << std::endl;
+// }
@@ -2873,6 +2888,15 @@ int main(int argc, char *argv[])
std::cout<< "q1_c: " << q1_c << std::endl;
std::cout<< "q2_c: " << q2_c << std::endl;
std::cout<< "q3_c: " << q3_c << std::endl;
+ std::cout<< "b1hat_c: " << B_hat[0] << std::endl;
+ std::cout<< "b2hat_c: " << B_hat[1] << std::endl;
+ std::cout<< "b3hat_c: " << B_hat[2] << std::endl;
+ std::cout<< "b1eff_c: " << Beff[0] << std::endl;
+ std::cout<< "b2eff_c: " << Beff[1] << std::endl;
+ std::cout<< "b3eff_c: " << Beff[2] << std::endl;
+
+
+ std::cout << "Theta: " << theta << std::endl;
std::cout << "Gamma: " << gamma << std::endl;
std::cout << "Number of Elements in each direction: " << nElements << std::endl;
log << "computed q1: " << q1_c << std::endl;
@@ -2894,35 +2918,56 @@ int main(int argc, char *argv[])
// double b1 = (mu1*p1/4)*(beta/(theta+(1-theta)*beta))*(1-theta*(1+alpha));
// double b2 = (mu1*p1/8)*(1-theta*(1+beta*alpha));
- double b1 = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
- double b2 = -(theta/4.0)*mu2;
- double b3 = 0.0;
+ double b1_hat = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
+ double b2_hat = -(theta/4.0)*mu2;
+ double b3_hat = 0.0;
+
+ double b1_eff = (-3.0/2.0)*theta;
+ double b2_eff = (-3.0*theta*mu2)/(mu1*(1-theta)+mu2*theta);
+ double b3_eff = 0.0;
+
+
+
+
double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
std::cout << " --- print analytic solutions --- " << std::endl;
- std::cout << "b1 : " << b1 << std::endl;
- std::cout << "b2 : " << b2 << std::endl;
- std::cout << "b3 : " << b3 << std::endl;
+ std::cout << "b1_hat : " << b1_hat << std::endl;
+ std::cout << "b2_hat : " << b2_hat << std::endl;
+ std::cout << "b3_hat : " << b3_hat << std::endl;
+ std::cout << "b1_eff : " << b1_eff << std::endl;
+ std::cout << "b2_eff : " << b2_eff << std::endl;
+ std::cout << "b3_eff : " << b3_eff << std::endl;
+
std::cout << "q1 : " << q1 << std::endl;
std::cout << "q2 : " << q2 << std::endl;
std::cout << "q3 should be between q1 and q2" << std::endl;
log << " --- analytic solutions: --- " << std::endl;
- log << "b1 : " << b1 << std::endl;
- log << "b2 : " << b2 << std::endl;
- log << "b3 : " << b3 << std::endl;
+ log << "b1_hat : " << b1_hat << std::endl;
+ log << "b2_hat : " << b2_hat << std::endl;
+ log << "b3_hat : " << b3_hat << std::endl;
+ log << "b1_eff : " << b1_eff << std::endl;
+ log << "b2_eff : " << b2_eff << std::endl;
+ log << "b3_eff : " << b3_eff << std::endl;
log << "q1 : " << q1 << std::endl;
log << "q2 : " << q2 << std::endl;
log << "q3 should be between q1 and q2" << std::endl;
+
+ Storage_Quantities.push_back(std::abs(q1 - q1_c));
+ Storage_Quantities.push_back(std::abs(q2 - q2_c));
+ Storage_Quantities.push_back(std::abs(q3 - q3_c));
+ Storage_Quantities.push_back(std::abs(b1_eff - b1eff_c));
+ Storage_Quantities.push_back(std::abs(b2_eff - b2eff_c));
+ Storage_Quantities.push_back(std::abs(b3_eff - b3eff_c));
+
+ auto symPhi_1_analytic = [mu1, mu2, theta, muTerm] (const Domain& x) {
+ return MatrixRT{{ (((mu1*mu2)/((theta*mu1 +(1-theta)*mu2)*muTerm(x))) - 1)*x[2] , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
-
- auto symPhi_1_analytic = [mu1, mu2, theta, muTerm] (const Domain& x) {
- return MatrixRT{{ (((mu1*mu2)/((theta*mu1 +(1-theta)*mu2)*muTerm(x))) - 1)*x[2] , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
-
- auto zeroFunction = [] (const Domain& x) {
- return MatrixRT{{0.0 , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
+ auto zeroFunction = [] (const Domain& x) {
+ return MatrixRT{{0.0 , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
if(write_L2Error)
@@ -2967,40 +3012,40 @@ int main(int argc, char *argv[])
double EOC = 0.0;
- Storage.push_back(L2SymError);
+ Storage_Error.push_back(L2SymError);
//Compute Experimental order of convergence (EOC)
if(levelCounter > 0)
{
- // Storage:: #1 level #2 L2SymError #3 L2SymErrorOrder #4 L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
- // TODO Storage as Vectors? (q1,q2,q3) ,(b1,b2,b3) ....
+ // Storage_Error:: #1 level #2 L2SymError #3 L2SymErrorOrder #4 L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
+ // TODO Storage_Error as Vectors? (q1,q2,q3) ,(b1,b2,b3) ....
- // Storage.push_back(5);
- // Storage.push_back("Second Entry");
+ // Storage_Error.push_back(5);
+ // Storage_Error.push_back("Second Entry");
std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter-1)*6)+1 << std::endl; // Besser std::map ???
- std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter)*6)+1 << std::endl; // muss Storage[idx] muss idx zur compile time feststehen?!
+ std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter)*6)+1 << std::endl; // muss Storage_Error[idx] muss idx zur compile time feststehen?!
- auto ErrorOld = std::get<double>(Storage[((levelCounter-1)*6)+1]);
- auto ErrorNew = std::get<double>(Storage[(levelCounter*6)+1]);
+ auto ErrorOld = std::get<double>(Storage_Error[((levelCounter-1)*6)+1]);
+ auto ErrorNew = std::get<double>(Storage_Error[(levelCounter*6)+1]);
//
EOC = std::log(ErrorNew/ErrorOld)/(-1*std::log(2));
- // std::cout << "Storage[0] :" << std::get<1>(Storage[0]) << std::endl;
- // std::cout << "output variant :" << std::get<std::string>(Storage[1]) << std::endl;
- // std::cout << "output variant :" << std::get<0>(Storage[1]) << std::endl;
+ // std::cout << "Storage_Error[0] :" << std::get<1>(Storage_Error[0]) << std::endl;
+ // std::cout << "output variant :" << std::get<std::string>(Storage_Error[1]) << std::endl;
+ // std::cout << "output variant :" << std::get<0>(Storage_Error[1]) << std::endl;
}
-// Storage.push_back(level);
- Storage.push_back(EOC);
- Storage.push_back(L2Norm_Symphi);
- Storage.push_back(L2Norm_SymAnalytic);
- Storage.push_back(L2Norm);
-// Storage.push_back();
-// Storage.push_back();
-// Storage.push_back();
-// Storage.push_back();
+// Storage_Error.push_back(level);
+ Storage_Error.push_back(EOC);
+ Storage_Error.push_back(L2Norm_Symphi);
+ Storage_Error.push_back(L2Norm_SymAnalytic);
+ Storage_Error.push_back(L2Norm);
+// Storage_Error.push_back();
+// Storage_Error.push_back();
+// Storage_Error.push_back();
+// Storage_Error.push_back();
}
@@ -3096,7 +3141,7 @@ int main(int argc, char *argv[])
std::vector<std::variant<const size_t , double>> tmpVec;
int StorageCount = 0;
- for(auto& v: Storage) {
+ for(auto& v: Storage_Error) {
// if(StorageCount % 6 == 0 )
// std::cout<< "Level: ";