Skip to content
Snippets Groups Projects
Commit 515759b3 authored by wenqing's avatar wenqing Committed by Dmitri Naumov
Browse files

[ProcessLib/LIE/HM] Roll back to use dNdx in the global coordinate system

parent 6a5a5de9
No related branches found
No related tags found
No related merge requests found
...@@ -196,8 +196,9 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement, ...@@ -196,8 +196,9 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
Eigen::MatrixXd const global2local_rotation = Eigen::MatrixXd const global2local_rotation =
R.template topLeftCorner<ShapeFunctionPressure::DIM, GlobalDim>(); R.template topLeftCorner<ShapeFunctionPressure::DIM, GlobalDim>();
DimVectorType const gravity_vec = GlobalDimVectorType const gravity_vec = global2local_rotation.transpose() *
global2local_rotation * _process_data.specific_body_force; global2local_rotation *
_process_data.specific_body_force;
ParameterLib::SpatialPosition x_position; ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID()); x_position.setElementID(_element.getID());
...@@ -260,7 +261,7 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement, ...@@ -260,7 +261,7 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
ip_data.permeability_state.get(), ip_data.aperture0, b_m); ip_data.permeability_state.get(), ip_data.aperture0, b_m);
// derivative of permeability respect to aperture // derivative of permeability respect to aperture
double const local_dk_db = double const dk_db =
frac_prop.permeability_model->dpermeability_daperture( frac_prop.permeability_model->dpermeability_daperture(
ip_data.permeability_state.get(), ip_data.aperture0, b_m); ip_data.permeability_state.get(), ip_data.aperture0, b_m);
...@@ -289,18 +290,17 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement, ...@@ -289,18 +290,17 @@ void HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
// //
// pressure equation, displacement jump part. // pressure equation, displacement jump part.
// //
auto const grad_head_over_mu = GlobalDimVectorType const grad_head_over_mu =
((dNdx_p * p + rho_fr * gravity_vec) / mu).eval(); (dNdx_p * p + rho_fr * gravity_vec) / mu;
Eigen::Matrix<double, 1, displacement_size> const mT_R_Hg = Eigen::Matrix<double, 1, displacement_size> const mT_R_Hg =
identity2.transpose() * R * H_g; identity2.transpose() * R * H_g;
// velocity in global coordinates // velocity in global coordinates
ip_data.darcy_velocity = ip_data.darcy_velocity = -k * grad_head_over_mu;
-global2local_rotation.transpose() * k * grad_head_over_mu;
J_pg.noalias() += N_p.transpose() * S * N_p * p_dot * mT_R_Hg * ip_w; J_pg.noalias() += N_p.transpose() * S * N_p * p_dot * mT_R_Hg * ip_w;
J_pg.noalias() += J_pg.noalias() +=
dNdx_p.transpose() * k * grad_head_over_mu * mT_R_Hg * ip_w; dNdx_p.transpose() * k * grad_head_over_mu * mT_R_Hg * ip_w;
J_pg.noalias() += dNdx_p.transpose() * b_m * local_dk_db * J_pg.noalias() += dNdx_p.transpose() * b_m * dk_db * grad_head_over_mu *
grad_head_over_mu * mT_R_Hg * ip_w; mT_R_Hg * ip_w;
} }
// displacement equation, pressure part // displacement equation, pressure part
...@@ -387,7 +387,7 @@ void HydroMechanicsLocalAssemblerFracture< ...@@ -387,7 +387,7 @@ void HydroMechanicsLocalAssemblerFracture<
typename HMatricesType::ForceVectorType ele_w = typename HMatricesType::ForceVectorType ele_w =
HMatricesType::ForceVectorType::Zero(GlobalDim); HMatricesType::ForceVectorType::Zero(GlobalDim);
double ele_Fs = -std::numeric_limits<double>::max(); double ele_Fs = -std::numeric_limits<double>::max();
GlobalDimVector ele_velocity = GlobalDimVector::Zero(); GlobalDimVectorType ele_velocity = GlobalDimVectorType::Zero();
for (auto const& ip : _ip_data) for (auto const& ip : _ip_data)
{ {
ele_b += ip.aperture; ele_b += ip.aperture;
......
...@@ -80,16 +80,14 @@ private: ...@@ -80,16 +80,14 @@ private:
// Types for displacement. // Types for displacement.
using ShapeMatricesTypeDisplacement = using ShapeMatricesTypeDisplacement =
ShapeMatrixPolicyType<ShapeFunctionDisplacement, ShapeMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>;
ShapeFunctionDisplacement::DIM>;
using HMatricesType = using HMatricesType =
HMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>; HMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>;
using HMatrixType = typename HMatricesType::HMatrixType; using HMatrixType = typename HMatricesType::HMatrixType;
// Types for pressure. // Types for pressure.
using ShapeMatricesTypePressure = using ShapeMatricesTypePressure =
ShapeMatrixPolicyType<ShapeFunctionPressure, ShapeMatrixPolicyType<ShapeFunctionPressure, GlobalDim>;
ShapeFunctionPressure::DIM>;
// Types for the integration point data // Types for the integration point data
using IntegrationPointDataType = using IntegrationPointDataType =
...@@ -98,8 +96,7 @@ private: ...@@ -98,8 +96,7 @@ private:
ShapeMatricesTypePressure, ShapeMatricesTypePressure,
GlobalDim>; GlobalDim>;
using GlobalDimVector = Eigen::Matrix<double, GlobalDim, 1>; using GlobalDimVectorType = Eigen::Matrix<double, GlobalDim, 1>;
using DimVectorType = typename ShapeMatricesTypePressure::DimVectorType;
HydroMechanicsProcessData<GlobalDim>& _process_data; HydroMechanicsProcessData<GlobalDim>& _process_data;
......
...@@ -29,11 +29,11 @@ struct IntegrationPointDataFracture final ...@@ -29,11 +29,11 @@ struct IntegrationPointDataFracture final
: fracture_material(fracture_material_), : fracture_material(fracture_material_),
material_state_variables( material_state_variables(
fracture_material.createMaterialStateVariables()), 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::HMatrixType H_u;
typename HMatricesType::ForceVectorType sigma_eff, sigma_eff_prev; typename HMatricesType::ForceVectorType sigma_eff, sigma_eff_prev;
...@@ -58,7 +58,7 @@ struct IntegrationPointDataFracture final ...@@ -58,7 +58,7 @@ struct IntegrationPointDataFracture final
Eigen::MatrixXd C; Eigen::MatrixXd C;
double integration_weight; double integration_weight;
GlobalDimVector darcy_velocity; GlobalDimVectorType darcy_velocity;
void pushBackState() void pushBackState()
{ {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment