diff --git a/src/elasticity/assemble/DGElasticity.cpp b/src/elasticity/assemble/DGElasticity.cpp
index 7c01614f2b9a105a0edf41bab9c44caeee7e24ef..a7c69061c2c222e0fa43c793ef0a1dcef4e21dbd 100644
--- a/src/elasticity/assemble/DGElasticity.cpp
+++ b/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
diff --git a/src/elasticity/assemble/DGElasticity.hpp b/src/elasticity/assemble/DGElasticity.hpp
index 9d501c10ff83d97964bf889f876ed6745cdaf463..3afb52479dc1bfde2c1b5bbbb34415344659e7fb 100644
--- a/src/elasticity/assemble/DGElasticity.hpp
+++ b/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