From 6a81f7e69f1bcff5f21b21b3044c2ee4580b7593 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Tue, 26 Jul 2022 14:22:07 +0200
Subject: [PATCH] Add more debugging tests
---
src/Cell-Problem.cc | 286 ++++++++++++++++++++++++++++++++++----------
1 file changed, 223 insertions(+), 63 deletions(-)
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index caf7ceb3..82b4a5ed 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -218,8 +218,8 @@ void computeElementStiffnessMatrix(const LocalView& localView,
///////////////////////////////////////////////
// Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant...
//////////////////////////////////////////////
- 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_1 {{1.0, 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.0}};
MatrixRT G_3 {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
@@ -254,17 +254,22 @@ void computeElementStiffnessMatrix(const LocalView& localView,
MatrixRT defGradientV(0);
defGradientV[l][0] = gradients[j][0]; // Y
defGradientV[l][1] = gradients[j][1]; //X2
- defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
+// defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
+ defGradientV[l][2] = gradients[j][2]; //X3
+
+ defGradientV = crossSectionDirectionScaling((1/gamma),defGradientV);
// "phi*phi"-part
for (size_t k=0; k < dimWorld; k++)
- for (size_t i=0; i < nSf; i++)
- {
+ for (size_t i=0; i < nSf; i++)
+ {
// (scaled) Deformation gradient of the ansatz basis function
MatrixRT defGradientU(0);
defGradientU[k][0] = gradients[i][0]; // Y
defGradientU[k][1] = gradients[i][1]; //X2
- defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
+// defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
+ defGradientU[k][2] = gradients[i][2]; //X3
// printmatrix(std::cout, defGradientU , "defGradientU", "--");
+ defGradientU = crossSectionDirectionScaling((1/gamma),defGradientU);
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
// double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works..
@@ -272,16 +277,16 @@ void computeElementStiffnessMatrix(const LocalView& localView,
size_t col = localView.tree().child(k).localIndex(i);
elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
- }
+ }
- // "m*phi" & "phi*m" - part
- for( size_t m=0; m<3; m++)
- {
- double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
- auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
- elementMatrix[row][localPhiOffset+m] += value;
- elementMatrix[localPhiOffset+m][row] += value;
- }
+ // "m*phi" & "phi*m" - part
+ for( size_t m=0; m<3; m++)
+ {
+ double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
+ auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
+ elementMatrix[row][localPhiOffset+m] += value;
+ elementMatrix[localPhiOffset+m][row] += value;
+ }
}
// "m*m"-part
@@ -397,8 +402,8 @@ void computeElementLoadVector( const LocalView& localView,
///////////////////////////////////////////////
// Basis for R_sym^{2x2}
//////////////////////////////////////////////
- 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_1 {{1.0, 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.0}};
MatrixRT G_3 {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
@@ -428,7 +433,10 @@ void computeElementLoadVector( const LocalView& localView,
MatrixRT defGradientV(0);
defGradientV[k][0] = gradients[i][0]; // Y
defGradientV[k][1] = gradients[i][1]; // X2
- defGradientV[k][2] = (1.0/gamma)*gradients[i][2]; // X3
+// defGradientV[k][2] = (1.0/gamma)*gradients[i][2]; //
+ defGradientV[k][2] = gradients[i][2]; // X3
+
+ defGradientV = crossSectionDirectionScaling((1/gamma),defGradientV);
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(quadPos), defGradientV );
size_t row = localView.tree().child(k).localIndex(i);
@@ -956,8 +964,8 @@ auto test_derivative(const Basis& basis,
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1133,8 +1141,8 @@ auto check_Orthogonality(const Basis& basis,
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1183,6 +1191,113 @@ auto check_Orthogonality(const Basis& basis,
+template<class Basis,class LocalFunction1, class LocalFunction2, class GVFunction , class MatrixFunction, class Matrix>
+auto computeFullQ(const Basis& basis,
+ LocalFunction1& mu,
+ LocalFunction2& lambda,
+ const double& gamma,
+ Matrix& M1,
+ Matrix& M2,
+ const GVFunction& DerPhi_1,
+ const GVFunction& DerPhi_2,
+// const GVFunction& matrixFieldA,
+ const MatrixFunction& matrixFieldFuncG1,
+ const MatrixFunction& matrixFieldFuncG2
+ )
+{
+
+// TEST HIGHER PRECISION
+// using float_50 = boost::multiprecision::cpp_dec_float_50;
+// float_50 higherPrecEnergy = 0.0;
+ double energy = 0.0;
+ constexpr int dim = Basis::LocalView::Element::dimension;
+ constexpr int dimWorld = dim;
+
+ auto localView = basis.localView();
+
+
+ auto DerPhi1 = localFunction(DerPhi_1);
+ auto DerPhi2 = localFunction(DerPhi_2);
+
+ auto matrixFieldG1GVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncG1, basis.gridView());
+ auto matrixFieldG1 = localFunction(matrixFieldG1GVF);
+ auto matrixFieldG2GVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncG2, basis.gridView());
+ auto matrixFieldG2 = localFunction(matrixFieldG2GVF);
+// auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
+// auto matrixFieldB = localFunction(matrixFieldBGVF);
+
+ using GridView = typename Basis::GridView;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+// TEST
+// FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
+// printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
+
+ for (const auto& e : elements(basis.gridView()))
+ {
+ localView.bind(e);
+ matrixFieldG1.bind(e);
+ matrixFieldG2.bind(e);
+ DerPhi1.bind(e);
+ DerPhi2.bind(e);
+ mu.bind(e);
+ lambda.bind(e);
+
+
+ double elementEnergy = 0.0;
+ //double elementEnergy_HP = 0.0;
+
+ auto geometry = e.geometry();
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
+
+ for (const auto& quadPoint : quad)
+ {
+ const auto& quadPos = quadPoint.position();
+ const double integrationElement = geometry.integrationElement(quadPos);
+
+ auto Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, DerPhi1(quadPos))) + *M1;
+ auto Chi2 = sym(crossSectionDirectionScaling(1.0/gamma, DerPhi2(quadPos))) + *M2;
+
+// auto strain1 = DerPhi1(quadPos);
+// // printmatrix(std::cout, strain1 , "strain1", "--");
+// //cale with GAMMA
+// strain1 = crossSectionDirectionScaling(1.0/gamma, strain1);
+// strain1 = sym(strain1);
+
+
+ // ADD M
+// auto test = strain1 + *M ;
+// std::cout << "test:" << test << std::endl;
+
+// for (size_t m=0; m<3; m++ )
+// for (size_t n=0; n<3; n++ )
+// strain1[m][n] += M[m][n];
+
+ auto G1 = matrixFieldG1(quadPos);
+ auto G2 = matrixFieldG2(quadPos);
+
+ auto X1 = G1 + Chi1;
+ auto X2 = G2 + Chi2;
+
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp1, tmp2);
+
+ elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
+
+
+// elementEnergy += strain1 * quadPoint.weight() * integrationElement;
+ //elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
+ }
+ energy += elementEnergy;
+ //higherPrecEnergy += elementEnergy_HP;
+ }
+// TEST
+// std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
+ return energy;
+}
@@ -1485,6 +1600,18 @@ int main(int argc, char *argv[])
Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {return -1.0*x3G_3(x);};
+
+
+
+ //TEST
+ std::cout << "Test crossSectionDirectionScaling:" << std::endl;
+
+ MatrixRT T {{1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}};
+ printmatrix(std::cout, T, "Matrix T", "--");
+
+ auto ST = crossSectionDirectionScaling((1.0/5.0),T);
+ printmatrix(std::cout, ST, "scaled Matrix T", "--");
+
//TEST
// auto QuadraticForm = [] (const double mu, const double lambda, const MatrixRT& M) {
//
@@ -1768,6 +1895,7 @@ int main(int argc, char *argv[])
+
//TEST
// auto local_cor1 = localFunction(correctorFunction_1);
@@ -1782,6 +1910,8 @@ int main(int argc, char *argv[])
auto Der2 = derivative(correctorFunction_2);
auto Der3 = derivative(correctorFunction_3);
+ const std::array<decltype(Der1)*,3> phiDerContainer = {&Der1, &Der2, &Der3};
+
// auto output_der = test_derivative(Basis_CE,Der1);
@@ -1795,6 +1925,8 @@ int main(int argc, char *argv[])
FieldVector<double,3> B_hat ;
+
+ //VARIANT 1
//Compute effective elastic law Q
for(size_t a = 0; a < 3; a++)
for(size_t b=0; b < 3; b++)
@@ -1807,21 +1939,33 @@ int main(int argc, char *argv[])
MGterm = energy_MG(Basis_CE, muLocal, lambdaLocal, mContainer[a], x3MatrixBasis[b]);
double tmp = 0.0;
- if(a==0)
- {
- tmp = test_derivative(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],Der1,x3MatrixBasis[b]);
- std::cout << "check_Orthogonality:" << check_Orthogonality(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[1],Der1,Der2,x3MatrixBasis[a]) << std::endl;
- }
- else if (a==1)
- {
- tmp = test_derivative(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],Der2,x3MatrixBasis[b]);
- std::cout << "check_Orthogonality:" << check_Orthogonality(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[1],Der2,Der2,x3MatrixBasis[a]) << std::endl;
- }
- else
- {
- tmp = test_derivative(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],Der3,x3MatrixBasis[b]);
- std::cout << "check_Orthogonality:" << check_Orthogonality(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[1],Der3,Der2,x3MatrixBasis[a]) << std::endl;
- }
+
+ tmp = test_derivative(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],*phiDerContainer[a],x3MatrixBasis[b]);
+
+
+
+ std::cout << "---- (" << a << "," << b << ") ---- " << std::endl;
+
+ std::cout << "check_Orthogonality:" << check_Orthogonality(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[b],*phiDerContainer[a],*phiDerContainer[b],x3MatrixBasis[a]) << std::endl;
+
+
+
+
+// if(a==0)
+// {
+// tmp = test_derivative(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],Der1,x3MatrixBasis[b]);
+// std::cout << "check_Orthogonality:" << check_Orthogonality(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[1],Der1,Der2,x3MatrixBasis[a]) << std::endl;
+// }
+// else if (a==1)
+// {
+// tmp = test_derivative(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],Der2,x3MatrixBasis[b]);
+// std::cout << "check_Orthogonality:" << check_Orthogonality(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[1],Der2,Der2,x3MatrixBasis[a]) << std::endl;
+// }
+// else
+// {
+// tmp = test_derivative(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],Der3,x3MatrixBasis[b]);
+// std::cout << "check_Orthogonality:" << check_Orthogonality(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[1],Der3,Der2,x3MatrixBasis[a]) << std::endl;
+// }
std::cout << "GGTerm:" << GGterm << std::endl;
@@ -1829,9 +1973,7 @@ int main(int argc, char *argv[])
std::cout << "tmp:" << tmp << std::endl;
std::cout << "(coeffContainer[a]*tmp1):" << (coeffContainer[a]*tmp1) << std::endl;
-
-
-
+
// TEST
// std::setprecision(std::numeric_limits<float>::digits10);
@@ -1845,31 +1987,49 @@ int main(int argc, char *argv[])
std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
}
}
-
-// // Compute effective elastic law Q
-// 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]); // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
-//
-// double GGterm = 0.0;
-// GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
-//
-// // TEST
-// // std::setprecision(std::numeric_limits<float>::digits10);
-//
-// Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric...check positiv definitness?
-//
-// if (print_debug)
-// {
-// std::cout << "analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
-// std::cout << "GGTerm:" << GGterm << std::endl;
-// std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
-// }
-// }
printmatrix(std::cout, Q, "Matrix Q", "--");
log << "Effective Matrix Q: " << std::endl;
log << Q << std::endl;
+
+
+
+ //VARIANT 2
+ //Compute effective elastic law Q
+ MatrixRT Q_2(0);
+ for(size_t a = 0; a < 3; a++)
+ for(size_t b=0; b < 3; b++)
+ {
+ Q_2[a][b] = computeFullQ(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[b],*phiDerContainer[a],*phiDerContainer[b],x3MatrixBasis[a],x3MatrixBasis[b]);
+ }
+ printmatrix(std::cout, Q_2, "Matrix Q_2", "--");
+
+ // VARIANT 3
+ // Compute effective elastic law Q
+ MatrixRT Q_3(0);
+ 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]); // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
+
+ double GGterm = 0.0;
+ GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
+
+ // TEST
+ // std::setprecision(std::numeric_limits<float>::digits10);
+
+ Q_3[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric...check positiv definitness?
+
+ if (print_debug)
+ {
+ std::cout << "analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
+ std::cout << "GGTerm:" << GGterm << std::endl;
+ std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
+ }
+ }
+ printmatrix(std::cout, Q_3, "Matrix Q_3", "--");
+
+
+
// compute B_hat
for(size_t a = 0; a<3; a++)
@@ -2347,5 +2507,5 @@ int main(int argc, char *argv[])
log.close();
-// std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
+ std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
}
--
GitLab