fixed prestress for DG and added needed (missing) methods

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