mpp @ e448af94

-Subproject commit cfd2f66c53e10c51415ddc429e56a002ad7b97b0
+Subproject commit e448af94f8abb4d860b9ce064dd363c663353690

src/elasticity/MainElasticity.cpp

@@ -96,9 +96,11 @@ Vector &MainElasticity::Run() {
       CalculatePrestress(*mechA, *displacement);
     }
     if (dyn == "Static") {
+      mout << "=== Running MainElasticity after Prestress: ===============" << endl;
       elasticityProblem->EvaluationResults(*displacement);
       // Generate IterativePressureSolver in src/elasticity/solvers/IterativePressureSolver.hpp
       IterativePressureSolver mSolver(*elasticityProblem, pressureSteps);
+      *displacement = 0;
       mSolver.Method(*mechA, *displacement);
     } else {
       // Generate ElastodynamicTimeIntegrator in src/elasticity/solvers/DynamicMechanicsSolver.hpp

src/elasticity/assemble/DGElasticity.cpp

@@ -568,6 +568,44 @@ void DGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
   }
 }
 
+//TODO stimmt das so?
+void DGElasticity::TractionStiffness(Matrix &matrix) const {
+  matrix = 0;
+  for (cell c = matrix.cells(); c != matrix.cells_end(); ++c) {
+    DGVectorFieldElement E(matrix, *c);
+    DGSizedRowEntries A_c(matrix, *c, E.shape_size(), *c, E.shape_size());
+
+    for (int f = 0; f < c.Faces(); ++f) {
+      DGVectorFieldFaceElement faceElem(matrix, *c, f); //passt das hier zb?
+      if (matrix.OnBoundary(*c, f)) {
+        int bc = matrix.GetMesh().BoundaryFacePart(c.Face(f));
+
+        if ((bc >= ROBIN_BC_EPI) && (bc <= ROBIN_BC_BASE)) {
+          for (int q = 0; q < faceElem.nQ(); ++q) {
+            VectorField n = faceElem.QNormal(q);
+            auto N = Product(n, n);
+
+            double w = faceElem.QWeight(q);
+            for (int i = 0; i < faceElem.shape_size(); ++i) {
+              for (int k = 0; k < E.Dim(); ++k) {
+                auto phi_ik = faceElem.VectorComponentValue(q, i, k);
+                auto Nphi_ik = N * faceElem.VectorComponentValue(q, i, k);
+                for (int j = 0; j < faceElem.shape_size(); ++j) {
+                  for (int l = 0; l < E.Dim(); ++l) {
+                    auto phi_jl = faceElem.VectorComponentValue(q, j, l);
+                    A_c(i, j, k, l)
+                        += w * eProblem.TractionStiffness(Time(), faceElem.QPoint(q),E.Bnd(f))
+                           * Nphi_ik * phi_jl;
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+}
 
 void DGElasticity::Project(const Vector &fine, Vector &u) const {
   int dlevel = fine.SpaceLevel() - u.SpaceLevel();
@@ -854,14 +892,11 @@ double DGElasticity::H1Error(const Vector &u, const Vector &reference) const {
 }
 
 void DGElasticity::prestress(const Vector &u, Vector &pStress) {
-  Prestress = std::make_unique<Vector>(0.0, u);
-
-  //TODO THis has to be updated
   for (cell c = u.cells(); c != u.cells_end(); ++c) {
     const auto &cellMat = eProblem.GetMaterial(c);
 
     DGVectorFieldElement E(u, *c);
-    DGSizedRowBndValues r_c(*Prestress, *c, E.shape_size());
+    DGSizedRowBndValues r_c(pStress, *c, E.shape_size());
 
     Tensor Q = OrientationMatrix(c.GetData());
 
@@ -880,10 +915,119 @@ void DGElasticity::prestress(const Vector &u, Vector &pStress) {
         }
       }
     }
-  }
+    for (int f = 0; f < c.Faces(); ++f) {
+      double s = cellMat.Isochoric().GetMaxParameter() *
+                 (pow(degree, 2) * penalty) / u.GetMesh().MaxMeshWidth();
+
+      DGVectorFieldFaceElement FE(u, *c, f);
+      if (E.Bnd(f) == 1) { // Dirichlet-RB
+        for (int q = 0; q < FE.nQ(); ++q) {
+          double w = FE.QWeight(q);
+          Point z = FE.QPoint(q);
+          VectorField N = FE.QNormal(q);
+          VectorField U = eProblem.Deformation(Time(), z);
+          VectorField u_diff = FE.VectorValue(q, u) - U;
+          Tensor F = DeformationGradient(FE, q, u);
+          auto Fiso = eProblem.IsochoricPart(F);
+
+          for (int i = 0; i < FE.shape_size(); ++i)
+            for (int k = 0; k < E.Dim(); ++k) {
+              auto u_ik = FE.VectorComponentValue(q, i, k);
+              auto Du_ik = FE.VectorRowGradient(q, i, k);
+              double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
+              double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik, Product(u_diff, N));
+
+              r_c(i, k) -= w * (SN + SN_ik * sign - s * u_diff * u_ik);
+            }
+        }
+
+      } else if (E.Bnd(f) == 199) { // Dirichlet-RB
+        for (int q = 0; q < FE.nQ(); ++q) {
+          double w = FE.QWeight(q);
+          Point z = FE.QPoint(q);
+          VectorField N = FE.QNormal(q);
+          VectorField u_diff = FE.VectorValue(q, u);
+          Tensor F = DeformationGradient(FE, q, u);
+          auto Fiso = eProblem.IsochoricPart(F);
+
+          for (int i = 0; i < FE.shape_size(); ++i)
+            for (int k = 0; k < E.Dim(); ++k) {
+              auto u_ik = FE.VectorComponentValue(q, i, k);
+              auto Du_ik = FE.VectorRowGradient(q, i, k);
+
+              double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
+              double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik, Product(u_diff, N));
+              r_c(i, k) -= w * (SN + SN_ik * sign - s * u_diff * u_ik);
+            }
+        }
+
+      } else if (E.Bnd(f) == -1) {
+        cell cf = pStress.find_neighbour_cell(c, f);
+        if (cf() < c()) continue;
+        int f1 = pStress.find_neighbour_face_id(c.Face(f), cf);
+        DGVectorFieldElement E_nghbr(u, *cf);
+        DGVectorFieldFaceElement FE_nghbr(pStress, *cf, f1);
+        DGSizedRowBndValues r_cf(pStress, *cf, E_nghbr.shape_size());
+
+        s = cellMat.Isochoric().GetMaxParameter() * (pow(degree, 2) * penalty) /
+            u.GetMesh().MaxMeshWidth();
+        for (int q = 0; q < FE.nQ(); ++q) {
+          double w = FE.QWeight(q);
+          Point z = FE.QPoint(q);
+          VectorField N = FE.QNormal(q);
+          // für erstes Face-Element
+          VectorField U = FE.VectorValue(q, u);
+
+          Tensor F = DeformationGradient(FE, q, u);
+          auto Fiso = eProblem.IsochoricPart(F);
+          Tensor inverseFA = mat::active_deformation(Q, FE.VectorValue(q, *gamma));
+          F = F * inverseFA;
+
+          // das gleiche für das zweite Face-Element
+          int q1 = FE.findQPointID(FE_nghbr, z);
+          VectorField u_nghbr = FE_nghbr.VectorValue(q1, u);
+
+          Tensor F_nghbr = DeformationGradient(FE_nghbr, q1, u);
+          auto Fiso_nghbr = eProblem.IsochoricPart(F_nghbr);
+          Tensor inverseFA_nghbr = mat::active_deformation(Q, FE_nghbr.VectorValue(q1, *gamma));
+          F_nghbr = F_nghbr * inverseFA_nghbr;
+
+          for (int i = 0; i < FE.shape_size(); ++i) {
+            for (int k = 0; k < E.Dim(); ++k) {
+              auto u_ik = FE.VectorComponentValue(q, i, k);
+              auto Du_ik = FE.VectorRowGradient(q, i, k);
 
-  Prestress->ClearDirichletValues();
-  Prestress->Collect();
+              double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
+              double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik, Product(U, N));
+
+              double SN_nghbr = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q, Product(u_ik, N));
+              double SN_ik_nghbr = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik,
+                                                                 Product(u_nghbr, N));
+              r_c(i, k) += w * (s * U * u_ik - 0.5 * SN - 0.5 * SN_ik * sign);
+              r_c(i, k) += w * (-s * u_nghbr * u_ik - 0.5 * SN_nghbr + 0.5 * SN_ik_nghbr * sign);
+            }
+          }
+
+          for (int j = 0; j < FE_nghbr.shape_size(); ++j)
+            for (int k = 0; k < E.Dim(); ++k) {
+              auto u_jk = FE_nghbr.VectorComponentValue(q1, j, k);
+              auto Du_jk = FE_nghbr.VectorRowGradient(q1, j, k);
+
+              double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_jk, N));
+              double SN_jk = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q, Du_jk,
+                                                           Product(U, N));
+
+              double SN_nghbr = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q, Product(u_jk, N));
+              double SN_jk_nghbr = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q, Du_jk,
+                                                                 Product(u_nghbr, N));
+
+              r_cf(j, k) += w * (s * u_nghbr * u_jk + 0.5 * SN_nghbr + 0.5 * SN_jk_nghbr * sign);
+              r_cf(j, k) += w * (-s * U * u_jk + 0.5 * SN - 0.5 * SN_jk * sign);
+            }
+        }
+      }
+    }
+  }
 }
 
 
@@ -1278,3 +1422,8 @@ std::pair<double, double> DGElasticity::detF(const Vector &u) const {
 
   return {Jmin, Jmax};
 }
+
+
+void DGElasticity::SetInitialValue(Vector &u) {
+  IElasticity::SetInitialValue(u);
+}
\ No newline at end of file

src/elasticity/assemble/DGElasticity.hpp

@@ -39,6 +39,8 @@ public:
 
   void Jacobi(const cell &c, const Vector &u, Matrix &A) const override;
 
+  void TractionStiffness(Matrix &matrix) const override;
+
   void MassMatrix(Matrix &massMatrix) const override;
 
   void SystemMatrix(Matrix &systemMatrix) const override;
@@ -129,6 +131,8 @@ public:
   }
 
   std::pair<double, double> detF(const Vector &u) const override;
+
+  void SetInitialValue(Vector &u) override;
 };
 
 #endif //DGELASTICITY_HPP

src/elasticity/assemble/EGElasticity.cpp

@@ -344,7 +344,7 @@ void EGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
     auto Fiso = eProblem.IsochoricPart(F);
     Tensor inverseFA = mat::active_deformation(Q, E.VectorValue(q, *gamma));
     F = F * inverseFA;
-    double T = eProblem.ActiveStress(Time(), E.QPoint(q));
+//    double T = eProblem.ActiveStress(Time(), E.QPoint(q));
 
     auto theta = E.PenaltyVectorValue(q);
     auto Dtheta = E.PenaltyVectorGradient(q) * transpose(inverseFA);
