[ComponentTransport] Vectorial velocity output.

ProcessLib/ComponentTransport/ComponentTransportFEM.h

@@ -16,6 +16,7 @@
 #include "ComponentTransportProcessData.h"
 #include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -51,23 +52,11 @@ class ComponentTransportLocalAssemblerInterface
       public NumLib::ExtrapolatableElement
 {
 public:
-    virtual std::vector<double> const& getIntPtDarcyVelocityX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const = 0;
-
-    virtual std::vector<double> const& getIntPtDarcyVelocityY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const = 0;
-
-    virtual std::vector<double> const& getIntPtDarcyVelocityZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const = 0;
+    virtual std::vector<double> const& getIntPtDarcyVelocity(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const = 0;
 };
 
 template <typename ShapeFunction, typename IntegrationMethod,
@@ -100,10 +89,7 @@ public:
                        ComponentTransportProcessData const& process_data)
         : _element(element),
           _process_data(process_data),
-          _integration_method(integration_order),
-          _darcy_velocities(
-              GlobalDim,
-              std::vector<double>(_integration_method.getNumberOfPoints()))
+          _integration_method(integration_order)
     {
         // This assertion is valid only if all nodal d.o.f. use the same shape
         // matrices.
@@ -267,24 +253,29 @@ public:
         }
     }
 
-    void computeSecondaryVariableConcrete(
-        double const t, std::vector<double> const& local_x) override
+    std::vector<double> const& getIntPtDarcyVelocity(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const override
     {
+        auto const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        auto const indices = NumLib::getIndices(_element.getID(), dof_table);
+        assert(!indices.empty());
+        auto const local_x = current_solution.get(indices);
+
+        cache.clear();
+        auto cache_mat = MathLib::createZeroedMatrix<
+            Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, GlobalDim, n_integration_points);
+
         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.fluid_properties->getValue(
-            MaterialLib::Fluid::FluidPropertyType::Viscosity, 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);
 
@@ -294,7 +285,14 @@ public:
             auto const& N = ip_data.N;
             auto const& dNdx = ip_data.dNdx;
 
-            GlobalDimVectorType velocity = -K_over_mu * dNdx * p_nodal_values;
+            auto const& K =
+                _process_data.porous_media_properties.getIntrinsicPermeability(
+                    t, pos);
+            auto const mu = _process_data.fluid_properties->getValue(
+                MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+            GlobalDimMatrixType const K_over_mu = K / mu;
+
+            cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
             if (_process_data.has_gravity)
             {
                 double C_int_pt = 0.0;
@@ -310,14 +308,11 @@ public:
                     MaterialLib::Fluid::FluidPropertyType::Density, 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];
+                cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
             }
         }
+
+        return cache;
     }
 
 
@@ -330,36 +325,6 @@ public:
         return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
     }
 
-    std::vector<double> const& getIntPtDarcyVelocityX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(_darcy_velocities.size() > 0);
-        return _darcy_velocities[0];
-    }
-
-    std::vector<double> const& getIntPtDarcyVelocityY(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(_darcy_velocities.size() > 1);
-        return _darcy_velocities[1];
-    }
-
-    std::vector<double> const& getIntPtDarcyVelocityZ(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(_darcy_velocities.size() > 2);
-        return _darcy_velocities[2];
-    }
-
 private:
     MeshLib::Element const& _element;
     ComponentTransportProcessData const& _process_data;
@@ -370,7 +335,6 @@ private:
         Eigen::aligned_allocator<
             IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>>
         _ip_data;
-    std::vector<std::vector<double>> _darcy_velocities;
 };
 
 }  // namespace ComponentTransport

ProcessLib/ComponentTransport/ComponentTransportProcess.cpp

@@ -45,27 +45,10 @@ void ComponentTransportProcess::initializeConcreteProcess(
         mesh.isAxiallySymmetric(), integration_order, _process_data);
 
     _secondary_variables.addSecondaryVariable(
-        "darcy_velocity_x",
-        makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                         &ComponentTransportLocalAssemblerInterface::
-                             getIntPtDarcyVelocityX));
-
-    if (mesh.getDimension() > 1)
-    {
-        _secondary_variables.addSecondaryVariable(
-            "darcy_velocity_y",
-            makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                             &ComponentTransportLocalAssemblerInterface::
-                                 getIntPtDarcyVelocityY));
-    }
-    if (mesh.getDimension() > 2)
-    {
-        _secondary_variables.addSecondaryVariable(
-            "darcy_velocity_z",
-            makeExtrapolator(1, getExtrapolator(), _local_assemblers,
-                             &ComponentTransportLocalAssemblerInterface::
-                                 getIntPtDarcyVelocityZ));
-    }
+        "darcy_velocity",
+        makeExtrapolator(
+            mesh.getDimension(), getExtrapolator(), _local_assemblers,
+            &ComponentTransportLocalAssemblerInterface::getIntPtDarcyVelocity));
 }
 
 void ComponentTransportProcess::assembleConcreteProcess(
@@ -99,16 +82,6 @@ void ComponentTransportProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac, coupling_term);
 }
 
-void ComponentTransportProcess::computeSecondaryVariableConcrete(
-    double const t, GlobalVector const& x,
-    StaggeredCouplingTerm const& coupling_term)
-{
-    DBUG("Compute the Darcy velocity for ComponentTransportProcess.");
-    GlobalExecutor::executeMemberOnDereferenced(
-        &ComponentTransportLocalAssemblerInterface::computeSecondaryVariable,
-        _local_assemblers, *_local_to_global_index_map, t, x, coupling_term);
-}
-
 }  // namespace ComponentTransport
 }  // namespace ProcessLib
 

ProcessLib/ComponentTransport/ComponentTransportProcess.h

@@ -95,10 +95,6 @@ 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,