grid.cpp

Go to the documentation of this file.

/**
 * @file grid.cpp
 * @brief Implementation of the uniform 1D grid and its validator.
 */
 
#include "grid.hpp"
#include <cmath>
#include <stdexcept>
#include <cassert>
 
Grid::Grid(double L_cm, double dx_cm) : L(L_cm), dx(dx_cm) {
    if (L <= 0.0 || dx <= 0.0) {
        throw std::invalid_argument("Grid: L and dx must be positive.");
    }
    // Nx = round(L/dx) + 1 to include both endpoints 0 and L
    Nx = static_cast<std::size_t>(std::llround(L / dx)) + 1;
    if (Nx < 2) {
        throw std::invalid_argument("Grid: Nx < 2 (check L and dx).");
    }
    x.resize(Nx);
    for (std::size_t i = 0; i < Nx; ++i) {
        x[i] = static_cast<double>(i) * dx;
    }
}
 
bool validate_grid(const Grid& g, double tol) {
    assert(g.Nx >= 2);
 
    // 1) Endpoints at 0 and ~L
    assert(std::fabs(g.x.front() - 0.0) <= tol);
    assert(std::fabs(g.x.back()  - g.L) <= tol);
 
    // 2) Monotonic non-decreasing
    for (std::size_t i = 1; i < g.Nx; ++i) {
        assert(g.x[i] >= g.x[i-1]);
    }
 
    // 3) Uniform spacing: x[i] - x[i-1] ≈ dx
    for (std::size_t i = 1; i < g.Nx; ++i) {
        const double step = g.x[i] - g.x[i-1];
        assert(std::fabs(step - g.dx) <= tol);
    }
 
    // 4) Consistency with construction rule: Nx == round(L/dx) + 1
    const std::size_t idxL = static_cast<std::size_t>(std::llround(g.L / g.dx));
    assert(idxL + 1 == g.Nx);
 
    return true;
}