Merge branch 'TH_fracture' into 'master'

[HT] Use lower dimensional fracture elements in HT process

See merge request ogs/ogs!4613

Documentation/ProjectFile/prj/processes/process/HT/aperture_size/i_aperture_size.md

@@ -0,0 +1 @@
+\copydoc ProcessLib::HT::HTProcessData::aperture_size

Documentation/ProjectFile/prj/processes/process/HT/aperture_size/t_parameter.md

@@ -0,0 +1,2 @@
+It stores the aperture size for homogeneous aperture, or the aperture size
+ of each fracture element for heterogeneous aperture.

ProcessLib/HT/CreateHTProcess.cpp

@@ -16,6 +16,8 @@
 #include "MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h"
 #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MeshLib/IO/readMeshFromFile.h"
+#include "MeshLib/Utils/GetElementRotationMatrices.h"
+#include "MeshLib/Utils/GetSpaceDimension.h"
 #include "NumLib/NumericalStability/CreateNumericalStabilization.h"
 #include "ParameterLib/ConstantParameter.h"
 #include "ParameterLib/Utils.h"
@@ -118,23 +120,25 @@ std::unique_ptr<Process> createHTProcess(
  }
 
  /// \section parametersht Process Parameters
- // Specific body force parameter.
- Eigen::VectorXd specific_body_force;
  std::vector<double> const b =
  //! \ogs_file_param{prj__processes__process__HT__specific_body_force}
  config.getConfigParameter<std::vector<double>>("specific_body_force");
  assert(!b.empty() && b.size() < 4);
- if (b.size() < mesh.getDimension())
+ int const mesh_space_dimension =
+ MeshLib::getSpaceDimension(mesh.getNodes());
+ if (static_cast<int>(b.size()) != mesh_space_dimension)
  {
  OGS_FATAL(
  "specific body force (gravity vector) has {:d} components, mesh "
  "dimension is {:d}",
- b.size(), mesh.getDimension());
+ b.size(), mesh_space_dimension);
  }
+
+ // Specific body force parameter.
+ Eigen::VectorXd specific_body_force(b.size());
  bool const has_gravity = MathLib::toVector(b).norm() > 0;
  if (has_gravity)
  {
- specific_body_force.resize(b.size());
  std::copy_n(b.data(), b.size(), specific_body_force.data());
  }
 
@@ -181,12 +185,40 @@ std::unique_ptr<Process> createHTProcess(
  DBUG("Media properties verified.");
 
  auto stabilizer = NumLib::createNumericalStabilization(mesh, config);
- HTProcessData process_data{
- std::move(media_map), has_fluid_thermal_expansion,
- *solid_thermal_expansion, *biot_constant,
- specific_body_force, has_gravity,
- heat_transport_process_id, hydraulic_process_id,
- std::move(stabilizer)};
+
+ auto const* aperture_size_parameter = &ParameterLib::findParameter<double>(
+ ProcessLib::Process::constant_one_parameter_name, parameters, 1);
+ auto const aperture_config =
+ //! \ogs_file_param{prj__processes__process__HT__aperture_size}
+ config.getConfigSubtreeOptional("aperture_size");
+ if (aperture_config)
+ {
+ aperture_size_parameter = &ParameterLib::findParameter<double>(
+ //! \ogs_file_param_special{prj__processes__process__HT__aperture_size__parameter}
+ *aperture_config, "parameter", parameters, 1);
+ }
+
+ auto const rotation_matrices = MeshLib::getElementRotationMatrices(
+ mesh_space_dimension, mesh.getDimension(), mesh.getElements());
+ std::vector<Eigen::VectorXd> projected_specific_body_force_vectors;
+ projected_specific_body_force_vectors.reserve(rotation_matrices.size());
+
+ std::transform(rotation_matrices.begin(), rotation_matrices.end(),
+ std::back_inserter(projected_specific_body_force_vectors),
+ [&specific_body_force](const auto& R)
+ { return R * R.transpose() * specific_body_force; });
+
+ HTProcessData process_data{std::move(media_map),
+ has_fluid_thermal_expansion,
+ *solid_thermal_expansion,
+ *biot_constant,
+ has_gravity,
+ heat_transport_process_id,
+ hydraulic_process_id,
+ std::move(stabilizer),
+ projected_specific_body_force_vectors,
+ mesh_space_dimension,
+ *aperture_size_parameter};
 
  SecondaryVariableCollection secondary_variables;
 

ProcessLib/HT/HTFEM.h

@@ -66,6 +66,11 @@ class HTFEM : public HTLocalAssemblerInterface
  _integration_method.getNumberOfPoints();
  _ip_data.reserve(n_integration_points);
 
+ ParameterLib::SpatialPosition pos;
+ pos.setElementID(_element.getID());
+
+ double const aperture_size = _process_data.aperture_size(0.0, pos)[0];
+
  auto const shape_matrices =
  NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
  GlobalDim>(element, is_axially_symmetric,
@@ -77,7 +82,7 @@ class HTFEM : public HTLocalAssemblerInterface
  shape_matrices[ip].N, shape_matrices[ip].dNdx,
  _integration_method.getWeightedPoint(ip).getWeight() *
  shape_matrices[ip].integralMeasure *
- shape_matrices[ip].detJ);
+ shape_matrices[ip].detJ * aperture_size);
  }
  }
 
