Merge pull request #1871 from endJunction/HMQuadraticMesh

Fixing HM specific error.

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;
 }
 
 }

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;
 }
-
 }
 

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() -=

ProcessLib/HydroMechanics/Tests.cmake

 # 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
 )

Data @ 106f66cb

-Subproject commit d52965530fe0564b5a3a7ad4ede851bf0a484014
+Subproject commit 106f66cbe93a67ed7139ffa3582f09570c29d4ad