From 313020fecc0e45751b42177bb89ab7783385cf53 Mon Sep 17 00:00:00 2001
From: "laura.stengel" <laura.stengel@kit.edu>
Date: Wed, 26 Mar 2025 11:59:35 +0100
Subject: [PATCH 1/9] Started with Prestress for EG and fixed error for CG
---
src/elasticity/assemble/EGElasticity.cpp | 67 +++++++++++++++++--
src/elasticity/assemble/EGElasticity.hpp | 11 ++-
src/elasticity/assemble/IElasticity.hpp | 2 +-
.../solvers/IterativePressureSolver.cpp | 14 +++-
test/elasticity/CMakeLists.txt | 4 +-
test/elasticity/TestPrestress.cpp | 13 ++--
6 files changed, 95 insertions(+), 16 deletions(-)
diff --git a/src/elasticity/assemble/EGElasticity.cpp b/src/elasticity/assemble/EGElasticity.cpp
index c4ec4b108..588fd82b3 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);//gleiche wie oben
+
+ 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 ((bc >= ROBIN_BC_EPI) && (bc <= ROBIN_BC_BASE)) {//auch hier gleiche wie oben
+ 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);
+ MixedRowBndValues 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);
@@ -769,8 +819,8 @@ void EGElasticity::prestress(const Vector &u, Vector &pStress) {
r_c(E.PenaltyIndex(), 0, 0) += w * cellMat.TotalDerivative(Fiso, F, Q, Dtheta);
}
}
- Prestress->ClearDirichletValues();
- Prestress->Collect();
+// pStress.ClearDirichletValues();
+// pStress.Collect();
}
double EGElasticity::PressureError(const Vector &u) const {
@@ -1244,3 +1294,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 24f70cd21..7b54bc53d 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 c8fb19966..3a4dde407 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 d3299f74e..6a0f164d4 100644
--- a/src/elasticity/solvers/IterativePressureSolver.cpp
+++ b/src/elasticity/solvers/IterativePressureSolver.cpp
@@ -52,7 +52,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();
@@ -98,3 +98,15 @@ void CalculatePrestress(IElasticity &assemble, Vector &u) {
assemble.InitializePrestress(u);
//u.Clear();
}
+
+//old version
+//void CalculatePrestress(IElasticity &assemble, Vector &u) {
+// int prestressSteps{1};
+// Config::Get("PrestressSteps", prestressSteps);
+//
+// IterativePressureSolver prestressSolver(assemble.GetElasticityProblem(), prestressSteps);
+// prestressSolver.Method(assemble, u);
+// mout << "u after preStressSolver.Method " << u << endl << endl;
+// assemble.InitializePrestress(u);
+// //u.Clear();
+//}
\ No newline at end of file
diff --git a/test/elasticity/CMakeLists.txt b/test/elasticity/CMakeLists.txt
index dc581d0f1..7fa4973f9 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 e0921f19a..14cf3dfc8 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"] = "73";
testConfig["PrestressSteps"] = "5";
testConfig["WithPrestress"] = "true";
testConfig["MechEpsilon"] = "1e-12";
testConfig["NewtonEpsilon"] = "1e-11";
testConfig["NewtonReduction"] = "1e-12";
- testConfig["WithPrestress"] = "true";
//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;
}
@@ -93,7 +98,7 @@ 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"}),
@@ -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?
/*
--
GitLab
From 4bacf9fe5cdd62ea4e50346291f12710dfcc6fb2 Mon Sep 17 00:00:00 2001
From: Christian Wieners <wieners@pde13.cluster.math.kit.edu>
Date: Wed, 26 Mar 2025 14:19:33 +0100
Subject: [PATCH 2/9] i.. Z
---
src/elasticity/MainElasticity.cpp | 1 +
src/elasticity/solvers/StaticSolver.cpp | 10 ++++++++++
test/elasticity/TestPrestress.cpp | 26 +++++++++++++------------
3 files changed, 25 insertions(+), 12 deletions(-)
diff --git a/src/elasticity/MainElasticity.cpp b/src/elasticity/MainElasticity.cpp
index bad17c051..e8f9aad63 100644
--- a/src/elasticity/MainElasticity.cpp
+++ b/src/elasticity/MainElasticity.cpp
@@ -99,6 +99,7 @@ Vector &MainElasticity::Run() {
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/solvers/StaticSolver.cpp b/src/elasticity/solvers/StaticSolver.cpp
index 3602e6b8d..96dc52f62 100644
--- a/src/elasticity/solvers/StaticSolver.cpp
+++ b/src/elasticity/solvers/StaticSolver.cpp
@@ -9,9 +9,19 @@ 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);
+
+ 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();
diff --git a/test/elasticity/TestPrestress.cpp b/test/elasticity/TestPrestress.cpp
index 14cf3dfc8..61e303849 100644
--- a/test/elasticity/TestPrestress.cpp
+++ b/test/elasticity/TestPrestress.cpp
@@ -29,7 +29,7 @@ protected:
testConfig["PressureSteps"] = "1";
testConfig["DGSign"] = "-1";
testConfig["DGPenalty"] = "73";
- testConfig["PrestressSteps"] = "5";
+ testConfig["PrestressSteps"] = "1";
testConfig["WithPrestress"] = "true";
testConfig["MechEpsilon"] = "1e-12";
@@ -83,11 +83,12 @@ 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(
@@ -101,11 +102,12 @@ INSTANTIATE_TEST_SUITE_P(DGBeam, PrestressTest, Values(
*/
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"})
));
@@ -122,4 +124,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
+}
--
GitLab
From 6d9ab165835254210fee61ffab926014494f7f30 Mon Sep 17 00:00:00 2001
From: "laura.stengel" <laura.stengel@kit.edu>
Date: Fri, 28 Mar 2025 10:22:46 +0100
Subject: [PATCH 3/9] Removed debugging output
---
src/elasticity/solvers/StaticSolver.cpp | 8 --------
1 file changed, 8 deletions(-)
diff --git a/src/elasticity/solvers/StaticSolver.cpp b/src/elasticity/solvers/StaticSolver.cpp
index 96dc52f62..cdaae5876 100644
--- a/src/elasticity/solvers/StaticSolver.cpp
+++ b/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();
--
GitLab
From 89c8d63abcbb2c871a08dc9d626d4454f3132a0f Mon Sep 17 00:00:00 2001
From: "laura.stengel" <laura.stengel@kit.edu>
Date: Fri, 28 Mar 2025 10:23:11 +0100
Subject: [PATCH 4/9] added output to see when starting calculations after
prestress
---
src/elasticity/MainElasticity.cpp | 1 +
1 file changed, 1 insertion(+)
diff --git a/src/elasticity/MainElasticity.cpp b/src/elasticity/MainElasticity.cpp
index e8f9aad63..9e190cd22 100644
--- a/src/elasticity/MainElasticity.cpp
+++ b/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);
--
GitLab
From ae1a2ef05681b55901263bf696fe4dd560e07e11 Mon Sep 17 00:00:00 2001
From: "laura.stengel" <laura.stengel@kit.edu>
Date: Fri, 28 Mar 2025 10:24:11 +0100
Subject: [PATCH 5/9] has to be updated again!!! (but first implement
prolongate also for EG and DG)
---
.../solvers/IterativePressureSolver.cpp | 41 +++++++++----------
1 file changed, 20 insertions(+), 21 deletions(-)
diff --git a/src/elasticity/solvers/IterativePressureSolver.cpp b/src/elasticity/solvers/IterativePressureSolver.cpp
index 6a0f164d4..c84a4facc 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);
@@ -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
--
GitLab
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 6/9] 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 588fd82b3..97ca62cae 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
From 3757389322faa11affaf268b8614b4df2370dc4f Mon Sep 17 00:00:00 2001
From: "laura.stengel" <laura.stengel@kit.edu>
Date: Fri, 28 Mar 2025 10:25:16 +0100
Subject: [PATCH 7/9] fixed prestress for DG and added needed (missing) methods
---
src/elasticity/assemble/DGElasticity.cpp | 163 ++++++++++++++++++++++-
src/elasticity/assemble/DGElasticity.hpp | 4 +
2 files changed, 160 insertions(+), 7 deletions(-)
diff --git a/src/elasticity/assemble/DGElasticity.cpp b/src/elasticity/assemble/DGElasticity.cpp
index 7c01614f2..a7c69061c 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 9d501c10f..3afb52479 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
--
GitLab
From ae36a207974b7f522a2ec9378efbd79a3a667434 Mon Sep 17 00:00:00 2001
From: "laura.stengel" <laura.stengel@kit.edu>
Date: Fri, 28 Mar 2025 10:30:03 +0100
Subject: [PATCH 8/9] fixed TestPrestress
---
test/elasticity/TestPrestress.cpp | 34 +++++++++++++++----------------
1 file changed, 16 insertions(+), 18 deletions(-)
diff --git a/test/elasticity/TestPrestress.cpp b/test/elasticity/TestPrestress.cpp
index 61e303849..c0c626b05 100644
--- a/test/elasticity/TestPrestress.cpp
+++ b/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"})
));
--
GitLab
From a7484144824644586fa4d12b85d8249a0f5ed4aa Mon Sep 17 00:00:00 2001
From: "laura.stengel" <laura.stengel@kit.edu>
Date: Fri, 28 Mar 2025 10:31:56 +0100
Subject: [PATCH 9/9] upd ref
---
mpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/mpp b/mpp
index cfd2f66c5..e448af94f 160000
--- a/mpp
+++ b/mpp
@@ -1 +1 @@
-Subproject commit cfd2f66c53e10c51415ddc429e56a002ad7b97b0
+Subproject commit e448af94f8abb4d860b9ce064dd363c663353690
--
GitLab