diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 2df43214ef36ef09800d67f5fb1ea236cf921b0c..19046161d0e9e64d1a88ff420e4ad1ee1fcd2c14 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -249,89 +249,14 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
"fracture_permeability", 1,
&LocalAssemblerInterface::getIntPtFracturePermeability);
- auto mesh_prop_sigma_xx = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "stress_xx",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_sigma_xx->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_stress_xx = mesh_prop_sigma_xx;
-
- auto mesh_prop_sigma_yy = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "stress_yy",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_sigma_yy->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_stress_yy = mesh_prop_sigma_yy;
-
- auto mesh_prop_sigma_zz = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "stress_zz",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_sigma_zz->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_stress_zz = mesh_prop_sigma_zz;
-
- auto mesh_prop_sigma_xy = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "stress_xy",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_sigma_xy->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_stress_xy = mesh_prop_sigma_xy;
-
- if (GlobalDim == 3)
- {
- auto mesh_prop_sigma_xz = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "stress_xz",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_sigma_xz->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_stress_xz = mesh_prop_sigma_xz;
-
- auto mesh_prop_sigma_yz = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "stress_yz",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_sigma_yz->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_stress_yz = mesh_prop_sigma_yz;
- }
-
- auto mesh_prop_epsilon_xx = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "strain_xx",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_epsilon_xx->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_strain_xx = mesh_prop_epsilon_xx;
-
- auto mesh_prop_epsilon_yy = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "strain_yy",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_epsilon_yy->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_strain_yy = mesh_prop_epsilon_yy;
-
- auto mesh_prop_epsilon_zz = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "strain_zz",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_epsilon_zz->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_strain_zz = mesh_prop_epsilon_zz;
-
- auto mesh_prop_epsilon_xy = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "strain_xy",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_epsilon_xy->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_strain_xy = mesh_prop_epsilon_xy;
-
- if (GlobalDim == 3)
- {
- auto mesh_prop_epsilon_xz = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "strain_xz",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_epsilon_xz->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_strain_xz = mesh_prop_epsilon_xz;
+ _process_data.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "sigma_avg",
+ MeshLib::MeshItemType::Cell,
+ MathLib::KelvinVector::KelvinVectorType<GlobalDim>::RowsAtCompileTime);
- auto mesh_prop_epsilon_yz = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "strain_yz",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_epsilon_yz->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_strain_yz = mesh_prop_epsilon_yz;
- }
-
- auto mesh_prop_velocity = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "element_velocity",
+ _process_data.element_velocities = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "velocity_avg",
MeshLib::MeshItemType::Cell, GlobalDim);
- mesh_prop_velocity->resize(mesh.getNumberOfElements() * GlobalDim);
- _process_data.mesh_prop_velocity = mesh_prop_velocity;
if (!_vec_fracture_elements.empty())
{
@@ -356,21 +281,26 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
(*mesh_prop_levelset)[e->getID()] = levelsets[0];
}
- auto mesh_prop_w_n = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "w_n",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_w_n->resize(mesh.getNumberOfElements());
- auto mesh_prop_w_s = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "w_s",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_w_s->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_w_n = mesh_prop_w_n;
- _process_data.mesh_prop_w_s = mesh_prop_w_s;
+ _process_data.element_local_jumps =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "local_jump_w_avg",
+ MeshLib::MeshItemType::Cell, GlobalDim);
+
+ _process_data.element_fracture_stresses =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "fracture_stress_avg",
+ MeshLib::MeshItemType::Cell, GlobalDim);
+
+ _process_data.element_fracture_velocities =
+ MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh), "fracture_velocity_avg",
+ MeshLib::MeshItemType::Cell, GlobalDim);
auto mesh_prop_b = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "aperture",
+ const_cast<MeshLib::Mesh&>(mesh), "fracture_aperture_avg",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_b->resize(mesh.getNumberOfElements());
+
auto const mesh_prop_matid = materialIDs(mesh);
if (!mesh_prop_matid)
{
@@ -397,27 +327,11 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
_process_data.mesh_prop_b = mesh_prop_b;
auto mesh_prop_k_f = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "k_f",
+ const_cast<MeshLib::Mesh&>(mesh), "fracture_permeability_avg",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_k_f->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_k_f = mesh_prop_k_f;
- auto mesh_prop_fracture_stress_shear =
- MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "f_stress_s",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_fracture_stress_shear->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_fracture_stress_shear =
- mesh_prop_fracture_stress_shear;
-
- auto mesh_prop_fracture_stress_normal =
- MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "f_stress_n",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_fracture_stress_normal->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_fracture_stress_normal =
- mesh_prop_fracture_stress_normal;
-
auto mesh_prop_fracture_shear_failure =
MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "f_shear_failure",
@@ -438,24 +352,6 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
mesh_prop_nodal_b->resize(mesh.getNumberOfNodes());
_process_data.mesh_prop_nodal_b = mesh_prop_nodal_b;
- if (GlobalDim == 3)
- {
- auto mesh_prop_w_s2 = MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "w_s2",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_w_s2->resize(mesh.getNumberOfElements());
- _process_data.mesh_prop_w_s2 = mesh_prop_w_s2;
-
- auto mesh_prop_fracture_stress_shear2 =
- MeshLib::getOrCreateMeshProperty<double>(
- const_cast<MeshLib::Mesh&>(mesh), "f_stress_s2",
- MeshLib::MeshItemType::Cell, 1);
- mesh_prop_fracture_stress_shear2->resize(
- mesh.getNumberOfElements());
- _process_data.mesh_prop_fracture_stress_shear2 =
- mesh_prop_fracture_stress_shear2;
- }
-
auto mesh_prop_nodal_p = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
MeshLib::MeshItemType::Node, 1);
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
index e334d368db2497be3681940c4b36bea9bbef40b6..ec6d974b10e4f6bd9adf1c359cd56328e0d9ceb7 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
@@ -110,27 +110,14 @@ struct HydroMechanicsProcessData
ParameterLib::Parameter<double> const* p0 = nullptr;
// mesh properties for output
- MeshLib::PropertyVector<double>* mesh_prop_stress_xx = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_stress_yy = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_stress_zz = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_stress_xy = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_stress_yz = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_stress_xz = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_strain_xx = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_strain_yy = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_strain_zz = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_strain_xy = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_strain_yz = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_strain_xz = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_velocity = nullptr;
+ MeshLib::PropertyVector<double>* element_stresses = nullptr;
+ MeshLib::PropertyVector<double>* element_velocities = nullptr;
+ MeshLib::PropertyVector<double>* element_local_jumps = nullptr;
+ MeshLib::PropertyVector<double>* element_fracture_stresses = nullptr;
+ MeshLib::PropertyVector<double>* element_fracture_velocities = nullptr;
+
MeshLib::PropertyVector<double>* mesh_prop_b = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_k_f = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_w_n = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_w_s = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_w_s2 = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_shear = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_shear2 = nullptr;
- MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_normal = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_fracture_shear_failure = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_nodal_w = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_nodal_b = nullptr;
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
index 1c4d00c8d9eb3f22e6d062ccaa3b27cf1c7e98e9..04145faea9fa2608de4b21528e66e55c58c2db50 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
@@ -341,7 +341,8 @@ void HydroMechanicsLocalAssemblerFracture<
auto constexpr index_normal = GlobalDim - 1;
ParameterLib::SpatialPosition x_position;
- x_position.setElementID(_element.getID());
+ auto const e_id = _element.getID();
+ x_position.setElementID(e_id);
unsigned const n_integration_points = _ip_data.size();
for (unsigned ip = 0; ip < n_integration_points; ip++)
@@ -391,17 +392,18 @@ void HydroMechanicsLocalAssemblerFracture<
HMatricesType::ForceVectorType::Zero(GlobalDim);
typename HMatricesType::ForceVectorType ele_w =
HMatricesType::ForceVectorType::Zero(GlobalDim);
- double ele_Fs = -std::numeric_limits<double>::max();
GlobalDimVectorType ele_velocity = GlobalDimVectorType::Zero();
+
+ double ele_Fs = -std::numeric_limits<double>::max();
for (auto const& ip : _ip_data)
{
ele_b += ip.aperture;
ele_k += ip.permeability;
ele_w += ip.w;
ele_sigma_eff += ip.sigma_eff;
+ ele_velocity += ip.darcy_velocity;
ele_Fs = std::max(
ele_Fs, ip.material_state_variables->getShearYieldFunctionValue());
- ele_velocity += ip.darcy_velocity;
}
ele_b /= static_cast<double>(n_integration_points);
ele_k /= static_cast<double>(n_integration_points);
@@ -411,20 +413,19 @@ void HydroMechanicsLocalAssemblerFracture<
auto const element_id = _element.getID();
(*_process_data.mesh_prop_b)[element_id] = ele_b;
(*_process_data.mesh_prop_k_f)[element_id] = ele_k;
- (*_process_data.mesh_prop_w_n)[element_id] = ele_w[index_normal];
- (*_process_data.mesh_prop_w_s)[element_id] = ele_w[0];
- (*_process_data.mesh_prop_fracture_stress_normal)[element_id] =
- ele_sigma_eff[index_normal];
- (*_process_data.mesh_prop_fracture_stress_shear)[element_id] =
- ele_sigma_eff[0];
- (*_process_data.mesh_prop_fracture_shear_failure)[element_id] = ele_Fs;
- if (GlobalDim == 3)
- {
- (*_process_data.mesh_prop_w_s2)[element_id] = ele_w[1];
- (*_process_data.mesh_prop_fracture_stress_shear2)[element_id] =
- ele_sigma_eff[1];
- }
+ Eigen::Map<GlobalDimVectorType>(
+ &(*_process_data.element_fracture_stresses)[e_id * GlobalDim]) =
+ ele_sigma_eff;
+
+ Eigen::Map<GlobalDimVectorType>(
+ &(*_process_data.element_fracture_velocities)[e_id * GlobalDim]) =
+ ele_velocity;
+
+ Eigen::Map<GlobalDimVectorType>(
+ &(*_process_data.element_local_jumps)[e_id * GlobalDim]) = ele_w;
+
+ (*_process_data.mesh_prop_fracture_shear_failure)[element_id] = ele_Fs;
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
index aec11079ddbdff9812473d2240ef8d6aa63ea351..6371fe136a7b1a7115354772f269884595084cad 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
@@ -331,7 +331,14 @@ void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
MPL::VariableArray variables;
MPL::VariableArray variables_prev;
ParameterLib::SpatialPosition x_position;
- x_position.setElementID(_element.getID());
+
+ auto const e_id = _element.getID();
+ x_position.setElementID(e_id);
+
+ using KV = MathLib::KelvinVector::KelvinVectorType<GlobalDim>;
+ KV sigma_avg = KV::Zero();
+ GlobalDimVector velocity_avg;
+ velocity_avg.setZero();
unsigned const n_integration_points = _ip_data.size();
for (unsigned ip = 0; ip < n_integration_points; ip++)
@@ -384,6 +391,8 @@ void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
MathLib::KelvinVector::KelvinMatrixType<GlobalDim> C;
std::tie(sigma_eff, state, C) = std::move(*solution);
+ sigma_avg += ip_data.sigma_eff;
+
if (!_process_data.deactivate_matrix_in_flow) // Only for hydraulically
// active matrix
{
@@ -396,51 +405,19 @@ void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
ip_data.darcy_velocity.head(GlobalDim).noalias() =
-k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
+ velocity_avg += ip_data.darcy_velocity.head(GlobalDim);
}
}
- int n = GlobalDim == 2 ? 4 : 6;
- Eigen::VectorXd ele_stress = Eigen::VectorXd::Zero(n);
- Eigen::VectorXd ele_strain = Eigen::VectorXd::Zero(n);
- GlobalDimVector ele_velocity = GlobalDimVector::Zero();
+ sigma_avg /= n_integration_points;
+ velocity_avg /= n_integration_points;
- for (auto const& ip_data : _ip_data)
- {
- ele_stress += ip_data.sigma_eff;
- ele_strain += ip_data.eps;
- ele_velocity += ip_data.darcy_velocity;
- }
+ Eigen::Map<KV>(
+ &(*_process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_avg);
- ele_stress /= static_cast<double>(n_integration_points);
- ele_strain /= static_cast<double>(n_integration_points);
- ele_velocity /= static_cast<double>(n_integration_points);
-
- auto const element_id = _element.getID();
- (*_process_data.mesh_prop_stress_xx)[_element.getID()] = ele_stress[0];
- (*_process_data.mesh_prop_stress_yy)[_element.getID()] = ele_stress[1];
- (*_process_data.mesh_prop_stress_zz)[_element.getID()] = ele_stress[2];
- (*_process_data.mesh_prop_stress_xy)[_element.getID()] = ele_stress[3];
- if (GlobalDim == 3)
- {
- (*_process_data.mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
- (*_process_data.mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
- }
-
- (*_process_data.mesh_prop_strain_xx)[_element.getID()] = ele_strain[0];
- (*_process_data.mesh_prop_strain_yy)[_element.getID()] = ele_strain[1];
- (*_process_data.mesh_prop_strain_zz)[_element.getID()] = ele_strain[2];
- (*_process_data.mesh_prop_strain_xy)[_element.getID()] = ele_strain[3];
- if (GlobalDim == 3)
- {
- (*_process_data.mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
- (*_process_data.mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
- }
-
- for (unsigned i = 0; i < GlobalDim; i++)
- {
- (*_process_data.mesh_prop_velocity)[element_id * GlobalDim + i] =
- ele_velocity[i];
- }
+ Eigen::Map<GlobalDimVector>(
+ &(*_process_data.element_velocities)[e_id * GlobalDim]) = velocity_avg;
NumLib::interpolateToHigherOrderNodes<
ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,