adding weights for C++ (#267)

CMakeLists.txt

@@ -59,8 +59,8 @@ if(POLICY CMP0135)
 endif()
 
 # Paths for downloading and building libraries
-set(ARMADILLO_VERSION "14.2.2")
-set(SUPERLU_VERSION "5.3.0")
+set(ARMADILLO_VERSION "15.2.2")
+set(SUPERLU_VERSION "7.0.1")
 set(BUILD_DIR "${CMAKE_BINARY_DIR}/third_party_build")
 set(INSTALL_DIR "${CMAKE_BINARY_DIR}/third_party_install")
 

src/cpp/divergence.cpp

@@ -27,8 +27,6 @@ Divergence::Divergence(u16 k, u32 m, Real dx) : sp_mat(m + 2, m + 1) {
       at(i, i - 1) = -1.0;
       at(i, i) = 1.0;
     }
-    // Weights
-    Q = { 1.0, 1.0, 1.0, 1.0, 1.0 };
     break;
   case 4:
     // A
@@ -50,9 +48,6 @@ Divergence::Divergence(u16 k, u32 m, Real dx) : sp_mat(m + 2, m + 1) {
       at(i, i) = 9.0 / 8.0;
       at(i, i + 1) = -1.0 / 24.0;
     }
-    Q  = { 2186.0 / 1943.0 , 2125.0 / 2828.0 , 1441.0 / 1240.0 , 648.0 / 673.0
-      , 349.0 / 350.0 , 648.0 / 673.0 , 1441.0 / 1240.0 , 2125.0 / 2828.0
-      , 2186.0 / 1943.0 };
     break;
   case 6:
     // A
@@ -94,11 +89,6 @@ Divergence::Divergence(u16 k, u32 m, Real dx) : sp_mat(m + 2, m + 1) {
       at(i, i + 1) = -25.0 / 384.0;
       at(i, i + 2) = 3.0 / 640.0;
     }
-    // Weights
-    Q  = { 2383.0 / 2005.0 , 929.0 / 2002.0 , 887.0 / 531.0 , 3124.0 / 5901.0
-      , 1706.0 / 1457.0 , 457.0 / 467.0 , 1057.0 / 1061.0 , 457.0 / 467.0
-      , 1706.0 / 1457.0 , 3124.0 / 5901.0 , 887.0 / 531.0 , 929.0 / 2002.0
-      , 2383.0 / 2005.0 };
     break;
   }
 
@@ -166,6 +156,3 @@ Divergence::Divergence(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz) {
         Utils::spkron(A1, D1) + Utils::spkron(A2, D2) + Utils::spkron(A3, D3);
   }
 }
-
-// Returns weights
-vec Divergence::getQ() { return Q; }

src/cpp/divergence.h

@@ -61,15 +61,6 @@ class Divergence : public sp_mat {
    */  
   Divergence(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz);
 
-  /**
-   * @brief Returns the weights used in the Mimeitc Divergence Operators.
-   *
-   * @note for informational purposes only, can be used in constructing new operators.
-   */    
-  vec getQ();
-
-private:
-  vec Q;
 };
 
 #endif // DIVERGENCE_H

src/cpp/gradient.cpp

@@ -37,8 +37,6 @@ Gradient::Gradient(u16 k, u32 m, Real dx) : sp_mat(m + 1, m + 2) {
       at(i, i) = -1.0;
       at(i, i + 1) = 1.0;
     }
-    // Weights
-    P = { 3.0 / 8.0 , 9.0 / 8.0 , 1.0 , 9.0 / 8.0 , 3.0 / 8.0 };
     break;
   case 4:
     // A
@@ -70,10 +68,6 @@ Gradient::Gradient(u16 k, u32 m, Real dx) : sp_mat(m + 1, m + 2) {
       at(i, i + 1) = 9.0 / 8.0;
       at(i, i + 2) = -1.0 / 24.0;
     }
-    // Weights
-    P = { 1606.0 / 4535.0 , 941.0 / 766.0 , 1384.0 / 1541.0 , 1371.0 / 1346.0
-      , 701.0 / 700.0 , 1371.0 / 1346.0 , 1384.0 / 1541.0 , 941.0 / 766.0
-      , 1606.0 / 4535.0 };
     break;
   case 6:
     // A
@@ -129,12 +123,6 @@ Gradient::Gradient(u16 k, u32 m, Real dx) : sp_mat(m + 1, m + 2) {
       at(i, i + 2) = -25.0 / 384.0;
       at(i, i + 3) = 3.0 / 640.0;
     }
