diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 2db0ee0d841185ceb6e75f8b7b7f4f9e15dd6878..95e6b07e2ee7c8dadb144823dd1258eb7e5abf04 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -54,6 +54,9 @@
+// #include <dune/fufem/discretizationerror.hh>
+
+
using namespace Dune;
@@ -103,12 +106,12 @@ void computeElementStiffnessMatrix(const LocalView& localView,
//////////////////////////////////////////////
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
- MatrixRT G_3 {{0.0, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
//print:
- printmatrix(std::cout, basisContainer[0] , "G_1", "--");
- printmatrix(std::cout, basisContainer[1] , "G_2", "--");
- printmatrix(std::cout, basisContainer[2] , "G_3", "--");
+// printmatrix(std::cout, basisContainer[0] , "G_1", "--");
+// printmatrix(std::cout, basisContainer[1] , "G_2", "--");
+// printmatrix(std::cout, basisContainer[2] , "G_3", "--");
////////////////////////////////////////////////////
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
@@ -497,7 +500,7 @@ void computeElementLoadVector( const LocalView& localView,
//////////////////////////////////////////////
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
- MatrixRT G_3 {{0.0, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
@@ -637,7 +640,7 @@ void assembleCellStiffness(const Basis& basis,
// Dune::Matrix<double> elementMatrix;
Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
- printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
+// printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
@@ -697,7 +700,7 @@ void assembleCellLoad(const Basis& basis,
)
{
- std::cout << "assemble load vector." << std::endl;
+// std::cout << "assemble load vector." << std::endl;
b.resize(basis.size()+3);
b = 0;
@@ -866,6 +869,8 @@ void setOneBaseFunctionToZero(const Basis& basis,
)
{
+ std::cout << "setting one Basis-function to zero" << std::endl;
+
constexpr int dim = 3;
auto localView = basis.localView();
@@ -884,7 +889,7 @@ void setOneBaseFunctionToZero(const Basis& basis,
{
auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
row[k] = localView.index(rowLocal);
- std::cout << "row[k]:" << row[k] << std::endl;
+// std::cout << "row[k]:" << row[k] << std::endl;
load_alpha1[row[k]] = 0.0; //verwende hier indexVector
load_alpha2[row[k]] = 0.0;
load_alpha3[row[k]] = 0.0;
@@ -957,7 +962,7 @@ void setIntegralZero (const Basis& basis,
{
auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
row[k] = localView.index(rowLocal);
- std::cout << "row[k]:" << row[k] << std::endl;
+// std::cout << "row[k]:" << row[k] << std::endl;
stiffnessMatrix[row[k]] = 0;
load_alpha1[row[k]] = 0.0;
load_alpha2[row[k]] = 0.0;
@@ -1202,19 +1207,31 @@ auto subtractIntegralMean(const Basis& basis,
-
-
-
// Check whether two points are equal on R/Z x R/Z x R
auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
{
return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
and (FloatCmp::eq(x[2],y[2])));
-// and (FloatCmp::eq(x[2],y[2]) or FloatCmp::eq(x[2]+1,y[2]) or FloatCmp::eq(x[2]-1,y[2])));
};
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
int main(int argc, char *argv[])
{
@@ -1240,18 +1257,19 @@ int main(int argc, char *argv[])
// Get Parameters/Data
///////////////////////////////////
- double gamma = parameterSet.get<double>("gamma",1.0); // ratio dimension reeduction to homogenization
+ double gamma = parameterSet.get<double>("gamma",5.0); // ratio dimension reduction to homogenization
// Material parameter laminat
- double E1 = parameterSet.get<double>("E1", 17e6); //material1
- double nu1 = parameterSet.get<double>("nu1", 0.3);
- double E2 = parameterSet.get<double>("E2", 35e6); //material2
- double nu2 = parameterSet.get<double>("nu2", 0.5);
+// double E1 = parameterSet.get<double>("E1", 17e6); //material1
+// double nu1 = parameterSet.get<double>("nu1", 0.0);
+// double nu1 = parameterSet.get<double>("nu1", 0.3);
+// double E2 = parameterSet.get<double>("E2", 35e6); //material2
+// double nu2 = parameterSet.get<double>("nu2", 0.5);
// // Plate geometry settings
-// double width = parameterSet.get<double>("width", 1.0); //geometry cell, cross section
+ double width = parameterSet.get<double>("width", 1.0); //geometry cell, cross section
// double len = parameterSet.get<double>("len", 10.0); //length
// double height = parameterSet.get<double>("h", 0.1); //rod thickness
// double eps = parameterSet.get<double>("eps", 0.1); //size of perioticity cell
@@ -1260,13 +1278,17 @@ int main(int argc, char *argv[])
///////////////////////////////////
// ---Get Prestrain ---
///////////////////////////////////
-// unsigned int prestraintype = parameterSet.get<unsigned int>("imp");
-// double p1 = parameterSet.get<double>("prestress_pressure1", 2.0); //0 x2 1 x3 2 y
-// double p2 = parameterSet.get<double>("prestress_pressure2", 0.5);
-// double a = parameterSet.get<double>("prestrainCoff",10.0);
-//
-// auto prestrainImp = PrestrainImp(p1, p2, a,width);
-// auto B_Term = prestrainImp.getPrestrain(prestraintype);
+ 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 = 0.5;
+ p1 = 1.0;
+ p2 = 1.0;
+ prestraintype = 2;
+
+ auto prestrainImp = PrestrainImp(p1, p2, theta, width);
+ auto B_Term = prestrainImp.getPrestrain(prestraintype);
@@ -1275,14 +1297,19 @@ int main(int argc, char *argv[])
///////////////////////////////////
using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
- FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
- FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
+ // Domain : [0,1)^2 x (-1/2, 1/2)
+// FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
+// FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
+
+ // (shifted) Domain (-1/2,1/2)^3
+ FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0});
+ FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0});
- int nE = 2;
+ int nE = 10;
std::array<int,dim> nElements={nE,nE,nE}; //#Elements in each direction
- CellGridType grid_CE(lower,upper,nElements); //Corrector problem Domain: [0,1)^2 x (-1/2, 1/2)
+ CellGridType grid_CE(lower,upper,nElements); //Corrector problem Domain:
using GridView = CellGridType::LeafGridView;
const GridView gridView_CE = grid_CE.leafGridView();
@@ -1309,18 +1336,45 @@ int main(int argc, char *argv[])
- double mu1 = 0.5 * 1.0/(1.0 + nu1) * E1;
- double lambda1 = nu1/(1.0-2.0*nu1) * 1.0/(1.0+nu1) * E1;
+
+ double alpha = 2;
+ double beta = 2;
+
+
+// double mu1 = 10;
+ double mu1 = 0.5*17e6;
+ double mu2 = beta*mu1;
+// //
+
+// double mu1 = 0.5 * 1.0/(1.0 + nu1) * E1;
+// double lambda1 = nu1/(1.0-2.0*nu1) * 1.0/(1.0+nu1) * E1;
+ double lambda1 = 0;
+
// Create Lambda-Functions for material Parameters
- auto muTerm = [mu1] (const Domain& z) {
- return mu1;
+// auto muTerm = [mu1] (const Domain& z) {
+// return mu1;
+// };
+
+ auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+ if (abs(z[0]) >= (theta/2))
+ return mu1;
+ else
+ return mu2;
};
+
+
+
+
auto lambdaTerm = [lambda1] (const Domain& z) {
return lambda1;
};
+
+
+
+
auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
auto muLocal = localFunction(muGridF);
auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
@@ -1356,14 +1410,14 @@ int main(int argc, char *argv[])
const int phiOffset = Basis_CE.size();
// TEST
- auto X1Basis = subspaceBasis(Basis_CE, 0);
- auto X2Basis = subspaceBasis(Basis_CE, 1);
- auto X3Basis = subspaceBasis(Basis_CE, 2);
-
- std::cout << "size X1Basis: " << X1Basis.size() << std::endl;
- std::cout << "X1Basis.dimension()" << X1Basis.dimension() <<std::endl;
- std::cout << "size X2Basis: " << X2Basis.size() << std::endl;
- std::cout << "size X3Basis: " << X3Basis.size() << std::endl;
+// auto X1Basis = subspaceBasis(Basis_CE, 0);
+// auto X2Basis = subspaceBasis(Basis_CE, 1);
+// auto X3Basis = subspaceBasis(Basis_CE, 2);
+//
+// std::cout << "size X1Basis: " << X1Basis.size() << std::endl;
+// std::cout << "X1Basis.dimension()" << X1Basis.dimension() <<std::endl;
+// std::cout << "size X2Basis: " << X2Basis.size() << std::endl;
+// std::cout << "size X3Basis: " << X3Basis.size() << std::endl;
@@ -1384,8 +1438,11 @@ MatrixCT stiffnessMatrix_CE;
-bool set_integral_zero = true;
-bool set_oneBasisFunction_Zero = true;
+// bool set_integral_zero = true;
+// bool set_oneBasisFunction_Zero = true;
+// bool substract_integralMean = true;
+bool set_oneBasisFunction_Zero = false;
+bool substract_integralMean = false;
@@ -1399,15 +1456,15 @@ Func2Tensor x3G_2 = [] (const Domain& x) {
return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*x[2], 0.0}, {0.0, 0.0, 0.0}}; };
Func2Tensor x3G_3 = [] (const Domain& x) {
- return MatrixRT{{0.0, 0.5*x[2], 0.0}, {0.5*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}}; };
+ return MatrixRT{{0.0, 1/sqrt(2)*x[2], 0.0}, {1/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}}; };
/////////////////////////////////////////////////////////////////////// TODO
// TODO : PrestrainImp.hh
-double theta = 0.5;
-double p1 = 1;
-double p2 = 2;
+// double theta = 0.5;
+// double p1 = 1;
+// double p2 = 2;
using std::abs;
@@ -1434,7 +1491,7 @@ Func2Tensor B = [] (const Domain& x) {
//////////////////////////////////////////////
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
-MatrixRT G_3 {{0.0, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
@@ -1458,7 +1515,7 @@ 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", "--");
+// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix before B.C", "--");
// set Integral to zero
@@ -1497,7 +1554,7 @@ FieldVector<int,3> row; //fill with Indices..
if(set_oneBasisFunction_Zero)
{
setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3,row);
- printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
+// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
}
@@ -1689,10 +1746,10 @@ for(size_t i=0; i<3; i++)
M_3 += m_3[i]*basisContainer[i];
}
-
+std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
printmatrix(std::cout, M_1, "Corrector-Matrix M_1", "--");
-printmatrix(std::cout, M_2, "Corrector-Matrix M_1", "--");
-printmatrix(std::cout, M_3, "Corrector-Matrix M_1", "--");
+printmatrix(std::cout, M_2, "Corrector-Matrix M_2", "--");
+printmatrix(std::cout, M_3, "Corrector-Matrix M_3", "--");
@@ -1701,10 +1758,19 @@ auto localPhi_1 = localFunction(correctorFunction_1);
auto localPhi_2 = localFunction(correctorFunction_2);
auto localPhi_3 = localFunction(correctorFunction_3);
+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;
+}
+
+if(substract_integralMean)
+{
+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>(
@@ -1718,6 +1784,30 @@ subtractIntegralMean(Basis_CE, localPhi_3 , phi_3);
// std::cout << "Integral-Mean (TEST) : " << A[i] << std::endl;
// }
////////////////////////////////////////////////////
+}
+
+
+
+
+// // PRINT COEFFICIENTS
+// std::cout << "print coefficient-vector phi_1" << std::endl;
+// for(size_t i=0; i<Basis_CE.size();i++)
+// {
+// std::cout << phi_1[i] << std::endl;
+// }
+// std::cout << "print coefficient-vector phi_2" << std::endl;
+// for(size_t i=0; i<Basis_CE.size();i++)
+// {
+// std::cout << phi_2[i] << std::endl;
+// }
+// std::cout << "print coefficient-vector phi_3" << std::endl;
+// for(size_t i=0; i<Basis_CE.size();i++)
+// {
+// std::cout << phi_3[i] << std::endl;
+// }
+//
+
+
@@ -1754,9 +1844,9 @@ for(size_t a = 0; a < 3; a++)
{
for(size_t b=0; b < 3; b++)
{
- assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, tmp1 ,x3MatrixBasis[b]); //
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, tmp1 ,x3MatrixBasis[b]); // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
- auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] );
+ auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
Q[a][b] = coeffContainer[a]*tmp1 + GGterm; // seems symmetric... check positiv definitness?
}
@@ -1767,9 +1857,9 @@ printmatrix(std::cout, Q, "Matrix Q", "--");
for(size_t a = 0; a < 3; a++)
{
- assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp2 ,B);
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp2 ,B_Term); // <L i(M_alpha) + sym(grad phi_alpha), B >
- auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B );
+ auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B_Term); // <L i(x3G_alpha) , B>
B_hat[a] = coeffContainer[a]*tmp2 + GGterm;
@@ -1796,7 +1886,13 @@ for(size_t i=0; i<3; i++)
VTKWriter<GridView> vtkWriter(gridView_CE);
vtkWriter.addVertexData(
correctorFunction_1,
- VTK::FieldInfo("solution", VTK::FieldInfo::Type::vector, dim));
+ VTK::FieldInfo("corrector phi_1", VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(
+ 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));
// vtkWriter.addVertexData(
// solutionFunction,
// VTK::FieldInfo("solution", VTK::FieldInfo::Type::scalar, 1));