Merge pull request #275 from landinjm/pinning_linear_solves

Pinning linear solves

include/matrixFreePDE.h

@@ -399,12 +399,21 @@ class MatrixFreePDE : public Subscriptor
   setPeriodicity();
   void
   setPeriodicityConstraints(AffineConstraints<double> *, const DoFHandler<dim> *) const;
+
+  /**
+   * \brief Set constraints to pin the solution to 0 at a certain vertex. This is
+   * automatically done at the origin if no value terms are detected in your dependencies
+   * in a time_independent or implicit solve.
+   *
+   * \param constraints The constraint set.
+   * \param dof_handler The list of the degrees of freedom.
+   * \param target_point The point where the solution is constrained. This is the origin
+   * by default.
+   */
   void
-  getComponentsWithRigidBodyModes(std::vector<int> &) const;
-  void
-  setRigidBodyModeConstraints(const std::vector<int> &,
-                              AffineConstraints<double> *,
-                              const DoFHandler<dim> *) const;
+  set_rigid_body_mode_constraints(AffineConstraints<double> *constraints,
+                                  const DoFHandler<dim>     *dof_handler,
+                                  const Point<dim> target_point = Point<dim>()) const;
 
   // methods to apply initial conditions
   /*Virtual method to apply initial conditions.  This is usually expected to be

include/userInputParameters.h

@@ -6,11 +6,13 @@
 
 #include <deal.II/base/conditional_ostream.h>
 #include <deal.II/base/exceptions.h>
+#include <deal.II/base/point.h>
 #include <deal.II/lac/la_parallel_vector.h>
 #include <deal.II/lac/vector.h>
 
 #include <boost/algorithm/string.hpp>
 #include <boost/algorithm/string/predicate.hpp>
+#include <boost/unordered_map.hpp>
 #include <boost/variant.hpp>
 
 #include "RefinementCriterion.h"
@@ -243,6 +245,9 @@ class userInputParameters
   // Nonlinear solver parameters
   NonlinearSolverParameters nonlinear_solver_parameters;
 
+  // Pinning point parameters
+  boost::unordered_map<unsigned int, dealii::Point<dim>> pinned_point;
+
   // Variable inputs (I might be able to leave some/all of these in
   // variable_attributes)
   unsigned int number_of_variables;

src/inputFileReader/inputFileReader.cc

@@ -560,7 +560,9 @@ inputFileReader::declare_parameters(dealii::ParameterHandler     &parameter_hand
     "The number of checkpoints (or number of checkpoints per decade for the "
     "N_PER_DECADE type).");
 
-  // Declare the boundary condition variables
+  /*----------------------
+  |  Boundary conditions
+  -----------------------*/
   for (unsigned int i = 0; i < var_types.size(); i++)
     {
       if (var_types[i] == SCALAR)
@@ -604,6 +606,31 @@ inputFileReader::declare_parameters(dealii::ParameterHandler     &parameter_hand
         }
     }
 
+  /*----------------------
+  |  Pinning point
+  -----------------------*/
+  std::string pinning_text = "Pinning point: ";
+  for (unsigned int i = 0; i < var_types.size(); i++)
+    {
+      pinning_text.append(var_names.at(i));
+      parameter_handler.enter_subsection(pinning_text);
+      {
+        parameter_handler.declare_entry("x",
+                                        "-1.0",
+                                        dealii::Patterns::Double(),
+                                        "X-coordinate of the point");
+        parameter_handler.declare_entry("y",
+                                        "0.0",
+                                        dealii::Patterns::Double(),
+                                        "Y-coordinate of the point");
+        parameter_handler.declare_entry("z",
+                                        "0.0",
+                                        dealii::Patterns::Double(),
+                                        "Z-coordinate of the point");
+      }
+      parameter_handler.leave_subsection();
+    }
+
   // Declare the nucleation parameters
   parameter_handler.declare_entry(
     "Enable evolution before nucleation",

src/matrixfree/boundaryConditions.cc

@@ -281,100 +281,35 @@ MatrixFreePDE<dim, degree>::setPeriodicityConstraints(
 #endif
 }
 
-// Determine which (if any) components of the current field have rigid body
-// modes (i.e no Dirichlet BCs) if the equation is elliptic
 template <int dim, int degree>
 void
-MatrixFreePDE<dim, degree>::getComponentsWithRigidBodyModes(
-  std::vector<int> &rigidBodyModeComponents) const
-{
-  // Rigid body modes only matter for elliptic equations
-  if (userInputs.var_eq_type[currentFieldIndex] == IMPLICIT_TIME_DEPENDENT ||
-      userInputs.var_eq_type[currentFieldIndex] == TIME_INDEPENDENT)
-    {
-      // First, get the variable index of the current field
-      unsigned int starting_BC_list_index = 0;
-      for (unsigned int i = 0; i < currentFieldIndex; i++)
-        {
-          if (userInputs.var_type[i] == SCALAR)
-            {
-              starting_BC_list_index++;
-            }
-          else
-            {
-              starting_BC_list_index += dim;
-            }
-        }
-
-      // Get number of components of the field
-      unsigned int num_components = 1;
-      if (userInputs.var_type[currentFieldIndex] == VECTOR)
-        {
-          num_components = dim;
-        }
-
-      // Loop over each component and determine if it has a rigid body mode
-      // (i.e. no Dirichlet BCs)
-      for (unsigned int component = 0; component < num_components; component++)
-        {
-          bool rigidBodyMode = true;
-          for (unsigned int direction = 0; direction < 2 * dim; direction++)
-            {
-              if (userInputs.BC_list[starting_BC_list_index + component]
-                    .var_BC_type[direction] == DIRICHLET)
-                {
-                  rigidBodyMode = false;
-                }
-            }
-          // If the component has a rigid body mode, add it to the list
-          if (rigidBodyMode == true)
-            {
-              rigidBodyModeComponents.push_back(component);
-            }
-        }
-    }
-}
-
-// Set constraints to pin the solution if there are no Dirichlet BCs for a
-// component of a variable in an elliptic equation
-template <int dim, int degree>
-void
-MatrixFreePDE<dim, degree>::setRigidBodyModeConstraints(
-  const std::vector<int>    &rigidBodyModeComponents,
+MatrixFreePDE<dim, degree>::set_rigid_body_mode_constraints(
   AffineConstraints<double> *constraints,
-  const DoFHandler<dim>     *dof_handler) const
+  const DoFHandler<dim>     *dof_handler,
+  const Point<dim>           target_point) const
 {
-  if (rigidBodyModeComponents.size() > 0)
-    {
-      // Choose the point where the constraint will be placed. Must be the
-      // coordinates of a vertex.
-      dealii::Point<dim> target_point; // default constructor places the point at the
-                                       // origin
+  // Determine the number of components in the field. For a scalar field this is 1, for a
+  // vector dim, etc.
+  unsigned int n_components = 0;
+  userInputs.var_type[currentFieldIndex] == VECTOR ? n_components = dim
+                                                   : n_components = 1;
 
-      unsigned int vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
-
-      // Loop over each locally owned cell
-      for (const auto &cell : dof_handler->active_cell_iterators())
+  // Loop over each locally owned cell
+  for (const auto &cell : dof_handler->active_cell_iterators())
+    {
+      if (cell->is_locally_owned())
         {
-          if (cell->is_locally_owned())
+          for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
             {
-              for (unsigned int i = 0; i < vertices_per_cell; ++i)
+              // Check if the vertex is the target vertex
+              if (target_point.distance(cell->vertex(i)) < 1.0e-2 * cell->diameter())
                 {
-                  // Check if the vertex is the target vertex
-                  if (target_point.distance(cell->vertex(i)) < 1e-2 * cell->diameter())
+                  // Loop through the number of components and add the constraint
+                  for (unsigned int component = 0; component < n_components; component++)
                     {
-                      // Loop through the list of components with rigid body
-                      // modes and add an inhomogeneous constraint for each
-                      for (unsigned int component_num = 0;
-                           component_num < rigidBodyModeComponents.size();
-                           component_num++)
-                        {
-                          // unsigned int nodeID = cell->vertex_dof_index(i,
-                          // component_num);
-                          //  Temporarily disabling the addition of inhomogeneous
-                          //  constraints constraints->add_line(nodeID);
-                          //  constraints->set_inhomogeneity(nodeID,0.0);
-                        }
+                      unsigned int nodeID = cell->vertex_dof_index(i, component);
+                      constraints->add_line(nodeID);
+                      constraints->set_inhomogeneity(nodeID, 0.0);
                     }
                 }
             }

src/matrixfree/init.cc

@@ -163,11 +163,14 @@ MatrixFreePDE<dim, degree>::init()
       // Get hanging node constraints
       DoFTools::make_hanging_node_constraints(*dof_handler, *constraintsOther);
 
-      // Add a constraint to fix the value at the origin to zero if all BCs are
-      // zero-derivative or periodic
-      std::vector<int> rigidBodyModeComponents;
-      // getComponentsWithRigidBodyModes(rigidBodyModeComponents);
-      // setRigidBodyModeConstraints(rigidBodyModeComponents,constraintsOther,dof_handler);
+      // Pin solution
+      if (userInputs.pinned_point.find(currentFieldIndex) !=
+          userInputs.pinned_point.end())
+        {
+          set_rigid_body_mode_constraints(constraintsOther,
+                                          dof_handler,
+                                          userInputs.pinned_point[currentFieldIndex]);
+        }
 
       // Get constraints for periodic BCs
       setPeriodicityConstraints(constraintsOther, dof_handler);

src/matrixfree/reinit.cc

@@ -50,11 +50,14 @@ MatrixFreePDE<dim, degree>::reinit()
       // Get hanging node constraints
       DoFTools::make_hanging_node_constraints(*dof_handler, *constraintsOther);
 
-      // Add a constraint to fix the value at the origin to zero if all BCs are
-      // zero-derivative or periodic
-      std::vector<int> rigidBodyModeComponents;
-      getComponentsWithRigidBodyModes(rigidBodyModeComponents);
-      setRigidBodyModeConstraints(rigidBodyModeComponents, constraintsOther, dof_handler);
+      // Pin solution
+      if (userInputs.pinned_point.find(currentFieldIndex) !=
+          userInputs.pinned_point.end())
+        {
+          set_rigid_body_mode_constraints(constraintsOther,
+                                          dof_handler,
+                                          userInputs.pinned_point[currentFieldIndex]);
+        }
 
       // Get constraints for periodic BCs
       setPeriodicityConstraints(constraintsOther, dof_handler);

src/userInputParameters/userInputParameters.cc

@@ -630,11 +630,45 @@ userInputParameters<dim>::userInputParameters(inputFileReader          &input_fi
           bc_text.append(", y component");
           list_of_BCs.push_back(parameter_handler.get(bc_text));
 
-          bc_text = "Boundary condition for variable ";
-          bc_text.append(var_name.at(i));
-          bc_text.append(", z component");
-          list_of_BCs.push_back(parameter_handler.get(bc_text));
+          if (dim > 2)
+            {
+              bc_text = "Boundary condition for variable ";
+              bc_text.append(var_name.at(i));
+              bc_text.append(", z component");
+              list_of_BCs.push_back(parameter_handler.get(bc_text));
+            }
+        }
+    }
+
+  /*----------------------
+  |  Pinning point
+  -----------------------*/
+  std::string pinning_text = "Pinning point: ";
+  for (unsigned int i = 0; i < number_of_variables; i++)
+    {
+      pinning_text.append(input_file_reader.var_names.at(i));
+      parameter_handler.enter_subsection(pinning_text);
+
+      // Skip if the default
+      if (parameter_handler.get_double("x") == -1.0)
+        {
+          parameter_handler.leave_subsection();
+          continue;
+        }
+
+      // Otherwise, fill out point
+      if (dim == 2)
+        {
+          pinned_point[i] = dealii::Point<dim>(parameter_handler.get_double("x"),
+                                               parameter_handler.get_double("y"));
         }
+      else
+        {
+          pinned_point[i] = dealii::Point<dim>(parameter_handler.get_double("x"),
+                                               parameter_handler.get_double("y"),
+                                               parameter_handler.get_double("z"));
+        }
+      parameter_handler.leave_subsection();
     }
 
   // Load the BC information from the strings into a varBCs object

tests/unit_tests/main.cc

@@ -107,34 +107,6 @@ main(int argc, char **argv)
   pass = computeStress_tester_3DT.test_computeStress();
   tests_passed += pass;
 
-  // Unit tests for the method "setRigidBodyModeConstraints"
-  total_tests++;
-  unitTest<2, double> setRigidBodyModeConstraints_tester_null;
-  std::vector<int>    rigidBodyModeComponents;
-  pass = setRigidBodyModeConstraints_tester_null
-           .test_setRigidBodyModeConstraints(rigidBodyModeComponents, userInputs);
-  tests_passed += pass;
-
-  total_tests++;
-  unitTest<2, double> setRigidBodyModeConstraints_tester_one;
-  rigidBodyModeComponents.clear();
-  rigidBodyModeComponents.push_back(0);
-  pass = setRigidBodyModeConstraints_tester_one
-           .test_setRigidBodyModeConstraints(rigidBodyModeComponents, userInputs);
-  tests_passed += pass;
-
-  // In debug mode this test has an exception because it is trying to access
-  // components higher than 0 on the mesh for a scalar. To get this test working
-  // again, I'll need a DoFHandler for a vector field. total_tests++;
-  // unitTest<2,double> setRigidBodyModeConstraints_tester_three;
-  // rigidBodyModeComponents.clear();
-  // rigidBodyModeComponents.push_back(0);
-  // rigidBodyModeComponents.push_back(1);
-  // rigidBodyModeComponents.push_back(2);
-  // pass =
-  // setRigidBodyModeConstraints_tester_three.test_setRigidBodyModeConstraints(rigidBodyModeComponents,
-  // userInputs); tests_passed += pass;
-
   // Unit tests for the "LinearSolverParameters" class
   total_tests++;
   unitTest<2, double> LinearSolverParameters_tester;