diff --git a/src/elasticity/assemble/EGElasticity.cpp b/src/elasticity/assemble/EGElasticity.cpp
index c4ec4b1086d5e82265d7ab564916186ea4c3fae1..588fd82b31baf0ac12a1219bfa1fd080e6cbf37f 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 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..6a0f164d45fcf52c873bf20638bc33cf546e10f3 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 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..14cf3dfc824f14383c4f79ea9b9f42b9a55d0392 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?
/*