diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 8cb35c60ee4b2aa87e9abb98d234abdb2371ddc1..1cd513ebd1bba927ad125b592c5084f6bfe58edb 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 2
#numLevels = 1 1 # 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
@@ -59,8 +59,8 @@ alpha=8.0
#--- Lame-Parameters
-mu1=8.0
-lambda1=5.0
+mu1=80.0
+lambda1=80.0
# better: pass material parameters as a vector
@@ -113,8 +113,8 @@ write_VTK = true
#############################################
# Assembly options
#############################################
-set_IntegralZero = true
-#set_IntegralZero = false
+#set_IntegralZero = true
+set_IntegralZero = false
#arbitraryLocalIndex = 7
#arbitraryElementNumber = 3
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index 6fb3b2acb23b1b72b0cbf74fdaa492ad2eff6856..223a5d3b17f5c655558899494dc599d636364353 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -55,6 +55,8 @@
#include <dune/solvers/solvers/umfpacksolver.hh> //TEST
+#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
+
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
#include <dune/fufem/dunepython.hh>
@@ -226,7 +228,10 @@ void computeElementStiffnessMatrix(const LocalView& localView,
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
// 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;
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
@@ -371,13 +376,13 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
// Compute the source term for a single element
// < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
-template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class Force>
+template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class LocalForce>
void computeElementLoadVector( const LocalView& localView,
LocalFunction1& mu,
LocalFunction2& lambda,
const double gamma,
Vector& elementRhs,
- const Force& forceTerm
+ const LocalForce& forceTerm
)
{
// | f*phi|
@@ -410,8 +415,14 @@ void computeElementLoadVector( const LocalView& localView,
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
+// int orderQR = 0;
+// int orderQR = 1;
+ int orderQR = 2;
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+
+// std::cout << "Quad-Rule order used: " << orderQR << std::endl;
for (const auto& quadPoint : quad)
{
@@ -426,6 +437,15 @@ void computeElementLoadVector( const LocalView& localView,
for (size_t i=0; i< gradients.size(); i++)
jacobian.mv(referenceGradients[i][0], gradients[i]);
+ //TEST
+// std::cout << "forceTerm(element.geometry().global(quadPos)):" << std::endl;
+// std::cout << forceTerm(element.geometry().global(quadPos)) << std::endl;
+// std::cout << "forceTerm(quadPos)" << std::endl;
+// std::cout << forceTerm(quadPos) << std::endl;
+
+ //TEST QUadrature
+ std::cout << "mu(quadPos) :" << mu(quadPos) << std::endl;
+ std::cout << "lambda(quadPos) :" << lambda(quadPos) << std::endl;
// "f*phi"-part
@@ -562,6 +582,7 @@ void assembleCellLoad(const Basis& basis,
BlockVector<FieldVector<double,1> > elementRhs;
computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
+// computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, forceTerm); //TEST
// printvector(std::cout, elementRhs, "elementRhs", "--");
// printvector(std::cout, elementRhs, "elementRhs", "--");
//////////////////////////////////////////////////////////////////////////////
@@ -624,7 +645,10 @@ 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;
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -971,7 +995,10 @@ auto test_derivative(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 = 0;
+// int orderQR = 1;
+ int orderQR = 2;
+// 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)
@@ -1148,7 +1175,10 @@ 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;
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1259,7 +1289,10 @@ 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;
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1402,6 +1435,11 @@ int main(int argc, char *argv[])
std::cout <<"mu: " <<parameterSet.get<std::array<double,3>>("mu", {1.0,3.0,2.0}) << std::endl;
std::cout <<"lambda: " << parameterSet.get<std::array<double,3>>("lambda", {1.0,3.0,2.0}) << std::endl;
}
+ else
+ {
+ std::cout <<"mu: " << parameterSet.get<double>("mu1",1.0) << std::endl;
+ std::cout <<"lambda: " << parameterSet.get<double>("lambda1",0.0) << std::endl;
+ }
@@ -1484,6 +1522,7 @@ int main(int argc, char *argv[])
CellGridType grid_CE(lower,upper,nElements);
using GridView = CellGridType::LeafGridView;
const GridView gridView_CE = grid_CE.leafGridView();
+ std::cout << "Host grid has " << gridView_CE.size(dim) << " vertices." << std::endl;
//TEST
@@ -1508,7 +1547,6 @@ int main(int argc, char *argv[])
///////////////////////////////////
// Create Lambda-Functions for material Parameters depending on microstructure
///////////////////////////////////
-
auto materialImp = IsotropicMaterialImp<dim>();
auto muTerm = materialImp.getMu(parameterSet);
auto lambdaTerm = materialImp.getLambda(parameterSet);
@@ -1527,7 +1565,6 @@ int main(int argc, char *argv[])
else
return lambda2;
};*/
-
auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
auto muLocal = localFunction(muGridF);
auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
@@ -1539,24 +1576,21 @@ int main(int argc, char *argv[])
// 2 : GMRES
// 3 : QR
///////////////////////////////////
- unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
+ unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 3);
unsigned int Solver_verbosity = parameterSet.get<unsigned int>("Solver_verbosity", 2);
// Print Options
bool print_debug = parameterSet.get<bool>("print_debug", false);
-
-
+
+ //VTK-Write
bool write_materialFunctions = parameterSet.get<bool>("write_materialFunctions", false);
bool write_prestrainFunctions = parameterSet.get<bool>("write_prestrainFunctions", false);
-
-
bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
bool write_corrector_phi2 = parameterSet.get<bool>("write_corrector_phi2", false);
bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
bool write_L2Error = parameterSet.get<bool>("write_L2Error", false);
bool write_IntegralMean = parameterSet.get<bool>("write_IntegralMean", false);
-
bool write_VTK = parameterSet.get<bool>("write_VTK", false);
/////////////////////////////////////////////////////////
@@ -1578,7 +1612,11 @@ int main(int argc, char *argv[])
gridView_CE,
power<dim>( // eig dimworld?!?!
Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
- flatLexicographic()));
+ flatLexicographic()
+// blockedInterleaved() // ERROR
+ ));
+
+ std::cout << "power<periodic> basis has " << Basis_CE.dimension() << " degrees of freedom" << std::endl;
log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
const int phiOffset = Basis_CE.size();
@@ -1589,7 +1627,7 @@ int main(int argc, char *argv[])
VectorCT load_alpha1, load_alpha2, load_alpha3;
MatrixCT stiffnessMatrix_CE;
- bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", true);
+ bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", false);
bool set_oneBasisFunction_Zero = false;
bool substract_integralMean = false;
if(set_IntegralZero)
@@ -1664,24 +1702,31 @@ int main(int argc, char *argv[])
int arbitraryLeafIndex = parameterSet.get<unsigned int>("arbitraryLeafIndex", 0); // localIdx..assert Number < #ShapeFcts
int arbitraryElementNumber = parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
- FieldVector<int,3> row;
- row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLeafIndex);
+
if(print_debug)
+ {
+ FieldVector<int,3> row;
+ row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLeafIndex);
printvector(std::cout, row, "row" , "--");
+ }
/////////////////////////////////////////////////////////
// Assemble the system
/////////////////////////////////////////////////////////
+ Dune::Timer StiffnessTimer;
assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE, parameterSet);
+ std::cout << "Stiffness assembly Timer: " << StiffnessTimer.elapsed() << std::endl;
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);
- // printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
- // printvector(std::cout, load_alpha1, "load_alpha1", "--");
+// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
+// printvector(std::cout, load_alpha1, "load_alpha1", "--");
//TODO Add Option for this
// CHECK SYMMETRY:
// checkSymmetry(Basis_CE,stiffnessMatrix_CE);
+
+
// set one basis-function to zero
@@ -1689,8 +1734,16 @@ int main(int argc, char *argv[])
{
setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha1, load_alpha2, load_alpha3, parameterSet);
// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix after setOneBasisFunctionToZero", "--");
- // printvector(std::cout, load_alpha1, "load_alpha1 after setOneBaseFunctionToZero", "--");
+ printvector(std::cout, load_alpha1, "load_alpha1 after setOneBaseFunctionToZero", "--");
}
+
+ //TEST: Compute Condition Number (Can be very expensive !)
+ const bool verbose = true;
+ const unsigned int arppp_a_verbosity_level = 2;
+ const unsigned int pia_verbosity_level = 1;
+ MatrixInfo<MatrixCT> matrixInfo(stiffnessMatrix_CE,verbose,arppp_a_verbosity_level,pia_verbosity_level);
+ std::cout << "Get condition number of Stiffness_CE: " << matrixInfo.getCond2(true) << std::endl;
+// std::cout << "Get condition number of Stiffness_CE: " << matrixInfo.getCond2(false) << std::endl;
////////////////////////////////////////////////////
// Compute solution