diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index 380efd9ecc5ae4bd05cd8b0142fb1c3f790282ba..98c1650242b85d062ab1325dac654065545728c8 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -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
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index 5fd9a296acea3f86fe3d9b37564ef56c14186d8f..caf7ceb383ef71af5b23c9df8e91ecf980a3bc03 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -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)
{