From 7d499416292aa0a38a20d5a65af08055812438f8 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Wed, 23 Feb 2022 23:10:53 +0100
Subject: [PATCH] Check Cell-Solver
---
.../prestrain_material_geometry.hh | 3 +-
src/Cell-Problem.cc | 106 +++++++++++++++---
2 files changed, 92 insertions(+), 17 deletions(-)
diff --git a/dune/microstructure/prestrain_material_geometry.hh b/dune/microstructure/prestrain_material_geometry.hh
index 4522baef..c023fe6a 100644
--- a/dune/microstructure/prestrain_material_geometry.hh
+++ b/dune/microstructure/prestrain_material_geometry.hh
@@ -457,7 +457,8 @@ public:
auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
-// std::cout << "Analytical-LAMBDA is used" << std::endl; //TEST
+// 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;
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index b64396b1..d05eb93b 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -861,7 +861,7 @@ int main(int argc, char *argv[])
// Output setter
// std::string outputPath = parameterSet.get("outputPath", "../../outputs/output.txt");
// std::string outputPath = parameterSet.get("outputPath", "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt");
- std::string outputPath = parameterSet.get("outputPath", " /home/klaus/Desktop/DUNE/dune-microstructure/outputs");
+ std::string outputPath = parameterSet.get("outputPath", "/home/klaus/Desktop/DUNE/dune-microstructure/outputs");
// std::string MatlabPath = parameterSet.get("MatlabPath", "/home/klaus/Desktop/DUNE/dune-microstructure/Matlab-Programs");
// std::string outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt";
std::fstream log;
@@ -887,7 +887,7 @@ int main(int argc, char *argv[])
std::string imp = parameterSet.get<std::string>("material_prestrain_imp", "analytical_Example");
log << "material_prestrain used: "<< imp << std::endl;
double beta = parameterSet.get<double>("beta",2.0);
- double mu1 = parameterSet.get<double>("mu1",10.0);;
+ double mu1 = parameterSet.get<double>("mu1",1.0);;
double mu2 = beta*mu1;
double lambda1 = parameterSet.get<double>("lambda1",0.0);;
double lambda2 = beta*lambda1;
@@ -1447,6 +1447,7 @@ int main(int argc, char *argv[])
std::cout << "q3 : " << std::fixed << q3 << std::endl;
// std::cout<< "q3 : " << q3 << std::endl;
std::cout<< std::fixed << std::setprecision(6) << "q_onetwo=" << Q[0][1] << std::endl;
+// std::cout<< std::fixed << std::setprecision(6) << "q_onetwo=" << Q[1][0] << std::endl; //TEST
printvector(std::cout, B_hat, "B_hat", "--");
printvector(std::cout, Beff, "Beff", "--");
@@ -1472,20 +1473,7 @@ int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////
// Define Analytic Solutions
//////////////////////////////////////////////////////////////
- // 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_hat_ana = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
- double b2_hat_ana = -(theta/4.0)*mu2;
- double b3_hat_ana = 0.0;
-
- double b1_eff_ana = (-3.0/2.0)*theta;
- double b2_eff_ana = (-3.0/2.0)*(theta*mu2)/(mu1*(1-theta)+mu2*theta);
- double b3_eff_ana = 0.0;
-
-// double q1_ana = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
-// double q2_ana = ((1.0-theta)*mu1+theta*mu2)/6.0;
- double q1_ana = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0); // 1/6 * harmonicMean
- double q2_ana = mu1*((1-theta)+theta*beta) *(1.0/6.0); // 1/6 * arithmeticMean
+
@@ -1495,12 +1483,76 @@ int main(int argc, char *argv[])
double p1 = parameterSet.get<double>("rho1", 1.0);
double alpha = parameterSet.get<double>("alpha", 2.0);
double p2 = alpha*p1;
+
+ if (imp == "parametrized_Laminate") // print Errors only for analytical_Example
+ {
+
+ double rho1 = parameterSet.get<double>("rho1", 1.0);
+// double b1_hat_ana = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
+// double b2_hat_ana = -(theta/4.0)*mu2;
+// double b3_hat_ana = 0.0;
+
+ double b1_eff_ana = (3.0/2.0)*rho1*(1-theta*(1+alpha));
+ double b2_eff_ana = (3.0/2.0)*rho1*((1-theta*(1+beta*alpha))/(1-theta+theta*beta));
+ double b3_eff_ana = 0.0;
+
+ // double q1_ana = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
+ // double q2_ana = ((1.0-theta)*mu1+theta*mu2)/6.0;
+ double q1_ana = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0); // 1/6 * harmonicMean
+ double q2_ana = mu1*((1-theta)+theta*beta) *(1.0/6.0); // 1/6 * arithmeticMean
+
+ std::cout << "----- print analytic solutions -----" << std::endl;
+// std::cout << "b1_hat_ana : " << b1_hat_ana << std::endl;
+// std::cout << "b2_hat_ana : " << b2_hat_ana << std::endl;
+// std::cout << "b3_hat_ana : " << b3_hat_ana << std::endl;
+ std::cout << "b1_eff_ana : " << b1_eff_ana << std::endl;
+ std::cout << "b2_eff_ana : " << b2_eff_ana << std::endl;
+ std::cout << "b3_eff_ana : " << b3_eff_ana << std::endl;
+
+ std::cout << "q1_ana : " << q1_ana << std::endl;
+ std::cout << "q2_ana : " << q2_ana << std::endl;
+ std::cout << "q3 should be between q1 and q2" << std::endl;
+ log << "----- print analytic solutions -----" << std::endl;
+// log << "b1_hat_ana : " << b1_hat_ana << std::endl;
+// log << "b2_hat_ana : " << b2_hat_ana << std::endl;
+// log << "b3_hat_ana : " << b3_hat_ana << std::endl;
+ log << "b1_eff_ana : " << b1_eff_ana << std::endl;
+ log << "b2_eff_ana : " << b2_eff_ana << std::endl;
+ log << "b3_eff_ana : " << b3_eff_ana << std::endl;
+ log << "q1_ana : " << q1_ana << std::endl;
+ log << "q2_ana : " << q2_ana << std::endl;
+ log << "q3 should be between q1 and q2" << std::endl;
+
+ Storage_Quantities.push_back(std::abs(q1_ana - q1));
+ Storage_Quantities.push_back(std::abs(q2_ana - q2));
+ Storage_Quantities.push_back(q3);
+ Storage_Quantities.push_back(std::abs(b1_eff_ana - Beff[0]));
+ Storage_Quantities.push_back(std::abs(b2_eff_ana - Beff[1]));
+ Storage_Quantities.push_back(std::abs(b3_eff_ana - Beff[2]));
+ }
+
+
if (imp == "analytical_Example") // print Errors only for analytical_Example
{
std::cout << ((3.0*p1)/2.0)*beta*(1-(theta*(1+alpha))) << std::endl; // TODO ERROR in paper ??
+
+ // 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_hat_ana = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
+ double b2_hat_ana = -(theta/4.0)*mu2;
+ double b3_hat_ana = 0.0;
+
+ double b1_eff_ana = (-3.0/2.0)*theta;
+ double b2_eff_ana = (-3.0/2.0)*(theta*mu2)/(mu1*(1-theta)+mu2*theta);
+ double b3_eff_ana = 0.0;
+
+ // double q1_ana = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
+ // double q2_ana = ((1.0-theta)*mu1+theta*mu2)/6.0;
+ double q1_ana = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0); // 1/6 * harmonicMean
+ double q2_ana = mu1*((1-theta)+theta*beta) *(1.0/6.0); // 1/6 * arithmeticMean
std::cout << "----- print analytic solutions -----" << std::endl;
std::cout << "b1_hat_ana : " << b1_hat_ana << std::endl;
@@ -1700,6 +1752,28 @@ int main(int argc, char *argv[])
// WORKS Too
VTKFiller.vtkProblemCell(gridView_VTK, B_Term, muLocal,"VTKProblemCell");;
+
+
+ // TEST
+ auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
+ auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
+
+ std::vector<double> B_CoeffP0;
+ Functions::interpolate(scalarP0FeBasis, B_CoeffP0, B_Term);
+ auto B_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, B_CoeffP0);
+
+
+
+ VTKWriter<GridView> PrestrainVtkWriter(gridView_VTK);
+
+ PrestrainVtkWriter.addCellData(
+ B_DGBF_P0,
+ VTK::FieldInfo("B_P0", VTK::FieldInfo::Type::scalar, 1));
+
+ PrestrainVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
+ std::cout << "wrote data to file:" + outputPath +"/PrestrainFunctions-level" + std::to_string(level) << std::endl;
+
+
}
--
GitLab