WIP vectorial extrapolation

ProcessLib/ComponentTransport/ComponentTransportFEM.h

@@ -285,6 +285,8 @@ public:
             auto const& N = ip_data.N;
             auto const& dNdx = ip_data.dNdx;
 
+            pos.setIntegrationPoint(ip);
+
             auto const& K =
                 _process_data.porous_media_properties.getIntrinsicPermeability(
                     t, pos);

ProcessLib/HydroMechanics/HydroMechanicsFEM.h

@@ -15,6 +15,7 @@
 #include "MaterialLib/SolidModels/KelvinVector.h"
 #include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
@@ -364,43 +365,6 @@ public:
         }
     }
 
-    void postTimestepConcrete(std::vector<double> const& local_x) override
-    {
-        double const& t = _process_data.t;
-
-        auto p =
-            Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
-                pressure_size> const>(local_x.data() + pressure_index,
-                                      pressure_size);
-        using GlobalDimVectorType =
-            typename ShapeMatricesTypePressure::GlobalDimVectorType;
-
-        unsigned const n_integration_points =
-            _integration_method.getNumberOfPoints();
-
-        SpatialPosition x_position;
-        x_position.setElementID(_element.getID());
-        for (unsigned ip = 0; ip < n_integration_points; ip++)
-        {
-            x_position.setIntegrationPoint(ip);
-            double const K_over_mu =
-                _process_data.intrinsic_permeability(t, x_position)[0] /
-                _process_data.fluid_viscosity(t, x_position)[0];
-
-            auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
-            auto const& b = _process_data.specific_body_force;
-
-            // Compute the velocity
-            auto const& dNdx_p = _ip_data[ip].dNdx_p;
-            GlobalDimVectorType const darcy_velocity =
-                -K_over_mu * dNdx_p * p - K_over_mu * rho_fr * b;
-            for (unsigned d = 0; d < DisplacementDim; ++d)
-            {
-                _darcy_velocities[d][ip] = darcy_velocity[d];
-            }
-        }
-    }
-
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {
@@ -522,34 +486,54 @@ public:
         return getIntPtEpsilon(cache, 5);
     }
 
-    std::vector<double> const& getIntPtDarcyVelocityX(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
+    std::vector<double> const& getIntPtDarcyVelocity(
+        const double t,
+        GlobalVector const& current_solution,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
+        std::vector<double>& cache) const override
     {
-        assert(!_darcy_velocities.empty());
-        return _darcy_velocities[0];
-    }
+        auto const num_intpts = _ip_data.size();
 
-    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];
-    }
+        auto const indices = NumLib::getIndices(_element.getID(), dof_table);
+        assert(!indices.empty());
+        auto const local_x = current_solution.get(indices);
 
-    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];
+        cache.clear();
+        auto cache_vec = MathLib::createZeroedMatrix<Eigen::Matrix<
+            double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, DisplacementDim, num_intpts);
+
+        SpatialPosition pos;
+        pos.setElementID(_element.getID());
+
+        auto p =
+            Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
+                pressure_size> const>(local_x.data() + pressure_index,
+                                      pressure_size);
+        using GlobalDimVectorType =
+            typename ShapeMatricesTypePressure::GlobalDimVectorType;
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        SpatialPosition x_position;
+        x_position.setElementID(_element.getID());
+        for (unsigned ip = 0; ip < n_integration_points; ip++) {
+            x_position.setIntegrationPoint(ip);
+            double const K_over_mu =
+                _process_data.intrinsic_permeability(t, x_position)[0] /
+                _process_data.fluid_viscosity(t, x_position)[0];
+
+            auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
+            auto const& b = _process_data.specific_body_force;
+
+            // Compute the velocity
+            auto const& dNdx_p = _ip_data[ip].dNdx_p;
+            cache_vec.col(ip).noalias() =
+                -K_over_mu * dNdx_p * p - K_over_mu * rho_fr * b;
+        }
+
+        return cache;
     }
 
 private:
@@ -602,11 +586,6 @@ private:
         typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType>
         _secondary_data;
 
-    std::vector<std::vector<double>> _darcy_velocities =
-        std::vector<std::vector<double>>(
-            DisplacementDim,
-            std::vector<double>(_integration_method.getNumberOfPoints()));
-
     static const int pressure_index = 0;
     static const int pressure_size = ShapeFunctionPressure::NPOINTS;
     static const int displacement_index = ShapeFunctionPressure::NPOINTS;

ProcessLib/HydroMechanics/HydroMechanicsProcess-impl.h

