tested sim with hint on cell index for location for constant uniform grid

nufi/fields.h

@@ -97,8 +97,7 @@ inline std::vector<double> eval(std::vector<double> &X,
     Points[i][0] = X[i];
   }
 
-  return eval_vector_grad(poisson.get_mapping(), poisson.get_dof_handler(),
+  return poisson.eval_vector_grad(solution, Points);
-                          solution, Points);
 }
 
 inline double integral_space_vector(const PoissonProblem<1> &poisson,

nufi/poisson_problem.h

@@ -22,6 +22,7 @@
 #include <deal.II/lac/vector.h>
 
 #include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
 #include <deal.II/grid/grid_tools.h>
 #include <deal.II/grid/tria.h>
 
@@ -41,7 +42,9 @@
 #include <deal.II/numerics/vector_tools_interpolate.h>
 #include <deal.II/numerics/vector_tools_point_gradient.h>
 #include <deal.II/numerics/vector_tools_point_value.h>
+#include <fstream>
 #include <functional>
+#include <iostream>
 #include <memory>
 #include <string>
 #include <utility>
@@ -72,6 +75,12 @@ public:
   std::vector<double> sample_electric_potential(double x_min, double x_max,
                                                 unsigned int Nx);
 
+  std::vector<double>
+  eval_vector_grad(const Vector<double> &solution,
+                   const std::vector<Point<dim>> &points) const;
+
+  void save_grid_to_file(const std::string &filename) const;
+
 private:
   void create_mesh();
   void setup_system();
@@ -93,6 +102,14 @@ private:
   std::function<double(const Point<dim> &)> rhs_function;
 
   MappingQ<dim> mapping;
+
+  // typename DoFHandler<dim>::active_cell_iterator
+  // find_cell(const DoFHandler<dim> &dof_handler, const Point<dim> x);
+
+  std::vector<typename DoFHandler<dim>::active_cell_iterator> active_cells;
+
+  mutable std::vector<double> local_solution_buffer;
+  mutable std::unique_ptr<FEPointEvaluation<dim, dim>> evaluator;
 };
 
 // Utilities
@@ -173,31 +190,57 @@ PoissonProblem<dim>::sample_electric_potential(double x_min, double x_max,
 }
 
 template <int dim>
