mpp @ e448af94

-Subproject commit cfd2f66c53e10c51415ddc429e56a002ad7b97b0
+Subproject commit e448af94f8abb4d860b9ce064dd363c663353690

src/elasticity/MainElasticity.cpp

@@ -96,6 +96,7 @@ 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);

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

@@ -590,14 +590,14 @@ 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);//gleiche wie oben
+    MixedRowEntries A_c(matrix, *c);
 
     for (int f = 0; f < c.Faces(); ++f) {
       MixedEGVectorFieldFaceElement faceElem(matrix, *c, f);
-      if (matrix.OnBoundary(*c, f)) {//auch hier gleiche wie oben
-        int bc = matrix.GetMesh().BoundaryFacePart(c.Face(f));//auch hier gleiche wie oben
+      if (matrix.OnBoundary(*c, f)) {
+        int bc = matrix.GetMesh().BoundaryFacePart(c.Face(f));
 
-        if ((bc >= ROBIN_BC_EPI) && (bc <= ROBIN_BC_BASE)) {//auch hier gleiche wie oben
+        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);
@@ -797,7 +797,7 @@ void EGElasticity::prestress(const Vector &u, Vector &pStress) {
     const auto &cellMat = eProblem.GetMaterial(c);
 
     MixedEGVectorFieldElement E(u, *c);
-    MixedRowBndValues r_c(pStress, *c);
+    MixedRowValues r_c(pStress, *c);
 
     Tensor Q = OrientationMatrix(c.GetData());
 
@@ -818,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);
+        }
+      }
+    }
   }
-//  pStress.ClearDirichletValues();
-//  pStress.Collect();
 }
 
 double EGElasticity::PressureError(const Vector &u) const {

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);
@@ -84,29 +83,29 @@ bool KlotzPressureSolver::Step(IElasticity &assemble, Vector &u) {
   return conv;
 }
 
-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
+//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, u);
-//  mout << "u after preStressSolver.Method " << u << endl << endl;
+//  prestressSolver.Method(assemble, uPrestressLevel);
+//  LagrangeTransfer lagrangeT(uPrestressLevel,u);
+//  lagrangeT.Prolongate(uPrestressLevel,u);
 //  assemble.InitializePrestress(u);
 //  //u.Clear();
-//}
\ No newline at end of file
+//}
+
+//old version
+void CalculatePrestress(IElasticity &assemble, Vector &u) {
+  int prestressSteps{1};
+  Config::Get("PrestressSteps", prestressSteps);
+
+  IterativePressureSolver prestressSolver(assemble.GetElasticityProblem(), prestressSteps);
+  prestressSolver.Method(assemble, u);
+  assemble.InitializePrestress(u);
+  //u.Clear();
+}
\ No newline at end of file

src/elasticity/solvers/StaticSolver.cpp

@@ -12,16 +12,8 @@ double StaticSolver::calculateResidualUpdate(const IElasticity &assemble, const
   defect = 0;
 
   assemble.Residual(u, defect);
-
-  mout << "Residual u " << endl << u << endl;
-  mout << "Residual d " << endl << defect << endl;
-  mout << "Residual d1 " << endl << defect.norm() << endl;
-  
   defect += *residualMatrix * u;
 
-  mout << "Residual d2 " << endl << defect.norm() << endl;
-
-  
   defect.ClearDirichletValues();
   defect.Collect();
 

test/elasticity/TestPrestress.cpp

@@ -28,13 +28,13 @@ protected:
     testConfig["MechDiscretization"] = GetParam().discType;
     testConfig["PressureSteps"] = "1";
     testConfig["DGSign"] = "-1";
-    testConfig["DGPenalty"] = "73";
-    testConfig["PrestressSteps"] = "1";
+    testConfig["DGPenalty"] = "10";
+    testConfig["PrestressSteps"] = "5";
     testConfig["WithPrestress"] = "true";
 
     testConfig["MechEpsilon"] = "1e-12";
-    testConfig["NewtonEpsilon"] = "1e-11";
-    testConfig["NewtonReduction"] = "1e-12";
+    testConfig["NewtonEpsilon"] = "1e-10";
+    testConfig["NewtonReduction"] = "1e-11";
 
 
     //testConfig["Mesh"] = GetParam().meshName;
@@ -83,14 +83,13 @@ TEST_P(PrestressTest, TestInitialValue) {
 
 INSTANTIATE_TEST_SUITE_P(ConformingBeam, 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", "Conforming"})
-    //,
-    //TestPrestressParameter({1, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"}),
-    //TestPrestressParameter({2, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"}),
-    //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"})
+    TestPrestressParameter({0, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"}),
+    TestPrestressParameter({1, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"}),
+    TestPrestressParameter({2, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "Conforming"}),
+    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"}),
@@ -99,15 +98,14 @@ INSTANTIATE_TEST_SUITE_P(DGBeam, PrestressTest, Values(
     TestPrestressParameter({0, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "DG"}),
     TestPrestressParameter({1, 2, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "DG"})
 ));
-*/
+
 INSTANTIATE_TEST_SUITE_P(EGBeam, 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", "EG"})
-    //,
-    //TestPrestressParameter({1, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"}),
-    //TestPrestressParameter({2, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"}),
-    //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"})
+    TestPrestressParameter({0, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"}),
+    TestPrestressParameter({1, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"}),
+    TestPrestressParameter({2, 1, {10.0, 0.0, 1.0, 0.0}, "PrestressBeam", "EG"}),
+    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"})
 ));