@@ -586,6 +586,56 @@ void EGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
   }
 }
 
+void EGElasticity::TractionStiffness(Matrix &matrix) const {
+  matrix = 0;
+  for (cell c = matrix.cells(); c != matrix.cells_end(); ++c) {
+    MixedEGVectorFieldElement E(matrix, *c);
+    MixedRowEntries A_c(matrix, *c);
+
+    for (int f = 0; f < c.Faces(); ++f) {
+      MixedEGVectorFieldFaceElement faceElem(matrix, *c, f);
+      if (matrix.OnBoundary(*c, f)) {
+        int bc = matrix.GetMesh().BoundaryFacePart(c.Face(f));
+
+        if ((bc >= ROBIN_BC_EPI) && (bc <= ROBIN_BC_BASE)) {
+          for (int q = 0; q < faceElem.nQ(); ++q) {
+            VectorField n = faceElem.QNormal(q);
+            auto N = Product(n, n);
+
+            double w = faceElem.QWeight(q);
+            auto theta = faceElem.PenaltyVectorValue(q);
+            auto Ntheta = N * faceElem.PenaltyVectorValue(q);
+            for (int i = 0; i < faceElem.ValueSize(); ++i) {
+              for (int k = 0; k < faceElem.Dim(); ++k) {
+                auto phi_ik = faceElem.VectorComponentValue(q, i, k);
+                auto Nphi_ik = N * faceElem.VectorComponentValue(q, i, k);
+                for (int j = 0; j < faceElem.ValueSize(); ++j) {
+                  for (int l = 0; l < faceElem.Dim(); ++l) {
+                    auto phi_jl = faceElem.VectorComponentValue(q, j, l);
+                    A_c(faceElem.ValueIndex(), faceElem.ValueIndex(),i, j, k, l)
+                        += w * eProblem.TractionStiffness(Time(), faceElem.QPoint(q),E.Bnd(f))
+                           * Nphi_ik * phi_jl;
+                  }
+                }
+                A_c(faceElem.ValueIndex(), faceElem.PenaltyIndex(),i, 0, k, 0)
+                    += w * eProblem.TractionStiffness(Time(), faceElem.QPoint(q),E.Bnd(f))
+                       * Ntheta * phi_ik;
+                A_c(faceElem.PenaltyIndex(), faceElem.ValueIndex(),0, i, 0, k)
+                    += w * eProblem.TractionStiffness(Time(), faceElem.QPoint(q),E.Bnd(f))
+                       * Nphi_ik * theta;
+              }
+            }
+            A_c(faceElem.PenaltyIndex(), faceElem.PenaltyIndex(),0, 0, 0, 0)
+                += w * eProblem.TractionStiffness(Time(), faceElem.QPoint(q),E.Bnd(f))
+                   * Ntheta * theta;
+          }
+
+        }
+      }
+    }
+  }
+}
+
 double EGElasticity::L2Norm(const Vector &u) const {
   double err = 0;
   for (cell c = u.cells(); c != u.cells_end(); ++c) {
@@ -743,14 +793,14 @@ double EGElasticity::EnergyError(const Vector &u) const {
 
 
 void EGElasticity::prestress(const Vector &u, Vector &pStress) {
-  Prestress = std::make_unique<Vector>(0.0, u);
-
   for (cell c = u.cells(); c != u.cells_end(); ++c) {
     const auto &cellMat = eProblem.GetMaterial(c);
 
     MixedEGVectorFieldElement E(u, *c);
-    MixedRowValues r_c(*Prestress, *c);
+    MixedRowValues r_c(pStress, *c);
+
     Tensor Q = OrientationMatrix(c.GetData());
+
     for (int q = 0; q < E.nQ(); ++q) {
       double w = E.QWeight(q);
 
@@ -759,7 +809,7 @@ void EGElasticity::prestress(const Vector &u, Vector &pStress) {
 
       for (int i = 0; i < E.ValueSize(); ++i) {
         for (int k = 0; k < E.Dim(); ++k) {
-          auto F_ik = Q * E.VectorRowGradient(q, i, k);
+          auto F_ik = E.VectorRowGradient(q, i, k);
 
           // Passive material response
           r_c(E.ValueIndex(), i, k) += w * cellMat.TotalDerivative(Fiso, F, Q, F_ik);
@@ -768,9 +818,153 @@ void EGElasticity::prestress(const Vector &u, Vector &pStress) {
       auto Dtheta = E.PenaltyVectorGradient(q);
       r_c(E.PenaltyIndex(), 0, 0) += w * cellMat.TotalDerivative(Fiso, F, Q, Dtheta);
     }
+    for (int f = 0; f < c.Faces(); ++f) {
+      double s = cellMat.Isochoric().GetMaxParameter() *
+                 (pow(degree, 2) * penalty) / u.GetMesh().MaxMeshWidth();
+      MixedEGVectorFieldFaceElement FE(u, *c, f);
+      if (E.Bnd(f) == 1) {
+        for (int q = 0; q < FE.nQ(); ++q) {
+          double w = FE.QWeight(q);
+          VectorField N = FE.QNormal(q);
+          VectorField u_diff = FE.VectorValue(q, u)- eProblem.Deformation(Time(), FE.QPoint(q));
+
+          Tensor F = DeformationGradient(FE, q, u);
+          auto Fiso = eProblem.IsochoricPart(F);
+          Tensor inverseFA = mat::active_deformation(Q, FE.VectorValue(q, *gamma));
+          F = F * inverseFA;
+
+          for (int i = 0; i < FE.ValueSize(); ++i) {
+            for (int k = 0; k < FE.Dim(); ++k) {
+              auto u_ik = FE.VectorComponentValue(q, i, k);
+              auto Du_ik = FE.VectorRowGradient(q, i, k);
+              double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
+              double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik,
+                                                           Product(u_diff, N));
+              r_c(FE.ValueIndex(), i, k) -= w * (SN_ik * sign + SN - s * u_diff * u_ik);
+            }
+          }
+          auto theta = FE.PenaltyVectorValue(q);
+          auto Dtheta = FE.PenaltyVectorGradient(q);
+          double SN_theta = cellMat.TotalDerivative(Fiso, F, Q, Product(theta, N));
+          double SN_ik_Dtheta = cellMat.TotalSecondDerivative(Fiso, F, Q, Dtheta,
+                                                              Product(u_diff, N));
+
+          r_c(FE.PenaltyIndex(), 0, 0) -= w * (SN_ik_Dtheta * sign + SN_theta - s * u_diff * theta);
+
+        }
+      } else if( E.Bnd(f) == 199) {
+        for (int q = 0; q < FE.nQ(); ++q) {
+          double w = FE.QWeight(q);
+          VectorField N = FE.QNormal(q);
+          VectorField u_diff = FE.VectorValue(q, u);
+
+          Tensor F = DeformationGradient(FE, q, u);
+          auto Fiso = eProblem.IsochoricPart(F);
+          Tensor inverseFA = mat::active_deformation(Q, FE.VectorValue(q, *gamma));
+          F = F * inverseFA;
+
+          for (int i = 0; i < FE.ValueSize(); ++i) {
+            for (int k = 0; k < FE.Dim(); ++k) {
+              auto u_ik = FE.VectorComponentValue(q, i, k);
+              auto Du_ik = FE.VectorRowGradient(q, i, k);
+              double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
+              double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik,
+                                                           Product(u_diff, N));
+              r_c(FE.ValueIndex(), i, k) -= w * ( SN_ik * sign  + SN - s * u_diff * u_ik);
+            }
+          }
+          auto theta = FE.PenaltyVectorValue(q);
+          auto Dtheta = FE.PenaltyVectorGradient(q);
+          double SN_theta = cellMat.TotalDerivative(Fiso, F, Q, Product(theta, N));
+          double SN_Dtheta = cellMat.TotalSecondDerivative(Fiso, F, Q, Dtheta,
+                                                           Product(u_diff, N));
+
+          r_c(FE.PenaltyIndex(), 0, 0) -= w *(SN_Dtheta * sign + SN_theta - s *  u_diff * theta);
+
+        }
+      } else if (E.Bnd(f) == -1) {
+        cell cf = pStress.find_neighbour_cell(c, f);
+        if (cf() < c()) continue;
+        int f1 = pStress.find_neighbour_face_id(c.Face(f), cf);
+        MixedEGVectorFieldElement E_nghbr(u, *cf);
+        MixedEGVectorFieldFaceElement FE_nghbr(pStress, *cf, f1);
+        MixedRowValues r_cf(pStress, *cf);
+
+        for (int q = 0; q < FE.nQ(); ++q) {
+          double w = FE.QWeight(q);
+          Point z = FE.QPoint(q);
+          VectorField N = FE.QNormal(q);
+          VectorField U = FE.VectorValue(q, u);
+          Tensor F = DeformationGradient(FE, q, u);
+          auto Fiso = eProblem.IsochoricPart(F);
+          int q1 = FE.FindQPointID(FE_nghbr, z);
+          VectorField u_nghbr = FE_nghbr.VectorValue(q1, u);
+          Tensor F_nghbr = DeformationGradient(FE_nghbr, q1, u);
+          auto Fiso_nghbr = eProblem.IsochoricPart(F_nghbr);
+
+          for (int i = 0; i < FE.ValueSize(); ++i) {
+            for (int k = 0; k < FE.Dim(); ++k) {
+              auto u_ik = FE.VectorComponentValue(q, i, k);
+              auto Du_ik = FE.VectorRowGradient(q, i, k);
+              double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_ik, N));
+              double SN_ik = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik, Product(U, N));
+              double SN_nghbr = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q, Product(u_ik, N));
+              double SN_ik_nghbr = cellMat.TotalSecondDerivative(Fiso, F, Q, Du_ik,
+                                                                 Product(u_nghbr, N));
+              r_c(FE.ValueIndex(), i, k) += w * (s * U * u_ik - 0.5 * SN - 0.5 * SN_ik * sign);
+              r_c(FE.ValueIndex(), i, k) +=
+                  w * (-s * u_nghbr * u_ik - 0.5 * SN_nghbr + 0.5 * SN_ik_nghbr * sign);
+            }
+          }
+          for (int j = 0; j < FE_nghbr.ValueSize(); ++j) {
+            for (int l = 0; l < FE_nghbr.Dim(); ++l) {
+              auto u_jl = FE_nghbr.VectorComponentValue(q1, j, l);
+              auto Du_jl = FE_nghbr.VectorRowGradient(q1, j, l);
+              double SN = cellMat.TotalDerivative(Fiso, F, Q, Product(u_jl, N));
+              double SN_jl = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q, Du_jl,
+                                                           Product(U, N));
+              double SN_nghbr = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q, Product(u_jl, N));
+              double SN_jl_nghbr = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q, Du_jl,
+                                                                 Product(u_nghbr, N));
+              r_cf(FE_nghbr.ValueIndex(), j, l) +=
+                  w * (s * u_nghbr * u_jl + 0.5 * SN_nghbr + 0.5 * SN_jl_nghbr * sign);
+              r_cf(FE_nghbr.ValueIndex(), j, l) +=
+                  w * (-s * U * u_jl + 0.5 * SN - 0.5 * SN_jl * sign);
+            }
+          }
+          auto theta = FE.PenaltyVectorValue(q);
+          auto Dtheta = FE.PenaltyVectorGradient(q);
+          double SN_theta = cellMat.TotalDerivative(Fiso, F, Q, Product(theta, N));
+          double SN_ik_Dtheta = cellMat.TotalSecondDerivative(Fiso, F, Q, Dtheta, Product(U, N));
+          double SN_nghbr_theta = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q,
+                                                          Product(theta, N));
+          double SN_ik_nghbr_Dtheta = cellMat.TotalSecondDerivative(Fiso, F, Q, Dtheta,
+                                                                    Product(u_nghbr, N));
+          r_c(FE.PenaltyIndex(), 0, 0) +=
+              w * (s * U * theta - 0.5 * SN_theta - 0.5 * SN_ik_Dtheta * sign);
+          r_c(FE.PenaltyIndex(), 0, 0) +=
+              w * (-s * u_nghbr * theta - 0.5 * SN_nghbr_theta + 0.5 * SN_ik_nghbr_Dtheta * sign);
+
+          auto theta_nghbr = FE_nghbr.PenaltyVectorValue(q1);
+          auto Dtheta_nghbr = FE_nghbr.PenaltyVectorGradient(q1);
+          double SN_theta_nghbr = cellMat.TotalDerivative(Fiso, F, Q, Product(theta_nghbr, N));
+          double SN_jl_Dtheta_nghbr = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q,
+                                                                    Dtheta_nghbr,
+                                                                    Product(U, N));
+          double SN_nghbr_theta_nghbr = cellMat.TotalDerivative(Fiso_nghbr, F_nghbr, Q,
+                                                                Product(theta_nghbr, N));
+          double SN_jl_nghbr_Dtheta_nghbr = cellMat.TotalSecondDerivative(Fiso_nghbr, F_nghbr, Q,
+                                                                          Dtheta_nghbr,
+                                                                          Product(u_nghbr, N));
+          r_cf(FE_nghbr.PenaltyIndex(), 0, 0) += w * (s * u_nghbr * theta_nghbr +
+                                                      0.5 * SN_nghbr_theta_nghbr +
+                                                      0.5 * SN_jl_nghbr_Dtheta_nghbr * sign);//?
+          r_cf(FE_nghbr.PenaltyIndex(), 0, 0) +=
+              w * (-s * U * theta_nghbr + 0.5 * SN_theta_nghbr - 0.5 * SN_jl_Dtheta_nghbr * sign);
+        }
+      }
+    }
   }
