[PL/LF] Init. impl. of LFLocalAssembler::getFlux.

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h

@@ -54,6 +54,47 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
     }
 }
 
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+Eigen::Vector3d
+LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::getFlux(
+    MathLib::Point3d const& p_local_coords, double const t,
+    std::vector<double> const& local_x) const
+{
+    // eval dNdx and invJ at p
+    auto const fe =
+        NumLib::createIsoparametricFiniteElement<ShapeFunction,
+                                                 ShapeMatricesType>(_element);
+
+    typename ShapeMatricesType::ShapeMatrices shape_matrices(
+        ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
+
+    // Note: Axial symmetry is set to false here, because we only need dNdx
+    // here, which is not affected by axial symmetry.
+    fe.computeShapeFunctions(p_local_coords.getCoords(), shape_matrices,
+                             GlobalDim, false);
+
+    // create pos object to access the correct media property
+    ParameterLib::SpatialPosition pos;
+    pos.setElementID(_element.getID());
+    const int material_id = _material_properties.getMaterialID(pos);
+
+    double pressure = 0.0;
+    NumLib::shapeFunctionInterpolate(local_x, shape_matrices.N, pressure);
+    const Eigen::MatrixXd& permeability = _material_properties.getPermeability(
+        material_id, t, pos, _element.getDimension(), pressure,
+        _reference_temperature);
+    const double mu =
+        _material_properties.getViscosity(pressure, _reference_temperature);
+
+    Eigen::Vector3d flux;
+    flux.head<GlobalDim>() =
+        -permeability / mu * shape_matrices.dNdx *
+        Eigen::Map<const NodalVectorType>(local_x.data(), local_x.size());
+
+    return flux;
+}
+
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 template <typename LaplacianGravityVelocityCalculator>

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h

@@ -128,6 +128,12 @@ public:
                   std::vector<double>& local_K_data,
                   std::vector<double>& local_b_data) override;
 
+    /// Computes the flux in the point \c p_local_coords that is given in local
+    /// coordinates using the values from \c local_x.
+    Eigen::Vector3d getFlux(MathLib::Point3d const& p_local_coords,
+                            double const t,
+                            std::vector<double> const& local_x) const override;
+
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {

ProcessLib/LiquidFlow/LiquidFlowProcess.cpp

@@ -125,5 +125,21 @@ void LiquidFlowProcess::computeSecondaryVariableConcrete(const double t,
         x, _coupled_solutions);
 }
 
+Eigen::Vector3d LiquidFlowProcess::getFlux(
+    std::size_t const element_id,
+    MathLib::Point3d const& p,
+    double const t,
+    std::vector<GlobalVector*> const& x) const
+{
+    // fetch local_x from primary variable
+    std::vector<GlobalIndexType> indices_cache;
+    auto const r_c_indices = NumLib::getRowColumnIndices(
+        element_id, *_local_to_global_index_map, indices_cache);
+    constexpr int process_id = 0;  // monolithic scheme.
+    std::vector<double> local_x(x[process_id]->get(r_c_indices.rows));
+
+    return _local_assemblers[element_id]->getFlux(p, t, local_x);
+}
+
 }  // namespace LiquidFlow
 }  // namespace ProcessLib

ProcessLib/LiquidFlow/LiquidFlowProcess.h

@@ -78,6 +78,12 @@ public:
                                           int const process_id) override;
 
     bool isLinear() const override { return true; }
+
+    Eigen::Vector3d getFlux(std::size_t const element_id,
+                            MathLib::Point3d const& p,
+                            double const t,
+                            std::vector<GlobalVector*> const& x) const override;
+
     int getGravitationalAxisID() const { return _gravitational_axis_id; }
     double getGravitationalAcceleration() const
     {