From 74509e93101b54648133ac2aa2cf98f7f1a4cf1e Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Sun, 31 Jul 2022 19:38:51 +0200
Subject: [PATCH] FIX BUG: (+)= Missing in ElementAssembly m*m-part
---
dune/microstructure/matrix_operations.hh | 25 +++---
inputs/cellsolver.parset | 7 +-
src/CMakeLists.txt | 2 +-
src/Cell-Problem.cc | 110 ++++++++++++++++-------
4 files changed, 96 insertions(+), 48 deletions(-)
diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index a4a563bc..1c89d36f 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -83,7 +83,7 @@ namespace MatrixOperations {
}
static double scalarProduct (MatrixRT M1, MatrixRT M2){
- double sum = 0;
+ double sum = 0.0;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
sum += M1[i][j] * M2[i][j];
@@ -99,7 +99,7 @@ namespace MatrixOperations {
// t2 *= lambda;
// return t1 + t2;
// }
-
+//
// static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2) // CHANGED
// {
// auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
@@ -111,15 +111,15 @@ namespace MatrixOperations {
auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
auto tmp1 = scalarProduct(t1,sym(E2));
- auto t2 = 2.0 * mu * sym(E2) + lambda * trace(sym(E2)) * Id();
- auto tmp2 = scalarProduct(t2,sym(E1));
-
- if(abs(tmp1-tmp2) > 1e-8 )
- {
- std::cout << "linearizedStVenantKirchhoffDensity NOT Symmetric!" << std::endl;
- std::cout << "tmp1: " << tmp1 << std::endl;
- std::cout << "tmp2: " << tmp2 << std::endl;
- }
+// auto t2 = 2.0 * mu * sym(E2) + lambda * trace(sym(E2)) * Id();
+// auto tmp2 = scalarProduct(t2,sym(E1));
+//
+// if(abs(tmp1-tmp2) > 1e-8 )
+// {
+// std::cout << "linearizedStVenantKirchhoffDensity NOT Symmetric!" << std::endl;
+// std::cout << "tmp1: " << tmp1 << std::endl;
+// std::cout << "tmp2: " << tmp2 << std::endl;
+// }
return tmp1;
}
@@ -130,7 +130,8 @@ namespace MatrixOperations {
auto tmp1 = sym(M);
double tmp2 = norm(tmp1);
- return lambda*std::pow(trace(M),2) + 2*mu*pow( tmp2 ,2);
+// double tmp2 = norm(M); //TEST
+ return lambda*std::pow(trace(M),2) + 2.0*mu*pow( tmp2 ,2);
// return lambda*std::pow(trace(M),2) + 2*mu*pow( norm( sym(M) ),2);
}
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index a16cca4f..afc14b79 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -27,7 +27,7 @@ cellDomain=1
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#----------------------------------------------------
-numLevels= 2 4
+numLevels= 2 4
#numLevels = 0 0 # computes all levels from first to second entry
#numLevels = 2 2 # computes all levels from first to second entry
#numLevels = 1 3 # computes all levels from first to second entry
@@ -113,8 +113,9 @@ write_VTK = true
#############################################
# Assembly options
#############################################
-set_IntegralZero = true
-#set_IntegralZero = false
+#set_IntegralZero = true
+set_IntegralZero = false
+#set_oneBasisFunction_Zero = true
#arbitraryLocalIndex = 7
#arbitraryElementNumber = 3
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 3264ac4d..f3f3b05f 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -16,4 +16,4 @@ foreach(_program ${programs})
endforeach()
-#set(CMAKE_BUILD_TYPE Debug)
+set(CMAKE_BUILD_TYPE Debug)
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index 3b81eed9..787eb9d0 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -232,14 +232,17 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// int orderQR = 0;
// int orderQR = 1;
// int orderQR = 2;
+// int orderQR = 3;
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
+// double elementContribution = 0.0;
// std::cout << "Print QuadratureOrder:" << orderQR << std::endl; //TEST`
+ int QPcounter= 0;
for (const auto& quadPoint : quad)
{
+ QPcounter++;
const auto& quadPos = quadPoint.position();
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
const auto integrationElement = geometry.integrationElement(quadPos);
@@ -254,8 +257,8 @@ void computeElementStiffnessMatrix(const LocalView& localView,
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
for (size_t l=0; l< dimWorld; l++)
- for (size_t j=0; j < nSf; j++ )
- {
+ for (size_t j=0; j < nSf; j++ )
+ {
size_t row = localView.tree().child(l).localIndex(j);
// (scaled) Deformation gradient of the test basis function
MatrixRT defGradientV(0);
@@ -264,7 +267,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
defGradientV[l][2] = gradients[j][2]; //X3
- defGradientV = crossSectionDirectionScaling((1/gamma),defGradientV);
+ defGradientV = crossSectionDirectionScaling((1.0/gamma),defGradientV);
// "phi*phi"-part
for (size_t k=0; k < dimWorld; k++)
for (size_t i=0; i < nSf; i++)
@@ -276,9 +279,10 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// 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);
+ defGradientU = crossSectionDirectionScaling((1.0/gamma),defGradientU);
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
+// double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), defGradientU, defGradientV); //TEST
// double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works..
size_t col = localView.tree().child(k).localIndex(i);
@@ -290,20 +294,41 @@ void computeElementStiffnessMatrix(const LocalView& localView,
for( size_t m=0; m<3; m++)
{
double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
+// double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], defGradientV); //TEST
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value;
elementMatrix[localPhiOffset+m][row] += value;
}
- }
+ }
// "m*m"-part
- for(size_t m=0; m<3; m++)
+ for(size_t m=0; m<3; m++) //TODO ist symmetric.. reicht die hälfte zu berechnen!!!
for(size_t n=0; n<3; n++)
{
+
+// std::cout << "QPcounter: " << QPcounter << std::endl;
+// std::cout << "m= " << m << " n = " << n << std::endl;
+// printmatrix(std::cout, basisContainer[m] , "basisContainer[m]", "--");
+// printmatrix(std::cout, basisContainer[n] , "basisContainer[n]", "--");
+// std::cout << "integrationElement:" << integrationElement << std::endl;
+// std::cout << "quadPoint.weight(): "<< quadPoint.weight() << std::endl;
+// std::cout << "mu(quadPos): " << mu(quadPos) << std::endl;
+// std::cout << "lambda(quadPos): " << lambda(quadPos) << std::endl;
+
+
double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
- elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
+// double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], basisContainer[n]); //TEST
+ elementMatrix[localPhiOffset+m][localPhiOffset+n] += energyDensityGG * quadPoint.weight() * integrationElement; // += !!!!! (Fixed-Bug)
+
+// std::cout << "energyDensityGG:" << energyDensityGG<< std::endl;
+// std::cout << "product " << energyDensityGG * quadPoint.weight() * integrationElement << std::endl;
+// printmatrix(std::cout, elementMatrix, "elementMatrix", "--");
+
}
}
+
+// std::cout << "Number of QuadPoints:" << QPcounter << std::endl;
+// printmatrix(std::cout, elementMatrix, "elementMatrix", "--");
}
@@ -411,7 +436,6 @@ 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.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};
@@ -419,8 +443,8 @@ void computeElementLoadVector( const LocalView& localView,
// int orderQR = 0;
// int orderQR = 1;
// int orderQR = 2;
- int orderQR = 3;
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+// int orderQR = 3;
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
// std::cout << "Quad-Rule order used: " << orderQR << std::endl;
@@ -445,11 +469,11 @@ void computeElementLoadVector( const LocalView& localView,
//
// //TEST QUadrature
// std::cout << "quadPos:" << quadPos << std::endl;
-// // std::cout << "element.geometry().global(quadPos):" << element.geometry().global(quadPos) << std::endl;
-// //
-// //
+// std::cout << "element.geometry().global(quadPos):" << element.geometry().global(quadPos) << std::endl;
+// // //
+// // //
// std::cout << "quadPoint.weight() :" << quadPoint.weight() << std::endl;
-// // std::cout << "integrationElement(quadPos):" << integrationElement << std::endl;
+// std::cout << "integrationElement(quadPos):" << integrationElement << std::endl;
// std::cout << "mu(quadPos) :" << mu(quadPos) << std::endl;
// std::cout << "lambda(quadPos) :" << lambda(quadPos) << std::endl;
@@ -465,9 +489,10 @@ void computeElementLoadVector( const LocalView& localView,
// defGradientV[k][2] = (1.0/gamma)*gradients[i][2]; //
defGradientV[k][2] = gradients[i][2]; // X3
- defGradientV = crossSectionDirectionScaling((1/gamma),defGradientV);
+ defGradientV = crossSectionDirectionScaling((1.0/gamma),defGradientV);
- double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(quadPos), defGradientV );
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(quadPos), defGradientV );
+// double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)),forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST
@@ -479,6 +504,7 @@ void computeElementLoadVector( const LocalView& localView,
for (size_t m=0; m<3; m++)
{
double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(quadPos),basisContainer[m] );
+// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(element.geometry().global(quadPos)),basisContainer[m] );//TEST
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
@@ -526,7 +552,7 @@ void assembleCellStiffness(const Basis& basis,
//std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
//TEST
- //Check Symmetry:
+ //Check Element-Stiffness-Symmetry:
for (size_t i=0; i<localPhiOffset; i++)
for (size_t j=0; j<localPhiOffset; j++ )
{
@@ -579,7 +605,7 @@ void assembleCellLoad(const Basis& basis,
auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
auto loadFunctional = localFunction(loadGVF);
- int counter = 1;
+// int counter = 1;
for (const auto& element : elements(basis.gridView()))
{
localView.bind(element);
@@ -592,7 +618,7 @@ void assembleCellLoad(const Basis& basis,
BlockVector<FieldVector<double,1> > elementRhs;
// std::cout << "----------------------------------Element : " << counter << std::endl; //TEST
- counter++;
+// counter++;
computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
// computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, forceTerm); //TEST
// printvector(std::cout, elementRhs, "elementRhs", "--");
@@ -657,11 +683,11 @@ auto energy(const Basis& basis,
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);
// int orderQR = 0;
// int orderQR = 1;
// int orderQR = 2;
- int orderQR = 3;
+// int orderQR = 3;
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1011,8 +1037,8 @@ auto test_derivative(const Basis& basis,
// int orderQR = 0;
// int orderQR = 1;
// int orderQR = 2;
- int orderQR = 3;
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+// int orderQR = 3;
+ 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)
@@ -1037,7 +1063,8 @@ auto test_derivative(const Basis& basis,
- double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), strain1, strain2);
+// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), strain1, strain2);
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), test, strain2);
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
@@ -1168,6 +1195,10 @@ auto check_Orthogonality(const Basis& basis,
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+
+// std::cout << "Press key.." << std::endl;
+// std::cin.get();
+
// TEST
// FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
// printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
@@ -1189,11 +1220,11 @@ auto check_Orthogonality(const Basis& basis,
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);
// int orderQR = 0;
// int orderQR = 1;
// int orderQR = 2;
- int orderQR = 3;
+// int orderQR = 3;
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1304,11 +1335,11 @@ auto computeFullQ(const Basis& basis,
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);
// int orderQR = 0;
// int orderQR = 1;
// int orderQR = 2;
- int orderQR = 3;
+// int orderQR = 3;
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1341,10 +1372,11 @@ auto computeFullQ(const Basis& basis,
auto X1 = G1 + Chi1;
auto X2 = G2 + Chi2;
-
+
+
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, X2);
- elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
+ elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ???
// elementEnergy += strain1 * quadPoint.weight() * integrationElement;
@@ -1644,7 +1676,8 @@ int main(int argc, char *argv[])
MatrixCT stiffnessMatrix_CE;
bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", false);
- bool set_oneBasisFunction_Zero = false;
+ bool set_oneBasisFunction_Zero = parameterSet.get<bool>("set_oneBasisFunction_Zero", false);
+// bool set_oneBasisFunction_Zero = false;
bool substract_integralMean = false;
if(set_IntegralZero)
{
@@ -1737,6 +1770,13 @@ int main(int argc, char *argv[])
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1neg);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2neg);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3neg);
+ //TEST
+// assembleCellStiffness(Basis_CE, muTerm, lambdaTerm, gamma, stiffnessMatrix_CE, parameterSet);
+// std::cout << "Stiffness assembly Timer: " << StiffnessTimer.elapsed() << std::endl;
+// assembleCellLoad(Basis_CE, muTerm, lambdaTerm,gamma, load_alpha1 ,x3G_1neg);
+// assembleCellLoad(Basis_CE, muTerm, lambdaTerm,gamma, load_alpha2 ,x3G_2neg);
+// assembleCellLoad(Basis_CE, muTerm, lambdaTerm,gamma, load_alpha3 ,x3G_3neg);
+
//TEST
// assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1);
// assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2);
@@ -1958,6 +1998,9 @@ int main(int argc, char *argv[])
log << " --------------------" << std::endl;
log << "Corrector-Matrix M_3: \n" << M_3 << std::endl;
log << " --------------------" << std::endl;
+
+
+
////////////////////////////////////////////////////////////////////////////
// Substract Integral-mean
@@ -2085,6 +2128,9 @@ int main(int argc, char *argv[])
// std::cout << "evaluate derivative " << output_der << std::endl;
+
+
+ // TODO : MOVE All of this into a separate class : 'computeEffectiveQuantities'
/////////////////////////////////////////////////////////
// Compute effective quantities: Elastic law & Prestrain
/////////////////////////////////////////////////////////
--
GitLab