From 515759b342a34aaef21a008753b9a48724e18314 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Thu, 25 Nov 2021 12:04:07 +0100
Subject: [PATCH] [ProcessLib/LIE/HM] Roll back to use dNdx in the global
 coordinate system

---
 ...ydroMechanicsLocalAssemblerFracture-impl.h | 20 +++++++++----------
 .../HydroMechanicsLocalAssemblerFracture.h    |  9 +++------
 .../IntegrationPointDataFracture.h            |  6 +++---
 3 files changed, 16 insertions(+), 19 deletions(-)

diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
index 3074dd4d01b..5e8cd289ac4 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
@@ -196,8 +196,9 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
     Eigen::MatrixXd const global2local_rotation =
         R.template topLeftCorner<ShapeFunctionPressure::DIM, GlobalDim>();
 
-    DimVectorType const gravity_vec =
-        global2local_rotation * _process_data.specific_body_force;
+    GlobalDimVectorType const gravity_vec = global2local_rotation.transpose() *
+                                            global2local_rotation *
+                                            _process_data.specific_body_force;
 
     ParameterLib::SpatialPosition x_position;
     x_position.setElementID(_element.getID());
@@ -260,7 +261,7 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
             ip_data.permeability_state.get(), ip_data.aperture0, b_m);
 
         // derivative of permeability respect to aperture
-        double const local_dk_db =
+        double const dk_db =
             frac_prop.permeability_model->dpermeability_daperture(
                 ip_data.permeability_state.get(), ip_data.aperture0, b_m);
 
@@ -289,18 +290,17 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
         //
         // pressure equation, displacement jump part.
         //
-        auto const grad_head_over_mu =
-            ((dNdx_p * p + rho_fr * gravity_vec) / mu).eval();
+        GlobalDimVectorType const grad_head_over_mu =
+            (dNdx_p * p + rho_fr * gravity_vec) / mu;
         Eigen::Matrix<double, 1, displacement_size> const mT_R_Hg =
             identity2.transpose() * R * H_g;
         // velocity in global coordinates
-        ip_data.darcy_velocity =
-            -global2local_rotation.transpose() * k * grad_head_over_mu;
+        ip_data.darcy_velocity = -k * grad_head_over_mu;
         J_pg.noalias() += N_p.transpose() * S * N_p * p_dot * mT_R_Hg * ip_w;
         J_pg.noalias() +=
             dNdx_p.transpose() * k * grad_head_over_mu * mT_R_Hg * ip_w;
-        J_pg.noalias() += dNdx_p.transpose() * b_m * local_dk_db *
-                          grad_head_over_mu * mT_R_Hg * ip_w;
+        J_pg.noalias() += dNdx_p.transpose() * b_m * dk_db * grad_head_over_mu *
+                          mT_R_Hg * ip_w;
     }
 
     // displacement equation, pressure part
@@ -387,7 +387,7 @@ void HydroMechanicsLocalAssemblerFracture<
     typename HMatricesType::ForceVectorType ele_w =
         HMatricesType::ForceVectorType::Zero(GlobalDim);
     double ele_Fs = -std::numeric_limits<double>::max();
-    GlobalDimVector ele_velocity = GlobalDimVector::Zero();
+    GlobalDimVectorType ele_velocity = GlobalDimVectorType::Zero();
     for (auto const& ip : _ip_data)
     {
         ele_b += ip.aperture;
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture.h
index 84e231dd342..35e60d70bea 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture.h
@@ -80,16 +80,14 @@ private:
 
     // Types for displacement.
     using ShapeMatricesTypeDisplacement =
-        ShapeMatrixPolicyType<ShapeFunctionDisplacement,
-                              ShapeFunctionDisplacement::DIM>;
+        ShapeMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>;
     using HMatricesType =
         HMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>;
     using HMatrixType = typename HMatricesType::HMatrixType;
 
     // Types for pressure.
     using ShapeMatricesTypePressure =
-        ShapeMatrixPolicyType<ShapeFunctionPressure,
-                              ShapeFunctionPressure::DIM>;
+        ShapeMatrixPolicyType<ShapeFunctionPressure, GlobalDim>;
 
     // Types for the integration point data
     using IntegrationPointDataType =
@@ -98,8 +96,7 @@ private:
                                      ShapeMatricesTypePressure,
                                      GlobalDim>;
 
-    using GlobalDimVector = Eigen::Matrix<double, GlobalDim, 1>;
-    using DimVectorType = typename ShapeMatricesTypePressure::DimVectorType;
+    using GlobalDimVectorType = Eigen::Matrix<double, GlobalDim, 1>;
 
     HydroMechanicsProcessData<GlobalDim>& _process_data;
 
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/IntegrationPointDataFracture.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/IntegrationPointDataFracture.h
index bdbb1bdcaef..81fea8ea7a4 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/IntegrationPointDataFracture.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/IntegrationPointDataFracture.h
@@ -29,11 +29,11 @@ struct IntegrationPointDataFracture final
         : fracture_material(fracture_material_),
           material_state_variables(
               fracture_material.createMaterialStateVariables()),
-          darcy_velocity(GlobalDimVector::Zero())
+          darcy_velocity(GlobalDimVectorType::Zero())
     {
     }
 
-    using GlobalDimVector = Eigen::Matrix<double, GlobalDim, 1>;
+    using GlobalDimVectorType = Eigen::Matrix<double, GlobalDim, 1>;
 
     typename HMatricesType::HMatrixType H_u;
     typename HMatricesType::ForceVectorType sigma_eff, sigma_eff_prev;
@@ -58,7 +58,7 @@ struct IntegrationPointDataFracture final
     Eigen::MatrixXd C;
     double integration_weight;
 
-    GlobalDimVector darcy_velocity;
+    GlobalDimVectorType darcy_velocity;
 
     void pushBackState()
     {
-- 
GitLab