From 3c22b7527c5f5160ed19ef165e43b0d9a20959f2 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Wed, 8 Mar 2017 17:54:55 +0100
Subject: [PATCH] [PL] Extend HM/LIE to 3D.
---
Applications/ApplicationsLib/ProjectData.cpp | 7 +++
.../CreateHydroMechanicsProcess.cpp | 8 ++++
.../HydroMechanicsProcess-fwd.h | 1 +
.../HydroMechanics/HydroMechanicsProcess.cpp | 38 ++++++++++++++++
.../HydroMechanics/HydroMechanicsProcess.h | 4 +-
.../HydroMechanicsProcessData.h | 6 +++
...ydroMechanicsLocalAssemblerFracture-impl.h | 6 ++-
.../HydroMechanicsLocalAssemblerMatrix-impl.h | 43 +++++++++++--------
.../LocalAssembler/LocalDataInitializer.h | 22 +++++-----
ProcessLib/LIE/HydroMechanics/Tests.cmake | 15 +++++++
10 files changed, 117 insertions(+), 33 deletions(-)
diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 4804da923c3..34f2d5a6588 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -371,6 +371,13 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
_process_variables, _parameters, integration_order,
process_config);
break;
+ case 3:
+ process = ProcessLib::LIE::HydroMechanics::
+ createHydroMechanicsProcess<3>(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config);
+ break;
default:
OGS_FATAL(
"HYDRO_MECHANICS_WITH_LIE process does not support "
diff --git a/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
index 618f338480c..76213f3603c 100644
--- a/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -32,6 +32,7 @@ template <unsigned GlobalDim>
class HydroMechanicsProcess;
extern template class HydroMechanicsProcess<2>;
+extern template class HydroMechanicsProcess<3>;
template <unsigned GlobalDim>
std::unique_ptr<Process> createHydroMechanicsProcess(
@@ -312,6 +313,13 @@ template std::unique_ptr<Process> createHydroMechanicsProcess<2>(
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
unsigned const integration_order,
BaseLib::ConfigTree const& config);
+template std::unique_ptr<Process> createHydroMechanicsProcess<3>(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config);
} // namespace HydroMechanics
} // namespace LIE
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess-fwd.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess-fwd.h
index 340890b30af..c82d43caa63 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess-fwd.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess-fwd.h
@@ -12,3 +12,4 @@
#include "HydroMechanicsProcess.h"
extern template class ProcessLib::LIE::HydroMechanics::HydroMechanicsProcess<2>;
+extern template class ProcessLib::LIE::HydroMechanics::HydroMechanicsProcess<3>;
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 814aa30eab8..ee45167a7e3 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -219,12 +219,31 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
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;
+
+ 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_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_epsilon_xx = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh),
"strain_xx", MeshLib::MeshItemType::Cell, 1);
@@ -237,12 +256,30 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
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;
+ 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_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;
+
auto mesh_prop_velocity = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh),
"velocity", MeshLib::MeshItemType::Cell, 3);
@@ -425,6 +462,7 @@ void HydroMechanicsProcess<GlobalDim>::computeSecondaryVariableConcrete(
// template instantiation
// ------------------------------------------------------------------------------------
template class HydroMechanicsProcess<2>;
+template class HydroMechanicsProcess<3>;
} // namespace HydroMechanics
} // namespace LIE
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
index a3ecb7442f5..91070d20779 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.h
@@ -29,9 +29,9 @@ class HydroMechanicsProcess final : public Process
{
using Base = Process;
- static_assert(GlobalDim==2,
+ static_assert(GlobalDim == 2 || GlobalDim == 3,
"Currently LIE::HydroMechanicsProcess "
- "supports only 2D.");
+ "supports only 2D or 3D.");
public:
HydroMechanicsProcess(
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
index 8d66a969365..a9bd8d1176a 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
@@ -128,10 +128,16 @@ struct HydroMechanicsProcessData
// 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>* mesh_prop_b = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_k_f = nullptr;
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
index 7001de80f5e..8297ccddadc 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
@@ -369,8 +369,10 @@ computeSecondaryVariableConcreteWithVector(
double ele_b = 0;
double ele_k = 0;
- Eigen::Vector2d ele_sigma_eff = Eigen::Vector2d::Zero();
- Eigen::Vector2d ele_w = Eigen::Vector2d::Zero();
+ typename HMatricesType::ForceVectorType ele_sigma_eff =
+ HMatricesType::ForceVectorType::Zero(GlobalDim);
+ typename HMatricesType::ForceVectorType ele_w =
+ HMatricesType::ForceVectorType::Zero(GlobalDim);
double ele_Fs = - std::numeric_limits<double>::max();
for (auto const& ip : _ip_data)
{
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
index fe74d623f9d..59af4ef4b65 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
@@ -367,22 +367,15 @@ computeSecondaryVariableConcreteWithBlockVectors(
q.noalias() = - k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
}
- Eigen::Vector3d ele_stress = Eigen::Vector3d::Zero();
- Eigen::Vector3d ele_strain = Eigen::Vector3d::Zero();
+ int n = GlobalDim == 2 ? 4 : 6;
+ Eigen::VectorXd ele_stress = Eigen::VectorXd::Zero(n);
+ Eigen::VectorXd ele_strain = Eigen::VectorXd::Zero(n);
Eigen::Vector3d ele_velocity = Eigen::Vector3d::Zero();
for (auto const& ip_data : _ip_data)
{
- ele_stress[0] += ip_data.sigma_eff[0];
- ele_stress[1] += ip_data.sigma_eff[1];
- // sigma_xy
- ele_stress[2] += ip_data.sigma_eff[3]; // sigma_eff[2] stores sigma_z
-
- ele_strain[0] += ip_data.eps[0];
- ele_strain[1] += ip_data.eps[1];
- // strain_xy
- ele_strain[2] += ip_data.eps[3];
-
+ ele_stress += ip_data.sigma_eff;
+ ele_strain += ip_data.eps;
ele_velocity += ip_data.darcy_velocity;
}
@@ -391,12 +384,26 @@ computeSecondaryVariableConcreteWithBlockVectors(
ele_velocity /= static_cast<double>(n_integration_points);
auto const element_id = _element.getID();
- (*_process_data.mesh_prop_stress_xx)[element_id] = ele_stress[0];
- (*_process_data.mesh_prop_stress_yy)[element_id] = ele_stress[1];
- (*_process_data.mesh_prop_stress_xy)[element_id] = ele_stress[2];
- (*_process_data.mesh_prop_strain_xx)[element_id] = ele_strain[0];
- (*_process_data.mesh_prop_strain_yy)[element_id] = ele_strain[1];
- (*_process_data.mesh_prop_strain_xy)[element_id] = ele_strain[2];
+ (*_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<3; i++)
(*_process_data.mesh_prop_velocity)[element_id*3 + i] = ele_velocity[i];
}
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/LocalDataInitializer.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/LocalDataInitializer.h
index 3827ae943ed..00f54e9f376 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/LocalDataInitializer.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/LocalDataInitializer.h
@@ -71,7 +71,6 @@ static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_CUBOID) | \
(ENABLED_ELEMENT_TYPE_PYRAMID) | (ENABLED_ELEMENT_TYPE_PRISM))
-
// Include only what is needed (Well, the conditions are not sharp).
#if OGS_ENABLED_ELEMENTS != 0
#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
@@ -115,7 +114,6 @@ namespace LIE
{
namespace HydroMechanics
{
-
/// The LocalDataInitializer is a functor creating a local assembler data with
/// corresponding to the mesh element type shape functions and calling
/// initialization of the new local assembler data.
@@ -125,12 +123,14 @@ namespace HydroMechanics
/// \attention This is modified version of the ProcessLib::LocalDataInitializer
/// class which does not include line elements, allows only shapefunction of
/// order 2, and provides additional template argument GlobalDim.
-template <
- typename LocalAssemblerInterface,
- template <typename, typename, typename, unsigned> class LocalAssemblerDataMatrix,
- template <typename, typename, typename, unsigned> class LocalAssemblerDataMatrixNearFracture,
- template <typename, typename, typename, unsigned> class LocalAssemblerDataFracture,
- unsigned GlobalDim, typename... ConstructorArgs>
+template <typename LocalAssemblerInterface,
+ template <typename, typename, typename, unsigned>
+ class LocalAssemblerDataMatrix,
+ template <typename, typename, typename, unsigned>
+ class LocalAssemblerDataMatrixNearFracture,
+ template <typename, typename, typename, unsigned>
+ class LocalAssemblerDataFracture,
+ unsigned GlobalDim, typename... ConstructorArgs>
class LocalDataInitializer final
{
public:
@@ -145,7 +145,7 @@ public:
"The given shape function order %d is not supported.\nOnly "
"shape functions of order 2 are supported.",
shapefunction_order);
- // /// Quads and Hexahedra ///////////////////////////////////
+// /// Quads and Hexahedra ///////////////////////////////////
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
@@ -155,7 +155,7 @@ public:
makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
#endif
- // /// Simplices ////////////////////////////////////////////////
+// /// Simplices ////////////////////////////////////////////////
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
@@ -163,7 +163,7 @@ public:
makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
#endif
- // /// Lines ///////////////////////////////////
+// /// Lines ///////////////////////////////////
#if OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Line3))] =
makeLocalAssemblerBuilder<NumLib::ShapeLine3>();
diff --git a/ProcessLib/LIE/HydroMechanics/Tests.cmake b/ProcessLib/LIE/HydroMechanics/Tests.cmake
index 7338e8f0eb7..10be51cf617 100644
--- a/ProcessLib/LIE/HydroMechanics/Tests.cmake
+++ b/ProcessLib/LIE/HydroMechanics/Tests.cmake
@@ -14,6 +14,21 @@ AddTest(
expected_single_fracture_pcs_0_ts_10_t_100.000000.vtu single_fracture_pcs_0_ts_10_t_100.000000.vtu displacement_jump1 displacement_jump1
)
+AddTest(
+ NAME LIE_HM_single_fracture_3D
+ PATH LIE/HydroMechanics
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS single_fracture_3D.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+ ABSTOL 1e-12 RELTOL 1e-12
+ DIFF_DATA
+ expected_single_fracture_3D_pcs_0_ts_10_t_100.000000.vtu single_fracture_3D_pcs_0_ts_10_t_100.000000.vtu pressure pressure
+ expected_single_fracture_3D_pcs_0_ts_10_t_100.000000.vtu single_fracture_3D_pcs_0_ts_10_t_100.000000.vtu displacement displacement
+ expected_single_fracture_3D_pcs_0_ts_10_t_100.000000.vtu single_fracture_3D_pcs_0_ts_10_t_100.000000.vtu displacement_jump1 displacement_jump1
+)
+
AddTest(
NAME LIE_HM_TaskB
PATH LIE/HydroMechanics
--
GitLab