From a0d8b828a71118dd84500da2a9ce129018b88c1a Mon Sep 17 00:00:00 2001
From: Klaus Boehnlein <klaus.boehnlein@tu-dresden.de>
Date: Mon, 12 Jul 2021 09:11:02 +0200
Subject: [PATCH] minor changes
---
src/dune-microstructure.cc | 136 ++++++++++++++++++++++++++-----------
1 file changed, 97 insertions(+), 39 deletions(-)
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 95e6b07e..f7d609b7 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -1257,7 +1257,7 @@ int main(int argc, char *argv[])
// Get Parameters/Data
///////////////////////////////////
- double gamma = parameterSet.get<double>("gamma",5.0); // ratio dimension reduction to homogenization
+ double gamma = parameterSet.get<double>("gamma",10000); // ratio dimension reduction to homogenization
// Material parameter laminat
// double E1 = parameterSet.get<double>("E1", 17e6); //material1
@@ -1281,11 +1281,13 @@ int main(int argc, char *argv[])
unsigned int prestraintype = parameterSet.get<unsigned int>("imp", 1);
double p1 = parameterSet.get<double>("prestress_pressure1", 1.0);
double p2 = parameterSet.get<double>("prestress_pressure2", 2.0);
- double theta = parameterSet.get<double>("theta",0.5);
+ double theta = parameterSet.get<double>("theta",0.3);
// double theta = 0.5;
p1 = 1.0;
- p2 = 1.0;
+ double alpha = 2;
+ p2 = alpha*1.0;
prestraintype = 2;
+// prestraintype = 1;
auto prestrainImp = PrestrainImp(p1, p2, theta, width);
auto B_Term = prestrainImp.getPrestrain(prestraintype);
@@ -1337,12 +1339,12 @@ int main(int argc, char *argv[])
- double alpha = 2;
- double beta = 2;
+ double beta = 2;
// double mu1 = 10;
double mu1 = 0.5*17e6;
+// double mu1 = 1000;
double mu2 = beta*mu1;
// //
@@ -1359,8 +1361,13 @@ int main(int argc, char *argv[])
auto muTerm = [mu1, mu2, theta] (const Domain& z) {
if (abs(z[0]) >= (theta/2))
return mu1;
+// if (abs(z[0]) < (theta/2) && z[2] < 0) // oder das?
+// return mu2;
else
return mu2;
+// return mu1;
+
+
};
@@ -1424,25 +1431,22 @@ int main(int argc, char *argv[])
/////////////////////////////////////////////////////////
// Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
-
-
-
VectorCT load_alpha1, load_alpha2, load_alpha3;
-
+MatrixCT stiffnessMatrix_CE;
// auto rhsBackend = Functions::istlVectorBackend(b);
// rhsBackend.resize(Basis_CE);
// b = 0;
-MatrixCT stiffnessMatrix_CE;
+
// bool set_integral_zero = true;
-// bool set_oneBasisFunction_Zero = true;
-// bool substract_integralMean = true;
-bool set_oneBasisFunction_Zero = false;
-bool substract_integralMean = false;
+bool set_oneBasisFunction_Zero = true;
+bool substract_integralMean = true;
+// bool set_oneBasisFunction_Zero = false;
+// bool substract_integralMean = false;
@@ -1502,19 +1506,16 @@ std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
-
/////////////////////////////////////////////////////////
// Assemble the system
/////////////////////////////////////////////////////////
-
-
assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE);
-
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3);
+
// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");
@@ -1710,8 +1711,6 @@ for(size_t i=0; i<3; i++)
////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
////////////////////////////////////////////////////////////////////////////
-
-
using SolutionRange = FieldVector<double,dim>;
// using SolutionRange = double;
@@ -1719,8 +1718,6 @@ using SolutionRange = FieldVector<double,dim>;
// Basis_CE,
// Dune::Functions::HierarchicVectorWrapper<VectorCT, double>(phi_1));
-
-
auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
phi_1);
@@ -1762,27 +1759,27 @@ std::cout << "check integralMean: " << std::endl;
auto A = integralMean(Basis_CE,localPhi_1);
for(size_t i=0; i<3; i++)
{
-std::cout << "Integral-Mean (TEST) : " << A[i] << std::endl;
+std::cout << "Integral-Mean : " << A[i] << std::endl;
}
if(substract_integralMean)
{
-std::cout << "substracting integralMean " << std::endl;
+std::cout << " --- substracting integralMean --- " << std::endl;
subtractIntegralMean(Basis_CE, localPhi_1 , phi_1);
subtractIntegralMean(Basis_CE, localPhi_2 , phi_2);
subtractIntegralMean(Basis_CE, localPhi_3 , phi_3);
////////////////////////////////////////////////////////////////////////
// CHECK INTEGRAL-MEAN:
-// auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
-// Basis_CE,
-// phi_1);
-//
-// auto AVPhi_1 = localFunction(AVFunction_1);
-// auto A = integralMean(Basis_CE, AVPhi_1);
-// for(size_t i=0; i<3; i++)
-// {
-// std::cout << "Integral-Mean (TEST) : " << A[i] << std::endl;
-// }
+auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+ Basis_CE,
+ phi_1);
+
+auto AVPhi_1 = localFunction(AVFunction_1);
+auto A = integralMean(Basis_CE, AVPhi_1);
+for(size_t i=0; i<3; i++)
+{
+std::cout << "Integral-Mean (CHECK) : " << A[i] << std::endl;
+}
////////////////////////////////////////////////////
}
@@ -1868,11 +1865,72 @@ for(size_t a = 0; a < 3; a++)
for(size_t i=0; i<3; i++)
{
- std::cout << B_hat[i] << std::endl;
+ std::cout <<"B_hat[i]: " << B_hat[i] << std::endl;
}
+// compute effective Prestrain
+FieldVector<double, 3> Beff;
+
+
+Q.solve(Beff,B_hat);
+
+std::cout << "Beff : " << std::endl;
+for(size_t i=0; i<3; i++)
+{
+ std::cout <<"Beff[i]: " << Beff[i] << std::endl;
+}
+
+
+
+// compute q1,q2,q3
+
+FieldVector<double ,3> g1 = {1.0 , 0 , 0};
+FieldVector<double ,3> g2 = {0 , 1.0 , 0};
+FieldVector<double ,3> g3 = {0 , 0, 1.0};
+FieldVector<double ,3> tmp_1;
+FieldVector<double ,3> tmp_2;
+FieldVector<double ,3> tmp_3;
+// auto tmp = g1- B_hat;
+// std::cout<< tmp << std::endl;
+
+Q.mv(g1,tmp_1);
+auto q1_c = tmp_1 * g1;
+std::cout<< "q1_c: " << q1_c << std::endl;
+
+Q.mv(g2,tmp_2);
+auto q2_c = tmp_2 * g2;
+std::cout<< "q2_c: " << q2_c << std::endl;
+
+Q.mv(g3,tmp_3);
+auto q3_c = tmp_3 * g3;
+std::cout<< "q3_c: " << q3_c << std::endl;
+
+
+std::cout << "Gamma: " << gamma << std::endl;
+
+
+// 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 = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
+double b2 = -(theta/4.0)*mu2;
+double b3 = 0.0;
+double q1 = (mu1/6.0)*mu2/(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 << "q1 : " << q1 << std::endl;
+std::cout << "q2 : " << q2 << std::endl;
+std::cout << "q3 should be between q1 and q2" << std::endl;
+
+
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
@@ -1888,11 +1946,11 @@ for(size_t i=0; i<3; i++)
correctorFunction_1,
VTK::FieldInfo("corrector phi_1", VTK::FieldInfo::Type::vector, dim));
vtkWriter.addVertexData(
- correctorFunction_1,
- VTK::FieldInfo("corrector phi_1", VTK::FieldInfo::Type::vector, dim));
+ correctorFunction_2,
+ VTK::FieldInfo("corrector phi_2", VTK::FieldInfo::Type::vector, dim));
vtkWriter.addVertexData(
- correctorFunction_1,
- VTK::FieldInfo("corrector phi_1", VTK::FieldInfo::Type::vector, dim));
+ correctorFunction_3,
+ VTK::FieldInfo("corrector phi_3", VTK::FieldInfo::Type::vector, dim));
// vtkWriter.addVertexData(
// solutionFunction,
// VTK::FieldInfo("solution", VTK::FieldInfo::Type::scalar, 1));
--
GitLab