From afe904eaede8257e7eaa68b112a72a5e3c069def Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue, 12 Dec 2023 13:50:02 +0100
Subject: [PATCH] [HM/LIE] Changed the output of element averaged tensor and
vector variables
to vectorial output, and removed the element averaged strain components.
In addition, assigned meaningful names to some element variables for
ouput.
---
.../HydroMechanics/HydroMechanicsProcess.cpp | 150 +++---------------
.../HydroMechanicsProcessData.h | 25 +--
...ydroMechanicsLocalAssemblerFracture-impl.h | 33 ++--
.../HydroMechanicsLocalAssemblerMatrix-impl.h | 59 +++----
4 files changed, 64 insertions(+), 203 deletions(-)
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 2df43214ef3..19046161d0e 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 e334d368db2..ec6d974b10e 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 1c4d00c8d9e..04145faea9f 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 aec11079ddb..6371fe136a7 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,
--
GitLab