From e5a083dbd9ea3a16d434d90a4ea7769b4cf0f52e Mon Sep 17 00:00:00 2001
From: "laura.stengel" <laura.stengel@kit.edu>
Date: Fri, 28 Mar 2025 10:24:33 +0100
Subject: [PATCH] fixed prestress for EG
---
src/elasticity/assemble/EGElasticity.cpp | 158 ++++++++++++++++++++++-
1 file changed, 151 insertions(+), 7 deletions(-)
diff --git a/src/elasticity/assemble/EGElasticity.cpp b/src/elasticity/assemble/EGElasticity.cpp
index 588fd82b..97ca62ca 100644
--- a/src/elasticity/assemble/EGElasticity.cpp
+++ b/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 {
--
GitLab