@@ -147,7 +152,9 @@ class HTFEM : public HTLocalAssemblerInterface
  liquid_phase
  .property(MaterialPropertyLib::PropertyType::density)
  .template value<double>(vars, pos, t, dt);
- auto const b = this->_process_data.specific_body_force;
+ auto const b =
+ this->_process_data.projected_specific_body_force_vectors
+ [this->_element.getID()];
  q += K_over_mu * rho_w * b;
  }
 
@@ -300,7 +307,9 @@ class HTFEM : public HTLocalAssemblerInterface
  liquid_phase
  .property(MaterialPropertyLib::PropertyType::density)
  .template value<double>(vars, pos, t, dt);
- auto const b = _process_data.specific_body_force;
+ auto const b =
+ _process_data.projected_specific_body_force_vectors
+ [_element.getID()];
  // here it is assumed that the vector b is directed 'downwards'
  cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
  }

ProcessLib/HT/HTProcess.cpp

@@ -48,24 +48,26 @@ void HTProcess::initializeConcreteProcess(
  MeshLib::Mesh const& mesh,
  unsigned const integration_order)
 {
+ int const mesh_space_dimension = _process_data.mesh_space_dimension;
+
  if (_use_monolithic_scheme)
  {
  ProcessLib::createLocalAssemblers<MonolithicHTFEM>(
- mesh.getDimension(), mesh.getElements(), dof_table,
+ mesh_space_dimension, mesh.getElements(), dof_table,
  _local_assemblers, NumLib::IntegrationOrder{integration_order},
  mesh.isAxiallySymmetric(), _process_data);
  }
  else
  {
  ProcessLib::createLocalAssemblers<StaggeredHTFEM>(
- mesh.getDimension(), mesh.getElements(), dof_table,
+ mesh_space_dimension, mesh.getElements(), dof_table,
  _local_assemblers, NumLib::IntegrationOrder{integration_order},
  mesh.isAxiallySymmetric(), _process_data);
  }
 
  _secondary_variables.addSecondaryVariable(
  "darcy_velocity",
- makeExtrapolator(mesh.getDimension(), getExtrapolator(),
+ makeExtrapolator(mesh_space_dimension, getExtrapolator(),
  _local_assemblers,
  &HTLocalAssemblerInterface::getIntPtDarcyVelocity));
 }

ProcessLib/HT/HTProcessData.h

@@ -16,6 +16,7 @@
 
 #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
 #include "NumLib/NumericalStability/NumericalStabilization.h"
+#include "ParameterLib/ConstantParameter.h"
 #include "ParameterLib/Parameter.h"
 
 namespace ProcessLib
@@ -29,12 +30,21 @@ struct HTProcessData final
  bool const has_fluid_thermal_expansion;
  ParameterLib::Parameter<double> const& solid_thermal_expansion;
  ParameterLib::Parameter<double> const& biot_constant;
- Eigen::VectorXd const specific_body_force;
+
  bool const has_gravity;
  int const heat_transport_process_id;
  int const hydraulic_process_id;
 
  NumLib::NumericalStabilization stabilizer;
+
+ /// Projected specific body force vector: R * R^T * b.
+ std::vector<Eigen::VectorXd> const projected_specific_body_force_vectors;
+ int const mesh_space_dimension;
+
+ /// The aperture size is the thickness of 2D element or the
+ /// cross section area of 1D element. For 3D element, the value is set to 1.
+ ParameterLib::Parameter<double> const& aperture_size =
+ ParameterLib::ConstantParameter<double>("constant_one", 1.0);
 };
 
 } // namespace HT

ProcessLib/HT/MonolithicHTFEM.h

@@ -108,7 +108,9 @@ class MonolithicHTFEM : public HTFEM<ShapeFunction, GlobalDim>
  auto const& liquid_phase = medium.phase("AqueousLiquid");
  auto const& solid_phase = medium.phase("Solid");
 
- auto const& b = process_data.specific_body_force;
+ auto const& b =
+ process_data
+ .projected_specific_body_force_vectors[this->_element.getID()];
 
  MaterialPropertyLib::VariableArray vars;
 

ProcessLib/HT/StaggeredHTFEM-impl.h

@@ -70,7 +70,9 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHydraulicEquation(
  auto const& liquid_phase = medium.phase("AqueousLiquid");
  auto const& solid_phase = medium.phase("Solid");
 
- auto const& b = process_data.specific_body_force;
+ auto const& b =
+ process_data
+ .projected_specific_body_force_vectors[this->_element.getID()];
 
  MaterialPropertyLib::VariableArray vars;
 
@@ -194,7 +196,9 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHeatTransportEquation(
  *process_data.media_map->getMedium(this->_element.getID());
  auto const& liquid_phase = medium.phase("AqueousLiquid");
 
- auto const& b = process_data.specific_body_force;
+ auto const& b =
+ process_data
+ .projected_specific_body_force_vectors[this->_element.getID()];
 
  MaterialPropertyLib::VariableArray vars;
 

ProcessLib/HT/Tests.cmake

@@ -380,6 +380,7 @@ endif()
 
 if (NOT (OGS_USE_MPI))
  OgsTest(PROJECTFILE Parabolic/HT/SimpleSynthetics/deactivated_subdomain/HT_DeactivatedSubdomain.prj)
+ OgsTest(PROJECTFILE Parabolic/HT/LowerDimensionalFracture/2D_single_fracture_HT.prj RUNTIME 55)
 endif()
 
 AddTest(