@@ -164,22 +164,10 @@ void HydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
             &LocalAssemblerInterface::getIntPtEpsilonXY));
 
     Base::_secondary_variables.addSecondaryVariable(
-        "velocity_x",
+        "velocity",
         makeExtrapolator(
             1, getExtrapolator(), _local_assemblers,
-            &LocalAssemblerInterface::getIntPtDarcyVelocityX));
-
-    Base::_secondary_variables.addSecondaryVariable(
-        "velocity_y",
-        makeExtrapolator(
-            1, getExtrapolator(), _local_assemblers,
-            &LocalAssemblerInterface::getIntPtDarcyVelocityY));
-
-    Base::_secondary_variables.addSecondaryVariable(
-        "velocity_z",
-        makeExtrapolator(
-            1, getExtrapolator(), _local_assemblers,
-            &LocalAssemblerInterface::getIntPtDarcyVelocityZ));
+            &LocalAssemblerInterface::getIntPtDarcyVelocity));
 }
 
 template <int DisplacementDim>

ProcessLib/HydroMechanics/LocalAssemblerInterface.h

@@ -91,19 +91,7 @@ struct LocalAssemblerInterface
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
         std::vector<double>& cache) const = 0;
 
-    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(
+    virtual std::vector<double> const& getIntPtDarcyVelocity(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h

@@ -19,6 +19,7 @@
 
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
@@ -39,19 +40,7 @@ class LiquidFlowLocalAssemblerInterface
       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(
+    virtual std::vector<double> const& getIntPtDarcyVelocity(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
@@ -122,34 +111,43 @@ 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.empty());
-        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
+    std::vector<double> const& getIntPtDarcyVelocity(
+        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];
+        // auto const num_nodes = ShapeFunction_::NPOINTS;
+        auto const num_intpts = _shape_matrices.size();
+
+        auto const indices = NumLib::getIndices(
+            _element.getID(), dof_table);
+        assert(!indices.empty());
+        auto const local_x = current_solution.get(indices);
+        auto const local_x_vec =
+            MathLib::toVector<Eigen::Matrix<double, ShapeFunction::NPOINTS, 1>>(
+                local_x, ShapeFunction::NPOINTS);
+
+        cache.clear();
+        auto cache_vec = MathLib::createZeroedMatrix<
+            Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, GlobalDim, num_intpts);
+
+        SpatialPosition pos;
+        pos.setElementID(_element.getID());
+
+        // TODO fix
+#if 0
+        for (unsigned i = 0; i < num_intpts; ++i) {
+            pos.setIntegrationPoint(i);
+            auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
+            // dimensions: (d x 1) = (d x n) * (n x 1)
+            cache_vec.col(i).noalias() =
+                -k * _shape_matrices[i].dNdx * local_x_vec;
+        }
+#endif
+
+        return cache;
     }
 
 private:

ProcessLib/LiquidFlow/LiquidFlowProcess.cpp

@@ -69,27 +69,10 @@ void LiquidFlowProcess::initializeConcreteProcess(
         *_material_properties);
 
     _secondary_variables.addSecondaryVariable(
-        "darcy_velocity_x",
-        makeExtrapolator(
-            1, getExtrapolator(), _local_assemblers,
-            &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
-
-    if (mesh.getDimension() > 1)
-    {
-        _secondary_variables.addSecondaryVariable(
-            "darcy_velocity_y",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityY));
-    }
-    if (mesh.getDimension() > 2)
-    {
-        _secondary_variables.addSecondaryVariable(
-            "darcy_velocity_z",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityZ));
-    }
+        "darcy_velocity",
+        makeExtrapolator(mesh.getDimension(),
+            getExtrapolator(), _local_assemblers,
+            &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocity));
 }
 
 void LiquidFlowProcess::assembleConcreteProcess(

ProcessLib/RichardsFlow/RichardsFlowFEM.h

@@ -14,6 +14,7 @@
 #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "MeshLib/CoordinateSystem.h"
+#include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -63,19 +64,7 @@ public:
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
         std::vector<double>& /*cache*/) const = 0;
 
