diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 1cd513ebd1bba927ad125b592c5084f6bfe58edb..a16cca4fed8e19837d285199b0f3b1be98b3e399 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -27,8 +27,8 @@ cellDomain=1
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#----------------------------------------------------
-numLevels= 2 2
-#numLevels = 1 1 # computes all levels from first to second entry
+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
#numLevels = 4 4 # computes all levels from first to second entry
@@ -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 223a5d3b17f5c655558899494dc599d636364353..3b81eed9220bbe1fef2adf8a9bea2012569a0bdd 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -228,15 +228,15 @@ 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;
+// int orderQR = 2;
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-// std::cout << "Print QuadratureOrder:" << orderQR << std::endl; //TEST
+// std::cout << "Print QuadratureOrder:" << orderQR << std::endl; //TEST`
for (const auto& quadPoint : quad)
{
@@ -411,17 +411,17 @@ 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};
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
// int orderQR = 0;
// int orderQR = 1;
- int orderQR = 2;
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
+// int orderQR = 2;
+ 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;
for (const auto& quadPoint : quad)
@@ -442,10 +442,16 @@ void computeElementLoadVector( const LocalView& localView,
// 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;
+//
+// //TEST QUadrature
+// std::cout << "quadPos:" << 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 << "mu(quadPos) :" << mu(quadPos) << std::endl;
+// std::cout << "lambda(quadPos) :" << lambda(quadPos) << std::endl;
// "f*phi"-part
@@ -462,6 +468,7 @@ void computeElementLoadVector( const LocalView& localView,
defGradientV = crossSectionDirectionScaling((1/gamma),defGradientV);
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(quadPos), defGradientV );
+// 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
size_t row = localView.tree().child(k).localIndex(i);
@@ -472,8 +479,10 @@ 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(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;
+// std::cout << "energyDensityfG * quadPoint.weight() * integrationElement: " << energyDensityfG * quadPoint.weight() * integrationElement << std::endl;
}
}
}
@@ -570,6 +579,7 @@ void assembleCellLoad(const Basis& basis,
auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
auto loadFunctional = localFunction(loadGVF);
+ int counter = 1;
for (const auto& element : elements(basis.gridView()))
{
localView.bind(element);
@@ -581,6 +591,8 @@ void assembleCellLoad(const Basis& basis,
// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
BlockVector<FieldVector<double,1> > elementRhs;
+// std::cout << "----------------------------------Element : " << counter << std::endl; //TEST
+ counter++;
computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
// computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, forceTerm); //TEST
// printvector(std::cout, elementRhs, "elementRhs", "--");
@@ -648,7 +660,8 @@ auto energy(const Basis& basis,
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
// int orderQR = 0;
// int orderQR = 1;
- int orderQR = 2;
+// int orderQR = 2;
+ int orderQR = 3;
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -997,7 +1010,8 @@ auto test_derivative(const Basis& basis,
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
// int orderQR = 0;
// int orderQR = 1;
- int orderQR = 2;
+// int orderQR = 2;
+ int orderQR = 3;
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
@@ -1178,7 +1192,8 @@ auto check_Orthogonality(const Basis& basis,
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
// int orderQR = 0;
// int orderQR = 1;
- int orderQR = 2;
+// int orderQR = 2;
+ int orderQR = 3;
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1292,7 +1307,8 @@ auto computeFullQ(const Basis& basis,
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
// int orderQR = 0;
// int orderQR = 1;
- int orderQR = 2;
+// int orderQR = 2;
+ int orderQR = 3;
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1655,6 +1671,8 @@ int main(int argc, char *argv[])
Func2Tensor x3G_2neg = [x3G_2] (const Domain& x) {return -1.0*x3G_2(x);};
Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {return -1.0*x3G_3(x);};
+ //TODO eigentlich funtkioniert es ja mit x3G_1 etc doch auch ?!
+
@@ -1719,6 +1737,10 @@ 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
+// assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1);
+// assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2);
+// assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3);
// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
// printvector(std::cout, load_alpha1, "load_alpha1", "--");
@@ -1734,7 +1756,7 @@ 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 !)