diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index 018e09a70312d4c32d1f88aa9d0b86216ca0b4d0..380efd9ecc5ae4bd05cd8b0142fb1c3f790282ba 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -83,14 +83,21 @@ namespace MatrixOperations {
return sum;
}
- static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2){
- E1= sym(E1);
- E2 = sym(E2);
- double t1 = scalarProduct(E1,E2);
- t1 *= 2* mu;
- double t2 = trace(E1)*trace(E2);
- t2 *= lambda;
- return t1 + t2;
+// static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2){ // ?? Check with Robert
+// E1= sym(E1);
+// E2 = sym(E2);
+// double t1 = scalarProduct(E1,E2);
+// t1 *= 2* mu;
+// double t2 = trace(E1)*trace(E2);
+// 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();
+ return scalarProduct(t1,E2);
+
}
// --- Generalization: Define Quadratic QuadraticForm
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 73e5b2f8f2b657342bffd258370bb8b884c96ac9..abee48768ad5ffe02a85ea78fb487a08a5ae1ebd 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -5,6 +5,7 @@
set(programs Cell-Problem
Cell-Problem_muGamma
+ Compute_MuGamma
)
diff --git a/src/Compute_MuGamma.cc b/src/Compute_MuGamma.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8482b7dbf91846cf29e968a0b42a06c72dcde7c6
--- /dev/null
+++ b/src/Compute_MuGamma.cc
@@ -0,0 +1,1043 @@
+#include <config.h>
+
+#include <vector>
+
+#include <dune/geometry/quadraturerules.hh>
+
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+#include <dune/istl/matrix.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/istl/io.hh>
+#include <dune/istl/matrix.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/matrixmarket.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh> // use for LagrangeBasis
+
+
+#include <dune/common/float_cmp.hh>
+#include <dune/common/math.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+
+#include <dune/functions/functionspacebases/periodicbasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+
+
+#include <dune/microstructure/prestrain_material_geometry.hh>
+// #include <dune/microstructure/matrix_operations.hh>
+
+
+// #include <math.h>
+#include <cmath>
+// #include <numbers>
+
+using namespace Dune;
+
+
+
+// ------------------------------------------------------------------------
+//
+// Solving the 2D "Poisson-Type" Equation ( (38) in the draft)
+// in order to compute mu_Gamma in a faster manner
+//
+// ------------------------------------------------------------------------
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+template<class LocalView, class Matrix, class LocalFunction>
+void assembleElementStiffnessMatrix(const LocalView& localView,
+ Matrix& elementMatrix,
+ const LocalFunction& mu,
+ const double gamma
+ )
+{
+ using Element = typename LocalView::Element;
+ constexpr int dim = Element::dimension;
+ const Element element = localView.element();
+ auto geometry = element.geometry();
+
+
+ const auto& localFiniteElement = localView.tree().finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
+
+
+// std::cout << "nSf:" << nSf << std::endl;
+
+ elementMatrix.setSize(localView.size(),localView.size());
+ elementMatrix = 0;
+
+// int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1) + 5; // doesnt change anything
+ int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
+// int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
+ const auto& quadRule = QuadratureRules<double, dim>::rule(element.type(),orderQR);
+
+// std::cout << "QuadRule-Order:" << orderQR << std::endl;
+
+// double muValue = 0.0;
+// std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
+
+ for (const auto& quadPoint : quadRule)
+ {
+
+
+ const auto& quadPos = quadPoint.position();
+
+// std::cout << " quadPOS: " << quadPos << std::endl;
+ // TEST
+// double muValue = mu(quadPos);
+// std::cout << "muValue:" << muValue << std::endl;
+
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+ const double integrationElement = geometry.integrationElement(quadPos);
+
+
+
+ std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceGradients);
+
+ // Compute the shape function gradients on the grid element
+ std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+
+ for (size_t i=0; i<gradients.size(); i++)
+ jacobian.mv(referenceGradients[i][0], gradients[i]); //TODO ? passt..
+
+
+// // print all gradients //TEST
+// FieldVector<double,dim> tmp = {0.0 , 0.0};
+// for (size_t i=0; i < nSf; i++)
+// {
+// printvector(std::cout, gradients[i], "gradients[i]", "--" );
+//
+//
+// tmp[0] += gradients[i][0];
+// tmp[1] += gradients[i][1];
+//
+// // tmp[0] += coeffVector[globalIdx]*gradients[i][0]; // 3D-Version
+// // tmp[1] += coeffVector[globalIdx]*gradients[i][2];
+// // printvector(std::cout, tmp, "tmp", "--" );
+//
+// }
+// printvector(std::cout, tmp, "sum of basisfunction Gradients", "--" );
+
+ for (size_t p=0; p<elementMatrix.N(); p++)
+ {
+ for (size_t q=0; q<elementMatrix.M(); q++)
+ {
+// auto localRow = localView.tree().localIndex(p); // VERTAUSCHT?!?!
+// auto localCol = localView.tree().localIndex(q);
+ auto localRow = localView.tree().localIndex(q);
+ auto localCol = localView.tree().localIndex(p);
+
+// auto value = mu(quadPos)*(2.0* gradients[p][0] * gradients[q][0])+ mu(quadPos)*((1.0/(std::pow(gamma,2)))*(gradients[p][1] * gradients[q][1]));
+// auto value = (mu(quadPos)*gradients[p][0] * gradients[q][0])+ ((1.0/gamma)*(gradients[p][1] * gradients[q][1]));
+
+
+
+// std::cout << "gradients[p][0]" << gradients[p][0] << std::endl;
+// std::cout << "gradients[q][0]" << gradients[q][0] << std::endl;
+// std::cout << "(1.0/std::pow(gamma,2))*gradients[p][1]" << (1.0/std::pow(gamma,2))*gradients[p][1]<< std::endl;
+
+
+
+// auto value3 = mu(quadPos)*(1.0/gamma)*gradients[p][2]; // 3D-Version
+// auto value4 = value3*gradients[q][2]; // 3D-Version
+
+ auto value1 = 2.0*mu(quadPos)* gradients[p][0] *gradients[q][0]; //#1
+// auto value1 = 2.0*mu(quadPos)*sqrt(2.0)* gradients[p][0] *gradients[q][0]; // TEST TODO warum passt es damit besser bei gamma = 1.0 ???
+// auto value2 = mu(quadPos)*(1.0/std::pow(gamma,2))*gradients[p][1] * gradients[q][1] ;
+
+
+
+
+
+ auto value2 = 2.0*mu(quadPos)*(1.0/std::pow(gamma,2))*gradients[p][1] * gradients[q][1] ; //#1 TEST FAKTOR 2 hat hier gefehlt?!?!
+
+ auto value = value1 + value2;
+
+ elementMatrix[localRow][localCol] += value * quadPoint.weight() * integrationElement;
+ }
+ }
+ }
+}
+
+
+
+
+// Compute the source term for a single element
+template<class LocalView,class Vector, class LocalFunction, class Force>
+void assembleElementVolumeTerm(const LocalView& localView,
+ Vector& elementRhs,
+ LocalFunction& mu,
+ const Force& forceTerm)
+{
+ using Element = typename LocalView::Element;
+ auto element = localView.element();
+ auto geometry = element.geometry();
+
+
+ constexpr int dim = Element::dimension;
+ // Set of shape functions for a single element
+ const auto& localFiniteElement = localView.tree().finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
+ // Set all entries to zero
+ elementRhs.resize(localFiniteElement.size());
+ elementRhs = 0;
+
+ // A quadrature rule
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;
+ int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
+// int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
+
+// std::cout << "elementRhs.size():" << elementRhs.size() << std::endl;
+
+
+
+ const auto& quadRule = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+ for (const auto& quadPoint : quadRule)
+ {
+ const FieldVector<double,dim>& quadPos = quadPoint.position();
+
+ const double integrationElement = element.geometry().integrationElement(quadPos);
+// double functionValue = forceTerm(element.geometry().global(quadPos));
+
+ // Evaluate all shape function values at this point
+// std::vector<FieldVector<double,1> > shapeFunctionValues;
+// localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
+
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+
+ std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceGradients);
+
+ // Compute the shape function gradients on the grid element
+ std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+ for (size_t i=0; i<gradients.size(); i++)
+ jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+
+
+ // Actually compute the vector entries
+ for (size_t p=0; p<nSf; p++)
+ {
+ auto localIndex = localView.tree().localIndex(p);
+// elementRhs[localIndex] += shapeFunctionValues[p] * forceTerm(quadPos) * quadPoint.weight() * integrationElement;
+// auto value = shapeFunctionValues[p]* (sqrt(2.0)*mu(quadPos) *forceTerm(quadPos));
+// auto value = -1.0*gradients[p][0] * (sqrt(2.0)*mu(quadPos) *forceTerm(quadPos)); //NEGATIVE!!! TODO
+// auto value = -2.0*mu(quadPos)*(sqrt(2.0)*forceTerm(quadPos))*gradients[p][0] ;
+// auto value = -1.0*mu(quadPos)*forceTerm(quadPos)*gradients[p][0] ;
+
+
+
+
+// auto value = -2.0*sqrt(2.0)*mu(quadPos)*forceTerm(quadPos)*gradients[p][0]; //#1
+
+
+ auto value = 2.0*sqrt(2.0)*mu(quadPos)*forceTerm(quadPos)*gradients[p][0]; // TEST
+
+
+
+
+
+// auto value = 2.0*mu(quadPos)*sqrt(2.0)*forceTerm(quadPos)*gradients[p][0];
+// auto value = -2.0*mu(quadPos)*sqrt(2.0)*quadPos[1]*gradients[p][0]; //TEST
+// std::cout << "value:" << value << std::endl;
+
+// std::cout << "forceTerm(quadPos):" << forceTerm(quadPos) << std::endl;
+// std::cout << "quadPos:" << quadPos << std::endl;
+// std::cout << "integrationElement:" << integrationElement << std::endl;
+
+// auto value = -1.0*sqrt(2.0)*forceTerm(quadPos)*gradients[p][0] ; //TEST
+
+// auto value = -1.0*mu(quadPos)*sqrt(2.0)*forceTerm(quadPos)*forceTerm(quadPos)*gradients[p][0]; //TEST
+
+ elementRhs[localIndex] += value * quadPoint.weight() * integrationElement;
+ }
+ }
+}
+// Get the occupation pattern of the stiffness matrix
+template<class Basis>
+void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
+{
+ nb.resize(basis.size(), basis.size());
+ auto gridView = basis.gridView();
+ // A loop over all elements of the grid
+ auto localView = basis.localView();
+ for (const auto& element : elements(gridView))
+ {
+ localView.bind(element);
+ for (size_t i=0; i<localView.size(); i++)
+ {
+ // The global index of the i-th vertex of the element
+ auto row = localView.index(i);
+ for (size_t j=0; j<localView.size(); j++ )
+ {
+ // The global index of the j-th vertex of the element
+ auto col = localView.index(j);
+ nb.add(row,col);
+ }
+ }
+ }
+}
+/** \brief Assemble the Laplace stiffness matrix on the given grid view */
+template<class Basis, class Matrix, class Vector, class LocalFunction, class Force>
+void assemblePoissonProblem(const Basis& basis,
+ Matrix& matrix,
+ Vector& b,
+ LocalFunction& muLocal,
+ const Force& forceTerm,
+ const double gamma)
+{
+ auto gridView = basis.gridView();
+
+
+ MatrixIndexSet occupationPattern;
+ getOccupationPattern(basis, occupationPattern);
+ occupationPattern.exportIdx(matrix);
+ matrix = 0;
+
+
+ auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
+ auto loadFunctional = localFunction(loadGVF);
+
+ b.resize(basis.dimension());
+ b = 0;
+
+ auto localView = basis.localView();
+
+ for (const auto& element : elements(gridView))
+ {
+
+ localView.bind(element);
+ muLocal.bind(element);
+ loadFunctional.bind(element);
+
+
+// Dune::Matrix<double> elementMatrix;
+ Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
+// Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
+ assembleElementStiffnessMatrix(localView, elementMatrix, muLocal, gamma);
+
+// std::cout << "elementMatrix.N() "<< elementMatrix.N() << std::endl;
+// std::cout << "elementMatrix.M() " << elementMatrix.M() << std::endl;
+
+
+ for(size_t p=0; p<elementMatrix.N(); p++)
+ {
+ auto row = localView.index(p);
+ for (size_t q=0; q<elementMatrix.M(); q++ )
+ {
+ auto col = localView.index(q);
+ matrix[row][col] += elementMatrix[p][q];
+ }
+ }
+
+
+// BlockVector<double> elementRhs;
+ BlockVector<FieldVector<double,1> > elementRhs;
+ assembleElementVolumeTerm(localView, elementRhs, muLocal, loadFunctional);
+ for (size_t p=0; p<elementRhs.size(); p++)
+ {
+ // The global index of the p-th vertex of the element
+ auto row = localView.index(p);
+ b[row] += elementRhs[p];
+ }
+ }
+}
+
+
+
+
+
+template<class Basis, class Vector, class LocalFunc1, class LocalFunc2>
+double computeMuGamma(const Basis& basis,
+ Vector& coeffVector,
+ const double gamma,
+ LocalFunc1& mu,
+ LocalFunc2& Func
+ )
+{
+
+
+// constexpr int dim = Basis::LocalView::Element::dimension;
+
+ double output = 0.0;
+ double Term1 = 0.0;
+ double Term2 = 0.0;
+ double Term11 = 0.0;
+ constexpr int dim = 2;
+// constexpr int dim = 3;
+ auto localView = basis.localView();
+
+// auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<double>(basis, coeffVector);
+// auto localSol = localFunction(solutionFunction);
+
+// auto loadGVF = Dune::Functions::makeGridViewFunction(x3Fun, basis.gridView());
+// auto x3Functional = localFunction(loadGVF);
+
+
+
+
+ double area = 0.0;
+
+ for(const auto& element : elements(basis.gridView()))
+ {
+
+// std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
+ double elementEnergy1 = 0.0;
+ double elementEnergy2 = 0.0;
+
+ localView.bind(element);
+ mu.bind(element);
+ // x3Functional.bind(element);
+ Func.bind(element);
+
+ auto geometry = element.geometry();
+ const auto& localFiniteElement = localView.tree().finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
+
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;
+ int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
+ // int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
+ const auto& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR); // TODO WARUM HIER KEINE COLOR ?? ERROR
+
+ // std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
+
+ for(const auto& quadPoint : quad)
+ {
+ const auto& quadPos = quadPoint.position();
+ // const FieldVector<double,dim>& quadPos = quadPoint.position();
+ // std::cout << " quadPOS: " << quadPos << std::endl;
+
+
+ const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+ const double integrationElement = geometry.integrationElement(quadPos);
+
+ area += quadPoint.weight() * integrationElement;
+
+ std::vector< FieldMatrix<double, 1, dim>> referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceGradients);
+
+ // Compute the shape function gradients on the grid element
+ std::vector<FieldVector<double,dim>> gradients(referenceGradients.size());
+
+ for (size_t i=0; i<gradients.size(); i++)
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+ FieldVector<double,dim> tmp(0.0);
+ // std::cout << "integrationElement :" << integrationElement << std::endl;
+ // std::cout << "quadPoint.weight() :" << quadPoint.weight() << std::endl;
+
+ for (size_t i=0; i < nSf; i++)
+ {
+ size_t localIdx = localView.tree().localIndex(i); // hier i:leafIdx
+ size_t globalIdx = localView.index(localIdx);
+ // printvector(std::cout, gradients[i], "gradients[i]", "--" );
+ // std::cout << "coeffVector[globalIdx]:" << coeffVector[globalIdx] << std::endl;
+ tmp[0] += coeffVector[globalIdx]*gradients[i][0];
+ tmp[1] += coeffVector[globalIdx]*gradients[i][1];
+ // tmp[0] += coeffVector[globalIdx]*gradients[i][0]; // 3D-Version
+ // tmp[1] += coeffVector[globalIdx]*gradients[i][2];
+ // printvector(std::cout, tmp, "tmp", "--" );
+ }
+ // printvector(std::cout, tmp, "gradient_w", "--" );
+
+ // auto value1 = std::pow( ((quadPos[1]/sqrt(2.0))+ (tmp[0]/2.0)) ,2); //TEST
+
+
+
+ auto q1 = (Func(quadPos)/sqrt(2.0)); //#1
+ auto q2 = (tmp[0]/2.0); //#1
+
+
+// auto q2 = (tmp[0]/sqrt(2.0)); // TEST !!!!!!!!!!
+
+
+
+
+
+ // CHECK : BEITRÄGE CHECKEN!!!!
+
+// std::cout << "Beitrag1: " << q2 << std::endl;
+// std::cout << "Beitrag2: " << (tmp[1]/(2.0*gamma)) << std::endl;
+//
+//
+//
+
+ auto q3 = q1 - q2; // TEST
+
+// auto q3 = q1 + q2; // #1
+ auto value1 = std::pow(q3,2);
+
+
+
+// auto value2 = std::pow( (tmp[1]/(2.0*gamma) ) , 2);
+ auto value2 = std::pow( (tmp[1]/(sqrt(2.0)*gamma) ) , 2);
+
+
+
+// auto value = 2.0*mu(quadPos)*(2.0*value1 + value2);
+
+
+
+ elementEnergy1 += 2.0*mu(quadPos)* 2.0*value1 * quadPoint.weight() * integrationElement;
+ elementEnergy2 += 2.0*mu(quadPos)* value2 * quadPoint.weight() * integrationElement;
+ // std::cout << "output:" << output << std::endl;
+ }
+// std::cout << "elementEnergy:" << elementEnergy << std::endl;
+
+ Term1 += elementEnergy1;
+ Term2 += elementEnergy2;
+
+// std::cout << "output: " << output << std::endl;
+ }
+ std::cout << "Term1: " << Term1 << std::endl;
+ std::cout << "Term2: " << Term2 << std::endl;
+ output = Term1 + Term2;
+ std::cout << "output: " << output << std::endl;
+ std::cout << "Domain-Area: " << area << std::endl;
+ return (1.0/area) * output;
+}
+
+
+// // -----------------------------------------------------------
+/*
+template<class Basis, class Vector, class LocalFunc1, class LocalFunc2>
+double computeMuGamma(const Basis& basis,
+ Vector& coeffVector,
+ const double gamma,
+ LocalFunc1& mu,
+ LocalFunc2& Func
+ )
+{
+
+
+// constexpr int dim = Basis::LocalView::Element::dimension;
+
+ double output = 0.0;
+ constexpr int dim = 2;
+// constexpr int dim = 3;
+ auto localView = basis.localView();
+
+// auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<double>(basis, coeffVector);
+// auto localSol = localFunction(solutionFunction);
+
+// auto loadGVF = Dune::Functions::makeGridViewFunction(x3Fun, basis.gridView());
+// auto x3Functional = localFunction(loadGVF);
+
+
+
+
+ double area = 0.0;
+
+ for(const auto& element : elements(basis.gridView()))
+ {
+
+// std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
+ double elementEnergy = 0.0;
+
+ localView.bind(element);
+ mu.bind(element);
+ // x3Functional.bind(element);
+ Func.bind(element);
+
+ auto geometry = element.geometry();
+ const auto& localFiniteElement = localView.tree().finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
+
+ // int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;
+ int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
+ // int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
+ const auto& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR); // TODO WARUM HIER KEINE COLOR ?? ERROR
+
+ // std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
+
+ for(const auto& quadPoint : quad)
+ {
+
+ const auto& quadPos = quadPoint.position();
+ // const FieldVector<double,dim>& quadPos = quadPoint.position();
+ // std::cout << " quadPOS: " << quadPos << std::endl;
+
+
+ const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+ const double integrationElement = geometry.integrationElement(quadPos);
+
+ area += quadPoint.weight() * integrationElement;
+
+ std::vector< FieldMatrix<double, 1, dim>> referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceGradients);
+
+ // Compute the shape function gradients on the grid element
+ std::vector<FieldVector<double,dim>> gradients(referenceGradients.size());
+
+ for (size_t i=0; i<gradients.size(); i++)
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+// FieldVector<double,dim> tmp = {0.0 , 0.0};
+ FieldVector<double,dim> tmp(0.0);
+ // FieldVector<double,dim> tmp = {0.0 ,0.0, 0.0}; //3D-Version
+
+ // std::cout << "integrationElement :" << integrationElement << std::endl;
+ // std::cout << "quadPoint.weight() :" << quadPoint.weight() << std::endl;
+
+ for (size_t i=0; i < nSf; i++)
+ {
+ size_t localIdx = localView.tree().localIndex(i); // hier i:leafIdx
+ size_t globalIdx = localView.index(localIdx);
+
+ // printvector(std::cout, gradients[i], "gradients[i]", "--" );
+ // std::cout << "coeffVector[globalIdx]:" << coeffVector[globalIdx] << std::endl;
+
+ tmp[0] += coeffVector[globalIdx]*gradients[i][0];
+ tmp[1] += coeffVector[globalIdx]*gradients[i][1];
+
+ // tmp[0] += coeffVector[globalIdx]*gradients[i][0]; // 3D-Version
+ // tmp[1] += coeffVector[globalIdx]*gradients[i][2];
+ // printvector(std::cout, tmp, "tmp", "--" );
+
+ }
+ // printvector(std::cout, tmp, "gradient_w", "--" );
+
+
+ // auto value1 = std::pow( ((quadPos[1]/sqrt(2.0))+ (tmp[0]/2.0)) ,2); //TEST
+
+
+
+ auto q1 = (Func(quadPos)/sqrt(2.0));
+ auto q2 = (tmp[0]/2.0);
+
+
+// auto q2 = (tmp[0]/sqrt(2.0)); // TEST !!!!!!!!!!
+
+
+
+
+
+ // CHECK : BEITRÄGE CHECKEN!!!!
+
+ std::cout << "Beitrag1: " << q2 << std::endl;
+ std::cout << "Beitrag2: " << (tmp[1]/(2.0*gamma)) << std::endl;
+
+
+
+
+ auto q3 = q1 + q2;
+ auto value1 = std::pow(q3,2);
+// auto value1 = std::pow( ((Func(quadPos)/sqrt(2.0)) + (tmp[0]/2.0)) , 2);
+
+
+ auto value2 = std::pow( (tmp[1]/(2.0*gamma) ) , 2);
+// auto value2 = std::pow( (tmp[1]/(sqrt(2.0)*gamma) ) , 2); //TEST
+
+
+ // auto value2 = (1.0/(std::pow(gamma,2)))* std::pow(tmp[1],2)/4.0 ; //TEST
+
+ auto value = 2.0*mu(quadPos)*(2.0*value1 + value2);
+
+ // std::cout << "quadPos[1]:" << quadPos[1]<< std::endl;
+ // std::cout << "Func(quadPos):" << Func(quadPos) << std::endl;
+ // std::cout << "sqrt(2.0):" << sqrt(2.0) << std::endl;
+ // std::cout << "tmp[0]:" << tmp[0] << std::endl;
+ // std::cout << "tmp[1]:" << tmp[1] << std::endl;
+ // std::cout << "value1:" << value1 << std::endl;
+ // std::cout << "value2:" << value2 << std::endl;
+ // std::cout << "value:" << value << std::endl;
+
+
+ // auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2) + std::pow( (tmp[1]/(2.0*gamma) ) ,2) );
+
+
+
+ // auto value = 2.0*mu(quadPos)*(2.0* (((x3Functional(quadPos)*x3Functional(quadPos))/2.0) + std::pow( (tmp[0]/2.0) ,2)) + std::pow( (tmp[1]/(2.0*gamma) ) ,2) ); //TEST
+
+ // auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2) ) + std::pow( (tmp[1]/(2.0*gamma) ) ,2) ; //TEST
+ // auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2)) + (1.0/gamma)*std::pow( (tmp[1]/2.0) ,2) ; //TEST
+ // auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2) + (1.0/gamma)*std::pow( (tmp[1]/2.0) ,2) ) ; //TEST
+
+ elementEnergy += value * quadPoint.weight() * integrationElement;
+ // std::cout << "output:" << output << std::endl;
+ }
+// std::cout << "elementEnergy:" << elementEnergy << std::endl;
+ output += elementEnergy;
+// std::cout << "output: " << output << std::endl;
+ }
+ std::cout << "Domain-Area: " << area << std::endl;
+ return (1.0/area) * output;
+}
+// -------------------------------------------------------------------------
+*/
+
+
+
+
+
+
+ // Check whether two points are equal on R/Z x R/Z
+ auto equivalent = [](const FieldVector<double,2>& x, const FieldVector<double,2>& y)
+ {
+ return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
+ and (FloatCmp::eq(x[1],y[1])) );
+ };
+
+
+// // Check whether two points are equal on R/Z x R/Z x R
+// auto equivalent3D = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
+// {
+// return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
+// and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
+// and (FloatCmp::eq(x[2],y[2]))
+// );
+// };
+//
+//
+
+
+
+
+int main(int argc, char *argv[])
+{
+
+ MPIHelper::instance(argc, argv);
+
+
+
+ ParameterTree parameterSet;
+ if (argc < 2)
+ ParameterTreeParser::readINITree("../../inputs/cellsolver.parset", parameterSet);
+ else
+ {
+ ParameterTreeParser::readINITree(argv[1], parameterSet);
+ ParameterTreeParser::readOptions(argc, argv, parameterSet);
+ }
+
+ //////////////////////////////////
+ // Generate the grid
+ //////////////////////////////////
+ constexpr int dim = 2;
+// constexpr int dim = 3; //TEST
+
+
+ // QUAD-GRID
+ FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0});
+ FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0});
+// std::array<int, dim> nElements = {128 ,128};
+ std::array<int, dim> nElements = {64,64};
+// std::array<int, dim> nElements = {32 ,32};
+// std::array<int, dim> nElements = {16,16};
+// std::array<int, dim> nElements = {8,8};
+
+ /*
+ FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0});
+ FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0});
+ std::array<int, dim> nElements = {16,16,16};*/
+
+ using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+ CellGridType grid_CE(lower,upper,nElements);
+ using GridView = CellGridType::LeafGridView;
+ const GridView gridView = grid_CE.leafGridView();
+
+
+ using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+
+
+
+ double gamma = parameterSet.get<double>("gamma", 1.0);
+ double theta = parameterSet.get<double>("theta", 1.0/4.0);
+ double mu1 = parameterSet.get<double>("mu1", 10.0);
+ double beta = parameterSet.get<double>("beta", 2.0);
+ double mu2 = beta*mu1;
+
+ std::cout << "Gamma:" << gamma << std::endl;
+ std::cout << "Theta:" << theta << std::endl;
+ std::cout << "mu1:" << mu1 << std::endl;
+ std::cout << "mu2:" << mu2 << std::endl;
+ std::cout << "beta:" << beta << std::endl;
+
+ auto muTerm = [mu1, mu2, theta] (const Domain& z) {
+
+// if (abs(z[0]) > (theta/2.0))
+ if (abs(z[0]) >= (theta/2.0))
+ return mu1;
+ else
+ return mu2;
+ };
+
+
+// auto materialImp = IsotropicMaterialImp();
+// auto muTermT = materialImp.getMu(parameterSet);
+
+ auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView);
+ auto muLocal = localFunction(muGridF);
+
+
+ /////////////////////////////////////////////////////////
+ // Stiffness matrix and right hand side vector
+ /////////////////////////////////////////////////////////
+// using Matrix = BCRSMatrix<double>;
+// using Vector = BlockVector<double>;
+ using Vector = BlockVector<FieldVector<double,1> >;
+ using Matrix = BCRSMatrix<FieldMatrix<double,1,1> >;
+
+ Matrix stiffnessMatrix;
+ Vector b;
+
+ /////////////////////////////////////////////////////////
+ // Assemble the system
+ /////////////////////////////////////////////////////////
+ using namespace Functions::BasisFactory;
+ Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
+
+ // Don't do the following in real life: It has quadratic run-time in the number of vertices.
+ for (const auto& v1 : vertices(gridView))
+ for (const auto& v2 : vertices(gridView))
+ if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
+ periodicIndices.unifyIndexPair({gridView.indexSet().index(v1)}, {gridView.indexSet().index(v2)});
+
+ auto basis = makeBasis(gridView, Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices)); // flatLexicographic()?
+
+
+
+// auto x3Func = [](const FieldVector<double,dim>& x){return x[1];}; //2D-Version
+// auto x3Func = [](const FieldVector<double,dim>& x){return x[2];}; // 3D-Version
+// auto forceTerm = [](const FieldVector<double,dim>& x){return 0.0;}; // gives q3= 2.08333 (close)...
+
+ auto forceTerm = [](const FieldVector<double,dim>& x){return x[1];}; //2D-Version
+// auto forceTermNEG = [](const FieldVector<double,dim>& x){return -1.0*x[1];}; //2D-Version
+// auto forceTerm = [](const FieldVector<double,dim>& x){return x[2];}; // 3D-Version
+
+ assemblePoissonProblem(basis, stiffnessMatrix, b, muLocal, forceTerm, gamma);
+// printmatrix(std::cout, stiffnessMatrix, "StiffnessMatrix", "--");
+// printvector(std::cout, b, "b", "--");
+
+
+
+ /////////////////////////////////////////////////////////////////////////////
+ // Write matrix and load vector to files, to be used in later examples
+ /////////////////////////////////////////////////////////////////////////////
+
+// // { matrix_rhs_writing_begin }
+// std::string baseName = "Dune-PeriodicPoisson";
+// //+ std::to_string(grid->maxLevel()) + "-refinements";
+// storeMatrixMarket(stiffnessMatrix, baseName + "-matrix.mtx");
+// storeMatrixMarket(b, baseName + "-rhs.mtx");
+// // { matrix_rhs_writing_end }
+//
+
+
+
+
+ ///////////////////////////
+ // Compute solution
+ ///////////////////////////
+
+ // { algebraic_solving_begin }
+ // Choose an initial iterate that fulfills the Dirichlet conditions
+ Vector x(basis.size());
+ x = b;
+
+ // Turn the matrix into a linear operator
+// MatrixAdapter<Matrix,Vector,Vector> linearOperator(stiffnessMatrix);
+//
+// // Sequential incomplete LU decomposition as the preconditioner
+// SeqILU<Matrix,Vector,Vector> preconditioner(stiffnessMatrix, 1.0); // Relaxation factor
+// // Preconditioned conjugate gradient solver
+// CGSolver<Vector> cg(linearOperator,
+// preconditioner,
+// 1e-5, // Desired residual reduction factor
+// 50,
+// // Maximum number of iterations
+// 2);
+// // Verbosity of the solver
+// // Object storing some statistics about the solving process
+// InverseOperatorResult statistics;
+
+
+ std::cout << "------------ CG - Solver ------------" << std::endl;
+ MatrixAdapter<Matrix, Vector, Vector> op(stiffnessMatrix);
+
+
+
+ // Sequential incomplete LU decomposition as the preconditioner
+ SeqILU<Matrix, Vector, Vector> ilu0(stiffnessMatrix,1.0);
+ int iter = parameterSet.get<double>("iterations_CG", 999);
+ // Preconditioned conjugate-gradient solver
+ CGSolver<Vector> solver(op,
+ ilu0, //NULL,
+ 1e-8, // desired residual reduction factorlack
+ iter, // maximum number of iterations
+ 2); // verbosity of the solver
+ InverseOperatorResult statistics;
+ // Solve!
+// cg.apply(x, b, statistics);
+ solver.apply(x, b, statistics);
+ // { algebraic_solving_end }
+
+
+// std::cout << "------------ GMRES - Solver ------------" << std::endl;
+// // Turn the matrix into a linear operator
+// MatrixAdapter<Matrix, Vector, Vector> op(stiffnessMatrix);
+//
+// // Fancy (but only) way to not have a preconditioner at all
+// Richardson<Vector,Vector> preconditioner(1.0);
+//
+// // Construct the iterative solver
+// RestartedGMResSolver<Vector> solver(
+// op, // Operator to invert
+// preconditioner, // Preconditioner
+// 1e-10, // Desired residual reduction factor
+// 500, // Number of iterations between restarts,
+// // here: no restarting
+// 500, // Maximum number of iterations
+// 2); // Verbosity of the solver
+//
+// // Object storing some statistics about the solving process
+// InverseOperatorResult statistics;
+//
+// // solve for different Correctors (alpha = 1,2,3)
+// solver.apply(x, b, statistics);
+//
+
+
+
+
+
+
+
+
+
+
+
+// using SolutionRange = FieldVector<double,dim>;
+// using SolutionRange = FieldVector<double,1>;
+ using SolutionRange = double;
+ auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(basis, x);
+
+
+
+// printvector(std::cout, x, "coeffVector", "--" );
+ std::cout << "Gamma:" << gamma << std::endl;
+
+ auto ForceGridF = Dune::Functions::makeGridViewFunction(forceTerm, gridView);
+ auto ForceLocal = localFunction(ForceGridF);
+
+
+
+// double mu_gamma = computeMuGamma(basis, x, gamma, muLocal, x3Func);
+// double mu_gamma = computeMuGamma(basis, x, gamma, muLocal, forceTerm);
+ double mu_gamma = computeMuGamma(basis, x, gamma, muLocal, ForceLocal);
+ std::cout << "mu_gamma:" << mu_gamma << std::endl;
+
+
+ std::cout.precision(10);
+ std::cout << "mu_gamma:" << std::fixed << mu_gamma << std::endl;
+
+
+
+ // TEST "computeMuGamma"
+// Vector zeroVec;
+// zeroVec.resize(Basis_CE.size());
+ Vector zeroVec(basis.size());
+ zeroVec = 0;
+
+ std::cout << "TEST computeMuGamma:" << computeMuGamma(basis, zeroVec, gamma, muLocal, ForceLocal)<< std::endl;
+
+ std::cout << " --- print analytic solutions --- " << std::endl;
+
+ double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
+ double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
+
+
+ double hm = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0);
+ double gm = mu1*((1-theta)+theta*beta) *(1.0/6.0);
+
+ std::cout << "q1 : " << q1 << std::endl;
+ std::cout << "q2 : " << q2 << std::endl;
+ std::cout << "q3 should be between q1 and q2" << std::endl;
+
+// std::cout << "hm : " << hm << std::endl;
+// std::cout << "gm : " << gm << std::endl;
+
+
+ if(mu_gamma > q2)
+ {
+ std::cout << "mu_gamma > q2!!.. (39) not satisfied" << std::endl;
+ }
+
+
+
+ std::cout << "beta : " << beta << std::endl;
+ std::cout << "beta : " << beta << std::endl;
+ std::cout << "theta : " << theta << std::endl;
+ std::cout << "Gamma : " << gamma << std::endl;
+ std::cout << "mu_gamma:" << mu_gamma << std::endl;
+
+
+
+ // Output result
+
+ VTKWriter<GridView> vtkWriter(gridView);
+// vtkWriter.addVertexData(x, "solution"); //--- Anpassen für P2
+ vtkWriter.addVertexData(
+ solutionFunction,
+ VTK::FieldInfo("solution", VTK::FieldInfo::Type::scalar, 1));
+// VTK::FieldInfo("solution", VTK::FieldInfo::Type::vector, dim));
+
+ vtkWriter.write("Dune-Poisson-Result");
+
+
+
+ using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+ VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0},{64,64});
+ using GridViewVTK = VTKGridType::LeafGridView;
+ const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
+
+ auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
+
+ std::vector<double> mu_CoeffP0;
+ Functions::interpolate(scalarP0FeBasis, mu_CoeffP0, muTerm);
+ auto mu_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, mu_CoeffP0);
+
+ VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
+
+ MaterialVtkWriter.addCellData(
+ mu_DGBF_P0,
+ VTK::FieldInfo("mu_P0", VTK::FieldInfo::Type::scalar, 1));
+
+ MaterialVtkWriter.write("MaterialFunctions-level" );
+ std::cout << "wrote data to file: MaterialFunctions" << std::endl;
+
+}