-    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(
+    virtual std::vector<double> const& getIntPtDarcyVelocity(
         const double /*t*/,
         GlobalVector const& /*current_solution*/,
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
@@ -110,9 +99,6 @@ public:
           _process_data(process_data),
           _integration_method(integration_order),
           _saturation(
-              std::vector<double>(_integration_method.getNumberOfPoints())),
-          _darcy_velocities(
-              GlobalDim,
               std::vector<double>(_integration_method.getNumberOfPoints()))
     {
         // This assertion is valid only if all nodal d.o.f. use the same shape
@@ -240,9 +226,44 @@ public:
         }  // end of mass lumping
     }
 
-    void computeSecondaryVariableConcrete(
-        double const t, std::vector<double> const& local_x) override
+    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+        const unsigned integration_point) const override
     {
+        auto const& N = _ip_data[integration_point].N;
+
+        // assumes N is stored contiguously in memory
+        return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+    }
+
+    std::vector<double> const& getIntPtSaturation(
+        const double /*t*/,
+        GlobalVector const& /*current_solution*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::vector<double>& /*cache*/) const override
+    {
+        assert(!_saturation.empty());
+        return _saturation;
+    }
+
+    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 num_nodes = ShapeFunction_::NPOINTS;
+        auto const num_intpts = _shape_matrices.size();
+
+        auto const indices = NumLib::getIndices(
+            _element.getID(), dof_table);
+        assert(!indices.empty());
+        auto const local_x = current_solution.get(indices);
+
+        cache.clear();
+        auto cache_vec = MathLib::createZeroedMatrix<
+            Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+            cache, GlobalDim, num_intpts);
+
         SpatialPosition pos;
         pos.setElementID(_element.getID());
         const int material_id =
@@ -278,7 +299,7 @@ public:
             auto const mu = _process_data.material->getFluidViscosity(
                 p_int_pt, temperature);
             auto const K_mat_coeff = permeability * (k_rel / mu);
-            GlobalDimVectorType velocity =
+            cache_vec.col(ip).noalias() =
                 -K_mat_coeff * _ip_data[ip].dNdx * p_nodal_values;
             if (_process_data.has_gravity)
             {
@@ -288,63 +309,11 @@ public:
                 assert(body_force.size() == GlobalDim);
                 // here it is assumed that the vector body_force is directed
                 // 'downwards'
-                velocity += K_mat_coeff * rho_w * body_force;
-            }
-
-            for (unsigned d = 0; d < GlobalDim; ++d)
-            {
-                _darcy_velocities[d][ip] = velocity[d];
+                cache_vec.col(ip).noalias() += K_mat_coeff * rho_w * body_force;
             }
         }
-    }
-
-    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
-        const unsigned integration_point) const override
-    {
-        auto const& N = _ip_data[integration_point].N;
-
-        // assumes N is stored contiguously in memory
-        return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
-    }
 
-    std::vector<double> const& getIntPtSaturation(
-        const double /*t*/,
-        GlobalVector const& /*current_solution*/,
-        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
-        std::vector<double>& /*cache*/) const override
-    {
-        assert(!_saturation.empty());
-        return _saturation;
-    }
-
-    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.empty());
-        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];
+        return cache;
     }
 
 private:
@@ -361,7 +330,6 @@ private:
             NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType>>>
         _ip_data;
     std::vector<double> _saturation;
-    std::vector<std::vector<double>> _darcy_velocities;
 };
 
 }  // namespace RichardsFlow

ProcessLib/RichardsFlow/RichardsFlowProcess.cpp

@@ -55,27 +55,10 @@ void RichardsFlowProcess::initializeConcreteProcess(
             &RichardsFlowLocalAssemblerInterface::getIntPtSaturation));
 
     _secondary_variables.addSecondaryVariable(
-        "darcy_velocity_x",
+        "darcy_velocity",
         makeExtrapolator(
-            1, getExtrapolator(), _local_assemblers,
-            &RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
-
-    if (mesh.getDimension() > 1)
-    {
-        _secondary_variables.addSecondaryVariable(
-            "darcy_velocity_y",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocityY));
-    }
-    if (mesh.getDimension() > 2)
-    {
-        _secondary_variables.addSecondaryVariable(
-            "darcy_velocity_z",
-            makeExtrapolator(
-                1, getExtrapolator(), _local_assemblers,
-                &RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocityZ));
-    }
+            mesh.getDimension(), getExtrapolator(), _local_assemblers,
+            &RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocity));
 }
 
 void RichardsFlowProcess::assembleConcreteProcess(
@@ -107,15 +90,5 @@ void RichardsFlowProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac, _coupling_term);
 }
 
-void RichardsFlowProcess::computeSecondaryVariableConcrete(
-    double const t, GlobalVector const& x)
-{
-    DBUG("Compute the Darcy velocity for RichardsFlowProcess");
-    GlobalExecutor::executeMemberOnDereferenced(
-        &RichardsFlowLocalAssemblerInterface::computeSecondaryVariable,
-        _local_assemblers, *_local_to_global_index_map, t, x,
-        _coupling_term);
-}
-
 }  // namespace RichardsFlow
 }  // namespace ProcessLib

ProcessLib/RichardsFlow/RichardsFlowProcess.h

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