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

Add Option to check Symmetry & Orthogonality(75)

parent a5f48f4a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/microstructure/matrix_operations.hh
+ 2
− 3
View file @ 27a4484e
......@@ -68,7 +68,7 @@ namespace MatrixOperations {
}
static MatrixRT crossSectionDirectionScaling(double w, MatrixRT M){
return {M[0], w*M[1], w*M[2]};
return {M[0], M[1], w*M[2]};
}
static double trace (MatrixRT M){
......@@ -96,8 +96,7 @@ namespace MatrixOperations {
static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2) // CHANGED
{
auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
return scalarProduct(t1,E2);
return scalarProduct(t1,sym(E2));
}
// --- Generalization: Define Quadratic QuadraticForm
......
src/Cell-Problem.cc
+ 293
− 13
View file @ 27a4484e
......@@ -104,6 +104,51 @@ std::string prd(const type x, const int decDigits, const int width) {
template<class Basis, class Matrix>
void checkSymmetry(const Basis& basis,
Matrix& matrix
)
{
std::cout << "--- Check Symmetry ---" << std::endl;
auto localView = basis.localView();
const int phiOffset = basis.dimension();
for (const auto& element : elements(basis.gridView()))
{
localView.bind(element);
const int localPhiOffset = localView.size();
for (size_t i=0; i<localPhiOffset; i++)
for (size_t j=0; j<localPhiOffset; j++ )
{
auto row = localView.index(i);
auto col = localView.index(j);
if(abs( matrix[row][col] - matrix[col][row]) > 1e-12 )
std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
for (size_t i=0; i<localPhiOffset; i++)
for(size_t m=0; m<3; m++)
{
auto row = localView.index(i);
if(abs( matrix[row][phiOffset+m] - matrix[phiOffset+m][row]) > 1e-12 )
std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
for (size_t m=0; m<3; m++ )
for (size_t n=0; n<3; n++ )
{
if(abs( matrix[phiOffset+m][phiOffset+n] - matrix[phiOffset+n][phiOffset+m]) > 1e-12 )
std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
}
std::cout << "--- Symmetry test passed ---" << std::endl;
}
template<class Basis>
auto arbitraryComponentwiseIndices(const Basis& basis,
......@@ -385,7 +430,7 @@ void computeElementLoadVector( const LocalView& localView,
defGradientV[k][1] = gradients[i][1]; // X2
defGradientV[k][2] = (1.0/gamma)*gradients[i][2]; // X3
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientV, forceTerm(quadPos));
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(quadPos), defGradientV );
size_t row = localView.tree().child(k).localIndex(i);
elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
}
......@@ -393,7 +438,7 @@ void computeElementLoadVector( const LocalView& localView,
// "f*m"-part
for (size_t m=0; m<3; m++)
{
double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], forceTerm(quadPos));
double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(quadPos),basisContainer[m] );
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
}
}
......@@ -436,7 +481,15 @@ void assembleCellStiffness(const Basis& basis,
//printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
//std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
//std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
//TEST
//Check Symmetry:
for (size_t i=0; i<localPhiOffset; i++)
for (size_t j=0; j<localPhiOffset; j++ )
{
if(abs( elementMatrix[i][j] - elementMatrix[j][i]) > 1e-12 )
std::cout << "ELEMENT-STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
//////////////////////////////////////////////////////////////////////////////
// GLOBAL STIFFNES ASSEMBLY
//////////////////////////////////////////////////////////////////////////////
......@@ -862,6 +915,7 @@ auto test_derivative(const Basis& basis,
const double& gamma,
Matrix& M,
const GVFunction& matrixFieldFuncA,
// const GVFunction& matrixFieldA,
const MatrixFunction& matrixFieldFuncB
)
{
......@@ -911,20 +965,214 @@ auto test_derivative(const Basis& basis,
const auto& quadPos = quadPoint.position();
const double integrationElement = geometry.integrationElement(quadPos);
// auto strain1 = matrixFieldA(quadPos);
// auto strain2 = matrixFieldB(quadPos);
auto strain1 = matrixFieldA(quadPos);
auto strain2 = matrixFieldB(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];
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), strain1, strain2);
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;
}
template<class Basis, class LocalFunction1, class LocalFunction2, class MatrixFunction, class Matrix>
auto energy_MG(const Basis& basis,
LocalFunction1& mu,
LocalFunction2& lambda,
Matrix& M,
const MatrixFunction& matrixFieldFuncB)
{
// 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 matrixFieldAGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncA, basis.gridView());
// auto matrixFieldA = localFunction(matrixFieldAGVF);
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);
// matrixFieldA.bind(e);
matrixFieldB.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 strain1 = matrixFieldA(quadPos);
auto strain2 = matrixFieldB(quadPos);
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), *M, strain2);
elementEnergy += energyDensity * 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;
}
template<class Basis,class LocalFunction1, class LocalFunction2, class GVFunction , class MatrixFunction, class Matrix>
auto check_Orthogonality(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& matrixFieldFuncG
)
{
// 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 matrixFieldGGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncG, basis.gridView());
auto matrixFieldG = localFunction(matrixFieldGGVF);
// 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);
matrixFieldG.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 Chi = sym(crossSectionDirectionScaling(1.0/gamma, DerPhi2(quadPos))) + *M2;
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), strain1, strain2);
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 G = matrixFieldG(quadPos);
auto tmp = G + *M1 + strain1;
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp, Chi);
// elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
// elementEnergy += strain1 * quadPoint.weight() * integrationElement;
//elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
}
// energy += elementEnergy;
energy += elementEnergy;
//higherPrecEnergy += elementEnergy_HP;
}
// TEST
......@@ -938,6 +1186,8 @@ auto test_derivative(const Basis& basis,
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
......@@ -992,6 +1242,9 @@ int main(int argc, char *argv[])
// Python::Module module = Python::import(parameterSet.get<std::string>("geometryFunction"));
// auto pythonInitialIterate = Python::make_function<double>(module.get("f"));
std::cout << "machine epsilon:" << std::numeric_limits<double>::epsilon() << std::endl;
constexpr int dim = 3;
constexpr int dimWorld = 3;
......@@ -1089,7 +1342,7 @@ int main(int argc, char *argv[])
Storage_Quantities.push_back(level);
std::array<int, dim> nElements = { (int)std::pow(2,level) , (int)std::pow(2,level) , (int)std::pow(2,level) };
std::array<unsigned int, dim> nElements_test = { (int)std::pow(2,level) , (int)std::pow(2,level) , (int)std::pow(2,level) };
// std::array<unsigned int, dim> nElements_test = { (int)std::pow(2,level) , (int)std::pow(2,level) , (int)std::pow(2,level) };
std::cout << "Number of Elements in each direction: " << nElements << std::endl;
log << "Number of Elements in each direction: " << nElements << std::endl;
......@@ -1224,7 +1477,7 @@ int main(int argc, char *argv[])
};
Func2Tensor x3G_3 = [] (const Domain& x) {
return MatrixRT{{0.0, 1.0/sqrt(2.0)*x[2], 0.0}, {1.0/sqrt(2.0)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
return MatrixRT{{0.0, (1.0/sqrt(2.0))*x[2], 0.0}, {(1.0/sqrt(2.0))*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
Func2Tensor x3G_1neg = [x3G_1] (const Domain& x) {return -1.0*x3G_1(x);};
......@@ -1281,6 +1534,9 @@ int main(int argc, char *argv[])
// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
// printvector(std::cout, load_alpha1, "load_alpha1", "--");
// CHECK SYMMETRY:
checkSymmetry(Basis_CE,stiffnessMatrix_CE);
// set one basis-function to zero
......@@ -1514,6 +1770,13 @@ int main(int argc, char *argv[])
//TEST
// auto local_cor1 = localFunction(correctorFunction_1);
// auto local_cor2 = localFunction(correctorFunction_2);
// auto local_cor3 = localFunction(correctorFunction_3);
//
// auto Der1 = derivative(local_cor1);
// auto Der2 = derivative(local_cor2);
// auto Der3 = derivative(local_cor3);
auto Der1 = derivative(correctorFunction_1);
auto Der2 = derivative(correctorFunction_2);
......@@ -1532,23 +1795,40 @@ int main(int argc, char *argv[])
FieldVector<double,3> B_hat ;
// Compute effective elastic law Q
//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;
double MGterm = 0.0;
GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
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;
}
std::cout << "GGTerm:" << GGterm << std::endl;
std::cout << "MGTerm:" << MGterm << std::endl;
std::cout << "tmp:" << tmp << std::endl;
std::cout << "(coeffContainer[a]*tmp1):" << (coeffContainer[a]*tmp1) << std::endl;
......@@ -1556,7 +1836,7 @@ int main(int argc, char *argv[])
// std::setprecision(std::numeric_limits<float>::digits10);
// Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric...check positiv definitness?
Q[a][b] = tmp + GGterm;
Q[a][b] = tmp + GGterm; // TODO : Zusammenfassen in einer Funktion ...
if (print_debug)
{
......
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