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

Run some tests

parent 9f5eab60
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
inputs/cellsolver.parset
+ 4
− 4
View file @ 9995a068
......@@ -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
......
src/Cell-Problem.cc
+ 38
− 16
View file @ 9995a068
......@@ -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 !)
......
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