diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index 25da2a88ce91f17897f0ad5d2f9532b1ce4f3138..5fd9a296acea3f86fe3d9b37564ef56c14186d8f 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -11,10 +11,14 @@
#include <dune/common/float_cmp.hh>
#include <dune/common/math.hh>
+
+
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh> //TEST
+
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/istl/matrix.hh>
@@ -845,6 +849,95 @@ double computeL2Norm(const Basis& basis,
}
return sqrt(error);
}
+
+
+
+
+
+
+template<class Basis,class LocalFunction1, class LocalFunction2, class GVFunction , class MatrixFunction, class Matrix>
+auto test_derivative(const Basis& basis,
+ LocalFunction1& mu,
+ LocalFunction2& lambda,
+ const double& gamma,
+ Matrix& M,
+ const GVFunction& matrixFieldFuncA,
+ 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(matrixFieldFuncA);
+ 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);
+// printmatrix(std::cout, strain1 , "strain1", "--");
+
+
+// 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;
+}
+
+
+
+
+
+
+
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -995,6 +1088,8 @@ int main(int argc, char *argv[])
Storage_Error.push_back(level);
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::cout << "Number of Elements in each direction: " << nElements << std::endl;
log << "Number of Elements in each direction: " << nElements << std::endl;
@@ -1003,6 +1098,20 @@ int main(int argc, char *argv[])
CellGridType grid_CE(lower,upper,nElements);
using GridView = CellGridType::LeafGridView;
const GridView gridView_CE = grid_CE.leafGridView();
+
+
+ //TEST
+// using CellGridTypeT = StructuredGridFactory<UGGrid<dim> >;
+// auto grid = StructuredGridFactory<UGGrid<dim> >::createCubeGrid(lower,upper,nElements_test);
+// auto gridView_CE = grid::leafGridView();
+
+// FieldVector<double,3> lowerLeft = {0.0, 0.0, 0.0};
+// FieldVector<double,3> upperRight = {1.0, 1.0, 1.0};
+// std::array<unsigned int,3> elements = {10, 10, 10};
+// auto grid = StructuredGridFactory<UGGrid<3> >::createCubeGrid(lowerLeft,
+// upperRight,
+// elements);
+//
using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
@@ -1401,6 +1510,19 @@ int main(int argc, char *argv[])
const std::array<VectorCT, 3> loadContainer = {load_alpha1, load_alpha2, load_alpha3};
const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
+
+
+ //TEST
+
+
+ auto Der1 = derivative(correctorFunction_1);
+ auto Der2 = derivative(correctorFunction_2);
+ auto Der3 = derivative(correctorFunction_3);
+
+// auto output_der = test_derivative(Basis_CE,Der1);
+
+
+// std::cout << "evaluate derivative " << output_der << std::endl;
/////////////////////////////////////////////////////////
// Compute effective quantities: Elastic law & Prestrain
@@ -1408,7 +1530,8 @@ int main(int argc, char *argv[])
MatrixRT Q(0);
VectorCT tmp1,tmp2;
FieldVector<double,3> B_hat ;
-
+
+
// Compute effective elastic law Q
for(size_t a = 0; a < 3; a++)
for(size_t b=0; b < 3; b++)
@@ -1417,11 +1540,23 @@ int main(int argc, char *argv[])
double GGterm = 0.0;
GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
+
+
+ double tmp = 0.0;
+ if(a==0)
+ tmp = test_derivative(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],Der1,x3MatrixBasis[b]);
+ else if (a==1)
+ tmp = test_derivative(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],Der2,x3MatrixBasis[b]);
+ else
+ tmp = test_derivative(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],Der3,x3MatrixBasis[b]);
+
+
// TEST
// std::setprecision(std::numeric_limits<float>::digits10);
- Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric...check positiv definitness?
+// Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric...check positiv definitness?
+ Q[a][b] = tmp + GGterm;
if (print_debug)
{
@@ -1430,6 +1565,28 @@ int main(int argc, char *argv[])
std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
}
}
+
+// // 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;
+// GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
+//
+// // TEST
+// // std::setprecision(std::numeric_limits<float>::digits10);
+//
+// Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric...check positiv definitness?
+//
+// if (print_debug)
+// {
+// std::cout << "analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
+// std::cout << "GGTerm:" << GGterm << std::endl;
+// std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
+// }
+// }
printmatrix(std::cout, Q, "Matrix Q", "--");
log << "Effective Matrix Q: " << std::endl;
log << Q << std::endl;