diff --git a/src/elasticity/assemble/DGElasticity.cpp b/src/elasticity/assemble/DGElasticity.cpp
index 4caa463dd1374cef9eac67dda563754052897b4b..73e39bb2a49f5351be54d9dbb1da2924d2f45e4b 100644
--- a/src/elasticity/assemble/DGElasticity.cpp
+++ b/src/elasticity/assemble/DGElasticity.cpp
@@ -132,6 +132,9 @@ void DGElasticity::Residual(const cell &c, const Vector &U, Vector &r) const {
Tensor F = DeformationGradient(FE, q, U);
VectorField N = mat::cofactor(F) * FE.QNormal(q);
+ if (contains(cellMat.Name(),"Linear") || contains(cellMat.Name(),"Laplace")) {
+ N = FE.QNormal(q);
+ }
double localPressure = -eProblem.Pressure(Time(), z, E.Bnd(f));
for (int i = 0; i < FE.shape_size(); ++i)
for (int k = 0; k < E.Dim(); ++k) {
@@ -458,6 +461,7 @@ void DGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
}
}
if (E.Bnd(face) >= PRESSURE_BC_LV && E.Bnd(face) <= PRESSURE_BC_RA) {
+ if (contains(cellMat.Name(),"Linear") || contains(cellMat.Name(),"Laplace")) continue;
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
VectorField N = FE.QNormal(q);
diff --git a/src/elasticity/assemble/EGElasticity.cpp b/src/elasticity/assemble/EGElasticity.cpp
index 6dafa3ed8925733904555f312185d347fd3d551c..ea8c1c4eecb70947ff7696f7b337b69c67da1056 100644
--- a/src/elasticity/assemble/EGElasticity.cpp
+++ b/src/elasticity/assemble/EGElasticity.cpp
@@ -123,6 +123,9 @@ void EGElasticity::Residual(const cell &c, const Vector &U, Vector &r) const {
Point z = FE.QPoint(q);
Tensor F = DeformationGradient(FE, q, U);
VectorField N = mat::cofactor(F) * FE.QNormal(q);
+ if (contains(cellMat.Name(),"Linear") || contains(cellMat.Name(),"Laplace")){
+ N = FE.QNormal(q);
+ }
double localPressure = -eProblem.Pressure(Time(), z, E.Bnd(f));
for (int i = 0; i < FE.ValueSize(); ++i) {
for (int k = 0; k < FE.Dim(); ++k) {
@@ -408,6 +411,7 @@ void EGElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const {
w * (SN_theta * sign + SN_Dtheta - s * (theta * theta));
}
} else if (E.Bnd(face) >= PRESSURE_BC_LV && E.Bnd(face) <= PRESSURE_BC_RA) {
+ if (contains(cellMat.Name(),"Linear") || contains(cellMat.Name(),"Laplace")) continue;
for (int q = 0; q < FE.nQ(); ++q) {
auto theta = FE.PenaltyVectorValue(q);
auto Dtheta = FE.PenaltyVectorGradient(q);
diff --git a/src/elasticity/assemble/LagrangeElasticity.cpp b/src/elasticity/assemble/LagrangeElasticity.cpp
index 10a37bfec23a5d19c07343c530dea1e9286d3b3d..c1e7aaba7778b1ef902180dc9ff58f33fcc598b5 100644
--- a/src/elasticity/assemble/LagrangeElasticity.cpp
+++ b/src/elasticity/assemble/LagrangeElasticity.cpp
@@ -260,11 +260,15 @@ void LagrangeElasticity::PressureBoundary(const cell &c, int face, int bc, const
VectorFieldFaceElement FE(U, *c, face);
RowValues r_c(r, FE);
+ const Material &cellMat = eProblem.GetMaterial(c);
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
Tensor F = DeformationGradient(FE, q, U);
VectorField N = mat::cofactor(F) * FE.QNormal(q);
+ if (contains(cellMat.Name(),"Linear") || contains(cellMat.Name(),"Laplace")) {
+ N = FE.QNormal(q);
+ }
double localPressure = -eProblem.Pressure(Time(), FE.QPoint(q), bc);
for (int i = 0; i < FE.NodalPoints(); ++i) {
@@ -340,7 +344,7 @@ void LagrangeElasticity::Jacobi(const cell &c, const Vector &U, Matrix &A) const
}
}
}
-
+ if (contains(cellMat.Name(),"Linear") || contains(cellMat.Name(),"Laplace")) return;
for (int face = 0; face < c.Faces(); ++face) {
if (E.Bnd(face) >= PRESSURE_BC_LV && E.Bnd(face) <= PRESSURE_BC_RA) {
VectorFieldFaceElement FE(U, *c, face);
diff --git a/src/elasticity/assemble/LagrangeElasticity.hpp b/src/elasticity/assemble/LagrangeElasticity.hpp
index a390c116ec24080a9c2563b82dac60e10950509e..cf527d2bf9ca7053973869a230e6903240d252ec 100644
--- a/src/elasticity/assemble/LagrangeElasticity.hpp
+++ b/src/elasticity/assemble/LagrangeElasticity.hpp
@@ -139,9 +139,9 @@ public:
void Scale (double s) const {
Material &cellMat = eProblem.GetMaterial();
- cellMat.ScalePenalty(s);
+ cellMat.ScalePenalty(s*s);
}
-
+
};
#endif //CARDMECH_LAGRANGEELASTICITY_HPP
diff --git a/src/elasticity/problem/BeamProblems.hpp b/src/elasticity/problem/BeamProblems.hpp
index 54eec3b54d00a0b0c1f2f25eb5c08ae461a8a1cc..0458df06eb54fe1e2f0540526b31b9e6dcf4ddd3 100644
--- a/src/elasticity/problem/BeamProblems.hpp
+++ b/src/elasticity/problem/BeamProblems.hpp
@@ -109,7 +109,7 @@ protected:
if (bndCondition == 230) {
return -(2 * mu + lambda);
} else if (bndCondition == 231 || bndCondition == 232) {
- return -lambda / 2;
+ return -lambda;
} else {
return 0;
}
@@ -145,7 +145,7 @@ protected:
if (bndCondition == 230) {
return (4 * mu + 2 * lambda);
} else if (bndCondition == 231 || bndCondition == 232) {
- return -2 * lambda;
+ return 2 * lambda;
} else {
return 0;
}
diff --git a/src/elasticity/problem/DynamicBoundaryProblems.cpp b/src/elasticity/problem/DynamicBoundaryProblems.cpp
index 5176b99f39779cd1eba3f13988972f4f6bfc1819..e10228d3eb5de7d3d20e56c05bb98f1712a41815 100644
--- a/src/elasticity/problem/DynamicBoundaryProblems.cpp
+++ b/src/elasticity/problem/DynamicBoundaryProblems.cpp
@@ -201,11 +201,11 @@ VectorField DPolynomialNeumannProblem::Load(double time, const Point &x, const C
double DPolynomialNeumannProblem::pressure(double t, const Point &x, int bndCondition) const {
if (bndCondition == 230) {
- return - ((t*t*x[0]) / ((1 + t*t*x[1]) * (1 + t*t*x[2])));
+ return - (t*t*x[0]);
} else if (bndCondition == 231) {
- return -((t*t*x[1]) / ((1 + t*t* x[0]) * (1 + t*t*x[2])));
+ return -(t*t*x[1]);
} else if (bndCondition == 232) {
- return - ((t*t*x[2]) / ((1 + t*t*x[0]) * (1 + t*t*x[1])));
+ return - (t*t*x[2]);
} else {
return 0;
}
@@ -246,14 +246,11 @@ VectorField DExponentialNeumannProblem::Load(double time, const Point &x, const
double DExponentialNeumannProblem::pressure(double t, const Point &x, int bndCondition) const {
if (bndCondition == 230) {
- return - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[0]) /((sin(0.5 * Pi * x[1]) * sin(0.5 * Pi * t) + 1) *
- (sin(0.5 * Pi * x[2]) * sin(0.5 * Pi * t) + 1));
+ return - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[0]) ;
} else if (bndCondition == 231) {
- return - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[1]) /((sin(0.5 * Pi * x[0]) * sin(0.5 * Pi * t) + 1) *
- (sin(0.5 * Pi * x[2]) * sin(0.5 * Pi * t) + 1));
+ return - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[1]);
} else if (bndCondition == 232) {
- return - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[2]) /((sin(0.5 * Pi * x[1]) * sin(0.5 * Pi * t) + 1) *
- (sin(0.5 * Pi * x[0]) * sin(0.5 * Pi * t) + 1));
+ return - sin(0.5 * Pi * t) * sin(0.5 * Pi * x[2]) ;
} else {
return 0;
}
diff --git a/src/elasticity/solvers/IterativePressureSolver.cpp b/src/elasticity/solvers/IterativePressureSolver.cpp
index 0a27fec408ca2b5cef3cb90341a72e2dd0d5e5d5..950b3bd76bf5bf74438a5dec7d94092841149269 100644
--- a/src/elasticity/solvers/IterativePressureSolver.cpp
+++ b/src/elasticity/solvers/IterativePressureSolver.cpp
@@ -23,8 +23,13 @@ bool IterativePressureSolver::Method(IElasticity &assemble, Vector &u) {
Vector uNew(u);
Initialize(assemble, uNew);
+ Vector d(uNew);
+ d = 0;
while (!iteration.IsFinished()) {
+ Vector uOld(uNew);
+ uNew += d;
bool converged = Step(assemble, uNew);
+ d = uNew - uOld;
if(converged) {
assemble.PrintPressureIteration(uNew);
assemble.PlotPressureIteration(uNew, iteration.Step());
diff --git a/test/elasticity/TestPrestress.cpp b/test/elasticity/TestPrestress.cpp
index 3f06cb63eb0c2e94a2aa6562ece59feab65483b7..2f06244d274f21f10c0aca4870ad4b50109d5fe3 100644
--- a/test/elasticity/TestPrestress.cpp
+++ b/test/elasticity/TestPrestress.cpp
@@ -30,6 +30,11 @@ protected:
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;
testConfig["Mechlevel"] = std::to_string(GetParam().level);
@@ -61,7 +66,8 @@ protected:
TEST_P(PrestressTest, TestInitialValue) {
Vector &solution = cmMain->Run();
- mout << solution << endl;
+ 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) {
EXPECT_NEAR(solution(r, i), 0.0, NEWTON_TOLERANCE);
@@ -109,6 +115,6 @@ INSTANTIATE_TEST_SUITE_P(OnEllipsoid, PrestressTest, Values(
*/
int main(int argc, char **argv) {
- MppTest mppTest = MppTestBuilder(argc, argv).WithPPM().WithoutDefaultConfig().WithScreenLogging();
+ MppTest mppTest = MppTestBuilder(argc, argv).WithPPM().WithoutDefaultConfig().WithScreenLogging().WithFileLogging();
return mppTest.RUN_ALL_MPP_TESTS();
}
\ No newline at end of file