removed position constraint of fixed gauge cell

nufi/grids.h

@@ -32,37 +32,17 @@ template <int dim> struct GridStructure {
   std::unique_ptr<DoFHandler<dim>> dof_handler;
   std::unique_ptr<MappingQ<dim>> mapping;
   std::unique_ptr<FE_Q<dim>> fe;
-  std::unique_ptr<AffineConstraints<double>> constraints;
-  std::unique_ptr<SparsityPattern> sparsity_pattern;
 
   std::unique_ptr<CellLocator<dim>> locator;
 
   unsigned int grid_version = 0;
 
-  // === // === //
-  // Evaluator  //
-  // === // === //
   std::vector<double>
   eval_vector_grad(const Vector<double> &solution,
                    const std::vector<Point<dim>> &points) const {
 
     std::vector<double> values(points.size());
 
-    // Uses point_gradient
-
-    // #pragma omp parallel for
-    //     for (unsigned int p = 0; p < points.size(); ++p) {
-    //
-    //       const Tensor<1, dim> grad_phi = VectorTools::point_gradient(
-    //           *mapping, *dof_handler, solution, points[p]);
-    //
-    //       values[p] = grad_phi[0];
-    //     }
-    //
-    //     return values;
-    //   }
-
-    // Uses cell locator
 #pragma omp parallel
     {
       std::vector<double> local_solution_buffer(fe->n_dofs_per_cell());
@@ -91,7 +71,6 @@ template <int dim> struct GridStructure {
 
 template <int dim>
 GridStructure<dim> make_grid_snapshot(PoissonProblem<dim> &poisson) {
-  bool PRINT_GAUGE_DOF_POSITION = true;
   GridStructure<dim> grid;
 
   grid.grid_version = 0;
@@ -106,61 +85,9 @@ GridStructure<dim> make_grid_snapshot(PoissonProblem<dim> &poisson) {
 
   grid.mapping = std::make_unique<MappingQ<dim>>(poisson.get_mapping());
 
-  grid.constraints = std::make_unique<AffineConstraints<double>>();
-
-  DoFTools::make_hanging_node_constraints(*grid.dof_handler, *grid.constraints);
-  DoFTools::make_periodicity_constraints(*grid.dof_handler, 0, 1, 0,
-                                         *grid.constraints);
-
-  const auto support_points =
-      DoFTools::map_dofs_to_support_points(*grid.mapping, *grid.dof_handler);
-
-  types::global_dof_index gauge_dof = numbers::invalid_dof_index;
-
-  // Search only inside the protected region
-  for (const auto &[dof, point] : support_points) {
-    if (grid.constraints->is_constrained(dof))
-      continue;
-    const double x = point[0];
-    if (x <= Parameters::X_DOMAIN_LEFT + .5) {
-      gauge_dof = dof;
-
-      if (PRINT_GAUGE_DOF_POSITION)
-        std::cout << " gauge_dof = " << gauge_dof
-                  << " gauge_point = " << point[0] << "\n";
-      break;
-    }
-  }
-
-  Assert(gauge_dof != numbers::invalid_dof_index,
-         ExcMessage("No gauge DoF found in protected gauge region."));
-
-  grid.constraints->add_line(gauge_dof);
-  grid.constraints->set_inhomogeneity(gauge_dof, 0.0);
-  grid.constraints->close();
-
   grid.locator = std::make_unique<CellLocator<dim>>();
   grid.locator->rebuild(*grid.dof_handler, *grid.triangulation);
 
-  // START: diagnostics
-  AssertThrow(
-      poisson.get_constraints().n_constraints() ==
-          grid.constraints->n_constraints(),
-      ExcMessage(
-          "PoissonProblem constraints doesn't match Snapshot constraints"));
-  // for (auto c1 = poisson.get_dof_handler().begin_active(),
-  //           c2 = grid.dof_handler->begin_active();
-  //      c1 != poisson.get_dof_handler().end(); ++c1, ++c2) {
-  //   std::vector<types::global_dof_index> d1(c1->get_fe().dofs_per_cell);
-  //   std::vector<types::global_dof_index> d2(c2->get_fe().dofs_per_cell);
-  //
-  //   c1->get_dof_indices(d1);
-  //   c2->get_dof_indices(d2);
-  //
-  //   AssertThrow(d1 == d2, ExcInternalError());
-  // }
-  // END: diagnostics
-
   return grid;
 }
 

nufi/poisson_problem.h

@@ -297,20 +297,17 @@ void PoissonProblem<dim>::setup_constraints(
 
   types::global_dof_index gauge_dof = numbers::invalid_dof_index;
 
-  // Search only inside the protected region
   for (const auto &[dof, point] : support_points) {
     if (constraints.is_constrained(dof))
       continue;
-    const double x = point[0];
+
-    if (x <= Parameters::X_DOMAIN_LEFT + .5) {
     gauge_dof = dof;
 
     if (PRINT_GAUGE_DOF_POSITION)
-        std::cout << " gauge_dof = " << gauge_dof
+      std::cout << " gauge_dof = " << gauge_dof << " gauge_point = " << point[0]
-                  << " gauge_point = " << point[0] << std::endl;
+                << std::endl;
     break;
   }
-  }
 
   Assert(gauge_dof != numbers::invalid_dof_index,
          ExcMessage("No gauge DoF found in protected gauge region."));

src/nufi_solver.cc

@@ -38,7 +38,6 @@ std::vector<double> NuFISolver::eval_ftilda(
   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;
   }
 
@@ -86,7 +85,6 @@ std::vector<double> NuFISolver::eval_ftilda(
   }
   for (size_t i = 0; i < x_size; ++i)
     results[i] = f0(X[i], U[i]);
-  // reset_x_eval(X);
   return results;
 }
 
@@ -102,7 +100,6 @@ 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;
   }
 
@@ -142,7 +139,6 @@ 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;
 }