-double
+std::vector<double> PoissonProblem<dim>::eval_vector_grad(
-eval_point_grad(const Mapping<dim> &mapping, const DoFHandler<dim> &dof_handler,
-                const Vector<double> &solution, const Point<dim> &point) {
-
-  Tensor<1, dim> grad =
-      VectorTools::point_gradient<dim>(mapping, dof_handler, solution, point);
-  double Ex = grad[0];
-
-  return Ex;
-}
-
-template <int dim>
-std::vector<double> eval_vector_grad(const Mapping<dim> &mapping,
-                                     const DoFHandler<dim> &dof_handler,
     const Vector<double> &solution,
-                                     const std::vector<Point<dim>> &points) {
+    const std::vector<Point<dim>> &points) const {
-  size_t p_size = points.size();
 
-  std::vector<double> Ex(p_size);
+  const unsigned int n_cells = active_cells.size();
-  for (size_t i = 0; i < p_size; ++i)
+  const double xmin = Parameters::X_DOMAIN_LEFT;
-    Ex[i] = eval_point_grad(mapping, dof_handler, solution, points[i]);
+  const double h = Parameters::LX / static_cast<double>(n_cells);
 
-  return Ex;
+  std::vector<double> values(points.size());
+
+  for (unsigned int p = 0; p < points.size(); ++p) {
+    const double x = points[p][0];
+
+    const unsigned int cell_index =
+        std::min(static_cast<unsigned int>(std::floor((x - xmin) / h)),
+                 n_cells - 1); // index is at most n_cells - 1
+
+    const auto cell = active_cells[cell_index];
+
+    const double xi = (x - cell->vertex(0)[0]) / h;
+
+    Point<dim> unit_point;
+    unit_point[0] = xi;
+
+    cell->get_dof_values(solution, local_solution_buffer.begin(),
+                         local_solution_buffer.end());
+
+    evaluator->reinit(cell, ArrayView<const Point<dim>>(&unit_point, 1));
+
+    evaluator->evaluate(local_solution_buffer, EvaluationFlags::gradients);
+
+    values[p] = evaluator->get_gradient(0)[0];
   }
 
+  return values;
+}
+
+// template <int dim>
+// std::vector<double> eval_vector_grad(const Mapping<dim> &mapping,
+//                                      const DoFHandler<dim> &dof_handler,
+//                                      const Vector<double> &solution,
+//                                      const std::vector<Point<dim>> &points) {
+//   size_t p_size = points.size();
+//
+//   std::vector<double> Ex(p_size);
+//   for (size_t i = 0; i < p_size; ++i)
+//     Ex[i] = eval_point_grad(mapping, dof_handler, solution, points[i]);
+//
+//   return Ex;
+// }
+
 template <int dim>
 double eval_point_value(const Mapping<dim> &mapping,
                         const DoFHandler<dim> &dof_handler,
@@ -206,6 +249,16 @@ double eval_point_value(const Mapping<dim> &mapping,
 
   return VectorTools::point_value<dim>(mapping, dof_handler, solution, point);
 }
+
+template <int dim>
+void PoissonProblem<dim>::save_grid_to_file(const std::string &filename) const {
+  std::ofstream out(filename);
+  GridOut grid_out;
+  grid_out.write_svg(triangulation, out);
+
+  std::cout << "Grid written to " << filename << "\n";
+}
+
 // dealii Poisson
 template <int dim> void PoissonProblem<dim>::create_mesh() {
 
@@ -262,6 +315,17 @@ template <int dim> void PoissonProblem<dim>::setup_system() {
 
   solution.reinit(dof_handler.n_dofs());
   system_rhs.reinit(dof_handler.n_dofs());
+
+  // used for evaluator to avoid running it anytime there is an eval
+  active_cells.clear();
+
+  for (auto cell = dof_handler.begin_active(); cell != dof_handler.end();
+       ++cell)
+    active_cells.push_back(cell);
+
+  local_solution_buffer.resize(fe.n_dofs_per_cell());
+  evaluator = std::make_unique<FEPointEvaluation<dim, dim>>(mapping, fe,
+                                                            update_gradients);
 }
 
 // Paul

src/nufi_solver.cc

@@ -9,6 +9,7 @@
 #include <deal.II/numerics/vector_tools.h>
 
 #include <cstddef>
+#include <fstream>
 #include <iostream>
 #include <memory>
 #include <ostream>
@@ -34,6 +35,7 @@ NuFISolver::eval_ftilda(unsigned int n, std::vector<double> &X, double u,
   if (n == 0) {
     for (size_t i = 0; i < x_size; ++i)
       results[i] = f0(X[i], U[i]);
+    reset_x_eval(X);
     return results;
   }
 
@@ -65,6 +67,7 @@ NuFISolver::eval_ftilda(unsigned int n, std::vector<double> &X, double u,
   }
   for (size_t i = 0; i < x_size; ++i)
     results[i] = f0(X[i], U[i]);
+  reset_x_eval(X);
   return results;
 }
 
@@ -80,6 +83,7 @@ NuFISolver::eval_f(unsigned int n, std::vector<double> &X, double u,
   if (n == 0) {
     for (size_t i = 0; i < x_size; ++i)
       results[i] = f0(X[i], U[i]);
+    reset_x_eval(X);
     return results;
   }
 
@@ -116,6 +120,7 @@ NuFISolver::eval_f(unsigned int n, std::vector<double> &X, double u,
 
   for (size_t i = 0; i < x_size; ++i)
     results[i] = f0(X[i], U[i]);
+  reset_x_eval(X);
   return results;
 }
 
@@ -134,14 +139,14 @@ NuFISolver::eval_rho(unsigned int n, std::vector<double> &X,
   std::vector<double> tmp_int(x_size);
 
   for (unsigned int i = 0; i < Nv; ++i) {
-    tmp_int = eval_ftilda(n, X, v_min + i * dv, poisson, phi_history);
+    tmp_int = eval_ftilda(n, X, v_min + i * dv, poisson,
+                          phi_history); // used eval_ftilda once per i
     for (size_t ii = 0; ii < x_size; ++ii)
       integral[ii] += tmp_int[ii];
   }
 
   for (size_t i = 0; i < x_size; ++i)
     integral[i] = 1 - integral[i] * dv;
-  reset_x_eval(X);
   return integral;
 }
 
@@ -155,20 +160,20 @@ void NuFISolver::run() {
       reinterpret_cast<double *>(std::aligned_alloc(64, sizeof(double) * Nx)),
       std::free};
 
-  if (rho == nullptr)
-    throw std::bad_alloc{};
-
   std::vector<double> int_E_squared;
   int_E_squared.reserve(Nt);
 
   std::vector<Vector<double>> phi_history;
 
-  size_t CALC_NX = Parameters::CALC_NX;
+  std::vector<double> x_eval(Parameters::CALC_NX);
-  std::vector<double> x_eval = make_x_eval(CALC_NX);
+
+  if (rho == nullptr)
+    throw std::bad_alloc{};
+
+  std::ofstream time_file("results/simulation_time.dat");
 
   double total_time = 0;
-  std::ofstream time_file("results/simulation_time.txt");
+  time_file << "it step_time total_time" << "\n";
-  time_file << "# it step_time total_time" << "\n";
 
   const double x_min = Parameters::X_DOMAIN_LEFT;
   double dx = Parameters::CALC_DX;
@@ -183,10 +188,11 @@ void NuFISolver::run() {
               << std::endl;
 
     // compute rho
+    std::vector<double> x_eval = make_x_eval(Nx);
     std::vector<double> tmp_rho =
         eval_rho(it, x_eval, poisson, phi_history, Parameters::NV);
 
-    for (size_t i = 0; i < CALC_NX; ++i) {
+    for (size_t i = 0; i < Nx; i++) {
       AssertThrow(std::isfinite(tmp_rho[i]), ExcMessage("NaN detected in rho"));
       rho.get()[i] = tmp_rho[i];
     }
@@ -223,12 +229,12 @@ void NuFISolver::run() {
       // save_Efield(it, coeffs.get(), 128, "results/field_" +
       // std::to_string(it) + ".dat");
 
-      std::vector<double> E_x(CALC_NX);
       std::vector<double> x_eval_Ex = make_x_eval(Parameters::PLOT_NX);
-      std::vector<double> tmp_Ex = eval(x_eval_Ex, poisson, phi_history[it]);
+      tmp_rho = eval(x_eval_Ex, poisson, phi_history[it]);
-      reset_x_eval(x_eval_Ex);
+
-      for (size_t i = 0; i < CALC_NX; ++i)
+      std::vector<double> E_x(Parameters::PLOT_NX);
-        E_x[i] = -tmp_Ex[i];
+      for (size_t i = 0; i < Parameters::PLOT_NX; ++i)
+        E_x[i] = -tmp_rho[i];
       save_space_vector(E_x, "field", it);
 
       double int_val =

src/save_results.cc

@@ -38,7 +38,7 @@ void save_f(const NuFISolver &solver, unsigned int n,
     for (unsigned int i = 0; i < Nx_out; ++i) {
       file << val[i];
 
-      if (j < Nx_out - 1)
+      if (i < Nx_out - 1)
         file << " ";
     }