-  Prestress->ClearDirichletValues();
-  Prestress->Collect();
 }
 
 double EGElasticity::PressureError(const Vector &u) const {
@@ -1244,3 +1438,6 @@ std::pair<double, double> EGElasticity::detF(const Vector &u) const {
   return {Jmin, Jmax};
 }
 
+void EGElasticity::SetInitialValue(Vector &u) {
+  IElasticity::SetInitialValue(u);
+}

src/elasticity/assemble/EGElasticity.hpp

@@ -54,6 +54,8 @@ public:
 
   void Jacobi(const cell &c, const Vector &u, Matrix &A) const override;
 
+  void TractionStiffness(Matrix &matrix) const override;
+
   void MassMatrix(Matrix &massMatrix) const override;
 
   void SystemMatrix(Matrix &systemMatrix) const override;
@@ -112,7 +114,14 @@ public:
 
   std::array<double, 4> ChamberVolume(const Vector &u) const override;
 
-  std::pair<double, double> detF(const Vector &u) const override; 
+  std::pair<double, double> detF(const Vector &u) const override;
+
+  void SetInitialValue(Vector &u) override;
+
+  void Scale (double s) const {
+    Material &cellMat = eProblem.GetMaterial();
+    //cellMat.ScalePenalty(s*s);
+  }
   
 };
 