-    // Weights
-    P = { 420249.0 / 1331069.0 , 2590978.0 / 1863105.0 , 882762.0 / 1402249.0
-      , 1677712.0 / 1359311.0 , 239985.0 / 261097.0 , 664189.0 / 657734.0
-      , 756049.0 / 754729.0 , 664189.0 / 657734.0 , 239985.0 / 261097.0
-      , 1677712.0 / 1359311.0 , 882762.0 / 1402249.0 , 2590978.0 / 1863105.0
-      , 420249.0 / 1331069.0 };
     break;
   case 8:
     // A
@@ -222,13 +210,6 @@ Gradient::Gradient(u16 k, u32 m, Real dx) : sp_mat(m + 1, m + 2) {
       at(i, i + 3) = 49.0 / 5120.0;
       at(i, i + 4) = -5.0 / 7168.0;
     }
-    // Weights
-    P = { 267425.0 / 904736.0 , 2307435.0 / 1517812.0 , 847667.0 / 3066027.0
-      , 4050911.0 / 2301238.0 , 498943.0 / 1084999.0 , 211042.0 / 170117.0
-      , 2065895.0 / 2191686.0 , 1262499.0 / 1258052.0 , 1314891.0 / 1312727.0
-      , 1262499.0 / 1258052.0 , 2065895.0 / 2191686.0 , 211042.0 / 170117.0
-      , 498943.0 / 1084999.0 , 4050911.0 / 2301238.0 , 847667.0 / 3066027.0
-      , 2307435.0 / 1517812.0 , 267425.0 / 904736.0 };
     break;
   }
 
@@ -296,6 +277,3 @@ Gradient::Gradient(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz) {
         Utils::spkron(A1, G1) + Utils::spkron(A2, G2) + Utils::spkron(A3, G3);
   }
 }
-
-// Returns weights
-vec Gradient::getP() { return P; }

src/cpp/gradient.h

@@ -62,16 +62,6 @@ class Gradient : public sp_mat {
    */  
   Gradient(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz);
 
-
-  /**
-   * @brief Returns the weights used in the Mimeitc Gradient Operators.
-   *
-   * @note for informational purposes only, can be used in constructing new operators.
-   */  
-  vec getP();
-
-private:
-  vec P;
 };
 
 #endif // GRADIENT_H

src/cpp/mole.h

@@ -24,5 +24,7 @@
 #include "operators.h"
 #include "robinbc.h"
 #include "utils.h"
+#include "weightsQ.h"
+#include "weightsP.h"
 
 #endif // MOLE_H

src/cpp/utils.cpp

@@ -18,8 +18,9 @@
 #include "utils.h"
 #include <cassert>
 
-#ifdef EIGEN
+//#ifdef EIGEN
 #include <eigen3/Eigen/SparseLU>
