From bb147190ae6182efa693bb01e9ed2a105d47d037 Mon Sep 17 00:00:00 2001
From: Klaus Boehnlein <klaus.boehnlein@tu-dresden.de>
Date: Tue, 15 Jun 2021 07:34:11 +0200
Subject: [PATCH] found error in local assembly

---
 .kdev4/dune-microstructure.kdev4       |   3 +
 build-cmake/CMakeFiles/CMakeError.log  |  22 ++
 build-cmake/CMakeFiles/CMakeOutput.log |   8 +
 build-cmake/DartConfiguration.tcl      |   4 +-
 src/.dune-microstructure.cc.kate-swp   | Bin 0 -> 203 bytes
 src/dune-microstructure.cc             | 488 +++++++++++++++++++++----
 6 files changed, 445 insertions(+), 80 deletions(-)
 create mode 100644 src/.dune-microstructure.cc.kate-swp

diff --git a/.kdev4/dune-microstructure.kdev4 b/.kdev4/dune-microstructure.kdev4
index 67c936e0..985d1fbf 100644
--- a/.kdev4/dune-microstructure.kdev4
+++ b/.kdev4/dune-microstructure.kdev4
@@ -35,3 +35,6 @@ Project Target=dune-microstructure,src,dune-microstructure
 Use External Terminal=false
 Working Directory=file:///home/klaus/Desktop/DUNE2.8/dune-microstructure/build-cmake/src
 isExecutable=false
+
+[Project]
+VersionControlSupport=kdevgit
diff --git a/build-cmake/CMakeFiles/CMakeError.log b/build-cmake/CMakeFiles/CMakeError.log
index daa832b4..c12aa03d 100644
--- a/build-cmake/CMakeFiles/CMakeError.log
+++ b/build-cmake/CMakeFiles/CMakeError.log
@@ -519,3 +519,25 @@ Determing location of SIONLIB failed:
 Include directory: 
 Library directory: 
 
+Determing location of ARPACK failed:
+Libraries to link against: 
+
+Determing location of ARPACK++ failed:
+Include directory: 
+Libraries to link against: 
+
+Determing location of SIONLIB failed:
+Include directory: 
+Library directory: 
+
+Determing location of ARPACK failed:
+Libraries to link against: 
+
+Determing location of ARPACK++ failed:
+Include directory: 
+Libraries to link against: 
+
+Determing location of SIONLIB failed:
+Include directory: 
+Library directory: 
+
diff --git a/build-cmake/CMakeFiles/CMakeOutput.log b/build-cmake/CMakeFiles/CMakeOutput.log
index 5aebca27..1504fbb9 100644
--- a/build-cmake/CMakeFiles/CMakeOutput.log
+++ b/build-cmake/CMakeFiles/CMakeOutput.log
@@ -862,3 +862,11 @@ Determining location of SuperLU succeeded:
 Include directory: /usr/include/superlu
 Library directory: /usr/lib/x86_64-linux-gnu/libsuperlu.so
 
+Determining location of SuperLU succeeded:
+Include directory: /usr/include/superlu
+Library directory: /usr/lib/x86_64-linux-gnu/libsuperlu.so
+
+Determining location of SuperLU succeeded:
+Include directory: /usr/include/superlu
+Library directory: /usr/lib/x86_64-linux-gnu/libsuperlu.so
+
diff --git a/build-cmake/DartConfiguration.tcl b/build-cmake/DartConfiguration.tcl
index 7ea68934..c442b5b7 100644
--- a/build-cmake/DartConfiguration.tcl
+++ b/build-cmake/DartConfiguration.tcl
@@ -57,9 +57,9 @@ P4UpdateOptions:
 P4UpdateCustom: 
 
 # Generic update command
-UpdateCommand: 
+UpdateCommand: /usr/bin/git
 UpdateOptions: 
-UpdateType: 
+UpdateType: git
 
 # Compiler info
 Compiler: /usr/bin/c++
diff --git a/src/.dune-microstructure.cc.kate-swp b/src/.dune-microstructure.cc.kate-swp
new file mode 100644
index 0000000000000000000000000000000000000000..ab134f664ac0f18ef9713f78c53389e793de0692
GIT binary patch
literal 203
zcmZQzU=Z?7EJ;-eE>A2_aLdd|RWQ;sU|?VniQlmAK&qs;W~iH@%?qa)0#o)VtP2ie
zU|`t800AMca4r)A1B0+@uqOip!(Ijk25$xi21W%0--m&Lfki<<0ZAwTQz!&eD1w23
jK~P^`L8%}kQ>!2&Q%ScVv8V)LFUTQGuEC)oSGfWJJ!d7U

