nonuniform dirichlet

applications/multigrid_test/CMakeLists.txt

@@ -15,7 +15,7 @@ include_directories(${CMAKE_SOURCE_DIR}/../../src)
 include_directories(${CMAKE_SOURCE_DIR})
 
 # Set the location of the main.cc file
-set(TARGET_SRC "${CMAKE_SOURCE_DIR}/main.cc" "${CMAKE_SOURCE_DIR}/equations.cc")
+set(TARGET_SRC "${CMAKE_SOURCE_DIR}/main.cc" "${CMAKE_SOURCE_DIR}/equations.cc" "${CMAKE_SOURCE_DIR}/ICs_and_BCs.cc")
 
 # Check if there has been updates to main library
 set(dir ${PROJECT_SOURCE_DIR}/../..)

applications/multigrid_test/ICs_and_BCs.cc

@@ -0,0 +1,24 @@
+#include <cmath>
+#include <config.h>
+#include <core/boundary_conditions/nonuniform_dirichlet.h>
+
+template <int dim>
+void
+customNonuniformDirichlet<dim>::set_nonuniform_dirichlet(
+  [[maybe_unused]] const uint               &index,
+  [[maybe_unused]] const uint               &boundary_id,
+  [[maybe_unused]] const uint               &component,
+  [[maybe_unused]] const dealii::Point<dim> &point,
+  [[maybe_unused]] double                   &scalar_value,
+  [[maybe_unused]] double                   &vector_component_value) const
+{
+  if (index == 0)
+    {
+      if (boundary_id == 0 || boundary_id == 1)
+        {
+          scalar_value = std::sin(point[1] * M_PI);
+        }
+    }
+}
+
+INSTANTIATE_UNI_TEMPLATE(customNonuniformDirichlet)

applications/multigrid_test/equations.cc

@@ -11,7 +11,7 @@ customAttributeLoader::loadVariableAttributes()
   set_variable_type(0, SCALAR);
   set_variable_equation_type(0, TIME_INDEPENDENT);
 
-  set_dependencies_value_term_RHS(0, "grad(phi)");
+  set_dependencies_value_term_RHS(0, "");
   set_dependencies_gradient_term_RHS(0, "grad(phi)");
   set_dependencies_value_term_LHS(0, "");
   set_dependencies_gradient_term_LHS(0, "grad(change(phi))");