+#include <eigen3/Eigen/SparseQR>
 
 vec Utils::spsolve_eigen(const sp_mat &A, const vec &b) {
   Eigen::SparseMatrix<Real> eigen_A(A.n_rows, A.n_cols);
@@ -47,7 +48,34 @@ vec Utils::spsolve_eigen(const sp_mat &A, const vec &b) {
 
   return vec(eigen_x.data(), eigen_x.size());
 }
-#endif
+
+vec Utils::spsolve_eigenQR(const sp_mat &A, const vec &b) {
+  Eigen::SparseMatrix<Real> eigen_A(A.n_rows, A.n_cols);
+  std::vector<Eigen::Triplet<Real>> triplets;
+  Eigen::SparseQR<Eigen::SparseMatrix<Real>, Eigen::COLAMDOrdering<int>> solver;
+
+  Eigen::VectorXd eigen_x(A.n_rows);
+  triplets.reserve(5 * A.n_rows);
+
+  auto it = A.begin();
+  while (it != A.end()) {
+    triplets.push_back(Eigen::Triplet<Real>(it.row(), it.col(), *it));
+    ++it;
+  }
+
+  eigen_A.setFromTriplets(triplets.begin(), triplets.end());
+  triplets.clear();
+
+  auto b_ = conv_to<std::vector<Real>>::from(b);
+  Eigen::Map<Eigen::VectorXd> eigen_b(b_.data(), b_.size());
+
+  solver.analyzePattern(eigen_A);
+  solver.factorize(eigen_A);
+  eigen_x = solver.solve(eigen_b);
+
+  return vec(eigen_x.data(), eigen_x.size());
+}
+//#endif
 
 // Basic implementation of Kronecker product
 /*

src/cpp/utils.h

@@ -68,6 +68,8 @@ class Utils {
   */
   static vec spsolve_eigen(const sp_mat &A, const vec &b);
 
+  static vec spsolve_eigenQR(const sp_mat &A, const vec &b);
+
   /**
   * @brief An analog to the MATLAB/Octave 2D meshgrid operation
   *

src/cpp/weightsP.cpp

@@ -0,0 +1,32 @@
+/*
+* SPDX-License-Identifier: GPL-3.0-or-later
+* © 2008-2024 San Diego State University Research Foundation (SDSURF).
+* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. 
+*/
+
+/*
+ * @file weights.cpp
+ * 
+ * @brief Mimetic Differences Q Weights
+ * 
+ * @date 2025/10/14
+ * 
+ */
+
+#include "utils.h"
+#include "divergence.h"
+#include "gradient.h"
+#include "weightsP.h"
+
+WeightsP::WeightsP(u16 k, u32 m, Real dx)
+{
+  Gradient G(k, m, dx);
+
+  arma::vec b(m+2);
+  b.at(0) = -1.0;
+  b.at(m+1) = 1.0;
+
+  sp_mat Gtranspose = G.t();
+  *this = Utils::spsolve_eigenQR(Gtranspose,b);
+
+}

src/cpp/weightsP.h

@@ -0,0 +1,42 @@
+/*
+* SPDX-License-Identifier: GPL-3.0-or-later
+* © 2008-2024 San Diego State University Research Foundation (SDSURF).
+* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. 
+*/
+
+ /*
+ * @file weights.h
+ * 
+ * @brief Mimetic Differences Weights
+ * 
+ * @date 2025/10/14
+ * 
+ */
+
+#ifndef WEIGHTSP_H
+#define WEIGHTSP_H
+
+#include "utils.h"
+#include <cassert>
+
+/**
+ * @brief Generate P Weights
+ *
+ */
+class WeightsP : public sp_vec {
+
+public:
+  using sp_vec::operator=;
+
+  /**
+   * @brief P Weights Constructor
+   *
+   * @param k Order of accuracy
+   * @param m Number of cells
+   * @param dx Spacing between cells
+   */  
+  WeightsP(u16 k, u32 m, Real dx);
+
+};
+
+#endif // WEIGHTSP_H

src/cpp/weightsQ.cpp

@@ -0,0 +1,43 @@
+/*
+* SPDX-License-Identifier: GPL-3.0-or-later
+* © 2008-2024 San Diego State University Research Foundation (SDSURF).
+* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. 
+*/
+
+/*
+ * @file weights.cpp
+ * 
+ * @brief Mimetic Differences Q Weights
+ * 
+ * @date 2025/10/14
+ * 
+ */
+
+#include "utils.h"
+#include "divergence.h"
+#include "gradient.h"
+#include "weightsQ.h"
+
+WeightsQ::WeightsQ(u16 k, u32 m, Real dx): sp_vec(m + 1) 
+{
+  Divergence D(k, m, dx);
+  D.shed_row(0);
+  D.shed_row(m);
+
+  vec b(m+1);
+  b.at(0) = -1.0;
+  b.at(m) = 1.0;
+
+  sp_mat Dtranspose = D.t();
+  sp_vec Q_prime = Utils::spsolve_eigenQR(Dtranspose, b);
+  sp_vec Q(m+2);
+  Q.at(0) = 1.0;
+  Q.at(m+1) = 1.0;
+  for (int i = 1; i < m+1; ++i) {
+    Q.at(i) = Q_prime.at(i-1);
+  }
+
+  *this = Q;
+
+}
+

src/cpp/weightsQ.h

@@ -0,0 +1,42 @@
+/*
+* SPDX-License-Identifier: GPL-3.0-or-later
+* © 2008-2024 San Diego State University Research Foundation (SDSURF).
+* See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. 
+*/
+
+ /*
+ * @file weights.h
+ * 
+ * @brief Mimetic Differences Weights
+ * 
+ * @date 2025/10/14
+ * 
+ */
+
+#ifndef WEIGHTSQ_H
+#define WEIGHTSQ_H
+
+#include "utils.h"
+#include <cassert>
+
+/**
+ * @brief Generate Q Weights
+ *
+ */
+class WeightsQ : public sp_vec {
+
+public:
+  using sp_vec::operator=;
+
+  /**
+   * @brief Q Weights Constructor
+   *
+   * @param k Order of accuracy
+   * @param m Number of cells
+   * @param dx Spacing between cells
+   */  
+  WeightsQ(u16 k, u32 m, Real dx);
+
+};
+
+#endif // WEIGHTSQ_H