diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index dde5d76f8718c72696906f8461f0165fbcaf5c2e..2db0ee0d841185ceb6e75f8b7b7f4f9e15dd6878 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -82,7 +82,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
// using Domain = typename LocalView::GridView::Codim<0>::Geometry::GlobalCoordinate;
// using FuncScalar = std::function< ScalarRT(const Domain&) >;
-// using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+// using Func2Tensor = std::function< Matload_alpha1 ,rixRT(const Domain&) >;
elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
@@ -638,6 +638,9 @@ void assembleCellStiffness(const Basis& basis,
Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
+
+
+
// std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
// std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
@@ -716,7 +719,7 @@ void assembleCellLoad(const Basis& basis,
loadFunctional.bind(element);
const int localPhiOffset = localView.size();
- std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
// BlockVector<double> elementRhs;
@@ -739,6 +742,116 @@ void assembleCellLoad(const Basis& basis,
}
}
+
+
+
+
+
+
+
+
+
+
+
+template<class Basis, class LocalFunction1, class LocalFunction2, class MatrixFunction>
+auto energy(const Basis& basis,
+ LocalFunction1& muLocal,
+ LocalFunction2& lambdaLocal,
+ const MatrixFunction& matrixFieldFuncA,
+ const MatrixFunction& matrixFieldFuncB)
+ {
+
+ auto energy = 0.0;
+ constexpr int dim = 3;
+ constexpr int nCompo = 3;
+
+ auto localView = basis.localView();
+
+ auto matrixFieldAGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncA, basis.gridView());
+ auto matrixFieldA = localFunction(matrixFieldAGVF);
+ 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, nCompo, nCompo>;
+
+
+ for (const auto& e : elements(basis.gridView()))
+ {
+ localView.bind(e);
+ matrixFieldA.bind(e);
+ matrixFieldB.bind(e);
+ muLocal.bind(e);
+ lambdaLocal.bind(e);
+
+
+ FieldVector<double,1> elementEnergy(0);
+
+ auto geometry = e.geometry();
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+ int order = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), order);
+
+ for (size_t pt=0; pt < quad.size(); pt++) {
+ const Domain& quadPos = quad[pt].position();
+ const double integrationElement = geometry.integrationElement(quadPos);
+
+
+ auto strain1 = matrixFieldA(quadPos);
+
+ // St.Venant-Kirchhoff stress
+ MatrixRT stressU(0);
+ stressU.axpy(2*muLocal(quadPos), strain1); //this += 2mu *strainU // eigentlich this += 2mu *strain1 ?
+
+ double trace = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ trace += strain1[ii][ii];
+
+ for (int ii=0; ii<nCompo; ii++)
+ stressU[ii][ii] += lambdaLocal(quadPos) * trace;
+
+ auto strain2 = matrixFieldB(quadPos);
+
+ // Local energy density: stress times strain
+ double energyDensity = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ energyDensity += stressU[ii][jj] * strain2[ii][jj];
+
+ elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
+ }
+ energy += elementEnergy;
+ }
+ return energy;
+ }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -1006,7 +1119,7 @@ auto integralMean(const Basis& basis,
area += quadPoint.weight() * integrationElement;
}
}
- std::cout << "Domain-Area: " << area << std::endl;
+// std::cout << "Domain-Area: " << area << std::endl;
for(const auto& element : elements(basis.gridView()))
@@ -1070,7 +1183,7 @@ auto subtractIntegralMean(const Basis& basis,
for(size_t k=0; k<dim; k++)
{
- std::cout << "Integral-Mean: " << IM[k] << std::endl;
+// std::cout << "Integral-Mean: " << IM[k] << std::endl;
auto idx = childToIndexMap(basis,k);
for ( int i : idx)
@@ -1290,6 +1403,30 @@ Func2Tensor x3G_3 = [] (const Domain& x) {
+/////////////////////////////////////////////////////////////////////// TODO
+// TODO : PrestrainImp.hh
+double theta = 0.5;
+double p1 = 1;
+double p2 = 2;
+
+using std::abs;
+
+
+Func2Tensor B = [] (const Domain& x) {
+ double theta = 0.5;
+ double p1 = 1;
+ double p2 = 2;
+ if (abs(x[0])>= (theta/2) && x[2] > 0)
+ return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
+ if (abs(x[0]) < (theta/2) && x[2] < 0)
+ return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
+ else
+ return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+
+};
+//////////////////////////////////////////////////////////////////////////////
+
+
///////////////////////////////////////////////
@@ -1561,46 +1698,88 @@ printmatrix(std::cout, M_3, "Corrector-Matrix M_1", "--");
auto localPhi_1 = localFunction(correctorFunction_1);
+auto localPhi_2 = localFunction(correctorFunction_2);
+auto localPhi_3 = localFunction(correctorFunction_3);
+subtractIntegralMean(Basis_CE, localPhi_1 , phi_1);
+subtractIntegralMean(Basis_CE, localPhi_2 , phi_2);
+subtractIntegralMean(Basis_CE, localPhi_3 , phi_3);
+////////////////////////////////////////////////////////////////////////
+// CHECK INTEGRAL-MEAN:
+// auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+// Basis_CE,
+// phi_1);
+//
+// auto AVPhi_1 = localFunction(AVFunction_1);
+// auto A = integralMean(Basis_CE, AVPhi_1);
+// for(size_t i=0; i<3; i++)
+// {
+// std::cout << "Integral-Mean (TEST) : " << A[i] << std::endl;
+// }
+////////////////////////////////////////////////////
-subtractIntegralMean(Basis_CE, localPhi_1 , phi_1);
+// auto df = derivative(correctorFunction_1); // does not work :(
-// INTEGRAL-MEAN TEST:
-auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
- Basis_CE,
- phi_1);
-auto AVPhi_1 = localFunction(AVFunction_1);
-auto A = integralMean(Basis_CE, AVPhi_1);
-for(size_t i=0; i<3; i++)
+
+
+
+const std::array<Func2Tensor, 3> x3MatrixBasis = {x3G_1, x3G_2, x3G_3};
+
+const std::array<VectorCT, 3> coeffContainer = {x_1, x_2, x_3};
+
+const std::array<decltype(localPhi_1)*, 3> phiContainer= {&localPhi_1, &localPhi_2, &localPhi_3};
+
+const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
+
+
+
+// compute Q
+
+MatrixRT Q(0);
+
+VectorCT tmp1,tmp2;
+
+FieldVector<double,3> B_hat ;
+
+
+
+for(size_t a = 0; a < 3; a++)
{
-std::cout << "Integral-Mean (TEST) : " << A[i] << std::endl;
+ for(size_t b=0; b < 3; b++)
+ {
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, tmp1 ,x3MatrixBasis[b]); //
+
+ auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] );
+
+ Q[a][b] = coeffContainer[a]*tmp1 + GGterm; // seems symmetric... check positiv definitness?
+ }
}
+printmatrix(std::cout, Q, "Matrix Q", "--");
+for(size_t a = 0; a < 3; a++)
+{
-// FieldVector<double,3> ones = {1.0, 1.0, 1.0};
-//
-// phi_1 - IM*;
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp2 ,B);
+ auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B );
-// // const-function
-// auto f1 = [](const auto& x)
-// {
-// return IM; // does not work..
-// };
-//
-//
-// VectorCT IM_1;
-// interpolate(Basis_CE, IM_1, f1);
+ B_hat[a] = coeffContainer[a]*tmp2 + GGterm;
+
+}
+for(size_t i=0; i<3; i++)
+{
+ std::cout << B_hat[i] << std::endl;
+}