@@ -36,10 +36,6 @@ customPDE<dim, degree, number>::compute_nonexplicit_RHS(
 {
   scalarGrad phix = variable_list.get_scalar_gradient(0);
 
-  scalarValue coefficient = 1.0 / (0.05 + 2.0 * q_point_loc.square());
-
-  variable_list.set_scalar_value_term(0,
-                                      make_vectorized_array<number>(1.0) / coefficient);
   variable_list.set_scalar_gradient_term(0, -phix);
 }
 

applications/multigrid_test/parameters.prm

@@ -8,9 +8,9 @@ set Subdivisions X = 1
 set Subdivisions Y = 1
 set Subdivisions Z = 1
 set Refine factor = 5
-set Element degree = 6
+set Element degree = 1
 
-set Boundary condition for variable phi = DIRICHLET: 0.0
+set Boundary condition for variable phi = NON_UNIFORM_DIRICHLET
 
 subsection Pinning point: phi
     set value = 2147483647

include/core/boundary_conditions/constraint_handler.h

@@ -7,6 +7,7 @@
 #include <deal.II/lac/affine_constraints.h>
 #include <deal.II/numerics/vector_tools_boundary.h>
 
+#include <core/boundary_conditions/nonuniform_dirichlet.h>
 #include <core/exceptions.h>
 #include <core/type_enums.h>
 #include <core/user_input_parameters.h>
@@ -24,8 +25,8 @@ class constraintHandler
    * \brief Constructor.
    */
   constraintHandler(const boundaryParameters<dim> &_boundary_parameters,
-                    const uint                    &global_variable_index,
-                    const fieldType               &field_type);
+                    const uint                    &_global_variable_index,
+                    const fieldType               &_field_type);
 
   /**
    * \brief Getter function for constraints (constant reference).
@@ -55,11 +56,11 @@ class constraintHandler
 
   std::unique_ptr<dealii::AffineConstraints<double>> constraints;
 
-  const boundaryParameters<dim> &boundary_parameters;
+  boundaryParameters<dim> boundary_parameters;
 
-  const uint &global_variable_index;
+  uint global_variable_index;
 
-  const fieldType &field_type;
+  fieldType field_type;
 };
 
 template <int dim>
@@ -104,6 +105,7 @@ constraintHandler<dim>::make_constraints(const dealii::MappingQ1<dim>  &mapping,
         {
           if (boundary_type == boundaryType::NATURAL)
             {
+              // Do nothing because they are naturally enforced.
               continue;
             }
           else if (boundary_type == boundaryType::DIRICHLET)
@@ -119,6 +121,7 @@ constraintHandler<dim>::make_constraints(const dealii::MappingQ1<dim>  &mapping,
             }
           else if (boundary_type == boundaryType::PERIODIC)
             {
+              // Do nothing because this is enforced in triangulation_handler.
               continue;
             }
           else if (boundary_type == boundaryType::NEUMANN)
@@ -127,8 +130,13 @@ constraintHandler<dim>::make_constraints(const dealii::MappingQ1<dim>  &mapping,
             }
           else if (boundary_type == boundaryType::NON_UNIFORM_DIRICHLET)
             {
-              Assert(false,
-                     FeatureNotImplemented("Nonuniform dirichlet boundary conditions"));
+              dealii::VectorTools::interpolate_boundary_values(
+                mapping,
+                dof_handler,
+                boundary_id,
+                nonuniformDirichlet<dim, fieldType::SCALAR>(global_variable_index,
+                                                            boundary_id),
+                *constraints);
               continue;
             }
           else if (boundary_type == boundaryType::NON_UNIFORM_NEUMANN)

include/core/boundary_conditions/nonuniform_dirichlet.h

@@ -0,0 +1,115 @@
+#ifndef nonuniform_dirichlet_h
+#define nonuniform_dirichlet_h
+
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+
+#include <core/type_enums.h>
+
+/**
+ * \brief Forward declaration of user-facing implementation
+ */
+template <int dim>
+class customNonuniformDirichlet;
+
+/**
+ * \brief Function for user-implemented nonuniform dirichlet boundary condition.
+ */
+template <int dim, fieldType field_type = fieldType::SCALAR>
+class nonuniformDirichlet : public dealii::Function<dim, double>
+{
+public:
+  /**
+   * \brief Constructor.
+   */
+  nonuniformDirichlet(const uint &_index, const uint &_boundary_id)
+    : dealii::Function<dim>((field_type == fieldType::VECTOR) ? dim : 1)
+    , index(_index)
+    , boundary_id(_boundary_id) {};
+
+  /**
+   * \brief Scalar value.
+   */
+  double
+  value(const dealii::Point<dim> &p, const unsigned int component = 0) const override;
+
+  /**
+   * \brief Vector value.
+   */
+  void
+  vector_value(const dealii::Point<dim> &p, dealii::Vector<double> &value) const override;
+
+private:
+  const uint index;
+
+  const uint boundary_id;
+
+  customNonuniformDirichlet<dim> custom_nonuniform_dirichlet;
+};
+
+template <int dim, fieldType field_type>
+inline double
+nonuniformDirichlet<dim, field_type>::value(
+  const dealii::Point<dim>           &p,
+  [[maybe_unused]] const unsigned int component) const
+{
+  // Initialize passed variables to zero
+  double                 scalar_value = 0.0;
+  dealii::Vector<double> vector_value(dim);
+
+  // Pass variables to user-facing function to evaluate
+  custom_nonuniform_dirichlet
+    .set_nonuniform_dirichlet(index, boundary_id, 0, p, scalar_value, vector_value(0));
+
+  return scalar_value;
+}
+
+template <int dim, fieldType field_type>
+inline void
+nonuniformDirichlet<dim, field_type>::vector_value(const dealii::Point<dim> &p,
+                                                   dealii::Vector<double>   &value) const
+{
+  // Initialize passed variables to zero
+  double                 scalar_value = 0.0;
+  dealii::Vector<double> vector_value(dim);
+
+  // Pass variables to user-facing function to evaluate
+  for (uint i = 0; i < dim; i++)
+    {
+      custom_nonuniform_dirichlet.set_nonuniform_dirichlet(index,
+                                                           boundary_id,
+                                                           i,
+                                                           p,
+                                                           scalar_value,
+                                                           vector_value(i));
+    }
+
+  value = vector_value;
+}
+
+/**
+ * \brief User-facing implementation of nonuniform boundary conditions
+ */
+template <int dim>
+class customNonuniformDirichlet
+{
+public:
+  /**
+   * \brief Constructor.
+   */
+  customNonuniformDirichlet() = default;
+
+  /**
+   * \brief Function that passes the value/vector and point that are set in the nonuniform
+   * dirichlet.
+   */
+  void
+  set_nonuniform_dirichlet(const uint               &index,
+                           const uint               &boundary_id,
+                           const uint               &component,
+                           const dealii::Point<dim> &point,
+                           double                   &scalar_value,
+                           double                   &vector_component_value) const;
+};
+
+#endif

include/core/pde_problem.h

@@ -176,18 +176,26 @@ PDEProblem<dim, degree>::init_system()
   system_matrix.initialize_dof_vector(solution);
   system_matrix.initialize_dof_vector(residual);
 
+  // Apply boundary conditions to the solution vector as an initial guess
+  constraint_handler.get_constraints().distribute(solution);
+
   const uint nlevels = triangulation_handler.get_n_global_levels();
   mg_matrices.resize(0, nlevels - 1, user_inputs);
 
-  const std::set<dealii::types::boundary_id> dirichlet_boundary_ids = {0};
+  // This should correspond to the userinput dirichlets
+  const std::set<dealii::types::boundary_id> dirichlet_boundary_ids = {0, 1, 2, 3};
+
   mg_constrained_dofs.initialize(dof_handler);
   mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, dirichlet_boundary_ids);
 
   for (uint level = 0; level < nlevels; ++level)
     {
+      // Setup the constraint set for the multigrid level
       dealii::AffineConstraints<double> level_constraints(
         dof_handler.locally_owned_mg_dofs(level),
         dealii::DoFTools::extract_locally_relevant_level_dofs(dof_handler, level));
+
+      // Loop through boundary indices and apply constraints
       for (const dealii::types::global_dof_index dof_index :
            mg_constrained_dofs.get_boundary_indices(level))
         {
@@ -282,7 +290,7 @@ PDEProblem<dim, degree>::solve_increment()
 
   system_matrix.compute_residual(residual, solution);
 
-  dealii::SolverControl solver_control(100, 1e-12 * residual.l2_norm());
+  dealii::SolverControl solver_control(1000, 1e-12 * residual.l2_norm());
   dealii::SolverCG<dealii::LinearAlgebra::distributed::Vector<double>> cg(solver_control);
 
   Timer time;
@@ -291,7 +299,7 @@ PDEProblem<dim, degree>::solve_increment()
   try
     {
       change_in_solution = 0.0;
-      cg.solve(system_matrix, change_in_solution, residual, preconditioner);
+      cg.solve(system_matrix, change_in_solution, residual, dealii::IdentityMatrix());
       constraint_handler.get_constraints().set_zero(change_in_solution);
     }
   catch (...)

include/core/triangulation_handler.h

@@ -2,10 +2,13 @@
 #define triangulation_handler_h
 
 #include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
 #include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
 
 #include <core/exceptions.h>
 #include <core/user_input_parameters.h>
+#include <fstream>
 
 /**
  * \brief This class handlers the generation and manipulation of triangulations.
@@ -42,7 +45,20 @@ class triangulationHandler
   void
   generate_mesh(const userInputParameters<dim> &user_inputs);
 
+  /**
+   * \brief Export triangulation to vtk.
+   */
+  void
+  export_triangulation_as_vtk(const std::string &filename) const;
+
 private:
+  /**
+   * \brief Mark the domain ids on the triangulation to get the proper mapping of
+   * specified boundary conditions.
+   */
+  void
+  mark_boundaries(const userInputParameters<dim> &user_inputs) const;
+
   std::unique_ptr<Triangulation> triangulation;
 };
 
@@ -88,8 +104,61 @@ triangulationHandler<dim>::generate_mesh(const userInputParameters<dim> &user_in
     dealii::Point<dim>(),
     dealii::Point<dim>(user_inputs.spatial_discretization.domain_size));
 
+  // Mark boundaries. This is done before global refinement to reduce the number of cells
+  // we have to loop through.
+  mark_boundaries(user_inputs);
+
+  // Output triangulation to vtk if in debug mode
+#ifdef DEBUG
+  export_triangulation_as_vtk("triangulation");
+#endif
+
   // Global refinement
   triangulation->refine_global(user_inputs.spatial_discretization.refine_factor);
 }
 
+template <int dim>
+void
+triangulationHandler<dim>::export_triangulation_as_vtk(const std::string &filename) const
+{
+  dealii::GridOut grid_out;
+  std::ofstream   out(filename + ".vtk");
+  grid_out.write_vtk(*triangulation, out);
+  std::cout << "Triangulation written to " << filename << ".vtk\n";
+}
+
+template <int dim>
+void
+triangulationHandler<dim>::mark_boundaries(
+  const userInputParameters<dim> &user_inputs) const
+{
+  double tolerance = 1e-12;
+
+  // Loop through the cells
+  for (const auto &cell : triangulation->active_cell_iterators())
+    {
+      // Mark the faces (faces_per_cell = 2*dim)
+      for (uint face_number = 0; face_number < dealii::GeometryInfo<dim>::faces_per_cell;
+           ++face_number)
+        {
+          // Direction for quad and hex cells
+          uint direction = std::floor(face_number / 2);
+
+          // Mark the boundary id for x=0, y=0, z=0
+          if (std::fabs(cell->face(face_number)->center()(direction) - 0) < tolerance)
+            {
+              cell->face(face_number)->set_boundary_id(face_number);
+            }
+          // Mark the boundary id for x=max, y=max, z=max
+          else if (std::fabs(
+                     cell->face(face_number)->center()(direction) -
+                     (user_inputs.spatial_discretization.domain_size[direction])) <
+                   tolerance)
+            {
+              cell->face(face_number)->set_boundary_id(face_number);
+            }
+        }
+    }
+}
+
 #endif