started infrastructure for f0 choice

nufi/fields.h

@@ -1,9 +1,9 @@
-#ifndef FIELDS_H
-#define FIELDS_H
+#ifndef NUFI_FIELDS_H_
+#define NUFI_FIELDS_H_
 
-#include "grids.h"
+#include "nufi/grids.h"
 #include "nufi/parameters.h"
-#include "nufi/poisson_problem.h"
+
 #include <cmath>
 #include <cstddef>
 #include <deal.II/base/function.h>
@@ -13,7 +13,6 @@
 using namespace dealii;
 
 inline std::vector<double> make_x_eval(size_t Nx) {
-
   std::vector<double> x_eval_E;
 
   double dx = (Parameters::X_DOMAIN_RIGHT - Parameters::X_DOMAIN_LEFT) / Nx;
@@ -36,20 +35,28 @@ inline void reset_x_eval(std::vector<double> &x_vals) {
 };
 
 inline double f0(const double x, const double v,
-                 const double eps = Parameters::EPS,
-                 const double k = Parameters::WAVE_NR) {
-  const double prefactor =
-      Parameters::F0_FACTOR * (1.0 + eps * std::cos(k * x));
-  const double gaussian = v * v * std::exp(-0.5 * v * v);
+                 const size_t f0_type = Parameters::f0_TYPE) {
+  const double eps = Parameters::EPS;
+  const double k = Parameters::WAVE_NR;
 
-  return prefactor * gaussian;
+  double result;
+
+  switch (f0_type) {
+  case 0:
+
+    const double prefactor =
+        Parameters::F0_FACTOR * (1.0 + eps * std::cos(k * x));
+    const double gaussian = v * v * std::exp(-0.5 * v * v);
+    result = prefactor * gaussian;
+  }
+
+  return result;
 }
 
 // wrapper for eval_point() { VectorTools::point_values() }
 inline std::vector<double> eval(std::vector<double> &X,
                                 const GridStructure<1> &grid,
                                 const Vector<double> &solution) noexcept {
-
   AssertThrow(grid.dof_handler->n_dofs() == solution.size(),
               ExcMessage("@ eval(...) grid's number of DoFs doesn't correspond "
                          "to solution's size"));
@@ -80,7 +87,7 @@ inline double integral_space_vector(const GridStructure<1> &grid,
   for (size_t i = 0; i < Nx; ++i)
     integral += tmp[i];
   return integral * dx;
-};
+}
 
 inline double integral_space_vector_squared(const GridStructure<1> &grid,
                                             const Vector<double> &solution,
@@ -96,7 +103,7 @@ inline double integral_space_vector_squared(const GridStructure<1> &grid,
   for (size_t i = 0; i < Nx; ++i)
     integral += tmp[i] * tmp[i];
   return integral * dx;
-};
+}
 
 inline std::vector<double>
 Point_vector_to_double_vector(const std::vector<Point<1>> &Points) {
@@ -108,4 +115,4 @@ Point_vector_to_double_vector(const std::vector<Point<1>> &Points) {
 
   return vector;
 }
-#endif
+#endif // NUFI_FIELDS_H_

nufi/nufi_solver.h

@@ -8,7 +8,7 @@
 #include <deal.II/numerics/vector_tools.h>
 #include <vector>
 
-#include "nufi/fields.h" //dont remove
+#include "nufi/fields.h" // dont remove
 #include "nufi/grids.h"
 #include "nufi/parameters.h"
 #include "nufi/poisson_problem.h"

nufi/parameters.h

@@ -2,6 +2,7 @@
 #define PARAMETERS_H
 
 #include <cmath>
+#include <cstddef>
 #include <cstdlib>
 #include <string>
 
@@ -22,6 +23,8 @@ constexpr double V_DOMAIN_RIGHT = 10.;
 constexpr unsigned int NV = 128;
 constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV;
 
+constexpr size_t f0_TYPE = 0; // for switch case
+
 // deal.ii options
 constexpr unsigned int GLOBAL_REFINEMENT = 6;
 constexpr unsigned int FE_DEGREE = 3;