diff --git a/NumLib/Fem/ShapeFunction/ShapeQuad8-impl.h b/NumLib/Fem/ShapeFunction/ShapeQuad8-impl.h
index 2ddb39bfa9a4a7b879f80b4161a22e71bba67661..dbbb49747764fe8ef3940083275877beb9f46cfd 100644
--- a/NumLib/Fem/ShapeFunction/ShapeQuad8-impl.h
+++ b/NumLib/Fem/ShapeFunction/ShapeQuad8-impl.h
@@ -13,15 +13,15 @@ namespace NumLib
template <class T_X, class T_N>
void ShapeQuad8::computeShapeFunction(const T_X &r, T_N &N)
{
- N[0] = -0.25 * (1.0 - r[0]) * (1.0 - r[1]) * (( 1.0 + r[0] + r[1]));
- N[1] = 0.25 * (1.0 + r[0]) * (1.0 - r[1]) * ((-1.0 + r[0] - r[1]));
- N[2] = 0.25 * (1.0 + r[0]) * (1.0 + r[1]) * ((-1.0 + r[0] + r[1]));
- N[3] = -0.25 * (1.0 - r[0]) * (1.0 + r[1]) * (( 1.0 + r[0] - r[1]));
+ N[0] = 0.25 * (1.0 + r[0]) * (1.0 + r[1]) * (-1.0 + r[0] + r[1]);
+ N[1] = -0.25 * (1.0 - r[0]) * (1.0 + r[1]) * (1.0 + r[0] - r[1]);
+ N[2] = -0.25 * (1.0 - r[0]) * (1.0 - r[1]) * (1.0 + r[0] + r[1]);
+ N[3] = 0.25 * (1.0 + r[0]) * (1.0 - r[1]) * (-1.0 + r[0] - r[1]);
//
- N[4] = 0.5 * (1.0 - r[0] * r[0]) * (1.0 - r[1]);
- N[5] = 0.5 * (1.0 - r[1] * r[1]) * (1.0 + r[0]);
- N[6] = 0.5 * (1.0 - r[0] * r[0]) * (1.0 + r[1]);
- N[7] = 0.5 * (1.0 - r[1] * r[1]) * (1.0 - r[0]);
+ N[4] = 0.5 * (1.0 - r[0] * r[0]) * (1.0 + r[1]);
+ N[5] = 0.5 * (1.0 - r[1] * r[1]) * (1.0 - r[0]);
+ N[6] = 0.5 * (1.0 - r[0] * r[0]) * (1.0 - r[1]);
+ N[7] = 0.5 * (1.0 - r[1] * r[1]) * (1.0 + r[0]);
}
template <class T_X, class T_N>
@@ -30,25 +30,27 @@ void ShapeQuad8::computeGradShapeFunction(const T_X &rs, T_N &dNdr)
const double r = rs[0];
const double s = rs[1];
- //dN/dr
- dNdr[0] = (1 - s) * (2 * r + s) * 0.25;
- dNdr[1] = (1 - s) * (2 * r - s) * 0.25;
- dNdr[2] = (1 + s) * (2 * r + s) * 0.25;
- dNdr[3] = (1 + s) * (2 * r - s) * 0.25;
- dNdr[4] = -r * (1 - s);
- dNdr[5] = (1 - s * s) * 0.5;
- dNdr[6] = -r * (1 + s);
- dNdr[7] = -dNdr[5];
-
- //dN/ds
- dNdr[8] = (1 - r) * (r + 2 * s) * 0.25;
- dNdr[9] = -(1 + r) * (r - 2 * s) * 0.25;
- dNdr[10] = (1 + r) * (r + 2 * s) * 0.25;
- dNdr[11] = -(1 - r) * (r - 2 * s) * 0.25;
- dNdr[12] = -(1 - r * r) * 0.5;
- dNdr[13] = -(1 + r) * s;
- dNdr[14] = -dNdr[12];
- dNdr[15] = -(1 - r) * s;
+ // dN/dr
+ dNdr[0] = (1 + s) * (2 * r + s) * 0.25;
+ dNdr[1] = (1 + s) * (2 * r - s) * 0.25;
+ dNdr[2] = (1 - s) * (2 * r + s) * 0.25;
+ dNdr[3] = (1 - s) * (2 * r - s) * 0.25;
+
+ dNdr[4] = -r * (1 + s);
+ dNdr[5] = -(1 - s * s) * 0.5;
+ dNdr[6] = -r * (1 - s);
+ dNdr[7] = (1 - s * s) * 0.5;
+
+ // dN/ds
+ dNdr[8] = (1 + r) * (r + 2 * s) * 0.25;
+ dNdr[9] = -(1 - r) * (r - 2 * s) * 0.25;
+ dNdr[10] = (1 - r) * (r + 2 * s) * 0.25;
+ dNdr[11] = -(1 + r) * (r - 2 * s) * 0.25;
+
+ dNdr[12] = (1 - r * r) * 0.5;
+ dNdr[13] = -(1 - r) * s;
+ dNdr[14] = -(1 - r * r) * 0.5;
+ dNdr[15] = -(1 + r) * s;
}
}
diff --git a/NumLib/Fem/ShapeFunction/ShapeQuad9-impl.h b/NumLib/Fem/ShapeFunction/ShapeQuad9-impl.h
index d6509dcba0e218102065ccdf33f07801db50906e..afcaed7a261729b4d490a5a1fb92bc4977e8f7e7 100644
--- a/NumLib/Fem/ShapeFunction/ShapeQuad9-impl.h
+++ b/NumLib/Fem/ShapeFunction/ShapeQuad9-impl.h
@@ -9,44 +9,42 @@
namespace NumLib
{
-
template <class T_X, class T_N>
-void ShapeQuad9::computeShapeFunction(const T_X &r, T_N &N)
+void ShapeQuad9::computeShapeFunction(const T_X& r, T_N& N)
{
- N[8] = (1.0 - r[0] * r[0]) * ( 1.0 - r[1] * r[1]);
- N[7] = 0.5 * (1.0 - r[1] * r[1]) * (1.0 + r[0]) - 0.5 * N[8];
- N[6] = 0.5 * (1.0 - r[0] * r[0]) * (1.0 - r[1]) - 0.5 * N[8];
- N[5] = 0.5 * (1.0 - r[1] * r[1]) * (1.0 - r[0]) - 0.5 * N[8];
- N[4] = 0.5 * (1.0 - r[0] * r[0]) * (1.0 + r[1]) - 0.5 * N[8];
- N[3] = 0.25 * (1.0 + r[0]) * (1.0 - r[1]) - 0.5 * N[6] - 0.5 * N[7] - 0.25 * N[8];
- N[2] = 0.25 * (1.0 - r[0]) * (1.0 - r[1]) - 0.5 * N[5] - 0.5 * N[6] - 0.25 * N[8];
- N[1] = 0.25 * (1.0 - r[0]) * (1.0 + r[1]) - 0.5 * N[4] - 0.5 * N[5] - 0.25 * N[8];
- N[0] = 0.25 * (1.0 + r[0]) * (1.0 + r[1]) - 0.5 * N[4] - 0.5 * N[7] - 0.25 * N[8];
+ N[0] = r[0] * (r[0] + 1) * r[1] * (r[1] + 1) / 4;
+ N[1] = r[0] * (r[0] - 1) * r[1] * (r[1] + 1) / 4;
+ N[2] = r[0] * (r[0] - 1) * r[1] * (r[1] - 1) / 4;
+ N[3] = r[0] * (r[0] + 1) * r[1] * (r[1] - 1) / 4;
+ N[4] = r[1] * (r[1] + 1) * (1 - r[0] * r[0]) / 2;
+ N[5] = r[0] * (r[0] - 1) * (1 - r[1] * r[1]) / 2;
+ N[6] = r[1] * (r[1] - 1) * (1 - r[0] * r[0]) / 2;
+ N[7] = r[0] * (r[0] + 1) * (1 - r[1] * r[1]) / 2;
+ N[8] = (1 - r[0] * r[0]) * (1 - r[1] * r[1]);
}
template <class T_X, class T_N>
-void ShapeQuad9::computeGradShapeFunction(const T_X &r, T_N &dNdr)
+void ShapeQuad9::computeGradShapeFunction(const T_X& r, T_N& dNdr)
{
- dNdr[8] = -2.0 * r[0] * (1.0 - r[1] * r[1]);
- dNdr[7] = +0.5 * (1.0 - r[1] * r[1]) - 0.5 * dNdr[8];
- dNdr[6] = -1.0 * r[0] * (1.0 - r[1]) - 0.5 * dNdr[8];
- dNdr[5] = -dNdr[7];
- dNdr[4] = -1.0 * r[0] * (1.0 + r[1]) - 0.5 * dNdr[8];
- dNdr[3] = +0.25 * (1 - r[1]) - 0.5 * dNdr[6] - 0.5 * dNdr[7] - 0.25 * dNdr[8];
- dNdr[2] = -0.25 * (1 - r[1]) - 0.5 * dNdr[5] - 0.5 * dNdr[6] - 0.25 * dNdr[8];
- dNdr[1] = -0.25 * (1 + r[1]) - 0.5 * dNdr[4] - 0.5 * dNdr[5] - 0.25 * dNdr[8];
- dNdr[0] = +0.25 * (1 + r[1]) - 0.5 * dNdr[4] - 0.5 * dNdr[7] - 0.25 * dNdr[8];
+ dNdr[0] = (r[0] + 0.5) * r[1] * (r[1] + 1) / 2;
+ dNdr[1] = (r[0] - 0.5) * r[1] * (r[1] + 1) / 2;
+ dNdr[2] = (r[0] - 0.5) * r[1] * (r[1] - 1) / 2;
+ dNdr[3] = (r[0] + 0.5) * r[1] * (r[1] - 1) / 2;
+ dNdr[4] = -r[0] * r[1] * (1 + r[1]);
+ dNdr[5] = (1 - r[1] * r[1]) * (r[0] - 0.5);
+ dNdr[6] = r[0] * r[1] * (1 - r[1]);
+ dNdr[7] = (1 - r[1] * r[1]) * (r[0] + 0.5);
+ dNdr[8] = 2 * r[0] * (r[1] * r[1] - 1);
- dNdr[17] = -2.0 * r[1] * (1.0 - r[0] * r[0]);
- dNdr[16] = -1.0 * r[1] * (1.0 + r[0]) - 0.5 * dNdr[17];
- dNdr[15] = -0.5 * (1.0 - r[0] * r[0]) - 0.5 * dNdr[17];
- dNdr[14] = -1.0 * r[1] * (1.0 - r[0]) - 0.5 * dNdr[17];
- dNdr[13] = -dNdr[15];
- dNdr[12] = -0.25 * (1 + r[0]) - 0.5 * dNdr[15] - 0.5 * dNdr[16] - 0.25 * dNdr[17];
- dNdr[11] = -0.25 * (1 - r[0]) - 0.5 * dNdr[14] - 0.5 * dNdr[15] - 0.25 * dNdr[17];
- dNdr[10] = +0.25 * (1 - r[0]) - 0.5 * dNdr[13] - 0.5 * dNdr[14] - 0.25 * dNdr[17];
- dNdr[9] = +0.25 * (1 + r[0]) - 0.5 * dNdr[13] - 0.5 * dNdr[16] - 0.25 * dNdr[17];
+ dNdr[10] = (r[1] + 0.5) * r[0] * (r[0] - 1) / 2;
+ dNdr[11] = (r[1] - 0.5) * r[0] * (r[0] - 1) / 2;
+ dNdr[12] = (r[1] - 0.5) * r[0] * (r[0] + 1) / 2;
+ dNdr[13] = (1 - r[0] * r[0]) * (r[1] + 0.5);
+ dNdr[14] = r[0] * r[1] * (1 - r[0]);
+ dNdr[15] = (1 - r[0] * r[0]) * (r[1] - 0.5);
+ dNdr[16] = -r[0] * r[1] * (1 + r[0]);
+ dNdr[17] = 2 * r[1] * (r[0] * r[0] - 1);
+ dNdr[9] = (r[1] + 0.5) * r[0] * (r[0] + 1) / 2;
}
-
}
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index c08978b8627c7b30845b20da3a11dd312b1f1cce..fb84183887b128d3dab10a69c2024a378ceb293f 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -379,19 +379,19 @@ public:
local_Jac
.template block<displacement_size, pressure_size>(
displacement_index, pressure_index)
- .noalias() -= Kup;
+ .noalias() = -Kup;
// pressure equation, pressure part.
local_Jac
.template block<pressure_size, pressure_size>(pressure_index,
pressure_index)
- .noalias() += laplace_p + storage_p / dt;
+ .noalias() = laplace_p + storage_p / dt;
// pressure equation, displacement part.
local_Jac
.template block<pressure_size, displacement_size>(
pressure_index, displacement_index)
- .noalias() += Kup.transpose() / dt;
+ .noalias() = Kup.transpose() / dt;
// pressure equation
local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
diff --git a/ProcessLib/HydroMechanics/Tests.cmake b/ProcessLib/HydroMechanics/Tests.cmake
index 7a1fa83ba73b5dbe21d90cac1258cb4177c106c7..2af2a31bed5ef6343fdda684acf708f49617faae 100644
--- a/ProcessLib/HydroMechanics/Tests.cmake
+++ b/ProcessLib/HydroMechanics/Tests.cmake
@@ -1,7 +1,7 @@
# HydroMechanics; Small deformations, linear poroelastic (HML)
AddTest(
- NAME HydroMechanics_HML_square_1e2_confined_compression
- PATH HydroMechanics/Linear
+ NAME HydroMechanics_HML_square_1e2_quad8_confined_compression
+ PATH HydroMechanics/Linear/Confined_Compression
EXECUTABLE ogs
EXECUTABLE_ARGS square_1e2.prj
WRAPPER time
@@ -26,6 +26,46 @@ AddTest(
expected_square_1e2_pcs_0_ts_120_t_1000.000000.vtu square_1e2_pcs_0_ts_120_t_1000.000000.vtu velocity_y velocity_y
expected_square_1e2_pcs_0_ts_420_t_4000.000000.vtu square_1e2_pcs_0_ts_420_t_4000.000000.vtu velocity_y velocity_y
)
+AddTest(
+ NAME HydroMechanics_HML_square_1e2_quad9_confined_compression
+ PATH HydroMechanics/Linear/Confined_Compression
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS square_1e2_quad9.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ ABSTOL 1e-15 RELTOL 1e-15
+ DIFF_DATA
+ expected_square_1e2_quad9_pcs_0_ts_1_t_5.000000.vtu square_1e2_quad9_pcs_0_ts_1_t_5.000000.vtu displacement displacement
+ expected_square_1e2_quad9_pcs_0_ts_1_t_5.000000.vtu square_1e2_quad9_pcs_0_ts_1_t_5.000000.vtu pressure pressure
+)
+AddTest(
+ NAME HydroMechanics_HML_square_1e2_tri6_confined_compression
+ PATH HydroMechanics/Linear/Confined_Compression
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS square_1e2_tri.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ ABSTOL 1e-15 RELTOL 1e-15
+ DIFF_DATA
+ expected_square_1e2_tri_pcs_0_ts_1_t_5.000000.vtu square_1e2_tri_pcs_0_ts_1_t_5.000000.vtu displacement displacement
+ expected_square_1e2_tri_pcs_0_ts_1_t_5.000000.vtu square_1e2_tri_pcs_0_ts_1_t_5.000000.vtu pressure pressure
+)
+AddTest(
+ NAME HydroMechanics_HML_cube_1e3_hex20_confined_compression
+ PATH HydroMechanics/Linear/Confined_Compression
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_1e3.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ ABSTOL 1e-15 RELTOL 1e-15
+ DIFF_DATA
+ expected_cube_1e3_pcs_0_ts_1_t_5.000000.vtu cube_1e3_pcs_0_ts_1_t_5.000000.vtu displacement displacement
+ expected_cube_1e3_pcs_0_ts_1_t_5.000000.vtu cube_1e3_pcs_0_ts_1_t_5.000000.vtu pressure pressure
+)
+
# HydroMechanics; Small deformation, linear poroelastic (unconfined compression early) The drainage process is ongoing and the displacement behaviour is related to water pressure and solid properties.
AddTest(
NAME HydroMechanics_HML_square_1e2_unconfined_compression_early
@@ -35,9 +75,10 @@ AddTest(
WRAPPER time
TESTER vtkdiff
REQUIREMENTS NOT OGS_USE_MPI
- ABSTOL 1e-3 RELTOL 1e-3
+ ABSTOL 1e-11 RELTOL 1e-16
DIFF_DATA
- UnconfinedCompressionAnalytical_1s.vtu square_1e2_UC_early_pcs_0_ts_10_t_1.000000.vtu displacement_ana displacement
+ expected_square_1e2_UC_early_pcs_0_ts_10_t_1.000000.vtu square_1e2_UC_early_pcs_0_ts_10_t_1.000000.vtu displacement displacement
+ expected_square_1e2_UC_early_pcs_0_ts_10_t_1.000000.vtu square_1e2_UC_early_pcs_0_ts_10_t_1.000000.vtu pressure pressure
)
# HydroMechanics; Small deformation, linear poroelastic (unconfined compression late) the drainage process is finished and the displacement of the porous media is only a result of solid properties.
AddTest(
@@ -48,7 +89,8 @@ AddTest(
WRAPPER time
TESTER vtkdiff
REQUIREMENTS NOT OGS_USE_MPI
- ABSTOL 1e-3 RELTOL 1e-3
+ ABSTOL 5e-14 RELTOL 1e-16
DIFF_DATA
- UnconfinedCompressionAnalytical_1000s.vtu square_1e2_UC_late_pcs_0_ts_10_t_1000.000000.vtu displacement_ana displacement
+ expected_square_1e2_UC_late_pcs_0_ts_10_t_1000.000000.vtu square_1e2_UC_late_pcs_0_ts_10_t_1000.000000.vtu displacement displacement
+ expected_square_1e2_UC_late_pcs_0_ts_10_t_1000.000000.vtu square_1e2_UC_late_pcs_0_ts_10_t_1000.000000.vtu pressure pressure
)
diff --git a/Tests/Data b/Tests/Data
index d52965530fe0564b5a3a7ad4ede851bf0a484014..106f66cbe93a67ed7139ffa3582f09570c29d4ad 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit d52965530fe0564b5a3a7ad4ede851bf0a484014
+Subproject commit 106f66cbe93a67ed7139ffa3582f09570c29d4ad