diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index 4a198995d3d0d73bec23d0535aac4d9af84d5246..eb7733650812342ab7c8c54edfe43ac91d16f94c 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -105,7 +105,7 @@ public:
     /// Computes the flux in the point \c p_local_coords that is given in local
     /// coordinates using the values from \c local_x.
     // TODO add time dependency
-    std::vector<double> getFlux(
+    Eigen::Vector3d getFlux(
         MathLib::Point3d const& p_local_coords,
         std::vector<double> const& local_x) const override
     {
@@ -123,8 +123,6 @@ public:
         // here, which is not affected by axial symmetry.
         fe.computeShapeFunctions(p_local_coords.getCoords(), shape_matrices,
                                  GlobalDim, false);
-        std::vector<double> flux;
-        flux.resize(3);
 
         // fetch hydraulic conductivity
         SpatialPosition pos;
@@ -133,9 +131,10 @@ public:
         double const t = 0.0;
         auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
 
-        Eigen::Map<Eigen::RowVectorXd>(flux.data(), flux.size()) =
+        Eigen::Vector3d flux;
+        flux.head<GlobalDim>() =
             -k * shape_matrices.dNdx *
-            Eigen::Map<const Eigen::VectorXd>(local_x.data(), local_x.size());
+            Eigen::Map<const NodalVectorType>(local_x.data(), local_x.size());
 
         return flux;
     }
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index eb58aa8e57d9fe0212775e4ee89d615130c0542c..5ce132d81af27f7de52b8dc6054649de298f10a9 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -47,7 +47,7 @@ public:
     bool isLinear() const override { return true; }
     //! @}
 
-    std::vector<double> getFlux(std::size_t element_id,
+    Eigen::Vector3d getFlux(std::size_t element_id,
                                 MathLib::Point3d const& p,
                                 GlobalVector const& x) const override
     {
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 3cc9cfca86f20e7f441f3af6ce7305f5a624dde7..ce156e71f35be02b9a46ab5d4f38a095ad7c70be 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -91,7 +91,7 @@ public:
     /// Computes the flux in the point \c pnt_local_coords that is given in
     /// local coordinates using the values from \c local_x.
     // TODO add time dependency
-    std::vector<double> getFlux(
+    Eigen::Vector3d getFlux(
         MathLib::Point3d const& pnt_local_coords,
         std::vector<double> const& local_x) const override
     {
@@ -149,10 +149,9 @@ public:
             auto const b = this->_material_properties.specific_body_force;
             q += K_over_mu * rho_w * b;
         }
-        std::vector<double> flux;
-        flux.resize(3);
-        Eigen::Map<GlobalDimVectorType>(flux.data(), flux.size()) = q;
 
+        Eigen::Vector3d flux;
+        flux.head<GlobalDim>() = q;
         return flux;
     }
 
diff --git a/ProcessLib/HT/HTLocalAssemblerInterface.h b/ProcessLib/HT/HTLocalAssemblerInterface.h
index 53117d6941206bbf139f851e14e9db2819614855..fbcd85a711d8bec0067a4eb744802b261e1068c8 100644
--- a/ProcessLib/HT/HTLocalAssemblerInterface.h
+++ b/ProcessLib/HT/HTLocalAssemblerInterface.h
@@ -59,7 +59,7 @@ public:
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
         std::vector<double>& /*cache*/) const = 0;
 
-    virtual std::vector<double> getFlux(
+    virtual Eigen::Vector3d getFlux(
         MathLib::Point3d const& pnt_local_coords,
         std::vector<double> const& local_x) const = 0;
 
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index 41159d516937cf4fb9a5220a9d258cd909573d1c..969b073abe99d6aba0dc079334fbbd5a4aeccdbf 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -240,9 +240,9 @@ void HTProcess::setCoupledSolutionsOfPreviousTimeStep()
     }
 }
 
-std::vector<double> HTProcess::getFlux(std::size_t element_id,
-                                       MathLib::Point3d const& p,
-                                       GlobalVector const& x) const
+Eigen::Vector3d HTProcess::getFlux(std::size_t element_id,
+                                   MathLib::Point3d const& p,
+                                   GlobalVector const& x) const
 {
     // fetch local_x from primary variable
     std::vector<GlobalIndexType> indices_cache;
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index 3972add55f2da571e2d7abf604bb7ed320d328bd..ab13614ba7b8cf025b5a5d48b14f83290b01d2e7 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -69,7 +69,7 @@ public:
     bool isLinear() const override { return false; }
     //! @}
 
-    std::vector<double> getFlux(std::size_t element_id,
+    Eigen::Vector3d getFlux(std::size_t element_id,
                                 MathLib::Point3d const& p,
                                 GlobalVector const& x) const override;
 
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index f750b6aa2a71b281b63bc99dabfb8f222f8a6ab8..aa24d8d390730ba6255e82ea771608fdd27db970 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -85,11 +85,11 @@ public:
 
     /// Computes the flux in the point \c p_local_coords that is given in local
     /// coordinates using the values from \c local_x.
-    virtual std::vector<double> getFlux(
+    virtual Eigen::Vector3d getFlux(
         MathLib::Point3d const& /*p_local_coords*/,
         std::vector<double> const& /*local_x*/) const
     {
-        return std::vector<double>();
+        return Eigen::Vector3d{};
     }
 
 private:
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 10bdeb8ba4c1a48214bcde2cbd9e0b8afd8b8939..17718362a435b56ecd34dc0890077b76560814f4 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -136,11 +136,11 @@ public:
 
     // Used as a call back for CalculateSurfaceFlux process.
 
-    virtual std::vector<double> getFlux(std::size_t /*element_id*/,
-                                        MathLib::Point3d const& /*p*/,
-                                        GlobalVector const& /*x*/) const
+    virtual Eigen::Vector3d getFlux(std::size_t /*element_id*/,
+                                    MathLib::Point3d const& /*p*/,
+                                    GlobalVector const& /*x*/) const
     {
-        return std::vector<double>{};
+        return Eigen::Vector3d{};
     }
 
 protected: