[PL/HC] Include computation of Darcy vel.

ProcessLib/HC/HCFEM.h

@@ -257,6 +257,58 @@ public:
         }
     }
 
+    void computeSecondaryVariableConcrete(
+        double const t, std::vector<double> const& local_x) override
+    {
+        SpatialPosition pos;
+        pos.setElementID(_element.getID());
+
+        auto const& K =
+            _process_data.porous_media_properties.getIntrinsicPermeability(t,
+                                                                           pos);
+        MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+        auto const mu = _process_data.viscosity_model->getValue(vars);
+        GlobalDimMatrixType const K_over_mu = K / mu;
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
+            &local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
+
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
+        {
+            auto const& ip_data = _ip_data[ip];
+            auto const& N = ip_data.N;
+            auto const& dNdx = ip_data.dNdx;
+
+            GlobalDimVectorType velocity = -K_over_mu * dNdx * p_nodal_values;
+            if (_process_data.has_gravity)
+            {
+                double C_int_pt = 0.0;
+                double p_int_pt = 0.0;
+                NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt,
+                                                 p_int_pt);
+                vars[static_cast<int>(
+                    MaterialLib::Fluid::PropertyVariableType::C)] = C_int_pt;
+                vars[static_cast<int>(
+                    MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
+
+                auto const rho_w = _process_data.fluid_density->getValue(vars);
+                auto const b = _process_data.specific_body_force;
+                // here it is assumed that the vector b is directed 'downwards'
+                velocity += K_over_mu * rho_w * b;
+            }
+
+            for (unsigned d = 0; d < GlobalDim; ++d)
+            {
+                _darcy_velocities[d][ip] = velocity[d];
+            }
+        }
+    }
+
+
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {

ProcessLib/HC/HCProcess.cpp

@@ -98,6 +98,16 @@ void HCProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac, coupling_term);
 }
 
+void HCProcess::computeSecondaryVariableConcrete(
+    double const t, GlobalVector const& x,
+    StaggeredCouplingTerm const& coupling_term)
+{
+    DBUG("Compute the Darcy velocity for HCProcess.");
+    GlobalExecutor::executeMemberOnDereferenced(
+        &HCLocalAssemblerInterface::computeSecondaryVariable,
+        _local_assemblers, *_local_to_global_index_map, t, x, coupling_term);
+}
+
 }  // namespace HC
 }  // namespace ProcessLib
 

ProcessLib/HC/HCProcess.h

@@ -57,6 +57,10 @@ public:
     bool isLinear() const override { return false; }
     //! @}
 
+    void computeSecondaryVariableConcrete(
+        double const t, GlobalVector const& x,
+        StaggeredCouplingTerm const& coupling_term) override;
+
 private:
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,