Validate dx/dy/dz spacing arguments across all C++ operator constructors (H2)

src/cpp/addscalarbc.cpp

@@ -145,6 +145,7 @@ void applyBCPairRhs(vec &b, const vec &dc, const vec &nc,
 void addScalarBClhs(u16 k, u32 m, Real dx,
                     const vec &dc, const vec &nc,
                     sp_mat &Al, sp_mat &Ar) {
+    mole::check_spacing(dx, "dx");
     Al = sp_mat(m+2, m+2);
     Ar = sp_mat(m+2, m+2);
 
@@ -171,6 +172,8 @@ void addScalarBClhs(u16 k, u32 m, Real dx,
 void addScalarBClhs(u16 k, u32 m, Real dx, u32 n, Real dy,
                     const vec &dc, const vec &nc,
                     sp_mat &Al, sp_mat &Ar, sp_mat &Ab, sp_mat &At) {
+    mole::check_spacing(dx, "dx");
+    mole::check_spacing(dy, "dy");
     Al.set_size(0, 0); Ar.set_size(0, 0);
     Ab.set_size(0, 0); At.set_size(0, 0);
 
@@ -202,6 +205,9 @@ void addScalarBClhs(u16 k, u32 m, Real dx, u32 n, Real dy, u32 o, Real dz,
                     const vec &dc, const vec &nc,
                     sp_mat &Al, sp_mat &Ar, sp_mat &Ab,
                     sp_mat &At, sp_mat &Af, sp_mat &Ak) {
+    mole::check_spacing(dx, "dx");
+    mole::check_spacing(dy, "dy");
+    mole::check_spacing(dz, "dz");
     Al.set_size(0,0); Ar.set_size(0,0);
     Ab.set_size(0,0); At.set_size(0,0);
     Af.set_size(0,0); Ak.set_size(0,0);
@@ -284,6 +290,7 @@ void addScalarBCrhs(vec &b, const vec &dc, const vec &nc,
 // ============================================================================
 
 void addScalarBC(sp_mat &A, vec &b, u16 k, u32 m, Real dx, const BC1D &bc) {
+    mole::check_spacing(dx, "dx");
     assert(bc.dc.n_elem == 2 && "dc must be a 2x1 vector");
     assert(bc.nc.n_elem == 2 && "nc must be a 2x1 vector");
     assert(A.n_rows == A.n_cols && "A must be square");
@@ -306,6 +313,8 @@ void addScalarBC(sp_mat &A, vec &b, u16 k, u32 m, Real dx, const BC1D &bc) {
 
 void addScalarBC(sp_mat &A, vec &b, u16 k, u32 m, Real dx,
                  u32 n, Real dy, const BC2D &bc) {
+    mole::check_spacing(dx, "dx");
+    mole::check_spacing(dy, "dy");
     assert(bc.dc.n_elem == 4 && "dc must be a 4x1 vector");
     assert(bc.nc.n_elem == 4 && "nc must be a 4x1 vector");
     assert(bc.v.size() == 4  && "v must have 4 vectors");
@@ -328,6 +337,9 @@ void addScalarBC(sp_mat &A, vec &b, u16 k, u32 m, Real dx,
 
 void addScalarBC(sp_mat &A, vec &b, u16 k, u32 m, Real dx,
                  u32 n, Real dy, u32 o, Real dz, const BC3D &bc) {
+    mole::check_spacing(dx, "dx");
+    mole::check_spacing(dy, "dy");
+    mole::check_spacing(dz, "dz");
     assert(bc.dc.n_elem == 6 && "dc must be a 6x1 vector");
     assert(bc.nc.n_elem == 6 && "nc must be a 6x1 vector");
     assert(bc.v.size() == 6  && "v must have 6 vectors");

src/cpp/divergence.cpp

@@ -32,6 +32,7 @@ int Divergence::isPeriodic(const ivec &dc, const ivec &nc) {
 }
 
 sp_mat Divergence::periodicDiv1D(u16 k, u32 m, Real dx) {
+  mole::check_spacing(dx, "dx");
   assert(!(k % 2));
   assert(k > 1 && k < 9);
   assert(m >= 2 * k);
@@ -99,6 +100,7 @@ sp_mat Divergence::periodicDiv1D(u16 k, u32 m, Real dx) {
 // ============================================================================
 
 Divergence::Divergence(u16 k, u32 m, Real dx) : sp_mat(m + 2, m + 1) {
+  mole::check_spacing(dx, "dx");
   assert(!(k % 2));
   assert(k > 1 && k < 9);
   assert(k > 1 && k < 9);
@@ -290,6 +292,8 @@ void Divergence::build_divergence(sp_mat &D_m, sp_mat &I, u16 k, u32 dim,
 // ============================================================================
 
 Divergence::Divergence(u16 k, u32 m, u32 n, Real dx, Real dy) {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
   Divergence Dx(k, m, dx);
   Divergence Dy(k, n, dy);
 
@@ -313,6 +317,9 @@ Divergence::Divergence(u16 k, u32 m, u32 n, Real dx, Real dy) {
 // ============================================================================
 
 Divergence::Divergence(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz) {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
+  mole::check_spacing(dz, "dz");
   Divergence Dx(k, m, dx);
   Divergence Dy(k, n, dy);
   Divergence Dz(k, o, dz);
@@ -347,6 +354,7 @@ Divergence::Divergence(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz) {
 
 Divergence::Divergence(u16 k, u32 m, Real dx, const ivec &dc, const ivec &nc)
     : sp_mat() {
+  mole::check_spacing(dx, "dx");
   assert(dc.n_elem == 2 && nc.n_elem == 2);
 
   if (isPeriodic(dc, nc)) {
@@ -367,6 +375,8 @@ Divergence::Divergence(u16 k, u32 m, Real dx, const ivec &dc, const ivec &nc)
 Divergence::Divergence(u16 k, u32 m, u32 n, Real dx, Real dy, const ivec &dc,
                        const ivec &nc)
     : sp_mat() {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
   assert(dc.n_elem == 4 && nc.n_elem == 4);
 
   // dc/nc index convention: entries [0,1] control the x-axis (left, right),
@@ -394,6 +404,9 @@ Divergence::Divergence(u16 k, u32 m, u32 n, Real dx, Real dy, const ivec &dc,
 Divergence::Divergence(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz,
                        const ivec &dc, const ivec &nc)
     : sp_mat() {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
+  mole::check_spacing(dz, "dz");
   assert(dc.n_elem == 6 && nc.n_elem == 6);
 
   // dc/nc index convention:

src/cpp/gradient.cpp

@@ -35,6 +35,7 @@ int Gradient::isPeriodic(const ivec &dc, const ivec &nc) {
 }
 
 sp_mat Gradient::periodicGrad1D(u16 k, u32 m, Real dx) {
+  mole::check_spacing(dx, "dx");
   assert(!(k % 2));
   assert(k > 1 && k < 9);
   assert(m >= 2 * k);
@@ -97,6 +98,7 @@ sp_mat Gradient::periodicGrad1D(u16 k, u32 m, Real dx) {
 // ============================================================================
 
 Gradient::Gradient(u16 k, u32 m, Real dx) : sp_mat(m + 1, m + 2) {
+  mole::check_spacing(dx, "dx");
   assert(!(k % 2));
   assert(k > 1 && k < 9);
   assert(m >= 2 * k);
@@ -325,6 +327,8 @@ void Gradient::build_gradient(sp_mat &G_m, sp_mat &I, u16 k, u32 dim,
 // ============================================================================
 
 Gradient::Gradient(u16 k, u32 m, u32 n, Real dx, Real dy) {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
   Gradient Gx(k, m, dx);
   Gradient Gy(k, n, dy);
 
@@ -348,6 +352,9 @@ Gradient::Gradient(u16 k, u32 m, u32 n, Real dx, Real dy) {
 // ============================================================================
 
 Gradient::Gradient(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz) {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
+  mole::check_spacing(dz, "dz");
   Gradient Gx(k, m, dx);
   Gradient Gy(k, n, dy);
   Gradient Gz(k, o, dz);
@@ -376,6 +383,7 @@ Gradient::Gradient(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz) {
 
 Gradient::Gradient(u16 k, u32 m, Real dx, const ivec &dc, const ivec &nc)
     : sp_mat() {
+  mole::check_spacing(dx, "dx");
   assert(dc.n_elem == 2 && nc.n_elem == 2);
 
   if (isPeriodic(dc, nc)) {
@@ -396,6 +404,8 @@ Gradient::Gradient(u16 k, u32 m, Real dx, const ivec &dc, const ivec &nc)
 Gradient::Gradient(u16 k, u32 m, u32 n, Real dx, Real dy, const ivec &dc,
                    const ivec &nc)
     : sp_mat() {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
   assert(dc.n_elem == 4 && nc.n_elem == 4);
   // dc/nc index convention: entries [0,1] control the x-axis (left, right),
   // entries [2,3] control the y-axis (bottom, top).
@@ -422,6 +432,9 @@ Gradient::Gradient(u16 k, u32 m, u32 n, Real dx, Real dy, const ivec &dc,
 Gradient::Gradient(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz,
                    const ivec &dc, const ivec &nc)
     : sp_mat() {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
+  mole::check_spacing(dz, "dz");
   assert((dc.n_elem == 6) && (nc.n_elem == 6));
   // dc/nc index convention:
   //   [0,1] = left, right  (x-axis)

src/cpp/laplacian.cpp

@@ -19,6 +19,7 @@
 
 // 1-D Constructor
 Laplacian::Laplacian(u16 k, u32 m, Real dx) {
+  mole::check_spacing(dx, "dx");
   Divergence div(k, m, dx);
   Gradient grad(k, m, dx);
 
@@ -28,6 +29,8 @@ Laplacian::Laplacian(u16 k, u32 m, Real dx) {
 
 // 2-D Constructor
 Laplacian::Laplacian(u16 k, u32 m, u32 n, Real dx, Real dy) {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
   Divergence div(k, m, n, dx, dy);
   Gradient grad(k, m, n, dx, dy);
 
@@ -37,6 +40,9 @@ Laplacian::Laplacian(u16 k, u32 m, u32 n, Real dx, Real dy) {
 
 // 3-D Constructor
 Laplacian::Laplacian(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz) {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
+  mole::check_spacing(dz, "dz");
   Divergence div(k, m, n, o, dx, dy, dz);
   Gradient grad(k, m, n, o, dx, dy, dz);
 

src/cpp/mixedbc.cpp

@@ -18,6 +18,7 @@
 MixedBC::MixedBC(u16 k, u32 m, Real dx, const std::string &left,
                  const std::vector<Real> &coeffs_left, const std::string &right,
                  const std::vector<Real> &coeffs_right) {
+  mole::check_spacing(dx, "dx");
   sp_mat A(m + 2, m + 2);
   sp_mat BG(m + 2, m + 2);
 
@@ -65,6 +66,8 @@ MixedBC::MixedBC(u16 k, u32 m, Real dx, u32 n, Real dy, const std::string &left,
                  const std::string &bottom,
                  const std::vector<Real> &coeffs_bottom, const std::string &top,
                  const std::vector<Real> &coeffs_top) {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
   MixedBC Bm(k, m, dx, left, coeffs_left, right, coeffs_right);
   MixedBC Bn(k, n, dy, bottom, coeffs_bottom, top, coeffs_top);
 
@@ -90,6 +93,9 @@ MixedBC::MixedBC(u16 k, u32 m, Real dx, u32 n, Real dy, u32 o, Real dz,
                  const std::vector<Real> &coeffs_top, const std::string &front,
                  const std::vector<Real> &coeffs_front, const std::string &back,
                  const std::vector<Real> &coeffs_back) {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
+  mole::check_spacing(dz, "dz");
   MixedBC Bm(k, m, dx, left, coeffs_left, right, coeffs_right);
   MixedBC Bn(k, n, dy, bottom, coeffs_bottom, top, coeffs_top);
   MixedBC Bo(k, o, dz, front, coeffs_front, back, coeffs_back);

src/cpp/robinbc.cpp

@@ -14,6 +14,7 @@
 #include "robinbc.h"
 
 RobinBC::RobinBC(u16 k, u32 m, Real dx, Real a, Real b) {
+  mole::check_spacing(dx, "dx");
   sp_mat A(m + 2, m + 2);
   sp_mat BG(m + 2, m + 2);
 
@@ -30,6 +31,8 @@ RobinBC::RobinBC(u16 k, u32 m, Real dx, Real a, Real b) {
 
 
 RobinBC::RobinBC(u16 k, u32 m, Real dx, u32 n, Real dy, Real a, Real b) {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
   RobinBC Bm(k, m, dx, a, b);
   RobinBC Bn(k, n, dy, a, b);
 
@@ -48,6 +51,9 @@ RobinBC::RobinBC(u16 k, u32 m, Real dx, u32 n, Real dy, Real a, Real b) {
 
 RobinBC::RobinBC(u16 k, u32 m, Real dx, u32 n, Real dy, u32 o, Real dz, Real a,
                  Real b) {
+  mole::check_spacing(dx, "dx");
+  mole::check_spacing(dy, "dy");
+  mole::check_spacing(dz, "dz");
   RobinBC Bm(k, m, dx, a, b);
   RobinBC Bn(k, n, dy, a, b);
   RobinBC Bo(k, o, dz, a, b);

src/cpp/utils.cpp

@@ -17,6 +17,9 @@
 
 #include "utils.h"
 #include <cassert>
+#include <cmath>
+#include <stdexcept>
+#include <string>
 
 #ifdef EIGEN
 #include <eigen3/Eigen/SparseLU>
@@ -244,3 +247,12 @@ double Utils::trapz(const vec &x, const vec &y) {
   return 0.5 * sum;
 }
 
+// Spacing validation shared across every operator entry point.
+void mole::check_spacing(Real h, const char* name) {
+  if (!std::isfinite(h) || h <= 0.0) {
+    throw std::invalid_argument(
+        std::string("MOLE: ") + name +
+        " must be a positive finite number, got " + std::to_string(h));
+  }
+}
+

src/cpp/utils.h

@@ -116,4 +116,23 @@ class Utils {
   static double trapz(const vec &x, const vec &y);
 };
 
+namespace mole {
+
+/**
+ * @brief Validate a cell-spacing argument.
+ *
+ * Throws std::invalid_argument if @p h is zero, negative, NaN, or Inf.
+ * Called at the top of every operator constructor and AddScalarBC free
+ * function that takes a dx / dy / dz argument, so validation is active
+ * in both Debug and Release builds (unlike assert(), which vanishes
+ * under NDEBUG).
+ *
+ * @param h    Spacing value to check.
+ * @param name Parameter name for the error message (e.g. "dx").
+ * @throws std::invalid_argument if h is not a positive finite number.
+ */
+void check_spacing(Real h, const char* name);
+
+} // namespace mole
+
 #endif // UTILS_H

tests/cpp/CMakeLists.txt

@@ -20,6 +20,7 @@ set(TEST_SOURCES
   test4.cpp
   test5.cpp
   test_addscalarbc.cpp
+  test_spacing_validation.cpp
 )
 
 set(TEST_EXECUTABLES "")