Fixing getMassMatrixFV with the new Foam2Eigen convention

Author: RomainTiphaigne

src/ITHACA_CORE/EigenFunctions/EigenFunctions.C

@@ -121,4 +121,25 @@ vectorTensorProduct(
     const Eigen::Matrix<float, Eigen::Dynamic, 1>& g,
     const Eigen::Tensor<float, 3 >& c,
     const Eigen::Matrix<float, Eigen::Dynamic, 1>& a);
+
+Eigen::VectorXd repeatElements(Eigen::VectorXd input, int n)
+{
+    if (!(n > 0))
+    {
+        Info << "The integer for the EigenFunctions::repeat method must be positive. Aborting" << endl;
+        abort();
+    }
+    else
+    {
+        Eigen::VectorXd output(input.size() * n);
+        for (int i = 0; i < input.size(); i++)
+        {
+            for (int r = 0; r < n; r++)
+            {
+                output(i * n + r) = input(i);
+            }
+        }
+        return output;
+    }
+}
 }

src/ITHACA_CORE/EigenFunctions/EigenFunctions.H

@@ -249,6 +249,17 @@ Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> vectorTensorProduct(
     const Eigen::Tensor<T, 3 >& c,
     const Eigen::Matrix<T, Eigen::Dynamic, 1>& a);
 
+
+//------------------------------------------------------------------------------
+/// @brief      Repeat each element n times after themselves
+///
+/// @details    Example: input = (a, b); n = 3
+///                      output = (a, a, a, b, b, b)
+///
+/// @param      input Eigen::VectorXd
+/// @param      n int 
+Eigen::VectorXd repeatElements(Eigen::VectorXd input, int n);
+
 };
 
 template <typename T>

src/ITHACA_CORE/ITHACAutilities/ITHACAcoeffsMass.C

@@ -321,7 +321,7 @@ Eigen::VectorXd getMassMatrixFV(
         Eigen::MatrixXd snapEigen = Foam2Eigen::field2Eigen(snapshot);
         label dim = std::nearbyint(snapEigen.rows() / (snapshot.mesh().V()).size());
         Eigen::VectorXd volumes = Foam2Eigen::field2Eigen(snapshot.mesh().V());
-        Eigen::VectorXd vol3 = volumes.replicate(dim, 1);
+        Eigen::VectorXd vol3 = EigenFunctions::repeatElements(volumes, dim);
         return vol3;
     }
     else if constexpr(std::is_same<pointMesh, GeoMesh>::value)
@@ -330,7 +330,7 @@ Eigen::VectorXd getMassMatrixFV(
         label dim = std::nearbyint(snapEigen.rows() / (snapshot.mesh()
                                    ().points()).size());
         Eigen::VectorXd pointsdata = Foam2Eigen::field2Eigen(snapshot);
-        Eigen::VectorXd points3 = pointsdata.replicate(dim, 1);
+        Eigen::VectorXd points3 = EigenFunctions::repeatElements(pointsdata, dim);
         return points3;
     }
 }