[LF] Changed the type of a member from unique_ptr reference to class reference

and Removed a member, _current_material_id for thread safe.

ProcessLib/LiquidFlow/CreateLiquidFlowMaterialProperties.h

@@ -10,7 +10,6 @@
  *   Created on December 14, 2016, 1:20 PM
  */
 
-
 #pragma once
 
 #include <memory>

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h

@@ -30,27 +30,30 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
 {
     SpatialPosition pos;
     pos.setElementID(_element.getID());
-    _material_properties->setMaterialID(pos);
+    const int material_id = _material_properties.getMaterialID(pos);
 
-    const Eigen::MatrixXd& perm =
-        _material_properties->getPermeability(t, pos, _element.getDimension());
+    const Eigen::MatrixXd& perm = _material_properties.getPermeability(
+        material_id, t, pos, _element.getDimension());
     // Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
     //  the assert must be changed to perm.rows() == _element->getDimension()
     assert(perm.rows() == GlobalDim || perm.rows() == 1);
 
     if (perm.size() == 1)  // isotropic or 1D problem.
-        local_assemble<IsotropicCalculator>(
-            t, local_x, local_M_data, local_K_data, local_b_data, pos, perm);
+        local_assemble<IsotropicCalculator>(material_id, t, local_x,
+                                            local_M_data, local_K_data,
+                                            local_b_data, pos, perm);
     else
-        local_assemble<AnisotropicCalculator>(
-            t, local_x, local_M_data, local_K_data, local_b_data, pos, perm);
+        local_assemble<AnisotropicCalculator>(material_id, t, local_x,
+                                              local_M_data, local_K_data,
+                                              local_b_data, pos, perm);
 }
 
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 template <typename LaplacianGravityVelocityCalculator>
 void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
-    local_assemble(double const t, std::vector<double> const& local_x,
+    local_assemble(const int material_id, double const t,
+                   std::vector<double> const& local_x,
                    std::vector<double>& local_M_data,
                    std::vector<double>& local_K_data,
                    std::vector<double>& local_b_data,
@@ -86,17 +89,17 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
             sm.integralMeasure * sm.detJ * wp.getWeight();
 
         // Assemble mass matrix, M
-        local_M.noalias() +=
-            _material_properties->getMassCoefficient(
-                t, pos, porosity_variable, storage_variable, p, _temperature) *
-            sm.N.transpose() * sm.N * integration_factor;
+        local_M.noalias() += _material_properties.getMassCoefficient(
+                                 material_id, t, pos, porosity_variable,
+                                 storage_variable, p, _temperature) *
+                             sm.N.transpose() * sm.N * integration_factor;
 
         // Compute density:
         const double rho_g =
-            _material_properties->getLiquidDensity(p, _temperature) *
+            _material_properties.getLiquidDensity(p, _temperature) *
             _gravitational_acceleration;
         // Compute viscosity:
-        const double mu = _material_properties->getViscosity(p, _temperature);
+        const double mu = _material_properties.getViscosity(p, _temperature);
 
         // Assemble Laplacian, K, and RHS by the gravitational term
         LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
@@ -113,9 +116,9 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
 {
     SpatialPosition pos;
     pos.setElementID(_element.getID());
-    _material_properties->setMaterialID(pos);
-    const Eigen::MatrixXd& perm =
-        _material_properties->getPermeability(t, pos, _element.getDimension());
+    const int material_id = _material_properties.getMaterialID(pos);
+    const Eigen::MatrixXd& perm = _material_properties.getPermeability(
+        material_id, t, pos, _element.getDimension());
 
     // Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
     //  the assert must be changed to perm.rows() == _element->getDimension()
@@ -155,10 +158,10 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
         // TODO : compute _temperature from the heat transport pcs
 
         const double rho_g =
-            _material_properties->getLiquidDensity(p, _temperature) *
+            _material_properties.getLiquidDensity(p, _temperature) *
             _gravitational_acceleration;
         // Compute viscosity:
-        const double mu = _material_properties->getViscosity(p, _temperature);
+        const double mu = _material_properties.getViscosity(p, _temperature);
 
         LaplacianGravityVelocityCalculator::calculateVelocity(
             _darcy_velocities, local_p_vec, sm, perm, ip, mu, rho_g,

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h

@@ -69,8 +69,7 @@ public:
         unsigned const integration_order,
         int const gravitational_axis_id,
         double const gravitational_acceleration,
-        std::unique_ptr<LiquidFlowMaterialProperties> const&
-            material_propertries)
+        LiquidFlowMaterialProperties const& material_propertries)
         : _element(element),
           _integration_method(integration_order),
           _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
@@ -174,7 +173,8 @@ private:
     };
 
     template <typename LaplacianGravityVelocityCalculator>
-    void local_assemble(double const t, std::vector<double> const& local_x,
+    void local_assemble(const int material_id, double const t,
+                        std::vector<double> const& local_x,
                         std::vector<double>& local_M_data,
                         std::vector<double>& local_K_data,
                         std::vector<double>& local_b_data,
@@ -189,7 +189,7 @@ private:
 
     const int _gravitational_axis_id;
     const double _gravitational_acceleration;
-    const std::unique_ptr<LiquidFlowMaterialProperties>& _material_properties;
+    const LiquidFlowMaterialProperties& _material_properties;
     double _temperature;
 };
 

ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp

@@ -33,16 +33,16 @@ namespace ProcessLib
 {
 namespace LiquidFlow
 {
-void LiquidFlowMaterialProperties::setMaterialID(const SpatialPosition& pos)
+int LiquidFlowMaterialProperties::getMaterialID(
+    const SpatialPosition& pos) const
 {
     if (!_has_material_ids)
     {
-        _current_material_id = 0;
-        return;
+        return 0;
     }
 
     assert(pos.getElementID().get() < _material_ids.size());
-    _current_material_id = _material_ids[pos.getElementID().get()];
+    return _material_ids[pos.getElementID().get()];
 }
 
 double LiquidFlowMaterialProperties::getLiquidDensity(const double p,
@@ -66,8 +66,8 @@ double LiquidFlowMaterialProperties::getViscosity(const double p,
 }
 
 double LiquidFlowMaterialProperties::getMassCoefficient(
-    const double /*t*/, const SpatialPosition& /*pos*/, const double p,
-    const double T, const double porosity_variable,
+    const int material_id, const double /*t*/, const SpatialPosition& /*pos*/,
+    const double p, const double T, const double porosity_variable,
     const double storage_variable) const
 {
     ArrayType vars;
@@ -81,16 +81,17 @@ double LiquidFlowMaterialProperties::getMassCoefficient(
     assert(rho > 0.);
 
     const double porosity =
-        _porosity_models[_current_material_id]->getValue(porosity_variable, T);
+        _porosity_models[material_id]->getValue(porosity_variable, T);
     const double storage =
-        _storage_models[_current_material_id]->getValue(storage_variable);
+        _storage_models[material_id]->getValue(storage_variable);
     return porosity * drho_dp / rho + storage;
 }
 
 Eigen::MatrixXd const& LiquidFlowMaterialProperties::getPermeability(
-    const double /*t*/, const SpatialPosition& /*pos*/, const int /*dim*/) const
+    const int material_id, const double /*t*/, const SpatialPosition& /*pos*/,
+    const int /*dim*/) const
 {
-    return _intrinsic_permeability_models[_current_material_id];
+    return _intrinsic_permeability_models[material_id];
 }
 
 }  // end of namespace

ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h

@@ -43,6 +43,9 @@ class PropertyVector;
 
 namespace ProcessLib
 {
+template <typename T>
+struct Parameter;
+
 class SpatialPosition;
 
 namespace LiquidFlow
@@ -75,7 +78,7 @@ public:
     {
     }
 
-    void setMaterialID(const SpatialPosition& pos);
+    int getMaterialID(const SpatialPosition& pos) const;
 
     /**
      * \brief Compute the coefficient of the mass term by
@@ -84,6 +87,7 @@ public:
      *      \f]
      *     where \f$n\f$ is the porosity, \f$rho_l\f$ is the liquid density,
      *     \f$bata_s\f$ is the storage.
+     * \param material_id        Material index.
      * \param t                  Time.
      * \param pos                Position of element.
      * \param p                  Pressure value.
@@ -93,12 +97,13 @@ public:
      *                           saturation, and invariant of stress or strain.
      * \param storage_variable   Variable for storage model.
      */
-    double getMassCoefficient(const double t, const SpatialPosition& pos,
-                              const double p, const double T,
-                              const double porosity_variable,
+    double getMassCoefficient(const int material_id, const double t,
+                              const SpatialPosition& pos, const double p,
+                              const double T, const double porosity_variable,
                               const double storage_variable) const;
 
-    Eigen::MatrixXd const& getPermeability(const double t,
+    Eigen::MatrixXd const& getPermeability(const int material_id,
+                                           const double t,
                                            const SpatialPosition& pos,
                                            const int dim) const;
 
@@ -106,11 +111,6 @@ public:
 
     double getViscosity(const double p, const double T) const;
 
-    MaterialLib::Fluid::FluidProperties* getFluidProperties() const
-    {
-        return _fluid_properties.get();
-    }
-
 private:
     /// A flag to indicate whether the reference member, _material_ids,
     /// is not assigned.
@@ -120,8 +120,6 @@ private:
      */
     MeshLib::PropertyVector<int> const& _material_ids;
 
-    int _current_material_id = 0;
-
     const std::unique_ptr<MaterialLib::Fluid::FluidProperties>
         _fluid_properties;
 

ProcessLib/LiquidFlow/LiquidFlowProcess.cpp

@@ -62,7 +62,7 @@ void LiquidFlowProcess::initializeConcreteProcess(
         mesh.getDimension(), mesh.getElements(), dof_table,
         pv.getShapeFunctionOrder(), _local_assemblers,
         mesh.isAxiallySymmetric(), integration_order, _gravitational_axis_id,
-        _gravitational_acceleration, _material_properties);
+        _gravitational_acceleration, *_material_properties);
 
     _secondary_variables.addSecondaryVariable(
         "darcy_velocity_x", 1,

ProcessLib/LiquidFlow/LiquidFlowProcess.h

@@ -76,12 +76,6 @@ public:
                                           GlobalVector const& x) override;
 
     bool isLinear() const override { return true; }
-
-    MaterialLib::Fluid::FluidProperties* getFluidProperties() const
-    {
-        return _material_properties->getFluidProperties();
-    }
-
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,

Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp

@@ -86,7 +86,7 @@ TEST(MaterialProcessLibLiquidFlow, checkLiquidFlowMaterialProperties)
     pos.setElementID(0);
 
     // Check permeability
-    const Eigen::MatrixXd& perm = lprop->getPermeability(0., pos, 1);
+    const Eigen::MatrixXd& perm = lprop->getPermeability(0, 0., pos, 1);
     ASSERT_EQ(2.e-10, perm(0, 0));
     ASSERT_EQ(0., perm(0, 1));
     ASSERT_EQ(0., perm(0, 2));
@@ -99,6 +99,7 @@ TEST(MaterialProcessLibLiquidFlow, checkLiquidFlowMaterialProperties)
 
     const double T = 273.15 + 60.0;
     const double p = 1.e+6;
-    const double mass_coef = lprop->getMassCoefficient(0., pos, p, T, 0., 0.);
+    const double mass_coef =
+        lprop->getMassCoefficient(0, 0., pos, p, T, 0., 0.);
     ASSERT_NEAR(0.000100000093, mass_coef, 1.e-10);
 }