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(