src/elasticity/assemble/IElasticity.hpp

@@ -81,8 +81,8 @@ public:
   using IAssemble::Residual;
 
   double Residual(const Vector &u, Vector &r) const override {
+    r=0;
     r = *Prestress;
-    r = 0;
     for (cell ce = u.cells(); ce != u.cells_end(); ++ce) {
       Residual(ce, u, r);
     }

src/elasticity/solvers/IterativePressureSolver.cpp

@@ -21,7 +21,6 @@ bool IterativePressureSolver::Method(IElasticity &assemble, Vector &u) {
           << endl;
 
   Vector uNew(u);
-
   Initialize(assemble, uNew);
   while (!iteration.IsFinished()) {
     bool converged = Step(assemble, uNew);
@@ -52,7 +51,7 @@ bool IterativePressureSolver::Step(IElasticity &assemble, Vector &u) {
   Config::Get("WithPrestress", prestress);
   std::string activeDeformation;
   Config::Get("activeDeformation", activeDeformation);
-  if (activeDeformation != "" && prestress == false) {
+  if (activeDeformation != "" && !prestress) {
     assemble.UpdateStretch(u);
   }
   iteration.NextTimeStep();
@@ -84,17 +83,29 @@ bool KlotzPressureSolver::Step(IElasticity &assemble, Vector &u) {
   return conv;
 }
 
+//TODO implement EG Prolongarte and use this newer version again!!!
+//void CalculatePrestress(IElasticity &assemble, Vector &u) {
+//  int prestressSteps{1};
+//  Config::Get("PrestressSteps", prestressSteps);
+//  int prestressLevel=0;
+//  Config::Get("PrestressLevel", prestressLevel);
+//  Vector uPrestressLevel(assemble.GetSharedDisc(),prestressLevel);
+//  uPrestressLevel.Accumulate();
+//  IterativePressureSolver prestressSolver(assemble.GetElasticityProblem(), prestressSteps);
+//  prestressSolver.Method(assemble, uPrestressLevel);
+//  LagrangeTransfer lagrangeT(uPrestressLevel,u);
+//  lagrangeT.Prolongate(uPrestressLevel,u);
+//  assemble.InitializePrestress(u);
+//  //u.Clear();
+//}
+
+//old version
 void CalculatePrestress(IElasticity &assemble, Vector &u) {
   int prestressSteps{1};
   Config::Get("PrestressSteps", prestressSteps);
-  int prestressLevel=0;
-  Config::Get("PrestressLevel", prestressLevel);
-  Vector uPrestressLevel(assemble.GetSharedDisc(),prestressLevel);
-  uPrestressLevel.Accumulate();
+
   IterativePressureSolver prestressSolver(assemble.GetElasticityProblem(), prestressSteps);
-  prestressSolver.Method(assemble, uPrestressLevel);
-  LagrangeTransfer lagrangeT(uPrestressLevel,u);
-  lagrangeT.Prolongate(uPrestressLevel,u);
+  prestressSolver.Method(assemble, u);
   assemble.InitializePrestress(u);
   //u.Clear();
-}
+}
\ No newline at end of file

src/elasticity/solvers/StaticSolver.cpp

@@ -9,6 +9,8 @@ void StaticSolver::Initialize(const IElasticity &assemble, Vector &u) {
 
 double StaticSolver::calculateResidualUpdate(const IElasticity &assemble, const Vector &u,
                                              Vector &defect) const {
+  defect = 0;
+
   assemble.Residual(u, defect);
   defect += *residualMatrix * u;
 

test/elasticity/CMakeLists.txt

@@ -19,7 +19,7 @@ add_mpp_test(TestLaplaceElasticity ELASTICITY)
 #add_mpp_test(TestDynamicBoundary ELASTICITY)
 #add_mpp_test(TestOscillation ELASTICITY)
 #add_mpp_test(TestCFLCondition ELASTICITY)
-#add_mpp_test(TestPrestress ELASTICITY)
+add_mpp_test(TestPrestress ELASTICITY)
 add_mpp_test(TestVolume ELASTICITY)
 #add_mpp_test(TestVolumePenalty ELASTICITY)
 
@@ -36,7 +36,7 @@ add_mpi_test(TestQuadraticBeamProblem ELASTICITY)
 #add_mpi_test(TestElasticityBlock ELASTICITY)
 add_mpi_test(TestDynamicBoundary ELASTICITY)
 #add_mpi_test(TestCFLCondition ELASTICITY)
-#add_mpi_test(TestPrestress ELASTICITY)
+add_mpi_test(TestPrestress ELASTICITY)
 add_mpi_test(TestVolume ELASTICITY)
 
 #add_mpi_test(TestOscillation ELASTICITY)

test/elasticity/TestPrestress.cpp

@@ -27,13 +27,14 @@ protected:
     testConfig["MechDynamics"] = "Static";
     testConfig["MechDiscretization"] = GetParam().discType;
     testConfig["PressureSteps"] = "1";
+    testConfig["DGSign"] = "-1";
+    testConfig["DGPenalty"] = "10";
     testConfig["PrestressSteps"] = "5";
     testConfig["WithPrestress"] = "true";
 
     testConfig["MechEpsilon"] = "1e-12";
-    testConfig["NewtonEpsilon"] = "1e-11";
-    testConfig["NewtonReduction"] = "1e-12";
-    testConfig["WithPrestress"] = "true";
+    testConfig["NewtonEpsilon"] = "1e-10";
+    testConfig["NewtonReduction"] = "1e-11";
 
 
     //testConfig["Mesh"] = GetParam().meshName;
@@ -69,10 +70,14 @@ TEST_P(PrestressTest, TestInitialValue) {
   mout << "solution " << solution << endl;
   mout << "num " << cmMain->EvaluateQuantity(solution, "L2") << endl << endl;
   for (row r = solution.rows(); r != solution.rows_end(); ++r) {
-    for (int i = 0; i < r.n(); ++i) {
+    for (int i = 0; i < r.NumberOfDofs(); ++i) {
       EXPECT_NEAR(solution(r, i), 0.0, NEWTON_TOLERANCE);
     }
   }
+  Vector exactSolution(cmMain->GetDisplacement());
+  cmMain->SetDisplacement(true, exactSolution);
+
+  mout << "exact " << exactSolution << endl << endl;
 
 }
 
@@ -84,7 +89,7 @@ INSTANTIATE_TEST_SUITE_P(ConformingBeam, PrestressTest, Values(
     TestPrestressParameter({0, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"}),
     TestPrestressParameter({1, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"})
 ));
-/*
+
 INSTANTIATE_TEST_SUITE_P(DGBeam, PrestressTest, Values(
     //TestPrestressParameter({0, 0, {0.0, 0.0,-0.0005,0.0}, "BenchmarkBeam3DTet"}),
     TestPrestressParameter({0, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "DG"}),
@@ -102,7 +107,7 @@ INSTANTIATE_TEST_SUITE_P(EGBeam, PrestressTest, Values(
     TestPrestressParameter({0, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"}),
     TestPrestressParameter({1, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"})
 ));
- */
+
 
 //TODO Why won't this work?
 /*
@@ -117,4 +122,4 @@ INSTANTIATE_TEST_SUITE_P(OnEllipsoid, PrestressTest, Values(
 int main(int argc, char **argv) {
   MppTest mppTest = MppTestBuilder(argc, argv).WithPPM().WithoutDefaultConfig().WithScreenLogging().WithFileLogging();
   return mppTest.RUN_ALL_MPP_TESTS();
-}
\ No newline at end of file
+}