literal 0
HcmV?d00001

diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 927bc607..5fb0bec0 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -54,40 +54,265 @@
 
 using namespace Dune;
 
+
+
 template<class LocalView, class Matrix, class localFunction1, class localFunction2>
-void computeElementStiffnessMatrix(const LocalView& localView,
-//                                     Matrix<FieldMatrix<double,1,1> >& elementMatrix,
+void computeElementStiffnessMatrixOLD(const LocalView& localView,
                                     Matrix& elementMatrix,
                                     const localFunction1& mu,
                                     const localFunction2& lambda,
-                                    const double gamma,
-                                    const double phiOffset
+                                    const double gamma
                                     )
 {
 
     // Local StiffnessMatrix of the form:
     // | phi*phi    m*phi |
     // | phi *m     m*m   |
+
     using Element = typename LocalView::Element;
     const Element element = localView.element();
+    auto geometry = element.geometry();
+
     constexpr int dim = Element::dimension;
     constexpr int nCompo = dim;
 
     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&) >;
 
+
+
+
+    elementMatrix.setSize(localView.size()+3, localView.size()+3);
+    elementMatrix = 0;
+
+
+    // LocalBasis-Offset
+    const int localPhiOffset = localView.size();
+
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();              // Unterscheidung children notwendig?
+    const auto nSf = localFiniteElement.localBasis().size();
+
+
+
+
+    ///////////////////////////////////////////////
+    // Basis for R_sym^{2x2}            // wird nicht als Funktion benötigt da konstant...
+    //////////////////////////////////////////////
+    MatrixRT G_1 {{1.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}};
+    MatrixRT G_3 {{0.0, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+
+
+    std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+    //print:
+    printmatrix(std::cout, basisContainer[0] , "G_1", "--");
+    printmatrix(std::cout, basisContainer[1] , "G_2", "--");
+    printmatrix(std::cout, basisContainer[2] , "G_3", "--");
+    ////////////////////////////////////////////////////
+
+    int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+    const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+
+    for (const auto& quadPoint : quad)
+    {
+
+
+        const auto& quadPos = quadPoint.position();
+        const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
+        const auto integrationElement = geometry.integrationElement(quadPoint.position());
+
+//         std::vector<FieldMatrix<double,1,dim> > referenceGradients;   // old
+        std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
+        localFiniteElement.localBasis().evaluateJacobian(
+                quadPoint.position(),
+                referenceGradients);
+
+        // Compute the shape function gradients on the grid element
+        std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+//         std::vector< VectorRT> gradients(referenceGradients.size());
+
+        for (size_t i=0; i<gradients.size(); i++)
+            jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+
+        for (size_t k=0; k < nCompo; k++)
+            for (size_t l=0; l< nCompo; l++)
+            {
+                for (size_t i=0; i < nSf; i++)
+                for (size_t j=0; j < nSf; j++ )
+                {
+
+                    // (scaled) Deformation gradient of the ansatz basis function
+                    MatrixRT defGradientU(0);
+                    defGradientU[k][0] = gradients[i][0];               // Y
+                    defGradientU[k][1] = gradients[i][1];               //X2
+                    defGradientU[k][2] = (1.0/gamma)*gradients[i][2];   //X3
+
+                    // (scaled)  Deformation gradient of the test basis function
+                    MatrixRT defGradientV(0);
+                    defGradientV[l][0] = gradients[j][0];               // Y
+                    defGradientV[l][1] = gradients[j][1];               //X2
+                    defGradientV[l][2] = (1.0/gamma)*gradients[j][2];   //X3
+
+
+                    // symmetric Gradient (Elastic Strains)
+                    MatrixRT strainU, strainV;
+                    for (int ii=0; ii<nCompo; ii++)
+                    for (int jj=0; jj<nCompo; jj++)
+                    {
+                        strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);         // symmetric gradient
+                        strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
+    //                     strainV[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);         // same ? genügt strainU
+
+                    }
+
+                    // St.Venant-Kirchhoff stress
+                    // < L sym[D_gamma*nabla phi_i], sym[D_gamma*nabla phi_j] >
+                    // stressU*strainV
+                    MatrixRT stressU(0);
+                    stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
+
+                    double trace = 0;
+                    for (int ii=0; ii<nCompo; ii++)
+                    trace += strainU[ii][ii];
+
+                    for (int ii=0; ii<nCompo; ii++)
+                    stressU[ii][ii] += lambda(quadPos) * trace;
+
+                    // 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] * strainU[ii][jj];        // "phi*phi"-part
+                        energyDensity += stressU[ii][jj] * strainV[ii][jj];
+
+               /*     size_t row = localView.tree().child(k).localIndex(i);             // kann auf Unterscheidung zwischen k & l verzichten?!
+                    size_t col = localView.tree().child(l).localIndex(j);   */          // siehe DUNE-book p.394
+    //                 size_t col = localView.tree().child(k).localIndex(j);          // Indizes mit k=l genügen .. Kroenecker-Delta_kl  NEIN???
+                    size_t col = localView.tree().child(k).localIndex(i);             // kann auf Unterscheidung zwischen k & l verzichten?!
+                    size_t row = localView.tree().child(l).localIndex(j);
+
+                    elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
+
+
+                    for( size_t m=0; m<3; m++)
+                    {
+
+                    // < L G_i, sym[D_gamma*nabla phi_j] >
+                    // L G_i* strainV
+
+                    // St.Venant-Kirchhoff stress
+                    MatrixRT stressG(0);
+                    stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
+
+                    double traceG = 0;
+                    for (int ii=0; ii<nCompo; ii++)
+                    traceG += basisContainer[m][ii][ii];
+
+                    for (int ii=0; ii<nCompo; ii++)
+                    stressG[ii][ii] += lambda(quadPos) * traceG;
+
+                    double energyDensityGphi = 0;
+                    for (int ii=0; ii<nCompo; ii++)
+                    for (int jj=0; jj<nCompo; jj++)
+                        energyDensityGphi += stressG[ii][jj] * strainV[ii][jj];             // "m*phi"-part
+
+
+                    auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
+
+                    elementMatrix[row][localPhiOffset+m] += value;
+                    elementMatrix[localPhiOffset+m][row] += value;                 // ---- reicht das ? --- TODO
+
+
+//                     // equivalent? :
+//                     // < L sym[D_gamma*nabla phi_i], G_j >
+//                     double energyDensityPhiG = 0;
+//                     for (int ii=0; ii<nCompo; ii++)
+//                     for (int jj=0; jj<nCompo; jj++)
+//                         energyDensityPhiG += stressU[ii][jj] * basisContainer[m][ii][jj];             // "phi*m"-part
+//
+//                     elementMatrix[localPhiOffset+m][row] += energyDensityPhiG * quadPoint.weight() * integrationElement;
+//
+
+
+                    // St.Venant-Kirchhoff stress
+                    // < L G_alpha, G_alpha >
+                    for(size_t n=0; n<3; n++)
+                    {
+                        double energyDensityGG = 0;
+                        for (int ii=0; ii<nCompo; ii++)
+                            for (int jj=0; jj<nCompo; jj++)
+                                energyDensityGG += stressG[ii][jj] * basisContainer[n][ii][jj];         // "m*m"-part
+
+                        elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
+                    }
+
+
+
+                    }
+
+
+
+                }
+
+
+            }
+
+
+
+
+    }
+
+
+}
+
+
+
+
+
+
+
+
+
+template<class LocalView, class Matrix, class localFunction1, class localFunction2>
+void computeElementStiffnessMatrix(const LocalView& localView,
+                                    Matrix& elementMatrix,
+                                    const localFunction1& mu,
+                                    const localFunction2& lambda,
+                                    const double gamma
+                                    )
+{
+
+    // Local StiffnessMatrix of the form:
+    // | phi*phi    m*phi |
+    // | phi *m     m*m   |
+
+    using Element = typename LocalView::Element;
+    const Element element = localView.element();
     auto geometry = element.geometry();
 
+    constexpr int dim = Element::dimension;
+    constexpr int nCompo = dim;
+
     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&) >;
+
+
+
 
     elementMatrix.setSize(localView.size()+3, localView.size()+3);
     elementMatrix = 0;
 
-    const auto& localFiniteElement = localView.tree().child(0).finiteElement();                                    // Unterscheidung children notwendig?
+
+    // LocalBasis-Offset
+    const int localPhiOffset = localView.size();
+
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();              // Unterscheidung children notwendig?
     const auto nSf = localFiniteElement.localBasis().size();
 
 
@@ -95,8 +320,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
 
     ///////////////////////////////////////////////
     // Basis for R_sym^{2x2}            // wird nicht als Funktion benötigt da konstant...
-    //
-
+    //////////////////////////////////////////////
     MatrixRT G_1 {{1.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}};
     MatrixRT G_3 {{0.0, 0.5, 0.0}, {0.5, 0.0, 0.0}, {0.0, 0.0, 0.0}};
@@ -107,7 +331,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
     printmatrix(std::cout, basisContainer[0] , "G_1", "--");
     printmatrix(std::cout, basisContainer[1] , "G_2", "--");
     printmatrix(std::cout, basisContainer[2] , "G_3", "--");
-
+    ////////////////////////////////////////////////////
 
     int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
     const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
@@ -135,112 +359,215 @@ void computeElementStiffnessMatrix(const LocalView& localView,
             jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
 
 
+        // "phi*phi"-part
+        for (size_t k=0; k < nCompo; k++)
+            for (size_t l=0; l< nCompo; l++)
+            {
+                for (size_t i=0; i < nSf; i++)
+                for (size_t j=0; j < nSf; j++ )
+                {
+
+                    // (scaled) Deformation gradient of the ansatz basis function
+                    MatrixRT defGradientU(0);
+                    defGradientU[k][0] = gradients[i][0];               // Y
+                    defGradientU[k][1] = gradients[i][1];               //X2
+                    defGradientU[k][2] = (1.0/gamma)*gradients[i][2];   //X3
+
+                    // (scaled)  Deformation gradient of the test basis function
+                    MatrixRT defGradientV(0);
+                    defGradientV[l][0] = gradients[j][0];               // Y
+                    defGradientV[l][1] = gradients[j][1];               //X2
+                    defGradientV[l][2] = (1.0/gamma)*gradients[j][2];   //X3
+
+
+                    // symmetric Gradient (Elastic Strains)
+                    MatrixRT strainU, strainV;
+                    for (int ii=0; ii<nCompo; ii++)
+                    for (int jj=0; jj<nCompo; jj++)
+                    {
+                        strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);         // symmetric gradient
+                        strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
+    //                     strainV[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);         // same ? genügt strainU
+
+                    }
+
+                    // St.Venant-Kirchhoff stress
+                    // < L sym[D_gamma*nabla phi_i], sym[D_gamma*nabla phi_j] >
+                    // stressU*strainV
+                    MatrixRT stressU(0);
+                    stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
+
+                    double trace = 0;
+                    for (int ii=0; ii<nCompo; ii++)
+                    trace += strainU[ii][ii];
+
+                    for (int ii=0; ii<nCompo; ii++)
+                    stressU[ii][ii] += lambda(quadPos) * trace;
 
-        for (size_t i=0; i < nSf; i++)
-	      for (size_t j=0; j < nSf; j++ )
-	        for (size_t k=0; k < nCompo; k++)
-// 	          for (size_t l=0; l< nCompo; l++)
-	          {
+                    // 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] * strainU[ii][jj];        // "phi*phi"-part
+                        energyDensity += stressU[ii][jj] * strainV[ii][jj];
 
-                // scale Gradient
-	            MatrixRT defGradientU(0);
-	            defGradientU[k][0] = gradients[i][0];               // Y
-	            defGradientU[k][1] = gradients[i][1];               //X2
-	            defGradientU[k][2] = (1.0/gamma)*gradients[i][2];   //X3
+               /*     size_t row = localView.tree().child(k).localIndex(i);             // kann auf Unterscheidung zwischen k & l verzichten?!
+                    size_t col = localView.tree().child(l).localIndex(j);   */          // siehe DUNE-book p.394
+    //                 size_t col = localView.tree().child(k).localIndex(j);          // Indizes mit k=l genügen .. Kroenecker-Delta_kl  NEIN???
+                    size_t col = localView.tree().child(k).localIndex(i);             // kann auf Unterscheidung zwischen k & l verzichten?!
+                    size_t row = localView.tree().child(l).localIndex(j);
 
+                    elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
 
-	            // symmetric Gradient (Elastic Strains)
-	            MatrixRT strainU, strainV;
-	            for (int ii=0; ii<nCompo; ii++)
-	              for (int jj=0; jj<nCompo; jj++)
-	              {
-	                strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);         // symmetric gradient
-                    strainV[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);         // same ? genügt strainU
-// 	                strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
-	              }
+                }
 
-	            // St.Venant-Kirchhoff stress
-	            // < L sym[D_gamma*nabla phi_i], sym[D_gamma*nabla phi_i] >
-	            MatrixRT stressU(0);
-	            stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
+            }
 
-	            double trace = 0;
-	            for (int ii=0; ii<nCompo; ii++)
-	              trace += strainU[ii][ii];
 
-	            for (int ii=0; ii<nCompo; ii++)
-	              stressU[ii][ii] += lambda(quadPos) * trace;
+        // "m*phi"-part
+        for (size_t l=0; l< nCompo; l++)
+            for (size_t j=0; j < nSf; j++ )
+            {
 
-	              // 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] * strainU[ii][jj];        // "phi*phi"-part
-// 	                energyDensity += stressU[ii][jj] * strainV[ii][jj];
 
-	            size_t row = localView.tree().child(k).localIndex(i);         // kann auf Unterscheidung zwischen k & l verzichten?!
-// 	            size_t col = localView.tree().child(l).localIndex(j);         // siehe DUNE-book p.394
-                size_t col = localView.tree().child(k).localIndex(j);         // Indizes mit k=l genügen .. Kroenecker-Delta_kl
+                // (scaled)  Deformation gradient of the test basis function
+                MatrixRT defGradientV(0);
+                defGradientV[l][0] = gradients[j][0];               // Y
+                defGradientV[l][1] = gradients[j][1];               //X2
+                defGradientV[l][2] = (1.0/gamma)*gradients[j][2];   //X3
 
+                 // symmetric Gradient (Elastic Strains)
+                MatrixRT strainV;
+                for (int ii=0; ii<nCompo; ii++)
+                for (int jj=0; jj<nCompo; jj++)
+                    {
+                        strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
+                    }
 
-	            elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
 
-                for( size_t j=0; j<3; j++)
+                for( size_t m=0; m<3; m++)
                 {
 
-                // < L G_alpha, sym[D_gamma*nabla phi_i] >
-                // St.Venant-Kirchhoff stress
-                MatrixRT stressG(0);
-	            stressG.axpy(2*mu(quadPos), basisContainer[j]); //this += 2mu *strainU
+                    // < L G_i, sym[D_gamma*nabla phi_j] >
+                    // L G_i* strainV
 
-                double traceG = 0;
-	            for (int ii=0; ii<nCompo; ii++)
-	              traceG += basisContainer[j][ii][ii];
+                    // St.Venant-Kirchhoff stress
+                    MatrixRT stressG(0);
+                    stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
 
-	            for (int ii=0; ii<nCompo; ii++)
-	              stressG[ii][ii] += lambda(quadPos) * traceG;
+                    double traceG = 0;
+                    for (int ii=0; ii<nCompo; ii++)
+                    traceG += basisContainer[m][ii][ii];
 
-                double energyDensityGphi = 0;
-	            for (int ii=0; ii<nCompo; ii++)
-	              for (int jj=0; jj<nCompo; jj++)
-                    energyDensityGphi += stressG[ii][jj] * strainU[ii][jj];             // "phi*m"-part
+                    for (int ii=0; ii<nCompo; ii++)
+                    stressG[ii][ii] += lambda(quadPos) * traceG;
 
+                    double energyDensityGphi = 0;
+                    for (int ii=0; ii<nCompo; ii++)
+                    for (int jj=0; jj<nCompo; jj++)
+                        energyDensityGphi += stressG[ii][jj] * strainV[ii][jj];             // "m*phi"-part
+
+                    size_t row = localView.tree().child(l).localIndex(j);
 
-                auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
+                    auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
 
-                elementMatrix[row][phiOffset+j] += value;
-                elementMatrix[phiOffset+j][row] += value;
+                    elementMatrix[row][localPhiOffset+m] += value;
+    //                 elementMatrix[localPhiOffset+m][row] += value;                 // ---- reicht das ? --- TODO
 
+                }
+
+            }
+
+
+
+        // "phi*m"-part
+        for (size_t k=0; k < nCompo; k++)
+            for (size_t i=0; i < nSf; i++)
+            {
+
+                // (scaled) Deformation gradient of the ansatz basis function
+                MatrixRT defGradientU(0);
+                defGradientU[k][0] = gradients[i][0];               // Y
+                defGradientU[k][1] = gradients[i][1];               //X2
+                defGradientU[k][2] = (1.0/gamma)*gradients[i][2];   //X3
+
+
+                // symmetric Gradient (Elastic Strains)
+                MatrixRT strainU;
+                for (int ii=0; ii<nCompo; ii++)
+                for (int jj=0; jj<nCompo; jj++)
+                {
+                    strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);         // symmetric gradient
+                }
 
                 // St.Venant-Kirchhoff stress
-	            // < L G_alpha, G_alpha >
-                for(size_t m=0; m<3; m++)
+                MatrixRT stressU(0);
+                stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
+
+                double trace = 0;
+                for (int ii=0; ii<nCompo; ii++)
+                trace += strainU[ii][ii];
+
+                for (int ii=0; ii<nCompo; ii++)
+                stressU[ii][ii] += lambda(quadPos) * trace;
+
+                for( size_t n=0; n<3; n++)
                 {
-                    double energyDensityGG = 0;
+
+                    // < L sym[D_gamma*nabla phi_i], G_j >
+                    double energyDensityPhiG = 0;
                     for (int ii=0; ii<nCompo; ii++)
-                        for (int jj=0; jj<nCompo; jj++)
-                            energyDensityGG += stressG[ii][jj] * basisContainer[m][ii][jj];         // "m*m"-part
+                    for (int jj=0; jj<nCompo; jj++)
+                        energyDensityPhiG += stressU[ii][jj] * basisContainer[n][ii][jj];             // "phi*m"-part
 
+                    size_t col = localView.tree().child(k).localIndex(i);
 
+                    elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
 
-                    elementMatrix[phiOffset+m][phiOffset+j]= energyDensityGG * quadPoint.weight() * integrationElement;
                 }
 
 
 
-                }
+            }
 
 
 
-	          }
 
+        // "m*m"-part
 
+        for(size_t m=0; m<3; m++)
+            for(size_t n=0; n<3; n++)
+            {
 
-    }
+                // St.Venant-Kirchhoff stress
+                MatrixRT stressG(0);
+                stressG.axpy(2*mu(quadPos), basisContainer[m]); //this += 2mu *strainU
 
-}
+                double traceG = 0;
+                for (int ii=0; ii<nCompo; ii++)
+                traceG += basisContainer[m][ii][ii];
 
+                for (int ii=0; ii<nCompo; ii++)
+                stressG[ii][ii] += lambda(quadPos) * traceG;
 
+                double energyDensityGG = 0;
+
+                for (int ii=0; ii<nCompo; ii++)
+                    for (int jj=0; jj<nCompo; jj++)
+                        energyDensityGG += stressG[ii][jj] * basisContainer[n][ii][jj];         // "m*m"-part
+
+                elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
+
+            }
+
+
+
+
+
+    }
+
+
+}
 
 
 
@@ -520,12 +847,11 @@ void assembleCellProblem(const Basis& basis,
     auto localView = basis.localView();
 //     std::cout << localView.maxSize() << std::endl;
 
-    const int phiOffset = basis.dimension();
+//     const int phiOffset = basis.dimension();
 
 
     // Transform G_alpha's to GridViewFunctions/LocalFunctions  -- why exactly?
     auto loadGVF  = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
-
     auto loadFunctional = localFunction(loadGVF);      // Dune::Functions:: notwendig?
 
 
@@ -539,9 +865,15 @@ void assembleCellProblem(const Basis& basis,
 
         Dune::Matrix<double> elementMatrix;  //eher das ??
 //         Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
-        computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma, phiOffset);
+        computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
         printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
 
+        Dune::Matrix<double> TestelementMatrix;
+        computeElementStiffnessMatrixOLD(localView, TestelementMatrix, muLocal, lambdaLocal, gamma);
+        printmatrix(std::cout, TestelementMatrix, "TESTElementMatrix", "--");
+
+
+
         for (size_t i=0; i<elementMatrix.N(); i++)
         {
             // The global index of the i-th local degree of freedom of the current element
-- 
GitLab