diff --git a/mpp b/mpp
index cfd2f66c53e10c51415ddc429e56a002ad7b97b0..e448af94f8abb4d860b9ce064dd363c663353690 160000
--- a/mpp
+++ b/mpp
@@ -1 +1 @@
-Subproject commit cfd2f66c53e10c51415ddc429e56a002ad7b97b0
+Subproject commit e448af94f8abb4d860b9ce064dd363c663353690
diff --git a/src/elasticity/MainElasticity.cpp b/src/elasticity/MainElasticity.cpp
index bad17c051dce7f649ebef055b39b3ee0fd7a640b..9e190cd2207c5ffda91b22786908f7e3b6cc04f9 100644
--- a/src/elasticity/MainElasticity.cpp
+++ b/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
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
diff --git a/src/elasticity/assemble/EGElasticity.cpp b/src/elasticity/assemble/EGElasticity.cpp
index c4ec4b1086d5e82265d7ab564916186ea4c3fae1..97ca62caef61d563f772bb690bd69560efbf7c08 100644
--- a/src/elasticity/assemble/EGElasticity.cpp
+++ b/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);
+}
diff --git a/src/elasticity/assemble/EGElasticity.hpp b/src/elasticity/assemble/EGElasticity.hpp
index 24f70cd21d9c7be64f6518fc1baf29feae70a826..7b54bc53de4b2f22a6b9b68cb78344f9f7a3a1cc 100644
--- a/src/elasticity/assemble/EGElasticity.hpp
+++ b/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);
+ }
};
diff --git a/src/elasticity/assemble/IElasticity.hpp b/src/elasticity/assemble/IElasticity.hpp
index c8fb19966fc7850c7afe8a112d1fa9c7d8a016fa..3a4dde40781612684bc611adc6fdfe574ce61b67 100644
--- a/src/elasticity/assemble/IElasticity.hpp
+++ b/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);
}
diff --git a/src/elasticity/solvers/IterativePressureSolver.cpp b/src/elasticity/solvers/IterativePressureSolver.cpp
index d3299f74e99b72b31a9ada6c92f4b4d95705f1b6..c84a4faccc97ba6d7456b0985593cd1fcafc471a 100644
--- a/src/elasticity/solvers/IterativePressureSolver.cpp
+++ b/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
diff --git a/src/elasticity/solvers/StaticSolver.cpp b/src/elasticity/solvers/StaticSolver.cpp
index 3602e6b8d6bc0765c170229598b9fad681d665e9..cdaae58767441316d52dc72430fd0eccef62ee2f 100644
--- a/src/elasticity/solvers/StaticSolver.cpp
+++ b/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;
diff --git a/test/elasticity/CMakeLists.txt b/test/elasticity/CMakeLists.txt
index dc581d0f1ed4032bf73eb1528ce7206c7a80ce59..7fa4973f96f096796154e4b1cb68bbb808733ef9 100644
--- a/test/elasticity/CMakeLists.txt
+++ b/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)
diff --git a/test/elasticity/TestPrestress.cpp b/test/elasticity/TestPrestress.cpp
index e0921f19a0debd1ab70a4ce39ae9e1ef5d8876d6..c0c626b057246da6bcc0940d58c99c1b6818ade9 100644
--- a/test/elasticity/TestPrestress.cpp
+++ b/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
+}