From 080d10851a44e1349deb06c4c5dd78b5b949bf5b Mon Sep 17 00:00:00 2001
From: "Vasco C. B. Ferreira" <vascocbf@gmail.com>
Date: Tue, 10 Mar 2026 14:58:07 +0100
Subject: [PATCH] nufi online, todo: plotting

---
 fields.hpp          |   1 +
 nufi_poisson.cc     |   3 +
 nufi_solver.hpp     | 146 +++++++++++++++++++----------------------
 parameters.hpp      |   6 +-
 poisson_problem.hpp |  64 ++++++++++++------
 spline_field.hpp    | 154 ++++++++++++++++++++++++++++++++++++++++++++
 6 files changed, 274 insertions(+), 100 deletions(-)
 create mode 100644 spline_field.hpp

diff --git a/fields.hpp b/fields.hpp
index 1e15714..fa78be3 100644
--- a/fields.hpp
+++ b/fields.hpp
@@ -33,6 +33,7 @@ inline double compute_rho(const double x,
     const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv;
     integral += f0(x, v) * dv;
   }
+
   return 1.0 - integral;
 }
 
diff --git a/nufi_poisson.cc b/nufi_poisson.cc
index 44e87bf..5e2575c 100644
--- a/nufi_poisson.cc
+++ b/nufi_poisson.cc
@@ -5,6 +5,9 @@ int main()
 {
   try
   {
+
+    //Used in the past to test the basic poisson problem
+    //
     // PoissonProblem<Parameters::DIMENSION> poisson_problem(Parameters::FE_DEGREE,
     //                                                       Parameters::NV);
     //
diff --git a/nufi_solver.hpp b/nufi_solver.hpp
index f8b7625..7318691 100644
--- a/nufi_solver.hpp
+++ b/nufi_solver.hpp
@@ -1,8 +1,3 @@
-/*
-Todo: 
-- update Nx between timesteps to account for adaptivity changes because Vector rho needs to be resized
- */
-
 #ifndef NUFI_SOLVER_HPP
 #define NUFI_SOLVER_HPP
 
@@ -12,12 +7,14 @@ Todo:
 #include <deal.II/base/tensor.h>
 #include <deal.II/numerics/fe_field_function.h>
 
+#include <iostream>
 #include <vector>
 #include <cstddef>
 
 #include "parameters.hpp"
 #include "poisson_problem.hpp"
 #include "fields.hpp" // holds f0(x,v), and compute_rho(x)
+#include "spline_field.hpp"
 
 using namespace dealii;
 
@@ -27,19 +24,16 @@ public:
   NuFISolver();
 
   void run();
-  double eval_rho(unsigned int n, double x, unsigned int Nv = Parameters::NV);
+  double eval_rho(unsigned int n, double x, const UniformSpline1D<double,4>& E_spline, unsigned int Nv = Parameters::NV);
 
 private:
 
-  double eval_ftilda(unsigned int n, double x, double u);
-  void solve_poisson(unsigned int n);
+  double eval_ftilda(unsigned int n, double x, double u, const UniformSpline1D<double, 4>& E_spline);
 
-  double evaluate_E(double x);
-  
   std::vector<double> rho;
 
   unsigned int Nt = std::floor(Parameters::TMAX/Parameters::DT);
-  unsigned int Nx;
+  unsigned int Nx = Parameters::SPLINE_NX;
 
   double Lx = Parameters::LX;
 
@@ -51,63 +45,29 @@ private:
 
 };
 
-inline double NuFISolver::evaluate_E(double x)
-{
-  // Wrap x into the periodic domain
-  double x_periodic = x - Lx * std::floor(x / Lx);
-  Point<1> p(x_periodic);
-
-  Functions::FEFieldFunction<1> E_field(
-        poisson.get_dof_handler(),
-        poisson.get_solution()
-    );
-
-  double E_val = 0.0;
-
-  try
-  {
-      // Evaluate the electric field at point p
-      // If your solution represents phi, take negative gradient
-      Tensor<1,1> grad = E_field.gradient(p);
-      E_val = -grad[0];  // -∂φ/∂x
-  }
-  catch (const VectorTools::ExcPointNotAvailableHere &)
-  {
-      // This happens if p lies in an artificial cell in parallel
-      AssertThrow(false, ExcMessage("Point not available on this process."));
-  }
-
-  return E_val;
-}
-
 inline double NuFISolver::eval_ftilda(unsigned int n,
                                       double x,
-                                      double u)
+                                      double u,
+                                      const UniformSpline1D<double, 4>& E_spline)
 {
   double Lu = std::abs(Parameters::V_DOMAIN_LEFT - Parameters::V_DOMAIN_RIGHT);
 
   if (n == 0)
       return f0(x, u);
- 
-  double Ex;
-
 
   // Initial half-step.
-  Ex = evaluate_E(x);
-  u += 0.5*dt*Ex;
-
+  u += 0.5*dt*E_spline.eval(x);
 
   while ( --n )
   {
       x -= dt*u;
-      Ex = evaluate_E(x);
-      u += dt*Ex;
+      u += dt*E_spline.eval(x);
   }
 
   // Final half-step.
   x -= dt*u;
-  Ex = evaluate_E(x);
-  u += 0.5*dt*Ex; // is this line useless ?
+  u += 0.5*dt*E_spline.eval(x);
+
   double x_periodic = x - Lx * std::floor(x / Lx);
   double u_periodic = u - Lu * std::floor(u / Lu);
 
@@ -116,6 +76,7 @@ inline double NuFISolver::eval_ftilda(unsigned int n,
 
 inline double NuFISolver::eval_rho(const unsigned int n,
                           const double x,
+                          const UniformSpline1D<double, 4>& E_spline,
                           const unsigned int Nv)
 {
   const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
@@ -124,62 +85,85 @@ inline double NuFISolver::eval_rho(const unsigned int n,
   for (unsigned int i = 0; i < Nv; ++i)
   {
     const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv;
-    integral += eval_ftilda(n, x, v) * dv;
+    AssertThrow(std::isfinite(E_spline.eval(x)), ExcMessage("NaN detected in E_spline.eval(x) inside NuFISolver::eval_rho integral loop"));
+    integral += eval_ftilda(n, x, v, E_spline) * dv;
   }
-
+  
   return 1.0 - integral;
 }
 
 class ChargeDensity_NuFI : public Function<1>
 {
 public:
-  ChargeDensity_NuFI(NuFISolver &solver, size_t n)
-      : solver(solver), n(n) {}
+  ChargeDensity_NuFI(NuFISolver &solver, size_t n, const UniformSpline1D<double,4> &E_spline)
+      : solver(solver), n(n), E_spline(E_spline) {}
 
   virtual double value(const Point<1> &p,
                        [[maybe_unused]] const unsigned int component = 0) const override
   {
       double x = p[0];
 
-      return solver.eval_rho(n, x); 
-      // u not used anymore
+      return solver.eval_rho(n, x, E_spline); 
   }
 
 private:
   NuFISolver &solver;
   size_t n;
+  const UniformSpline1D<double, 4> &E_spline;
 };
 
 
-inline void NuFISolver::solve_poisson(unsigned int n)
-{
-  ChargeDensity_NuFI rho_function(*this, n);
-
-  poisson.set_rhs_function(rho_function);
-
-  poisson.solve_step();
-}
-
 inline void NuFISolver::run()
 {
-  std::cout << "Starting NuFI solver\n";
+  std::cout << "Start of NuFISolver::run()\n";
+  
+  // init E_spline
 
-  for (unsigned int n = 0; n < Nt; ++n)
+  unsigned int Nx = Parameters::SPLINE_NX;
+  // Nx grid points
+  double dx = Lx / (Nx-1);
+
+  std::vector<double> E_grid(Nx);
+
+  //set initial E points
+  for (unsigned int i=0; i<Nx; ++i)
   {
-    std::cout << "Timestep " << n << " / " << Nt << std::endl;
+    [[maybe_unused]] double x = Parameters::X_DOMAIN_LEFT + i*dx;
+    E_grid[i] = 1; 
+  }
+
+  UniformSpline1D<double,4> E_spline(E_grid, Parameters::X_DOMAIN_LEFT, Parameters::X_DOMAIN_RIGHT);
+
+  for (unsigned int it = 0; it < Nt; ++it)
+  {
+      std::cout << "Timestep " << it << " / " << Nt << std::endl;
+
+      // Step 1: Evaluate rho^n(x) using current E_spline
+      std::cout << "Start of eval_rho loop\n";
+      std::vector<double> rho(Nx);
+      for (unsigned int i = 0; i < (Nx); ++i)
+      {
+          double x = (i + 0.5) * dx;
+          rho[i] = eval_rho(it, x, E_spline, Parameters::NV);
+      }
+      std::cout << "End of eval_rho loop\n";
+
+      ChargeDensity_NuFI rho_function(*this, it, E_spline);
+
+      poisson.set_rhs_function(rho_function);
+
+      for (unsigned int i=0; i< rho.size(); ++i) // check for bad rho[i]
+      {
+        AssertThrow(std::isfinite(rho[i]), ExcMessage("NaN detected in rho"));
+      }
+
+      poisson.solve_step();
 
 
-    double dx = Lx / Nx;
-    
-    std::cout << "Start of eval_rho step with Nx = "<< Nx<< "\n";
-    for (unsigned int i = 0; i < Nx; ++i)
-    {
-      double x = (i + 0.5) * dx;
-      rho[i] = eval_rho(n, x);
-    }
-    std::cout << "End of eval_rho step\n";
-        
-    solve_poisson(n);
+      E_grid = poisson.sample_electric_field(poisson, Nx, 0.0, Lx);
+
+      // Step 4: Build spline for E^{n+1} (used in next time step)
+      E_spline = UniformSpline1D<double, 4>(E_grid, 0.0, Lx);
   }
 
   std::cout << "NuFI simulation finished.\n";
@@ -187,7 +171,7 @@ inline void NuFISolver::run()
 
 inline NuFISolver::NuFISolver()
   : order(Parameters::FE_DEGREE),
-    poisson(order, Parameters::NV)
+    poisson(order)
 {
   std::cout << "Initializing Poisson\n";
   poisson.initialize();
diff --git a/parameters.hpp b/parameters.hpp
index 97d9fe9..672e0cf 100644
--- a/parameters.hpp
+++ b/parameters.hpp
@@ -21,8 +21,12 @@ namespace Parameters
   constexpr double EPS = 0.01;
   constexpr double WAVE_NR = 0.5;
 
+  // NUFI options
   constexpr double DT=0.05;
-  constexpr unsigned int TMAX = 10;
+  constexpr unsigned int TMAX = 2;
+
+  //spline options
+  constexpr int SPLINE_NX = 256;
 }
 
 #endif
diff --git a/poisson_problem.hpp b/poisson_problem.hpp
index 472a987..f98cac1 100644
--- a/poisson_problem.hpp
+++ b/poisson_problem.hpp
@@ -31,6 +31,9 @@
 
 #include <deal.II/numerics/data_out.h>
 #include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/fe_field_function.h>
+
+#include <vector>
 
 #include "parameters.hpp"
 #include "fields.hpp"
@@ -43,18 +46,22 @@ template <int dim>
 class PoissonProblem
 {
 public:
-  PoissonProblem(unsigned int degree, unsigned int Nv);
+  PoissonProblem(unsigned int degree);
 
   void initialize();
   void solve_step();
   void run();
 
-  void set_Nv(unsigned int new_Nv);
   void set_rhs_function(const Function<dim> &rhs);
 
   const Vector<double> &get_solution() const { return solution; }
   const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
 
+  std::vector<double> sample_electric_field(const PoissonProblem<dim> &problem,  // sampling to save as spline
+                                            unsigned int Nx,
+                                            double x_min,
+                                            double x_max);
+
 private:
   void create_mesh();
   void setup_system();
@@ -77,16 +84,9 @@ private:
   const Function<dim> *rhs_function;
 
   MappingQ<dim> mapping;
-
-  unsigned int Nv;
 };
 
-
-template <int dim>
-void PoissonProblem<dim>::set_Nv(unsigned int new_Nv)
-{
-  Nv = new_Nv;
-}
+// Utilities
 
 template <int dim>
 void PoissonProblem<dim>::set_rhs_function(const Function<dim> &rhs)
@@ -95,16 +95,48 @@ void PoissonProblem<dim>::set_rhs_function(const Function<dim> &rhs)
 }
 
 template <int dim>
-PoissonProblem<dim>::PoissonProblem(unsigned int degree, unsigned int Nv)
+PoissonProblem<dim>::PoissonProblem(unsigned int degree)
   : fe(degree)
   , dof_handler(triangulation)
   , mapping(degree)
-  , Nv(Nv)
 {}
 
+template <int dim>
+std::vector<double> PoissonProblem<dim>::sample_electric_field(
+    const PoissonProblem<dim> &problem,
+    unsigned int Nx,
+    double x_min,
+    double x_max)
+{
 
-// =-=-=-=-= Make Grid =-=-=-=-=
+    const auto &dof_handler = problem.get_dof_handler();
+    const auto &solution = problem.get_solution();
 
+    Functions::FEFieldFunction<dim, Vector<double>>
+      field_function(dof_handler, solution, mapping);
+
+    std::vector<double> values(Nx);
+
+    double Lx = x_max - x_min;
+    double dx = Lx / Nx;
+
+    for (unsigned int i = 0; i < Nx; ++i)
+    {
+        double x = x_min + i * dx;
+
+        Point<dim> p;
+        p[0] = x;
+
+        Tensor<1, dim> grad = field_function.gradient(p);
+
+        values[i] = -grad[0];  // E = -dφ/dx
+    }
+
+    return values;
+}
+
+// dealii Poisson
+ 
 template<int dim>
 void PoissonProblem<dim>::create_mesh()
 {
@@ -273,7 +305,7 @@ void PoissonProblem<dim>::output_results() const
     x_coordinate[i] = support_points[i][0];  // x-component in 1D
                                              
   //---- Output density ----
-  ChargeDensity<dim> rho(Parameters::EPS, Parameters::WAVE_NR, Nv);
+  ChargeDensity<dim> rho(Parameters::EPS, Parameters::WAVE_NR, Parameters::NV);
 
   DataOut<dim> data_out_rho;
   data_out_rho.attach_dof_handler(dof_handler);
@@ -309,8 +341,6 @@ void PoissonProblem<dim>::output_results() const
 template <int dim>
 void PoissonProblem<dim>::initialize()
 {
-  set_Nv(Parameters::NV);
-
   create_mesh();      // build grid
   setup_system();     // distribute DoFs and matrices
 }
@@ -320,7 +350,6 @@ void PoissonProblem<dim>::solve_step()
 {
   system_matrix = 0;
   system_rhs = 0;
-  std::cout << "Calling PoissonProblem::solve_step()\n";
   assemble_system();
   solve();
 }
@@ -330,7 +359,6 @@ void PoissonProblem<dim>::solve_step()
 template <int dim>
 void PoissonProblem<dim>::run()
 {
-  set_Nv(Parameters::NV); // dont use anywhere else! Other functions still use Parameters::NV.
   create_mesh();
   setup_system();
   assemble_system();
diff --git a/spline_field.hpp b/spline_field.hpp
new file mode 100644
index 0000000..fba9d05
--- /dev/null
+++ b/spline_field.hpp
@@ -0,0 +1,154 @@
+#ifndef SPLINE_FIELD_HPP
+#define SPLINE_FIELD_HPP
+
+#include <iostream>
+#include <stdexcept>
+#include <cstddef>
+#include <vector>
+#include <cmath>
+
+template<typename Real, size_t Order>
+class UniformSpline1D
+{
+public:
+
+    struct Config
+    {
+        size_t Nx;
+        Real x_min;
+        Real Lx;
+
+        Real dx;
+        Real dx_inv;
+        Real Lx_inv;
+    };
+
+    Config config;
+    std::vector<Real> coeffs;
+
+public:
+
+    UniformSpline1D(const std::vector<Real>& values,
+                    Real x_min,
+                    Real x_max)
+    {
+        config.Nx = values.size();
+        config.x_min = x_min;
+        config.Lx = x_max - x_min;
+
+        config.dx = config.Lx / (config.Nx-1);
+        config.dx_inv = 1.0 / config.dx;
+        config.Lx_inv = 1.0 / config.Lx;
+
+        coeffs.resize(config.Nx + Order - 1);
+
+        interpolate(values);
+    }
+
+    Real eval(Real x) const
+    {
+      x -= config.x_min;
+      x = std::fmod(x, config.Lx);
+      if (x<0) x+= config.Lx;
+
+      Real x_cell = x * config.dx_inv;
+      size_t i = std::min(static_cast<size_t>(std::floor(x_cell)), config.Nx - 1);
+
+      Real local_x = x_cell - i;
+
+      Real basis[Order];
+      basisFunctions(local_x, basis);
+
+      Real result = 0;
+
+      for(size_t j = 0; j < Order; ++j)
+      {
+        size_t idx = (i+j) % coeffs.size();
+        result += basis[j] * coeffs[idx];
+      }
+
+      return result;
+    }
+
+private:
+
+    static void basisFunctions(Real x, Real* N)
+    {
+        Real v[Order];
+        v[Order-1] = 1;
+
+        for(size_t k = 1; k < Order; ++k)
+        {
+            v[Order-k-1] = (1-x)*v[Order-k];
+
+            for(int i = 1-k; i < 0; ++i)
+                v[Order-1+i] = (x-i)*v[Order-1+i] + (k+1+i-x)*v[Order+i];
+
+            v[Order-1] *= x;
+        }
+
+        Real factor = 1;
+        for(size_t i=2;i<Order;i++) factor*=i;
+
+        for(size_t i=0;i<Order;i++)
+            N[i] = v[i]/factor;
+    }
+
+    void interpolate(const std::vector<Real>& values)
+    {
+
+      std::cout << "interpolating";
+        size_t N = config.Nx;
+
+        std::vector<Real> rhs(values);
+        std::vector<std::vector<Real>> A(N, std::vector<Real>(N,0));
+
+        Real Nbasis[Order];
+        basisFunctions(0, Nbasis);
+
+        for(size_t i=0;i<N;i++)
+        {
+            for(size_t j=0;j<Order;j++)
+            {
+                size_t col = (i+j)%N;
+                A[i][col] = Nbasis[j];
+            }
+        }
+
+        // naive Gaussian elimination
+        std::vector<Real> x = rhs;
+
+        for(size_t k=0;k<N;k++)
+        {
+            Real pivot = A[k][k];
+            for(size_t j=k;j<N;j++)
+                A[k][j] /= pivot;
+
+            x[k] /= pivot;
+
+            for(size_t i=k+1;i<N;i++)
+            {
+                Real f = A[i][k];
+
+                for(size_t j=k;j<N;j++)
+                    A[i][j] -= f*A[k][j];
+
+                x[i] -= f*x[k];
+            }
+        }
+
+        for(int i=N-1;i>=0;i--)
+        {
+            for(size_t j=i+1;j<N;j++)
+                x[i] -= A[i][j]*x[j];
+        }
+
+        for(size_t i=0;i<N;i++)
+            coeffs[i] = x[i];
+
+        for(size_t i=0;i<Order-1;i++)
+            coeffs[N+i] = coeffs[i];
+    }
+};
+
+#endif // !SPLINE_FIELD_HPP