Skip to content
Snippets Groups Projects
Commit d7b03048 authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

Add more debugging tests

parent 270c411b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 223 additions and 63 deletions
src/Cell-Problem.cc
+ 223
− 63
View file @ d7b03048
......@@ -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;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment