there wont be any changes if you dont save your field....

src/blas.cc

@@ -0,0 +1,95 @@
+#include "nufi/blas.h"
+#include <cblas.h>
+
+namespace blas
+{
+
+double dot( const size_t n, const double *x, size_t incx,
+                            const double *y, size_t incy )
+{
+    return cblas_ddot(n,x,incx,y,incy);
+}
+
+float dot( const size_t n, const float *x, size_t incx,
+                           const float *y, size_t incy )
+{
+    return cblas_sdot(n,x,incx,y,incy);
+}
+
+void axpy( size_t n, double alpha, const double *x, size_t incx,
+                                         double *y, size_t incy )
+{
+    cblas_daxpy(n,alpha,x,incx,y,incy);
+}
+
+void axpy( size_t n, float  alpha, const float  *x, size_t incx,
+                                         float  *y, size_t incy )
+{
+    cblas_saxpy(n,alpha,x,incx,y,incy);
+}
+
+void scal( size_t n, double alpha, double *x, size_t incx )
+{
+    cblas_dscal(n,alpha,x,incx);
+}
+
+void scal( size_t n, float  alpha, float  *x, size_t incx )
+{
+    cblas_sscal(n,alpha,x,incx);
+}
+
+void copy( size_t n, const double *x, size_t incx, double *y, size_t incy )
+{
+    cblas_dcopy(n,x,incx,y,incy);
+}
+
+void copy( size_t n, const float  *x, size_t incx, float  *y, size_t incy )
+{
+    cblas_scopy(n,x,incx,y,incy);
+}
+
+void ger( const size_t M, const size_t N, const double alpha,
+          const double *X, const size_t incX, const double *Y, const size_t incY,
+          double *A, const size_t lda)
+{
+    cblas_dger( CblasColMajor, M, N, alpha, X, incX, Y, incY, A, lda );
+}
+
+void ger( const size_t M, const size_t N, const float alpha,
+          const float *X, const size_t incX, const float *Y, const size_t incY,
+          float *A, const size_t lda)
+{
+    cblas_sger( CblasColMajor, M, N, alpha, X, incX, Y, incY, A, lda );
+}
+
+void gemv( const char trans, size_t m, size_t n,
+		   double alpha, const double *a, size_t lda,
+           const double *x, size_t incx, double beta,
+           double *y, size_t incy )
+{
+    if ( trans == 'T' || trans == 'Y' )
+    {
+        cblas_dgemv( CblasColMajor, CblasTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
+    }
+    else
+    {
+        cblas_dgemv( CblasColMajor, CblasNoTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
+    }
+}
+
+void gemv( const char trans, size_t m, size_t n,
+		   float alpha, const float *a, size_t lda,
+           const float *x, size_t incx, float beta,
+           float *y, size_t incy )
+{
+    if ( trans == 'T' || trans == 'Y' )
+    {
+        cblas_sgemv( CblasColMajor, CblasTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
+    }
+    else
+    {
+        cblas_sgemv( CblasColMajor, CblasNoTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
+    }
+}
+
+}

src/nufi_solver.cc

@@ -1,5 +1,6 @@
 #include "nufi/nufi_solver.h"
 
+#include <algorithm>
 #include <cmath>
 #include <cstdlib>
 #include <deal.II/base/point.h>
@@ -10,6 +11,7 @@
 #include <memory>
 #include <ostream>
 #include <cstddef>
+#include <vector>
 
 #include "nufi/parameters.h"
 #include "nufi/save_results.h"
@@ -90,32 +92,19 @@ void NuFISolver::run()
 
     	for(size_t i = 0; i<Nx; i++, x+=dx)
     	{
-        double ith_rho = eval_rho(it, x, coeffs.get(), Parameters::NV); 
+        double ith_rho = eval_rho(it, x, coeffs.get(),Parameters::NV); 
         AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho"));
     		rho.get()[i] =  ith_rho;
     	}
 
-      // bool rho_written = false;
-      //
-      // if (!rho_written)
-      // {
-      //     std::ofstream rho_file("rho_initial.dat");
-      //
-      //     if (!rho_file)
-      //         throw std::runtime_error("Could not open rho_initial.dat");
-      //
-      //     rho_file.precision(16);
-      //     rho_file << std::scientific;
-      //
-      //     for (size_t i = 0; i < Nx; ++i)
-      //         rho_file << rho.get()[i] << "\n";
-      //
-      //     rho_written = true;
-      // }
-
       poisson.set_rhs_function(std::make_unique<ChargeDensity_NuFI<1>>(rho.get(), Parameters::SPLINE_NX));
       poisson.solve_step();
 
+      std::vector<double> E_vals = poisson.sample_electric_field(Parameters::SPLINE_NX, Parameters::X_DOMAIN_LEFT, Parameters::X_DOMAIN_RIGHT);
+
+      double* current_coeffs = coeffs.get() + it*stride_t;
+      interpolate<double, Parameters::SPLINE_ORDER>(current_coeffs, E_vals.data());
+
       if (it % Parameters::PLOT_FREQUENCY == 0)
       {
         std::cout << "Saving results... \n\n";