From d09fbef0a16087e3cd0b52879f64ddd0346f9565 Mon Sep 17 00:00:00 2001
From: Audrey Collard-Daigneault
 <71884806+acdaigneault@users.noreply.github.com>
Date: Sat, 10 Aug 2024 20:05:57 -0400
Subject: [PATCH] Modify some classes to functions and struct in DEM (#1232)

Description
Some classes should not has been classes but functions or struct.
I modified those classes to a more coherent architechture.
I also gathered all the particle-object contact info structures in one file because it is more convenient that having 3 files for small structs.

Testing
Unit tests have been change according to that.

Former-commit-id: ce0736bbbe7a475fc29303051914693d21842378
---
 ...cle_wall_contact_info.h => contact_info.h} |  78 +++-
 include/dem/data_containers.h                 |  27 +-
 include/dem/dem.h                             |   5 -
 include/dem/dem_contact_manager.h             |  23 +-
 include/dem/find_boundary_cells_information.h |  23 +-
 include/dem/find_cell_neighbors.h             | 245 ++++++-----
 include/dem/find_contact_detection_step.h     |  17 +-
 include/dem/force_chains_visualization.h      |   2 +-
 include/dem/grid_motion.h                     |   6 +-
 include/dem/integrator.h                      |   2 +-
 include/dem/lagrangian_post_processing.h      | 156 +++----
 include/dem/output_force_torque_calculation.h |  50 +--
 include/dem/particle_particle_broad_search.h  | 237 +++++------
 include/dem/particle_particle_contact_force.h |   2 +-
 include/dem/particle_particle_contact_info.h  |  65 ---
 include/dem/particle_particle_fine_search.h   |  84 ++--
 .../dem/particle_point_line_broad_search.h    | 153 ++++---
 .../dem/particle_point_line_contact_force.h   |  14 +-
 .../particle_point_line_contact_info_struct.h |  40 --
 include/dem/particle_point_line_fine_search.h | 123 +++---
 include/dem/particle_wall_broad_search.h      | 386 ++++++++----------
 include/dem/particle_wall_contact_force.h     |   2 +-
 include/dem/particle_wall_dmt_force.h         |   1 -
 include/dem/particle_wall_fine_search.h       | 144 +++----
 include/dem/particle_wall_jkr_force.h         |   1 -
 include/dem/particle_wall_linear_force.h      |   1 -
 include/dem/particle_wall_nonlinear_force.h   |   1 -
 include/dem/update_fine_search_candidates.h   |   3 +-
 .../dem/update_local_particle_containers.h    |   2 +-
 include/fem-dem/cfd_dem_coupling.h            |   9 +-
 include/fem-dem/ib_particles_dem.h            |   8 +-
 source/dem/CMakeLists.txt                     |   4 +-
 source/dem/dem.cc                             |  26 +-
 source/dem/dem_contact_manager.cc             | 247 ++++++-----
 source/dem/find_boundary_cells_information.cc |  45 +-
 source/dem/find_cell_neighbors.cc             | 133 +++---
 source/dem/find_contact_detection_step.cc     |   8 +-
 source/dem/grid_motion.cc                     |   1 -
 source/dem/input_parameter_inspection.cc      |  21 +-
 source/dem/lagrangian_post_processing.cc      | 120 +++---
 source/dem/output_force_torque_calculation.cc |  54 +++
 source/dem/particle_particle_broad_search.cc  | 187 ++++++---
 source/dem/particle_particle_fine_search.cc   |  41 +-
 .../dem/particle_point_line_broad_search.cc   | 272 ++++++------
 .../dem/particle_point_line_contact_force.cc  |   6 +-
 source/dem/particle_point_line_fine_search.cc | 195 +++++----
 source/dem/particle_wall_broad_search.cc      | 185 ++++++++-
 source/dem/particle_wall_fine_search.cc       |  77 +++-
 source/dem/update_fine_search_candidates.cc   |   2 +-
 .../dem/update_local_particle_containers.cc   |  16 +-
 source/fem-dem/cfd_dem_coupling.cc            |   3 +-
 tests/dem/find_cell_neighbors.cc              |   7 +-
 tests/dem/find_contact_pairs.cc               |  22 +-
 tests/dem/find_full_cell_neighbors.cc         |   4 +-
 tests/dem/normal_force.cc                     |   8 +-
 .../particle_particle_contact_force_linear.cc |  27 +-
 ...rticle_particle_contact_force_nonlinear.cc |  23 +-
 ...icle_particle_contact_on_two_processors.cc |  30 +-
 tests/dem/particle_particle_fine_search.cc    |  23 +-
 tests/dem/particle_point_contact.cc           |  25 +-
 .../dem/particle_wall_contact_force_linear.cc |  13 +-
 .../particle_wall_contact_force_nonlinear.cc  |  12 +-
 tests/dem/particle_wall_contact_pairs.cc      |   7 +-
 tests/dem/particle_wall_fine_search.cc        |   9 +-
 tests/dem/post_collision_velocity.cc          |   8 +-
 ...wo_particles_multiple_contacts_parallel.cc |  30 +-
 66 files changed, 1952 insertions(+), 1849 deletions(-)
 rename include/dem/{particle_wall_contact_info.h => contact_info.h} (60%)
 delete mode 100644 include/dem/particle_particle_contact_info.h
 delete mode 100644 include/dem/particle_point_line_contact_info_struct.h

diff --git a/include/dem/particle_wall_contact_info.h b/include/dem/contact_info.h
similarity index 60%
rename from include/dem/particle_wall_contact_info.h
rename to include/dem/contact_info.h
index 8a5a8d3aaa..df7ea28072 100644
--- a/include/dem/particle_wall_contact_info.h
+++ b/include/dem/contact_info.h
@@ -14,20 +14,29 @@
  * ---------------------------------------------------------------------
  */
 
-#ifndef lethe_particle_wall_contact_info_h
-#define lethe_particle_wall_contact_info_h
+#ifndef lethe_contact_info_h
+#define lethe_contact_info_h
 
 #include <deal.II/base/point.h>
 #include <deal.II/base/tensor.h>
 
 #include <deal.II/particles/particle_iterator.h>
 
+using namespace dealii;
+
 /**
- * @brief This struct handles the information related to the calculation of the
- * particle-wall contact force
+ * @brief Handle the information related to the calculation of the
+ * particle-particle contact force. Notably it is responsible for storing
+ * information that has to be preserved over multiple iterations of a contact,
+ * namely everything related to tangential overlaps
  */
-
-using namespace dealii;
+template <int dim>
+struct particle_particle_contact_info
+{
+  Particles::ParticleIterator<dim> particle_one;
+  Particles::ParticleIterator<dim> particle_two;
+  Tensor<1, 3>                     tangential_overlap;
+};
 
 template <int dim>
 class particle_wall_contact_info
@@ -39,13 +48,16 @@ class particle_wall_contact_info
    * boundary id. This is the commonly used constructor since it houses all the
    * information required to perform the contact calculation.
    *
-   * @param particle The iterator to the particle in contact with the wall
-   * @param normal_vector The outward pointing normal vector on the wall
-   * @param point_on_boundary A point that lies on the face
-   * @param boundary_id The boundary id. This id corresponds to the number attributed to the boundary condition.
+   * @param particle The iterator to the particle in contact with the wall.
+   * @param normal_vector The outward pointing normal vector on the wall.
+   * @param point_on_boundary A point that lies on the face.
+   * @param boundary_id The boundary id. This id corresponds to the number
+   * attributed to the boundary condition.
    *
+   * TODO: This should be a struct and normal_overlap, normal_relative_velocity
+   * and tangential_relative_velocity should be removed and be calculated on the
+   * fly as done for pp forces with particle_particle_contact_info
    */
-
   particle_wall_contact_info(const Particles::ParticleIterator<dim> &particle,
                              const Tensor<1, 3>       &normal_vector,
                              const Point<3>           &point_on_boundary,
@@ -89,4 +101,48 @@ class particle_wall_contact_info
   Tensor<1, 3>                     tangential_relative_velocity;
 };
 
+/**
+ * @brief Handle information related to the calculation of the particle-line
+ * contact forces.
+ */
+template <int dim>
+struct particle_line_contact_info
+{
+  Particles::ParticleIterator<dim> particle;
+  Point<3>                         point_one;
+  Point<3>                         point_two;
+};
+
+/**
+ * @brief Handle information related to the calculation of the particle-point
+ * contact forces.
+ */
+template <int dim>
+struct particle_point_contact_info
+{
+  Particles::ParticleIterator<dim> particle;
+  Point<3>                         point;
+};
+
+/**
+ * @brief Handle information related to the cell-line matching.
+ */
+template <int dim>
+struct cell_line_info
+{
+  typename Triangulation<dim>::active_cell_iterator cell;
+  Point<3>                                          point_one;
+  Point<3>                                          point_two;
+};
+
+/**
+ * @brief Handle information related to the cell-point matching.
+ */
+template <int dim>
+struct cell_point_info
+{
+  typename Triangulation<dim>::active_cell_iterator cell;
+  Point<3>                                          point;
+};
+
 #endif
diff --git a/include/dem/data_containers.h b/include/dem/data_containers.h
index 274aa17539..c1a2f722be 100644
--- a/include/dem/data_containers.h
+++ b/include/dem/data_containers.h
@@ -18,9 +18,7 @@
 #define lethe_data_containers_h
 
 #include <dem/boundary_cells_info_struct.h>
-#include <dem/particle_particle_contact_info.h>
-#include <dem/particle_point_line_contact_info_struct.h>
-#include <dem/particle_wall_contact_info.h>
+#include <dem/contact_info.h>
 
 #include <deal.II/base/tensor.h>
 
@@ -121,21 +119,24 @@ namespace DEM
                                          Particles::ParticleIterator<dim>>
       particle_index_iterator_map;
 
-    // <particle id, point-line info>
+    // <particle id, line info>
     typedef std::unordered_map<types::particle_index,
-                               particle_point_line_contact_info_struct<dim>>
-      particle_point_line_contact_info;
+                               particle_line_contact_info<dim>>
+      particle_line_in_contact;
 
-    // <particle id, (particle iterator, point, point)>
-    typedef std::unordered_map<
-      types::particle_index,
-      std::tuple<Particles::ParticleIterator<dim>, Point<dim>, Point<dim>>>
+    // <particle id, point info>
+    typedef std::unordered_map<types::particle_index,
+                               particle_point_contact_info<dim>>
+      particle_point_in_contact;
+
+    // <particle id, line info>
+    typedef std::unordered_map<types::particle_index,
+                               particle_line_contact_info<dim>>
       particle_line_candidates;
 
     // <particle id, (particle iterator, point)>
-    typedef std::unordered_map<
-      types::particle_index,
-      std::pair<Particles::ParticleIterator<dim>, Point<dim>>>
+    typedef std::unordered_map<types::particle_index,
+                               particle_point_contact_info<dim>>
       particle_point_candidates;
 
     // <particle id, <global face id, particle iterator>>
diff --git a/include/dem/dem.h b/include/dem/dem.h
index 1973738030..b68530784f 100644
--- a/include/dem/dem.h
+++ b/include/dem/dem.h
@@ -376,11 +376,6 @@ class DEMSolver
    */
   std::shared_ptr<ParticlesForceChainsBase<dim>> particles_force_chains_object;
 
-  /**
-   * @brief The post-processing object.
-   */
-  LagrangianPostProcessing<dim> post_processing_object;
-
   /**
    * @brief The integrator object.
    */
diff --git a/include/dem/dem_contact_manager.h b/include/dem/dem_contact_manager.h
index 9091240d94..bb3f057767 100644
--- a/include/dem/dem_contact_manager.h
+++ b/include/dem/dem_contact_manager.h
@@ -279,10 +279,11 @@ class DEMContactManager
     particle_floating_wall_in_contact;
   typename dem_data_structures<dim>::particle_wall_in_contact
     particle_wall_in_contact;
-  typename dem_data_structures<dim>::particle_point_line_contact_info
-    particle_points_in_contact;
-  typename dem_data_structures<dim>::particle_point_line_contact_info
+  typename dem_data_structures<dim>::particle_line_in_contact
     particle_lines_in_contact;
+  typename dem_data_structures<dim>::particle_point_in_contact
+    particle_points_in_contact;
+
 
   // Container with all the contact information of adjacent
   // local/ghost-local for pairwise contact force calculation
@@ -299,20 +300,6 @@ class DEMContactManager
 
   // Containers with other information
   typename DEM::dem_data_structures<dim>::cell_vector periodic_cells_container;
-
-private:
-  // Broad search objects
-  ParticleParticleBroadSearch<dim>  particle_particle_broad_search_object;
-  ParticleWallBroadSearch<dim>      particle_wall_broad_search_object;
-  ParticlePointLineBroadSearch<dim> particle_point_line_broad_search_object;
-
-  // Fine search objects
-  ParticleParticleFineSearch<dim>  particle_particle_fine_search_object;
-  ParticleWallFineSearch<dim>      particle_wall_fine_search_object;
-  ParticlePointLineFineSearch<dim> particle_point_line_fine_search_object;
-
-  // Other relevant objects
-  FindCellNeighbors<dim> cell_neighbors_object;
 };
 
-#endif // lethe_dem_contact_manager_h
+#endif
diff --git a/include/dem/find_boundary_cells_information.h b/include/dem/find_boundary_cells_information.h
index c09d949677..3433865cc6 100644
--- a/include/dem/find_boundary_cells_information.h
+++ b/include/dem/find_boundary_cells_information.h
@@ -18,6 +18,7 @@
 #define lethe_find_boundary_cells_information_h
 
 #include <dem/boundary_cells_info_struct.h>
+#include <dem/contact_info.h>
 #include <dem/data_containers.h>
 #include <dem/dem_solver_parameters.h>
 
@@ -87,19 +88,13 @@ class BoundaryCellsInformation
     return boundary_cells_information;
   }
 
-  std::unordered_map<
-    std::string,
-    std::pair<typename Triangulation<dim>::active_cell_iterator, Point<dim>>> &
+  std::unordered_map<std::string, cell_point_info<dim>> &
   get_boundary_cells_with_points()
   {
     return boundary_cells_with_points;
   }
 
-  std::unordered_map<
-    std::string,
-    std::tuple<typename Triangulation<dim>::active_cell_iterator,
-               Point<dim>,
-               Point<dim>>> &
+  std::unordered_map<std::string, cell_line_info<dim>> &
   get_boundary_cells_with_lines()
   {
     return boundary_cells_with_lines;
@@ -242,11 +237,7 @@ class BoundaryCellsInformation
     global_boundary_cells_information;
 
   // Structure that contains the boundary cells which have a line
-  std::unordered_map<
-    std::string,
-    std::tuple<typename Triangulation<dim>::active_cell_iterator,
-               Point<dim>,
-               Point<dim>>>
+  std::unordered_map<std::string, cell_line_info<dim>>
     boundary_cells_with_lines;
 
   // A vector of all the local cells with boundary lines. This vector is used in
@@ -255,9 +246,7 @@ class BoundaryCellsInformation
     local_cells_with_boundary_lines;
 
   // Structure that contains the boundary cells which have a point
-  std::unordered_map<
-    std::string,
-    std::pair<typename Triangulation<dim>::active_cell_iterator, Point<dim>>>
+  std::unordered_map<std::string, cell_point_info<dim>>
     boundary_cells_with_points;
 
   // Structure that contains the boundary cells with floating walls
@@ -286,4 +275,4 @@ class BoundaryCellsInformation
   bool display_pw_contact_expansion_warning = true;
 };
 
-#endif /* find_boundary_cells_information_h */
+#endif
diff --git a/include/dem/find_cell_neighbors.h b/include/dem/find_cell_neighbors.h
index 7665439a38..225f5042ce 100644
--- a/include/dem/find_cell_neighbors.h
+++ b/include/dem/find_cell_neighbors.h
@@ -26,139 +26,126 @@
 using namespace dealii;
 
 /**
- * @brief Finds the neighbors lists of all the active cells in the input triangulation.
- * find_cell_neighbors() is written to avoid any repetition, for instance if
- * cell B is recognized as the neighbor of cell A once, cell A will not appear
- * in the neighbor list of cell B again. On the other hand,
- * find_full_cell_neighbors() function finds the neighbors list with this
- * repetition.
+ * @brief Finds the neighbor list (without repetition) of all the active
+ * cells in the triangulation. It gets the vertices of the cells to get lists
+ * of the neighbor for each cells. There is some check to prevent repetition
+ * of a cell in a list (up to 8 vertices can have the same cell in common in
+ * 3D).
+ * 2 types of container are used for cell neighbors : local-local cells
+ * and local-ghost cells.
+ *
+ * @param triangulation Triangulation to access the information of the cells
+ * @param cells_local_neighbor_list A vector (with size of the local cell
+ * number) of vectors (local adjacent cells of each local cell). First element
+ * of each set shows the main cell itself
+ * @param cells_ghost_neighbor_list A vector (with size of the local cell
+ * number) of vectors (ghost adjacent cells of each local cell). First element
+ * of each set shows the main cell itself
  */
 template <int dim>
-class FindCellNeighbors
-{
-public:
-  FindCellNeighbors();
-
-  /**
-   * @brief Finds the neighbor list (without repetition) of all the active
-   * cells in the triangulation. It gets the vertices of the cells to get lists
-   * of the neighbor for each cells. There is some check to prevent repetition
-   * of a cell in a list (up to 4 vertices can have the same cell in common in
-   * 3D). 2 types of container are used for cell neighbors : local-local cells
-   * and local-ghost cells.
-   *
-   * @param triangulation Triangulation to access the information of the cells
-   * @param cells_local_neighbor_list A vector (with size of the local cell
-   * number) of vectors (local adjacent cells of each local cell). First element
-   * of each set shows the main cell itself
-   * @param cells_ghost_neighbor_list A vector (with size of the local cell
-   * number) of vectors (ghost adjacent cells of each local cell). First element
-   * of each set shows the main cell itself
-   */
+void
+find_cell_neighbors(
+  const parallel::distributed::Triangulation<dim> &triangulation,
+  typename DEM::dem_data_structures<dim>::cells_neighbor_list
+    &cells_local_neighbor_list,
+  typename DEM::dem_data_structures<dim>::cells_neighbor_list
+    &cells_ghost_neighbor_list);
 
-  void
-  find_cell_neighbors(
-    const parallel::distributed::Triangulation<dim> &triangulation,
-    typename DEM::dem_data_structures<dim>::cells_neighbor_list
-      &cells_local_neighbor_list,
-    typename DEM::dem_data_structures<dim>::cells_neighbor_list
-      &cells_ghost_neighbor_list);
-
-  /**
-   * @brief Finds the periodic neighbor list (without repetition) of all the
-   * active cells in the triangulation. It gets the coinciding vertices of the
-   * cells at the periodic boundary 0 to get lists of the periodic neighbor
-   * cells on the periodic boundary 1 for each cells. There is some check to
-   * prevent repetition of a cell in a list (up to 4 vertices can have the
-   * same cell in common in 3D). Also, 3 types of container are used for
-   * periodic mapping of cell neighbors : local-local cells, local-ghost cells
-   * and ghost-local cells. The last container is necessary since the mapping
-   * are only from periodic boundary 0 to periodic boundary 1 and the
-   * ghost-local particle pairs need distinction for proper handling of search
-   * of particle pairs and contact forces.
-   *
-   * @param triangulation Triangulation to access the information of the cells
-   * @param periodic_boundaries_cells_information A container of information
-   * related to the pairs of cell at periodic boundaries, used to get the cells
-   * on periodic boundary 0
-   * @param cells_local_periodic_neighbor_list A vector (with size of the local
-   * cell number) of vectors (local adjacent cells of each local cell). First
-   * element of each set shows the main cell itself
-   * @param cells_ghost_periodic_neighbor_list A vector (with size of the local
-   * cell number) of vectors (ghost adjacent cells of each local cell). First
-   * element of each set shows the main cell itself
-   * @param cells_ghost_local_periodic_neighbor_list A vector (with size of the
-   * ghost cell number) of vectors (local adjacent cells of each ghost cell).
-   * First element of each set shows the main ghost cell itself
-   */
-  void
-  find_cell_periodic_neighbors(
-    const parallel::distributed::Triangulation<dim> &triangulation,
-    const typename DEM::dem_data_structures<dim>::periodic_boundaries_cells_info
-      &periodic_boundaries_cells_information,
-    typename DEM::dem_data_structures<dim>::cells_neighbor_list
-      &cells_local_periodic_neighbor_list,
-    typename DEM::dem_data_structures<dim>::cells_neighbor_list
-      &cells_ghost_periodic_neighbor_list,
-    typename DEM::dem_data_structures<dim>::cells_neighbor_list
-      &cells_ghost_local_periodic_neighbor_list);
+/**
+ * @brief Finds the periodic neighbor list (without repetition) of all the
+ * active cells in the triangulation. It gets the coinciding vertices of the
+ * cells at the periodic boundary 0 to get lists of the periodic neighbor
+ * cells on the periodic boundary 1 for each cells. There is some check to
+ * prevent repetition of a cell in a list (up to 8 vertices can have the
+ * same cell in common in 3D).
+ * 3 types of container are used for
+ * periodic mapping of cell neighbors : local-local cells, local-ghost cells
+ * and ghost-local cells. The last container is necessary since the mapping
+ * are only from periodic boundary 0 to periodic boundary 1 and the
+ * ghost-local particle pairs need distinction for proper handling of search
+ * of particle pairs and contact forces.
+ *
+ * @param triangulation Triangulation to access the information of the cells
+ * @param periodic_boundaries_cells_information A container of information
+ * related to the pairs of cell at periodic boundaries, used to get the cells
+ * on periodic boundary 0
+ * @param cells_local_periodic_neighbor_list A vector (with size of the local
+ * cell number) of vectors (local adjacent cells of each local cell). First
+ * element of each set shows the main cell itself
+ * @param cells_ghost_periodic_neighbor_list A vector (with size of the local
+ * cell number) of vectors (ghost adjacent cells of each local cell). First
+ * element of each set shows the main cell itself
+ * @param cells_ghost_local_periodic_neighbor_list A vector (with size of the
+ * ghost cell number) of vectors (local adjacent cells of each ghost cell).
+ * First element of each set shows the main ghost cell itself
+ */
+template <int dim>
+void
+find_cell_periodic_neighbors(
+  const parallel::distributed::Triangulation<dim> &triangulation,
+  const typename DEM::dem_data_structures<dim>::periodic_boundaries_cells_info
+    &periodic_boundaries_cells_information,
+  typename DEM::dem_data_structures<dim>::cells_neighbor_list
+    &cells_local_periodic_neighbor_list,
+  typename DEM::dem_data_structures<dim>::cells_neighbor_list
+    &cells_ghost_periodic_neighbor_list,
+  typename DEM::dem_data_structures<dim>::cells_neighbor_list
+    &cells_ghost_local_periodic_neighbor_list);
 
-  /**
-   * @brief Finds the full neighbor list (with repetition) of all the active
-   * cells in the triangulation. This function is used in particle-floating mesh
-   * contacts, as in this situation, we need to define all the particles located
-   * in the neighbor cells of the background cell (cut by the floating mesh) as
-   * contact candidates
-   *
-   * @param triangulation Triangulation to access the information of the cells
-   * @param cells_total_neighbor_list An unordered_map (with size of the local
-   * cell number) of vectors (all adjacent cells of each local cell)
-   */
-  void
-  find_full_cell_neighbors(
-    const parallel::distributed::Triangulation<dim> &triangulation,
-    typename DEM::dem_data_structures<dim>::cells_total_neighbor_list
-      &cells_total_neighbor_list);
+/**
+ * @brief Finds the full neighbor list (with repetition) of all the active
+ * cells in the triangulation. This function is used in particle-floating mesh
+ * contacts, as in this situation, we need to define all the particles located
+ * in the neighbor cells of the background cell (cut by the floating mesh) as
+ * contact candidates
+ *
+ * @param triangulation Triangulation to access the information of the cells
+ * @param cells_total_neighbor_list An unordered_map (with size of the local
+ * cell number) of vectors (all adjacent cells of each local cell)
+ */
+template <int dim>
+void
+find_full_cell_neighbors(
+  const parallel::distributed::Triangulation<dim> &triangulation,
+  typename DEM::dem_data_structures<dim>::cells_total_neighbor_list
+    &cells_total_neighbor_list);
 
-private:
-  /**
-   * @brief Generate a periodic neighbor cells list of the cell. With the
-   * coinciding_vertex_groups and the vertex_to_coinciding_vertex_group, we can
-   * get the coinciding vertices on a periodic boundary and with v_to_c we can
-   * get the list of the periodic neighbor cells to the main cell.
-   *
-   * Here is an example to understand how the search works :
-   * Cell 8 have vertices 0, 1, 2, 3 on the periodic boundary. First,
-   * it checks in the coinciding_vertex_groups if the vertex 0 is a key in the
-   * map, and if it is, it gets a label, let's say label is 10. In the
-   * vertex_to_coinciding_vertex_group, it can finds the the coinciding vertices
-   * with the label. Using the label 10, it gets the vertices 0 and 34, which
-   * means that vertices 0 and 34 are periodic. We do not want the cells
-   * attached to the vertex 0 since they are already found with the regular
-   * find_cell_neighbor function, so it skips vertex 0, but gets the cells 21,
-   * 22, 23, 24 attached to the vertex 34. The same checks are done for vertices
-   * 1, 2 and 3 and all periodic cells are return as a vector.
-   *
-   * @param cell The cell that needs the periodic neighbor list
-   * @param coinciding_vertex_groups A map of coinciding vertices labeled by an
-   * arbitrary element from them
-   * @param vertex_to_coinciding_vertex_group Map of a vertex to the label of a
-   * group of coinciding vertices
-   * @param v_to_c A vector of set with adjacent cells of all the vertices
-   * @param periodic_neighbor_list A vector which is the list of periodic cell
-   * neighbors
-   */
-  void
-  get_periodic_neighbor_list(
-    const typename Triangulation<dim>::active_cell_iterator &cell,
-    const std::map<unsigned int, std::vector<unsigned int>>
-      &coinciding_vertex_groups,
-    const std::map<unsigned int, unsigned int>
-      &vertex_to_coinciding_vertex_group,
-    const std::vector<
-      std::set<typename Triangulation<dim>::active_cell_iterator>> &v_to_c,
-    typename DEM::dem_data_structures<dim>::cell_vector
-      &periodic_neighbor_list);
-};
+/**
+ * @brief Generate a periodic neighbor cells list of the cell. With the
+ * coinciding_vertex_groups and the vertex_to_coinciding_vertex_group, we can
+ * get the coinciding vertices on a periodic boundary and with v_to_c we can
+ * get the list of the periodic neighbor cells to the main cell.
+ *
+ * Here is an example to understand how the search works :
+ * Cell 8 have vertices 0, 1, 2, 3 on the periodic boundary. First,
+ * it checks in the coinciding_vertex_groups if the vertex 0 is a key in the
+ * map, and if it is, it gets a label, let's say label is 10. In the
+ * vertex_to_coinciding_vertex_group, it can finds the the coinciding vertices
+ * with the label. Using the label 10, it gets the vertices 0 and 34, which
+ * means that vertices 0 and 34 are periodic. We do not want the cells
+ * attached to the vertex 0 since they are already found with the regular
+ * find_cell_neighbor function, so it skips vertex 0, but gets the cells 21,
+ * 22, 23, 24 attached to the vertex 34. The same checks are done for vertices
+ * 1, 2 and 3 and all periodic cells are return as a vector.
+ *
+ * @param cell The cell that needs the periodic neighbor list
+ * @param coinciding_vertex_groups A map of coinciding vertices labeled by an
+ * arbitrary element from them
+ * @param vertex_to_coinciding_vertex_group Map of a vertex to the label of a
+ * group of coinciding vertices
+ * @param v_to_c A vector of set with adjacent cells of all the vertices
+ * @param periodic_neighbor_list A vector which is the list of periodic cell
+ * neighbors
+ */
+template <int dim>
+void
+get_periodic_neighbor_list(
+  const typename Triangulation<dim>::active_cell_iterator &cell,
+  const std::map<unsigned int, std::vector<unsigned int>>
+                                             &coinciding_vertex_groups,
+  const std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group,
+  const std::vector<std::set<typename Triangulation<dim>::active_cell_iterator>>
+                                                      &v_to_c,
+  typename DEM::dem_data_structures<dim>::cell_vector &periodic_neighbor_list);
 
 #endif
diff --git a/include/dem/find_contact_detection_step.h b/include/dem/find_contact_detection_step.h
index 74a6236814..3a82a0cea5 100644
--- a/include/dem/find_contact_detection_step.h
+++ b/include/dem/find_contact_detection_step.h
@@ -35,18 +35,17 @@ using namespace dealii;
  * dynamic contact search steps. This value is defined as the minimum of
  * particle-particle and particle-wall displacement threshold values
  * @param mpi_communicator
- * @param sorting_in_subdomains_step True if it is insertion, load-
- * balance or contact detection step
  * @param displacement Displacement of particles since last sorting step
- * @param parallel_update Update the identification of the contact detection step in parallel.
- * If this parameter is set to false, the distance will be calculated but the
+ * @param parallel_update Update the identification of the contact detection
+ * step in parallel. If this parameter is set to false, the distance will be
+ * calculated but the
  * logical OR statement won't be called and a false value will be returned. In
  * essence, this will only update the displacement.
- * @return Returns true if the maximum cumulative
- * displacement of particles exceeds the threshold and false otherwise
  *
- */
+ * @return Returns true if the maximum cumulative displacement of particles
+ * exceeds the threshold and false otherwise
 
+ */
 template <int dim>
 void
 find_particle_contact_detection_step(
@@ -58,16 +57,14 @@ find_particle_contact_detection_step(
   const bool                       parallel_update = true);
 
 /**
- * @brief Carries out finding steps for dynamic contact search in particle-floating
+ * @brief Find steps for dynamic contact search in particle-floating
  * mesh contacts
  *
  * @param smallest_contact_search_criterion A criterion which defines the maximal displacement that a solid face may have displaced
  * @param solids All solid objects used in the simulation
  * @return Returns true if the maximum cumulative
  * displacement of particles exceeds the threshold and false otherwise
- *
  */
-
 template <int dim>
 void
 find_floating_mesh_mapping_step(
diff --git a/include/dem/force_chains_visualization.h b/include/dem/force_chains_visualization.h
index 98fa3d0242..fe88bca021 100644
--- a/include/dem/force_chains_visualization.h
+++ b/include/dem/force_chains_visualization.h
@@ -19,12 +19,12 @@
 
 #include <core/auxiliary_math_functions.h>
 #include <core/dem_properties.h>
+#include <core/pvd_handler.h>
 
 #include <dem/data_containers.h>
 #include <dem/dem_contact_manager.h>
 #include <dem/dem_solver_parameters.h>
 #include <dem/particle_particle_contact_force.h>
-#include <dem/particle_particle_contact_info.h>
 #include <dem/rolling_resistance_torque_models.h>
 
 #include <deal.II/particles/particle_handler.h>
diff --git a/include/dem/grid_motion.h b/include/dem/grid_motion.h
index f7411da884..ce8f1a30ac 100644
--- a/include/dem/grid_motion.h
+++ b/include/dem/grid_motion.h
@@ -17,9 +17,9 @@
 #ifndef lethe_grid_motion_h
 #define lethe_grid_motion_h
 
+#include <dem/contact_info.h>
 #include <dem/data_containers.h>
 #include <dem/dem_solver_parameters.h>
-#include <dem/particle_wall_contact_info.h>
 
 #include <deal.II/distributed/tria.h>
 
@@ -98,7 +98,7 @@ class GridMotion
   move_grid_rotational(Triangulation<dim, spacedim> &triangulation);
 
   /**
-   * Carries out translational motion of the triangulation
+   * @brief Carries out translational motion of the triangulation
    *
    * @param triangulation Triangulation
    */
@@ -128,4 +128,4 @@ class GridMotion
 };
 
 
-#endif /* grid_motion_h */
+#endif
diff --git a/include/dem/integrator.h b/include/dem/integrator.h
index 4d4c596bd7..b5f38121d2 100644
--- a/include/dem/integrator.h
+++ b/include/dem/integrator.h
@@ -96,4 +96,4 @@ class Integrator
             AdaptiveSparseContacts<dim> &sparse_contacts_object) = 0;
 };
 
-#endif /* integration_h */
+#endif
diff --git a/include/dem/lagrangian_post_processing.h b/include/dem/lagrangian_post_processing.h
index 3f97c5957b..94808a7aaf 100644
--- a/include/dem/lagrangian_post_processing.h
+++ b/include/dem/lagrangian_post_processing.h
@@ -41,87 +41,87 @@ using namespace dealii;
  * @todo Add tracer particles
  */
 
+/**
+ * @brief Carries out writing the cell data of the domain.
+ *
+ * @param triangulation Triangulation of the domain.
+ * @param grid_pvdhandler PVD handler for grid.
+ * @param background_dh Background DoF handler.
+ * @param particle_handler Particle handler.
+ * @param dem_parameters DEM parameters.
+ * @param current_time Simulation time.
+ * @param step_number DEM step number.
+ * @param mpi_communicator MPI communicator.
+ * @param sparse_contacts_object Àdaptive sparse contacts object.
+ */
 template <int dim>
-class LagrangianPostProcessing
-{
-public:
-  LagrangianPostProcessing();
-
-  /**
-   * Carries out writing the grid of the domain
-   *
-   * @param triangulation Triangulation
-   * @param grid_pvdhandler
-   * @param particle_handler
-   * @param dem_parameters
-   * @param background_dh
-   * @param current_time Simulation time
-   * @param step_number DEM step number
-   * @param mpi_communicator
-   */
-  void
-  write_post_processing_results(
-    const parallel::distributed::Triangulation<dim> &triangulation,
-    PVDHandler                                      &grid_pvdhandler,
-    const DoFHandler<dim>                           &background_dh,
-    const Particles::ParticleHandler<dim>           &particle_handler,
-    const DEMSolverParameters<dim>                  &dem_parameters,
-    const double                                     current_time,
-    const unsigned int                               step_number,
-    const MPI_Comm                                  &mpi_communicator,
-    AdaptiveSparseContacts<dim>                     &sparse_contacts_object);
-
-private:
-  /**
-   * Carries out the calculation of the average particles velocity in each local
-   * cell. These values are summed up during the post-processing steps (imposed
-   * by post-processing frequency), and finally divided by the sampling counter
-   * to calculate time-averaged particles velocity distribution.
-   *
-   * @param triangulation Triangulation
-   * @param particle_handler Particle handler
-   */
-  void
-  calculate_average_particles_velocity(
-    const parallel::distributed::Triangulation<dim> &triangulation,
-    const Particles::ParticleHandler<dim>           &particle_handler);
+void
+write_post_processing_results(
+  const parallel::distributed::Triangulation<dim> &triangulation,
+  PVDHandler                                      &grid_pvdhandler,
+  const DoFHandler<dim>                           &background_dh,
+  const Particles::ParticleHandler<dim>           &particle_handler,
+  const DEMSolverParameters<dim>                  &dem_parameters,
+  const double                                     current_time,
+  const unsigned int                               step_number,
+  const MPI_Comm                                  &mpi_communicator,
+  AdaptiveSparseContacts<dim>                     &sparse_contacts_object);
 
-  /**
-   * Carries out the calculation of the granular temperature in each local cell.
-   * These values are summed up during the post-processing steps (imposed by
-   * post-processing frequency), and finally divided by the sampling counter to
-   * calculate time-averaged granular temperature distribution.
-   *
-   * @param triangulation Triangulation
-   * @param particle_handler Particle handler
-   */
-  void
-  calculate_average_granular_temperature(
-    const parallel::distributed::Triangulation<dim> &triangulation,
-    const Particles::ParticleHandler<dim>           &particle_handler);
-
-  /**
-   * Carries out the calculation of average particles velocity for a single
-   * cell. It is defined as a separate private function since the
-   * calculate_average_granular_temperature function also needs to obtain the
-   * average velocity in the cell.
-   *
-   * @param cell A single cell in the triangulation
-   * @param particle_handler
-   * @return A tensor which stores the average particles velocity in the input cell
-   */
-  Tensor<1, dim>
-  calculate_cell_average_particles_velocity(
-    const typename parallel::distributed::Triangulation<dim>::cell_iterator
-                                          &cell,
-    const Particles::ParticleHandler<dim> &particle_handler);
+/**
+ * @brief Carries out the calculation of the average particles velocity in each local
+ * cell. These values are summed up during the post-processing steps (imposed
+ * by post-processing frequency), and finally divided by the sampling counter
+ * to calculate time-averaged particles velocity distribution.
+ *
+ * @param triangulation Triangulation.
+ * @param particle_handler Particle handler.
+ * @param velocity_average_x Average velocity in x-direction.
+ * @param velocity_average_y Average velocity in y-direction.
+ * @param velocity_average_z Average velocity in z-direction.
+ * @param velocity_average_magnitude Average velocity magnitude.
+ */
+template <int dim>
+void
+calculate_average_particles_velocity(
+  const parallel::distributed::Triangulation<dim> &triangulation,
+  const Particles::ParticleHandler<dim>           &particle_handler,
+  Vector<double>                                  &velocity_average_x,
+  Vector<double>                                  &velocity_average_y,
+  Vector<double>                                  &velocity_average_z,
+  Vector<double>                                  &velocity_average_magnitude);
 
-  Vector<double> velocity_average_x;
-  Vector<double> velocity_average_y;
-  Vector<double> velocity_average_z;
-  Vector<double> velocity_average_magnitude;
+/**
+ * @brief Carries out the calculation of the granular temperature in each local cell.
+ * These values are summed up during the post-processing steps (imposed by
+ * post-processing frequency), and finally divided by the sampling counter to
+ * calculate time-averaged granular temperature distribution.
+ *
+ * @param triangulation Triangulation.
+ * @param particle_handler Particle handler.
+ * @param granular_temperature_average Average granular temperature.
+ */
+template <int dim>
+void
+calculate_average_granular_temperature(
+  const parallel::distributed::Triangulation<dim> &triangulation,
+  const Particles::ParticleHandler<dim>           &particle_handler,
+  Vector<double> &granular_temperature_average);
 
-  Vector<double> granular_temperature_average;
-};
+/**
+ * @brief Carries out the calculation of average particles velocity for a single
+ * cell. It is defined as a separate private function since the
+ * calculate_average_granular_temperature function also needs to obtain the
+ * average velocity in the cell.
+ *
+ * @param cell A single cell in the triangulation
+ * @param particle_handler
+ *
+ * @return A tensor which stores the average particles velocity in the input cell
+ */
+template <int dim>
+Tensor<1, dim>
+calculate_cell_average_particles_velocity(
+  const typename parallel::distributed::Triangulation<dim>::cell_iterator &cell,
+  const Particles::ParticleHandler<dim> &particle_handler);
 
 #endif
diff --git a/include/dem/output_force_torque_calculation.h b/include/dem/output_force_torque_calculation.h
index 0a1e2de202..65e86400f2 100644
--- a/include/dem/output_force_torque_calculation.h
+++ b/include/dem/output_force_torque_calculation.h
@@ -37,7 +37,6 @@ write_forces_torques_output_locally(
  * Writes the results of force and torque calculations in a file, and it depends
  * on the verbosity, in the terminal
  */
-template <int dim>
 void
 write_forces_torques_output_results(
   const std::string                                filename,
@@ -45,51 +44,6 @@ write_forces_torques_output_results(
   const std::vector<unsigned int>                  boundary_index,
   const double                                     time_step,
   DEM::dem_data_structures<3>::vector_on_boundary &forces_boundary_information,
-  DEM::dem_data_structures<3>::vector_on_boundary &torques_boundary_information)
-{
-  unsigned int this_mpi_process =
-    Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
-  unsigned int n_mpi_processes =
-    Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
-  for (unsigned int i = 0; i < boundary_index.size(); i++)
-    {
-      std::string filename_boundary_description =
-        std::to_string(boundary_index[i]);
-      std::string filename_output =
-        filename + "_boundary_" + filename_boundary_description;
-      if (this_mpi_process == i || (i + 1) > n_mpi_processes)
-        {
-          TableHandler table;
-
-          for (unsigned int j = 0; j < forces_boundary_information.size();
-               j += output_frequency)
-            {
-              table.add_value("Force x", forces_boundary_information[j][i][0]);
-              table.add_value("Force y", forces_boundary_information[j][i][1]);
-              table.add_value("Force z", forces_boundary_information[j][i][2]);
-              table.add_value("Torque x",
-                              torques_boundary_information[j][i][0]);
-              table.add_value("Torque y",
-                              torques_boundary_information[j][i][1]);
-              table.add_value("Torque z",
-                              torques_boundary_information[j][i][2]);
-              table.add_value("Current time", j * time_step);
-            }
-
-          table.set_precision("Force x", 9);
-          table.set_precision("Force y", 9);
-          table.set_precision("Force z", 9);
-          table.set_precision("Torque x", 9);
-          table.set_precision("Torque y", 9);
-          table.set_precision("Torque z", 9);
-          table.set_precision("Current time", 9);
-
-          // output
-          std::ofstream out_file(filename_output);
-          table.write_text(out_file);
-          out_file.close();
-        }
-    }
-}
-
+  DEM::dem_data_structures<3>::vector_on_boundary
+    &torques_boundary_information);
 #endif
diff --git a/include/dem/particle_particle_broad_search.h b/include/dem/particle_particle_broad_search.h
index 2835e6886b..fb380311d1 100644
--- a/include/dem/particle_particle_broad_search.h
+++ b/include/dem/particle_particle_broad_search.h
@@ -18,9 +18,9 @@
 #define lethe_particle_particle_broad_search_h
 
 #include <dem/adaptive_sparse_contacts.h>
+#include <dem/contact_info.h>
 #include <dem/data_containers.h>
 #include <dem/find_boundary_cells_information.h>
-#include <dem/particle_particle_contact_info.h>
 
 #include <deal.II/distributed/tria.h>
 
@@ -28,154 +28,111 @@
 #include <deal.II/particles/particle_handler.h>
 #include <deal.II/particles/particle_iterator.h>
 
-using namespace dealii;
-
 template <int dim>
 class DEMContactManager;
 
 /**
- * @brief This class is used for broad particle-particle contact search.
- * Broad search is used to obtain all the particle-particle pairs in adjacent
- * cells.
+ * @brief Finds a vector of pairs (particle_particle_candidates) which shows the
+ * candidate particle-particle collision pairs. These collision pairs will be
+ * used in the fine search to investigate if they are in contact or not.
+ *
+ * @param particle_handler The particle handler of particles in the broad
+ * search
+ * @param container_manager The container manager object that contains
+ * containers to modify of contact pair candidates with other
+ * containers with neighbors lists
  */
-
 template <int dim>
-class ParticleParticleBroadSearch
-{
-public:
-  ParticleParticleBroadSearch();
-
-  /**
-   * @brief Finds a vector of pairs (particle_particle_candidates) which shows the
-   * candidate particle-particle collision pairs. These collision pairs will be
-   * used in the fine search to investigate if they are in contact or not.
-   *
-   * @param particle_handler The particle handler of particles in the broad
-   * search
-   * @param container_manager The container manager object that contains
-   * containers to modify of contact pair candidates with other
-   * containers with neighbors lists
-   */
-
-  void
-  find_particle_particle_contact_pairs(
-    dealii::Particles::ParticleHandler<dim> &particle_handler,
-    DEMContactManager<dim>                  &container_manager);
-
-  /**
-   * @brief Finds a vector of pairs (particle_particle_candidates) which shows the
-   * candidate particle-particle collision pairs. These collision pairs will be
-   * used in the fine search to investigate if they are in contact or not.
-   * This version of the function is used when adaptive sparse contacts regards
-   * mobility is enable.
-   *
-   * @param particle_handler The particle handler of particles in the broad
-   * search
-   * @param container_manager The container manager object that contains
-   * containers to modify of contact pair candidates with other
-   * containers with neighbors lists
-   * @param sparse_contacts_object The object that contains the
-   * information about the mobility status of cells
-   */
-
-  void
-  find_particle_particle_contact_pairs(
-    dealii::Particles::ParticleHandler<dim> &particle_handler,
-    DEMContactManager<dim>                  &container_manager,
-    const AdaptiveSparseContacts<dim>       &sparse_contacts_object);
+void
+find_particle_particle_contact_pairs(
+  dealii::Particles::ParticleHandler<dim> &particle_handler,
+  DEMContactManager<dim>                  &container_manager);
 
-  /**
-   * @brief Finds a vector of pairs (particle_particle_candidates) which contains the
-   * candidate particle-particle collision pairs. These collision pairs will be
-   * used in the fine search to investigate if they are in contact or not.
-   *
-   * @param particle_handler The particle handler of particles in the broad
-   * search
-   * @param container_manager The container manager object that contains
-   * containers to modify of contact pair periodic candidates with other
-   * containers with periodic neighbors lists
-   */
-
-  void
-  find_particle_particle_periodic_contact_pairs(
-    dealii::Particles::ParticleHandler<dim> &particle_handler,
-    DEMContactManager<dim>                  &container_manager);
-
-  /**
-   * @brief Finds a vector of pairs (particle_particle_candidates) which contains the
-   * candidate particle-particle collision pairs. These collision pairs will be
-   * used in the fine search to investigate if they are in contact or not.
-   * This version of the function is used when adaptive sparse contacts regards
-   * mobility is enable.
-   *
-   * @param particle_handler The particle handler of particles in the broad
-   * search
-   * @param container_manager The container manager object that contains
-   * containers to modify of contact pair periodic candidates with other
-   * containers with periodic neighbors lists
-   * @param sparse_contacts_object The object that contains the
-   * information about the mobility status of cells
-   */
-
-  void
-  find_particle_particle_periodic_contact_pairs(
-    dealii::Particles::ParticleHandler<dim> &particle_handler,
-    DEMContactManager<dim>                  &container_manager,
-    const AdaptiveSparseContacts<dim>       &sparse_contacts_object);
-
-private:
-  /**
-   * @brief Stores the candidate particle-particle collision pairs with a given
-   * particle iterator. particle_begin iterator is useful to skip storage of the
-   * first particle in main cell (particle_begin will be the iterator after the
-   * particles_to_evaluate.begin() in that case). When particle_begin is
-   * particles_to_evaluate.begin(), it stores all the particle id in
-   * contact_pair_candidates.
-   *
-   * @param particle_begin The particle iterator to start storing particle ids.
-   * @param particles_to_evaluate The particle range to evaluate with the
-   * particle iterator prior storing ids.
-   * @param contact_pair_candidates A map which will contain all the particle
-   * pairs candidate.
-   */
-  inline void
-  store_candidates(
-    const types::particle_index &main_particle_id,
-    const typename Particles::ParticleHandler<
-      dim>::particle_iterator_range::iterator &particle_begin,
-    const typename Particles::ParticleHandler<dim>::particle_iterator_range
-      &particles_to_evaluate,
-    typename DEM::dem_data_structures<dim>::particle_particle_candidates
-      &contact_pair_candidates)
-  {
-    // Find the contact candidate container of the main particle
-    auto candidates_container_it =
-      contact_pair_candidates.find(main_particle_id);
+/**
+ * @brief Finds a vector of pairs (particle_particle_candidates) which shows the
+ * candidate particle-particle collision pairs. These collision pairs will be
+ * used in the fine search to investigate if they are in contact or not.
+ * This version of the function is used when adaptive sparse contacts is
+ * enabled.
+ *
+ * @param particle_handler The particle handler of particles in the broad
+ * search
+ * @param container_manager The container manager object that contains
+ * containers to modify of contact pair candidates with other
+ * containers with neighbors lists
+ * @param sparse_contacts_object The object that contains the
+ * information about the mobility status of cells
+ */
+template <int dim>
+void
+find_particle_particle_contact_pairs(
+  dealii::Particles::ParticleHandler<dim> &particle_handler,
+  DEMContactManager<dim>                  &container_manager,
+  const AdaptiveSparseContacts<dim>       &sparse_contacts_object);
 
-    // Reserve arbitrary vector capacity and store if the particle does not have
-    // contact candidate yet
-    if (candidates_container_it == contact_pair_candidates.end())
-      {
-        std::vector<types::particle_index> candidates_container;
-        candidates_container.reserve(40);
+/**
+ * @brief Finds a vector of pairs (particle_particle_candidates) which contains the
+ * candidate particle-particle collision pairs. These collision pairs will be
+ * used in the fine search to investigate if they are in contact or not.
+ *
+ * @param particle_handler The particle handler of particles in the broad
+ * search
+ * @param container_manager The container manager object that contains
+ * containers to modify of contact pair periodic candidates with other
+ * containers with periodic neighbors lists
+ */
+template <int dim>
+void
+find_particle_particle_periodic_contact_pairs(
+  dealii::Particles::ParticleHandler<dim> &particle_handler,
+  DEMContactManager<dim>                  &container_manager);
 
-        // Insert the empty vector and get the iterator to the inserted element
-        // prior storing the particle ids
-        auto pair_it_bool =
-          contact_pair_candidates.emplace(main_particle_id,
-                                          candidates_container);
-        candidates_container_it = pair_it_bool.first;
-      }
+/**
+ * @brief Finds a vector of pairs (particle_particle_candidates) which contains the
+ * candidate particle-particle collision pairs. These collision pairs will be
+ * used in the fine search to investigate if they are in contact or not.
+ * This version of the function is used when adaptive sparse contacts is
+ * enabled.
+ *
+ * @param particle_handler The particle handler of particles in the broad
+ * search
+ * @param container_manager The container manager object that contains
+ * containers to modify of contact pair periodic candidates with other
+ * containers with periodic neighbors lists
+ * @param sparse_contacts_object The object that contains the
+ * information about the mobility status of cells
+ */
+template <int dim>
+void
+find_particle_particle_periodic_contact_pairs(
+  dealii::Particles::ParticleHandler<dim> &particle_handler,
+  DEMContactManager<dim>                  &container_manager,
+  const AdaptiveSparseContacts<dim>       &sparse_contacts_object);
 
-    // Store particle ids from the selected particle iterator
-    for (auto particle_iterator = particle_begin;
-         particle_iterator != particles_to_evaluate.end();
-         ++particle_iterator)
-      {
-        candidates_container_it->second.emplace_back(
-          particle_iterator->get_id());
-      }
-  }
-};
+/**
+ * @brief Stores the candidate particle-particle collision pairs with a given
+ * particle iterator. particle_begin iterator is useful to skip storage of the
+ * first particle in main cell (particle_begin will be the iterator after the
+ * particles_to_evaluate.begin() in that case). When particle_begin is
+ * particles_to_evaluate.begin(), it stores all the particle id in
+ * contact_pair_candidates.
+ *
+ * @param main_particle_id The id of the main particle to store the candidate.
+ * @param particle_begin The particle iterator to start storing particle ids.
+ * @param particles_to_evaluate The particle range to evaluate with the
+ * particle iterator prior storing ids.
+ * @param contact_pair_candidates A map which will contain all the particle
+ * pairs candidate.
+ */
+template <int dim>
+inline void
+store_candidates(
+  const types::particle_index &main_particle_id,
+  const typename dealii::Particles::ParticleHandler<
+    dim>::particle_iterator_range::iterator &particle_begin,
+  const typename dealii::Particles::ParticleHandler<
+    dim>::particle_iterator_range &particles_to_evaluate,
+  typename DEM::dem_data_structures<dim>::particle_particle_candidates
+    &contact_pair_candidates);
 
 #endif
diff --git a/include/dem/particle_particle_contact_force.h b/include/dem/particle_particle_contact_force.h
index 35fa57239e..c34f212b94 100644
--- a/include/dem/particle_particle_contact_force.h
+++ b/include/dem/particle_particle_contact_force.h
@@ -20,10 +20,10 @@
 #include <core/auxiliary_math_functions.h>
 #include <core/dem_properties.h>
 
+#include <dem/contact_info.h>
 #include <dem/data_containers.h>
 #include <dem/dem_contact_manager.h>
 #include <dem/dem_solver_parameters.h>
-#include <dem/particle_particle_contact_info.h>
 #include <dem/rolling_resistance_torque_models.h>
 
 #include <deal.II/particles/particle_handler.h>
diff --git a/include/dem/particle_particle_contact_info.h b/include/dem/particle_particle_contact_info.h
deleted file mode 100644
index 1991f80744..0000000000
--- a/include/dem/particle_particle_contact_info.h
+++ /dev/null
@@ -1,65 +0,0 @@
-/* ---------------------------------------------------------------------
- *
- * Copyright (C) 2019 - 2024 by the Lethe authors
- *
- * This file is part of the Lethe library
- *
- * The Lethe library is free software; you can use it, redistribute
- * it, and/or modify it under the terms of the GNU Lesser General
- * Public License as published by the Free Software Foundation; either
- * version 2.1 of the License, or (at your option) any later version.
- * The full text of the license can be found in the file LICENSE at
- * the top level of the Lethe distribution.
- *
- * ---------------------------------------------------------------------
- */
-
-#ifndef lethe_particle_particle_contact_info_struct_h
-#define lethe_particle_particle_contact_info_struct_h
-
-#include <deal.II/base/tensor.h>
-
-#include <deal.II/particles/particle_iterator.h>
-
-using namespace dealii;
-
-/**
- * @brief Handle the information related to the calculation of the
- * particle-particle contact force. Notably it is responsible for storing
- * information that has to be preserved over multiple iterations of a contact,
- * namely everything related to tangential overlaps
- */
-template <int dim>
-class particle_particle_contact_info
-{
-public:
-  /** @brief Constructor for cases where the two particle iterators are known.
-   *  This constructor is used everywhere in the regular DEM and unresolved
-   * CFD-DEM code
-   *
-   *  @param particle_one The iterator of particle one
-   *  @param particle_two The iterator of particle two
-   */
-
-  inline particle_particle_contact_info(
-    Particles::ParticleIterator<dim> &particle_one,
-    Particles::ParticleIterator<dim> &particle_two)
-    : particle_one(particle_one)
-    , particle_two(particle_two)
-  {}
-
-  /** @brief Dummy constructor for cases where the two particle iterators are unknown.
-   *  This constructor is only used in the resolved CFD-DEM and should be
-   * deprecated eventually.
-   */
-
-  inline particle_particle_contact_info()
-  {}
-
-
-  Tensor<1, 3>                     tangential_overlap;
-  Particles::ParticleIterator<dim> particle_one;
-  Particles::ParticleIterator<dim> particle_two;
-};
-
-#endif /* particle_particle_contact_info_struct_h */
diff --git a/include/dem/particle_particle_fine_search.h b/include/dem/particle_particle_fine_search.h
index 71e7051aad..78378b04ce 100644
--- a/include/dem/particle_particle_fine_search.h
+++ b/include/dem/particle_particle_fine_search.h
@@ -17,17 +17,10 @@
 #ifndef lethe_particle_particle_fine_search_h
 #define lethe_particle_particle_fine_search_h
 
-#include <core/dem_properties.h>
-
 #include <dem/data_containers.h>
-#include <dem/particle_particle_contact_info.h>
 
 #include <deal.II/base/tensor.h>
 
-#include <deal.II/particles/particle.h>
-#include <deal.II/particles/particle_handler.h>
-
-#include <boost/range/adaptor/map.hpp>
 
 template <int dim>
 class DEMContactManager;
@@ -35,52 +28,39 @@ class DEMContactManager;
 using namespace dealii;
 
 /**
- * @brief This class is used for local-local and local-ghost fine
- * particle-particle contact search. Fine search is used to find all the
- * particle pairs which are physically in contact and obtain all the required
- * information for calculation of the contact force.
+ * @brief Iterates over a vector of maps (pairs_in_contact) to see if the
+ * particles which were in contact in the last time step, are still in contact
+ * or not. If they are still in contact it will update the collision info,
+ * including tangential overlap, based on new properties of the particle pair,
+ * if they are not in contact anymore it will delete the pair from the
+ * pairs_in_contact and also its information from pairs_in_contact_info.
+ * Then it iterates over the contact candidates from broad search to see if
+ * they already exist in the pairs_in_contact or not, if they are not in the
+ * pairs_in_contact and they have overlap, the pair will be added to the
+ * pairs_in_contact and its contact information will be stored in the
+ * corresponding element of the pairs_in_contact_info
+ *
+ * @param particle_container A container that is used to obtain iterators to
+ * particles using their ids
+ * @param adjacent_particles A map of maps which stores all the required
+ * information for calculation of the contact force of particle pairs
+ * @param contact_pair_candidates The output of broad search which shows
+ * contact pair candidates
+ * @param neighborhood_threshold A value which defines the neighbor particles
+ * @param periodic_offset A tensor of the periodic offset to change the
+ * particle location of the particles on the periodic boundary 1 side,
+ * the tensor as 0.0 values by default
  */
-
 template <int dim>
-class ParticleParticleFineSearch
-{
-public:
-  ParticleParticleFineSearch();
-
-  /**
-   * @brief Iterates over a vector of maps (pairs_in_contact) to see if the
-   * particles which were in contact in the last time step, are still in contact
-   * or not. If they are still in contact it will update the collision info,
-   * including tangential overlap, based on new properties of the particle pair,
-   * if they are not in contact anymore it will delete the pair from the
-   * pairs_in_contact and also its information from pairs_in_contact_info.
-   * Then it iterates over the contact candidates from broad search to see if
-   * they already exist in the pairs_in_contact or not, if they are not in the
-   * pairs_in_contact and they have overlap, the pair will be added to the
-   * pairs_in_contact and its contact information will be stored in the
-   * corresponding element of the pairs_in_contact_info
-   *
-   * @param particle_container A container that is used to obtain iterators to
-   * particles using their ids
-   * @param adjacent_particles A map of maps which stores all the required
-   * information for calculation of the contact force of particle pairs
-   * @param contact_pair_candidates The output of broad search which shows
-   * contact pair candidates
-   * @param neighborhood_threshold A value which defines the neighbor particles
-   * @param periodic_offset A tensor of the periodic offset to change the
-   * particle location of the particles on the periodic boundary 1 side,
-   * the tensor as 0.0 values by default
-   */
-  void
-  particle_particle_fine_search(
-    typename DEM::dem_data_structures<dim>::particle_index_iterator_map const
-      &particle_container,
-    typename DEM::dem_data_structures<dim>::adjacent_particle_pairs
-      &adjacent_particles,
-    const typename DEM::dem_data_structures<dim>::particle_particle_candidates
-                        &contact_pair_candidates,
-    const double         neighborhood_threshold,
-    const Tensor<1, dim> periodic_offset = Tensor<1, dim>());
-};
+void
+particle_particle_fine_search(
+  typename DEM::dem_data_structures<dim>::particle_index_iterator_map const
+    &particle_container,
+  typename DEM::dem_data_structures<dim>::adjacent_particle_pairs
+    &adjacent_particles,
+  const typename DEM::dem_data_structures<dim>::particle_particle_candidates
+                      &contact_pair_candidates,
+  const double         neighborhood_threshold,
+  const Tensor<1, dim> periodic_offset = Tensor<1, dim>());
 
 #endif
diff --git a/include/dem/particle_point_line_broad_search.h b/include/dem/particle_point_line_broad_search.h
index f88e6044bd..8d097dfb96 100644
--- a/include/dem/particle_point_line_broad_search.h
+++ b/include/dem/particle_point_line_broad_search.h
@@ -18,6 +18,7 @@
 #define lethe_particle_point_line_broad_search_h
 
 #include <dem/adaptive_sparse_contacts.h>
+#include <dem/contact_info.h>
 #include <dem/data_containers.h>
 
 #include <deal.II/distributed/tria.h>
@@ -28,79 +29,101 @@
 using namespace dealii;
 
 /**
- * @brief This class is used for broad particle-line and particle-point contact
- * search.This broad search is used to obtain all the particles located at
- * boundary cells with lines and points.
+ * @brief Find a map of pairs (pair of particle and the boundary vertex
+ * location) which shows the candidate particle-point collision pairs. These
+ * collision pairs will be investigated in the fine search to check if they
+ * are in contact or not.
+ *
+ * @param particle_handler Particle handler of particles located in boundary
+ * cells.
+ * @param boundary_cells_with_points A container of cells which are located at
+ * boundaries with only one vertex.
+ * @param particle_point_contact_candidates The map of particle-point pairs.
+ * Each element of map (pair) contains a contact pair particle located near
+ * boundaries with vertices and the vertex location).
  */
 template <int dim>
-class ParticlePointLineBroadSearch
-{
-public:
-  ParticlePointLineBroadSearch();
-
-  /**
-   * @brief Find a map of pairs (pair of particle and the boundary vertex
-   * location) which shows the candidate particle-point collision pairs. These
-   * collision pairs will be investigated in the fine search to check if they
-   * are in contact or not.
-   *
-   * @param particle_handler Particle handler of particles located in boundary
-   * cells
-   * @param boundary_cells_with_points A container of cells which are located at
-   * boundaries with only one vertex
-   *
-   * @return A map of pairs. Each element of map (pair) contains a contact pair
-   * (particle located near boundaries with vertices and the vertex location)
-   */
+void
+find_particle_point_contact_pairs(
+  const Particles::ParticleHandler<dim> &particle_handler,
+  const std::unordered_map<std::string, cell_point_info<dim>>
+    &boundary_cells_with_points,
   typename DEM::dem_data_structures<dim>::particle_point_candidates
-  find_particle_point_contact_pairs(
-    const Particles::ParticleHandler<dim> &particle_handler,
-    const std::unordered_map<
-      std::string,
-      std::pair<typename Triangulation<dim>::active_cell_iterator, Point<dim>>>
-      &boundary_cells_with_points);
+    &particle_point_contact_candidates);
 
+/**
+ * @brief Find a map of pairs (pair of particle and the boundary vertex
+ * location) which shows the candidate particle-point collision pairs. These
+ * collision pairs will be investigated in the fine search to check if they
+ * are in contact or not.
+ *
+ * @param particle_handler Particle handler of particles located in boundary
+ * cells.
+ * @param boundary_cells_with_points A container of cells which are located at
+ * boundaries with only one vertex.
+ * @param particle_point_contact_candidates The map of particle-point pairs.
+ * Each element of map (pair) contains a contact pair particle located near
+ * boundaries with vertices and the vertex location).
+ * @param sparse_contacts_object The Adaptive Sparse Contacts for mobility
+ * status checks.
+ */
+template <int dim>
+void
+find_particle_point_contact_pairs(
+  const Particles::ParticleHandler<dim> &particle_handler,
+  const std::unordered_map<std::string, cell_point_info<dim>>
+    &boundary_cells_with_points,
   typename DEM::dem_data_structures<dim>::particle_point_candidates
-  find_particle_point_contact_pairs(
-    const Particles::ParticleHandler<dim> &particle_handler,
-    const std::unordered_map<
-      std::string,
-      std::pair<typename Triangulation<dim>::active_cell_iterator, Point<dim>>>
-                                      &boundary_cells_with_points,
-    const AdaptiveSparseContacts<dim> &sparse_contacts_object);
-
-  /**
-   * @brief Find a map of tuples (tuple of particle and the locations of
-   * beginning and ending vertices of the boundary lines) which shows the
-   * candidate particle-line collision pairs. These collision pairs will be
-   * investigated in the fine search to check if they are in contact or not.
-   *
-   * @param particle_handler Particle handler of particles located in boundary
-   * cells
-   * @param boundary_cells_with_lines A container of cells which are located at
-   * boundaries with only one line
-   * @return A map of tuples. Each element of map (tuple) contains a particle
-   * and the locations of beginning and ending vertices of the boundary lines
-   */
+                                    &particle_point_contact_candidates,
+  const AdaptiveSparseContacts<dim> &sparse_contacts_object);
 
+/**
+ * @brief Find a map of tuples (tuple of particle and the locations of
+ * beginning and ending vertices of the boundary lines) which shows the
+ * candidate particle-line collision pairs. These collision pairs will be
+ * investigated in the fine search to check if they are in contact or not.
+ *
+ * @param particle_handler Particle handler of particles located in boundary
+ * cells.
+ * @param boundary_cells_with_lines A container of cells which are located at
+ * boundaries with only one line.
+ * @param particle_line_contact_candidates Each element of map (tuple) contains
+ * a particle and the locations of beginning and ending vertices of the boundary
+ * lines.
+ */
+template <int dim>
+void
+find_particle_line_contact_pairs(
+  const Particles::ParticleHandler<dim> &particle_handler,
+  const std::unordered_map<std::string, cell_line_info<dim>>
+    &boundary_cells_with_lines,
   typename DEM::dem_data_structures<dim>::particle_line_candidates
-  find_particle_line_contact_pairs(
-    const Particles::ParticleHandler<dim> &particle_handler,
-    const std::unordered_map<
-      std::string,
-      std::tuple<typename Triangulation<dim>::active_cell_iterator,
-                 Point<dim>,
-                 Point<dim>>> &boundary_cells_with_lines);
+    &particle_line_contact_candidates);
 
+/**
+ * @brief Find a map of tuples (tuple of particle and the locations of
+ * beginning and ending vertices of the boundary lines) which shows the
+ * candidate particle-line collision pairs. These collision pairs will be
+ * investigated in the fine search to check if they are in contact or not.
+ *
+ * @param particle_handler Particle handler of particles located in boundary
+ * cells.
+ * @param boundary_cells_with_lines A container of cells which are located at
+ * boundaries with only one line.
+ * @param particle_line_contact_candidates Each element of map (tuple) contains
+ * a particle and the locations of beginning and ending vertices of the boundary
+ * lines.
+ * @param sparse_contacts_object The Adaptive Sparse Contacts for mobility
+ * status checks.
+ */
+template <int dim>
+void
+find_particle_line_contact_pairs(
+  const Particles::ParticleHandler<dim> &particle_handler,
+  const std::unordered_map<std::string, cell_line_info<dim>>
+    &boundary_cells_with_lines,
   typename DEM::dem_data_structures<dim>::particle_line_candidates
-  find_particle_line_contact_pairs(
-    const Particles::ParticleHandler<dim> &particle_handler,
-    const std::unordered_map<
-      std::string,
-      std::tuple<typename Triangulation<dim>::active_cell_iterator,
-                 Point<dim>,
-                 Point<dim>>>         &boundary_cells_with_lines,
-    const AdaptiveSparseContacts<dim> &sparse_contacts_object);
-};
+                                    &particle_line_contact_candidates,
+  const AdaptiveSparseContacts<dim> &sparse_contacts_object);
 
 #endif
diff --git a/include/dem/particle_point_line_contact_force.h b/include/dem/particle_point_line_contact_force.h
index e2925ed09e..bfe0b5837e 100644
--- a/include/dem/particle_point_line_contact_force.h
+++ b/include/dem/particle_point_line_contact_force.h
@@ -19,9 +19,9 @@
 
 #include <core/dem_properties.h>
 
+#include <dem/contact_info.h>
 #include <dem/data_containers.h>
 #include <dem/dem_solver_parameters.h>
-#include <dem/particle_point_line_contact_info_struct.h>
 
 #include <iostream>
 #include <vector>
@@ -53,10 +53,10 @@ class ParticlePointLineForce
    */
   void
   calculate_particle_point_contact_force(
-    const typename DEM::dem_data_structures<dim>::
-      particle_point_line_contact_info *particle_point_line_pairs_in_contact,
+    const typename DEM::dem_data_structures<dim>::particle_point_in_contact
+      *particle_point_pairs_in_contact,
     const Parameters::Lagrangian::LagrangianPhysicalProperties
-                              &lagrangian_physical_properties,
+                              &physical_properties,
     std::vector<Tensor<1, 3>> &force);
 
   /**
@@ -71,10 +71,10 @@ class ParticlePointLineForce
    */
   void
   calculate_particle_line_contact_force(
-    const typename DEM::dem_data_structures<
-      dim>::particle_point_line_contact_info *particle_line_pairs_in_contact,
+    const typename DEM::dem_data_structures<dim>::particle_line_in_contact
+      *particle_line_pairs_in_contact,
     const Parameters::Lagrangian::LagrangianPhysicalProperties
-                              &lagrangian_physical_properties,
+                              &physical_properties,
     std::vector<Tensor<1, 3>> &force);
 
 private:
diff --git a/include/dem/particle_point_line_contact_info_struct.h b/include/dem/particle_point_line_contact_info_struct.h
deleted file mode 100644
index fd13e2e081..0000000000
--- a/include/dem/particle_point_line_contact_info_struct.h
+++ /dev/null
@@ -1,40 +0,0 @@
-/* ---------------------------------------------------------------------
- *
- * Copyright (C) 2019 - 2024 by the Lethe authors
- *
- * This file is part of the Lethe library
- *
- * The Lethe library is free software; you can use it, redistribute
- * it, and/or modify it under the terms of the GNU Lesser General
- * Public License as published by the Free Software Foundation; either
- * version 2.1 of the License, or (at your option) any later version.
- * The full text of the license can be found in the file LICENSE at
- * the top level of the Lethe distribution.
- *
- * ---------------------------------------------------------------------
- */
-
-#ifndef lethe_particle_point_line_contact_info_struct_h
-#define lethe_particle_point_line_contact_info_struct_h
-
-#include <deal.II/base/point.h>
-#include <deal.II/base/tensor.h>
-
-#include <deal.II/particles/particle_iterator.h>
-
-/**
- * @brief Handle information related to the calculation of the particle-point
- * and particle-line contact forces.
- */
-
-using namespace dealii;
-
-template <int dim>
-struct particle_point_line_contact_info_struct
-{
-  Particles::ParticleIterator<dim> particle;
-  Point<3>                         point_one;
-  Point<3>                         point_two;
-};
-
-#endif
diff --git a/include/dem/particle_point_line_fine_search.h b/include/dem/particle_point_line_fine_search.h
index 31e8d1c676..a33fe79e5b 100644
--- a/include/dem/particle_point_line_fine_search.h
+++ b/include/dem/particle_point_line_fine_search.h
@@ -17,83 +17,68 @@
 #ifndef lethe_particle_point_line_fine_search_h
 #define lethe_particle_point_line_fine_search_h
 
-#include <core/dem_properties.h>
-
 #include <dem/data_containers.h>
-#include <dem/particle_point_line_contact_info_struct.h>
-
-#include <deal.II/particles/particle_handler.h>
-
-#include <iostream>
-#include <vector>
 
 using namespace dealii;
 
 /**
- * @brief Handle the fine particle-point and particle-line contact search.
- * Fine search is used to find all the particles which are physically in contact
- * with boundary lines and points, and obtain all the required information for
- * calculation of the corresponding contact force.
+ * @brief Iterate over a map of pairs (particle_point_contact_candidates)
+ * to see if the particle-point pairs are in contact or not. If they are in
+ * contact, the normal overlap, normal vector of contact and contact normal
+ * relative velocity are stored in a map which is the output of this function.
+ *
+ * @param particle_point_contact_candidates The output of particle-point broad
+ * search which shows contact pair candidates.
+ * @param neighborhood_threshold A value which defines the neighbor particles.
+ * @param particle_point_pairs_in_contact A map which contains all the required
+ * information (normal overlap, normal vector and contact normal relative
+ * velocity) for calculation of particle-point contact force.
  */
 template <int dim>
-class ParticlePointLineFineSearch
-{
-public:
-  ParticlePointLineFineSearch();
+void
+particle_point_fine_search(
+  const typename DEM::dem_data_structures<dim>::particle_point_candidates
+              &particle_point_contact_candidates,
+  const double neighborhood_threshold,
+  typename DEM::dem_data_structures<dim>::particle_point_in_contact
+    &particle_point_pairs_in_contact);
 
-  /**
-   * @brief Iterate over a map of pairs (particle_point_contact_candidates)
-   * to see if the particle-point pairs are in contact or not. If they are in
-   * contact, the normal overlap, normal vector of contact and contact normal
-   * relative velocity are stored in a map which is the output of this function
-   *
-   * @param particle_point_contact_candidates The output of particle-point broad
-   * search which shows contact pair candidates
-   * @param neighborhood_threshold A value which defines the neighbor particles
-   * @return A map which contains all the required information (normal overlap,
-   * normal vector and contact normal relative velocity) for calculation of
-   * particle-point contact force
-   */
-  typename DEM::dem_data_structures<dim>::particle_point_line_contact_info
-  particle_point_fine_search(
-    const typename DEM::dem_data_structures<dim>::particle_point_candidates
-                &particle_point_contact_candidates,
-    const double neighborhood_threshold);
-
-  /**
-   * @brief Iterate over a map of tuples (particle_line_contact_candidates) to
-   * see if the particle-line pairs are in contact or not. If they are in
-   * contact, the normal overlap, normal vector of contact and contact normal
-   * relative velocity are stored in a map which is the output of this function
-   *
-   * @param particle_line_contact_candidates The output of particle-line broad
-   * search which shows contact pair candidates
-   * @param neighborhood_threshold A value which defines the neighbor particles
-   * @return A map which contains all the required information (normal overlap,
-   * normal vector and contact normal relative velocity) for calculation of
-   * particle-line contact force
-   */
-  typename DEM::dem_data_structures<dim>::particle_point_line_contact_info
-  particle_line_fine_search(
-    const typename DEM::dem_data_structures<dim>::particle_line_candidates
-                &particle_line_contact_candidates,
-    const double neighborhood_threshold);
-
-private:
-  /**
-   * @brief This private function is used to find the projection of point_p on
-   * a line with beginning and ending vertices of point_a and point_b,
-   * respectively
-   * @param point_p A point which is going to be projected on the line
-   * @param point_a Beginning point of the line
-   * @param point_b Ending point of the line
-   * @return The projection of point_p on the line (from point_a to point_b)
-   */
+/**
+ * @brief Iterate over a map of tuples (particle_line_contact_candidates) to
+ * see if the particle-line pairs are in contact or not. If they are in
+ * contact, the normal overlap, normal vector of contact and contact normal
+ * relative velocity are stored in a map which is the output of this function.
+ *
+ * @param particle_line_contact_candidates The output of particle-line broad
+ * search which shows contact pair candidates.
+ * @param neighborhood_threshold A value which defines the neighbor particles.
+ * @param particle_line_pairs_in_contact A map which contains all the required
+ * information (normal overlap, normal vector and contact normal relative
+ * velocity) for calculation of particle-line contact force.
+ */
+template <int dim>
+void
+particle_line_fine_search(
+  const typename DEM::dem_data_structures<dim>::particle_line_candidates
+              &particle_line_contact_candidates,
+  const double neighborhood_threshold,
+  typename DEM::dem_data_structures<dim>::particle_line_in_contact
+    &particle_line_pairs_in_contact);
 
-  Point<3>
-  find_projection_point(const Point<3> &point_p,
-                        const Point<3> &point_a,
-                        const Point<3> &point_b);
-};
+/**
+ * @brief This private function is used to find the projection of point_p on
+ * a line with beginning and ending vertices of point_a and point_b,
+ * respectively
+ *
+ * @param point_p A point which is going to be projected on the line
+ * @param point_a Beginning point of the line
+ * @param point_b Ending point of the line
+ *
+ * @return The projection of point_p on the line (from point_a to point_b)
+ */
+inline Point<3>
+find_projection_point(const Point<3> &point_p,
+                      const Point<3> &point_a,
+                      const Point<3> &point_b);
 
 #endif
diff --git a/include/dem/particle_wall_broad_search.h b/include/dem/particle_wall_broad_search.h
index 77c9c18fb5..54f0c0fb5d 100644
--- a/include/dem/particle_wall_broad_search.h
+++ b/include/dem/particle_wall_broad_search.h
@@ -20,229 +20,201 @@
 
 #include <dem/adaptive_sparse_contacts.h>
 #include <dem/boundary_cells_info_struct.h>
+#include <dem/contact_info.h>
 #include <dem/data_containers.h>
 #include <dem/dem_solver_parameters.h>
 
 #include <deal.II/distributed/tria.h>
 
-#include <deal.II/dofs/dof_handler.h>
-
-#include <deal.II/particles/particle.h>
 #include <deal.II/particles/particle_handler.h>
-#include <deal.II/particles/particle_iterator.h>
 
 #include <iostream>
+#include <map>
 #include <vector>
 
 using namespace dealii;
 
 /**
- * @brief Broad search is used to obtain all the particles located at boundary
+ * @brief Finds unordered map of tuples (tuple of particle located in
+ * boundary cells, normal vector of the boundary face, a
+ * point on the face and the corresponding boundary cell) which shows the
+ * candidate particle-wall collision pairs. These collision candidates will be
+ * investigated in the fine search to check if they are in contact or not.
+ *
+ * @param boundary_cells_information Information of the boundary cells and
+ * faces. This is the output of the FindBoundaryCellsInformation class.
+ * @param particle_handler Particle handler of particles located in boundary
  * cells.
+ * @param particle_wall_contact_candidates A two-layered unordered map of tuples. Each
+ * tuple contains a particle located near boundaries, the normal vector of
+ * the corresponding face boundary, a point on the boundary and the boundary
+ * cell. The contact pair is used in the fine search.
+ */
+template <int dim>
+void
+find_particle_wall_contact_pairs(
+  const std::map<int, boundary_cells_info_struct<dim>>
+                                        &boundary_cells_information,
+  const Particles::ParticleHandler<dim> &particle_handler,
+  typename DEM::dem_data_structures<dim>::particle_wall_candidates
+    &particle_wall_contact_candidates);
+
+template <int dim>
+void
+find_particle_wall_contact_pairs(
+  const std::map<int, boundary_cells_info_struct<dim>>
+                                        &boundary_cells_information,
+  const Particles::ParticleHandler<dim> &particle_handler,
+  typename DEM::dem_data_structures<dim>::particle_wall_candidates
+                                    &particle_wall_contact_candidates,
+  const AdaptiveSparseContacts<dim> &sparse_contacts_object);
+
+/**
+ * @brief Find a two-layered unordered map of particle iterators which shows the
+ * candidate particle-floating wall collision candidates. These collision
+ * pairs will be investigated in the fine search to check if they are in
+ * contact or not
+ *
+ * @param boundary_cells_for_floating_walls Boundary cells located adjacent to
+ * floating walls
+ * @param particle_handler Particle handler of particles located in boundary
+ * cells
+ * @param floating_wall_properties Properties of the floating walls specified
+ * in the parameter handler file
+ * @param simulation_time Simulation time
+ * @param particle_floating_wall_candidates Output of particle-floating wall
+ * broad search which contains all the particle-floating wall collision
+ * candidates
+ */
+template <int dim>
+void
+find_particle_floating_wall_contact_pairs(
+  const std::unordered_map<
+    types::global_dof_index,
+    std::set<typename Triangulation<dim>::active_cell_iterator>>
+                                        &boundary_cells_for_floating_walls,
+  const Particles::ParticleHandler<dim> &particle_handler,
+  const Parameters::Lagrangian::FloatingWalls<dim> &floating_wall_properties,
+  const double                                      simulation_time,
+  typename DEM::dem_data_structures<dim>::particle_floating_wall_candidates
+    &particle_floating_wall_candidates);
+template <int dim>
+void
+find_particle_floating_wall_contact_pairs(
+  const std::unordered_map<
+    types::global_dof_index,
+    std::set<typename Triangulation<dim>::active_cell_iterator>>
+                                        &boundary_cells_for_floating_walls,
+  const Particles::ParticleHandler<dim> &particle_handler,
+  const Parameters::Lagrangian::FloatingWalls<dim> &floating_wall_properties,
+  const double                                      simulation_time,
+  typename DEM::dem_data_structures<dim>::particle_floating_wall_candidates
+                                    &particle_floating_wall_candidates,
+  const AdaptiveSparseContacts<dim> &sparse_contacts_object);
+
+/**
+ * @brief Find a two-layered unordered map
+ * (particle_floating_mesh_contact_candidates) of particle iterators that
+ * shows the candidate particle-floating mesh collision candidates. These
+ * collision pairs will be investigated in the fine search to check if they
+ * are in contact or not
+ *
+ * @param solid_surfaces_mesh_information Information of the solid surfaces mapped
+ * in the background triangulation.
+ * @param particle_handler
+ * @param particle_floating_mesh_contact_candidates Particle-floating mesh contact
+ * candidates
+ * @param cells_total_neighbor_list A container in which all the neighbor cells
+ * of the local cells are stored
+ * @param sparse_contacts_object The object that contains the
+ * information about the mobility status of cells
+ */
+template <int dim>
+void
+particle_solid_surfaces_contact_search(
+  const typename DEM::dem_data_structures<dim>::solid_surfaces_mesh_information
+                                        &solid_surfaces_mesh_information,
+  const Particles::ParticleHandler<dim> &particle_handler,
+  typename DEM::dem_data_structures<dim>::particle_floating_mesh_candidates
+    &particle_floating_mesh_contact_candidates,
+  typename DEM::dem_data_structures<dim>::cells_total_neighbor_list
+    &cells_total_neighbor_list);
+
+template <int dim>
+void
+particle_solid_surfaces_contact_search(
+  const typename DEM::dem_data_structures<dim>::solid_surfaces_mesh_information
+                                        &solid_surfaces_mesh_information,
+  const Particles::ParticleHandler<dim> &particle_handler,
+  typename DEM::dem_data_structures<dim>::particle_floating_mesh_candidates
+    &particle_floating_mesh_contact_candidates,
+  typename DEM::dem_data_structures<dim>::cells_total_neighbor_list
+                                    &cells_total_neighbor_list,
+  const AdaptiveSparseContacts<dim> &sparse_contacts_object);
+
+
+/**
+ * @brief Store the candidate particle-wall contact pairs.
+ *
+ * @param particle_iterator The particle iterator to start boundary id.
+ * @param boundary_cells_content The info of the current boundary cell to
+ * store.
+ * @param contact_pair_candidates A map which will contain all the
+ * particle/wall pairs candidate.
  */
 template <int dim>
-class ParticleWallBroadSearch
+inline void
+store_candidates(
+  const typename Particles::ParticleHandler<
+    dim>::particle_iterator_range::iterator &particle_iterator,
+  const boundary_cells_info_struct<dim>     &boundary_cells_content,
+  typename DEM::dem_data_structures<dim>::particle_wall_candidates
+    &contact_pair_candidates)
 {
-public:
-  ParticleWallBroadSearch();
-
-  /**
-   * @brief Finds unordered map of tuples (tuple of particle located in
-   * boundary cells, normal vector of the boundary face, a
-   * point on the face and the corresponding boundary cell) which shows the
-   * candidate particle-wall collision pairs. These collision candidates will be
-   * investigated in the fine search to check if they are in contact or not
-   *
-   * @param boundary_cells_information Information of the boundary cells and
-   * faces. This is the output of the FindBoundaryCellsInformation class
-   * @param particle_handler Particle handler of particles located in boundary
-   * cells
-   * @param particle_wall_contact_candidates A two-layered unordered map of tuples. Each
-   * tuple contains a particle located near boundaries, the normal vector of
-   * the corresponding face boundary, a point on the boundary and the boundary
-   * cell. The contact pair is used in the fine search
-   */
-  void
-  find_particle_wall_contact_pairs(
-    const std::map<int, boundary_cells_info_struct<dim>>
-                                          &boundary_cells_information,
-    const Particles::ParticleHandler<dim> &particle_handler,
-    typename DEM::dem_data_structures<dim>::particle_wall_candidates
-      &particle_wall_contact_candidates);
-
-  void
-  find_particle_wall_contact_pairs(
-    const std::map<int, boundary_cells_info_struct<dim>>
-                                          &boundary_cells_information,
-    const Particles::ParticleHandler<dim> &particle_handler,
-    typename DEM::dem_data_structures<dim>::particle_wall_candidates
-                                      &particle_wall_contact_candidates,
-    const AdaptiveSparseContacts<dim> &sparse_contacts_object);
-
-  /**
-   * @brief Find a two-layered unordered map of particle iterators which shows the
-   * candidate particle-floating wall collision candidates. These collision
-   * pairs will be investigated in the fine search to check if they are in
-   * contact or not
-   *
-   * @param boundary_cells_for_floating_walls Boundary cells located adjacent to
-   * floating walls
-   * @param particle_handler Particle handler of particles located in boundary
-   * cells
-   * @param floating_wall_properties Properties of the floating walls specified
-   * in the parameter handler file
-   * @param simulation_time Simulation time
-   * @param particle_floating_wall_candidates Output of particle-floating wall
-   * broad search which contains all the particle-floating wall collision
-   * candidates
-   */
-
-  void
-  find_particle_floating_wall_contact_pairs(
-    const std::unordered_map<
-      types::global_dof_index,
-      std::set<typename Triangulation<dim>::active_cell_iterator>>
-                                          &boundary_cells_for_floating_walls,
-    const Particles::ParticleHandler<dim> &particle_handler,
-    const Parameters::Lagrangian::FloatingWalls<dim> &floating_wall_properties,
-    const double                                      simulation_time,
-    typename DEM::dem_data_structures<dim>::particle_floating_wall_candidates
-      &particle_floating_wall_candidates);
-
-  void
-  find_particle_floating_wall_contact_pairs(
-    const std::unordered_map<
-      types::global_dof_index,
-      std::set<typename Triangulation<dim>::active_cell_iterator>>
-                                          &boundary_cells_for_floating_walls,
-    const Particles::ParticleHandler<dim> &particle_handler,
-    const Parameters::Lagrangian::FloatingWalls<dim> &floating_wall_properties,
-    const double                                      simulation_time,
-    typename DEM::dem_data_structures<dim>::particle_floating_wall_candidates
-                                      &particle_floating_wall_candidates,
-    const AdaptiveSparseContacts<dim> &sparse_contacts_object);
-
-  /**
-   * @brief Find a two-layered unordered map
-   * (particle_floating_mesh_contact_candidates) of particle iterators that
-   * shows the candidate particle-floating mesh collision candidates. These
-   * collision pairs will be investigated in the fine search to check if they
-   * are in contact or not
-   *
-   * @param solid_surfaces_mesh_information Information of the solid surfaces mapped
-   * in the background triangulation.
-   * @param particle_handler
-   * @param particle_floating_mesh_contact_candidates Particle-floating mesh contact
-   * candidates
-   * @param cells_total_neighbor_list A container in which all the neighbor cells
-   * of the local cells are stored
-   * @param sparse_contacts_object The object that contains the
-   * information about the mobility status of cells
-   */
-  void
-  particle_solid_surfaces_contact_search(
-    const typename DEM::dem_data_structures<
-      dim>::solid_surfaces_mesh_information &solid_surfaces_mesh_information,
-    const Particles::ParticleHandler<dim>   &particle_handler,
-    typename DEM::dem_data_structures<dim>::particle_floating_mesh_candidates
-      &particle_floating_mesh_contact_candidates,
-    typename DEM::dem_data_structures<dim>::cells_total_neighbor_list
-      &cells_total_neighbor_list);
-
-  void
-  particle_solid_surfaces_contact_search(
-    const typename DEM::dem_data_structures<
-      dim>::solid_surfaces_mesh_information &solid_surfaces_mesh_information,
-    const Particles::ParticleHandler<dim>   &particle_handler,
-    typename DEM::dem_data_structures<dim>::particle_floating_mesh_candidates
-      &particle_floating_mesh_contact_candidates,
-    typename DEM::dem_data_structures<dim>::cells_total_neighbor_list
-                                      &cells_total_neighbor_list,
-    const AdaptiveSparseContacts<dim> &sparse_contacts_object);
-
-private:
-  /**
-   * @brief Store the candidate particle-wall contact pairs.
-   *
-   * @param particle_iterator The particle iterator to start boundary id.
-   * @param boundary_cells_content The info of the current boundary cell to
-   * store.
-   * @param contact_pair_candidates A map which will contain all the
-   * particle/wall pairs candidate.
-   */
-  inline void
-  store_candidates(
-    const typename Particles::ParticleHandler<
-      dim>::particle_iterator_range::iterator &particle_iterator,
-    const boundary_cells_info_struct<dim>     &boundary_cells_content,
-    typename DEM::dem_data_structures<dim>::particle_wall_candidates
-      &contact_pair_candidates)
-  {
-    // Find the contact candidate container of the particle
-    const types::particle_index particle_id = particle_iterator->get_id();
-    auto candidates_container_it = contact_pair_candidates.find(particle_id);
-
-    // Reserve arbitrary vector capacity and store if the particle does not have
-    // contact candidate yet
-    if (candidates_container_it == contact_pair_candidates.end())
-      {
-        auto pair_it_bool = contact_pair_candidates.emplace(
-          particle_id,
-          ankerl::unordered_dense::map<
-            DEM::global_face_id,
-            std::tuple<Particles::ParticleIterator<dim>,
-                       Tensor<1, dim>,
-                       Point<dim>,
-                       DEM::global_face_id>>());
-
-        candidates_container_it = pair_it_bool.first;
-      }
-
-    // Store particle ids from the selected particle iterator
-    candidates_container_it->second.emplace(
-      boundary_cells_content.global_face_id,
-      std::make_tuple(particle_iterator,
-                      boundary_cells_content.normal_vector,
-                      boundary_cells_content.point_on_face,
-                      boundary_cells_content.boundary_id));
-  }
-
-  /**
-   * @brief Store the candidate particle-floating wall contact pairs.
-   *
-   * @param particle_iterator The particle iterator to start boundary id.
-   * @param floating_wall_id The floating wall id.
-   * @param contact_pair_candidates A map which will contain all the
-   * particle/wall pairs candidate.
-   */
-  inline void
-  store_candidates(
-    const typename Particles::ParticleHandler<
-      dim>::particle_iterator_range::iterator &particle_iterator,
-    const DEM::global_face_id                 &floating_wall_id,
-    typename DEM::dem_data_structures<dim>::particle_floating_wall_candidates
-      &contact_pair_candidates)
-  {
-    // Find the contact candidate container of the particle
-    const types::particle_index particle_id = particle_iterator->get_id();
-    auto candidates_container_it = contact_pair_candidates.find(particle_id);
-
-    // Reserve arbitrary vector capacity and store if the particle does not have
-    // contact candidate yet
-    if (candidates_container_it == contact_pair_candidates.end())
-      {
-        auto pair_it_bool = contact_pair_candidates.emplace(
-          particle_id,
-          std::unordered_map<DEM::global_face_id,
-                             Particles::ParticleIterator<dim>>());
-
-        candidates_container_it = pair_it_bool.first;
-      }
-
-    // Store particle ids from the selected particle iterator
-    candidates_container_it->second.emplace(floating_wall_id,
-                                            particle_iterator);
-  }
-};
+  // Find the contact candidate container of the particle
+  const types::particle_index particle_id = particle_iterator->get_id();
+  auto candidates_container_it = contact_pair_candidates.find(particle_id);
+
+  // Reserve arbitrary vector capacity and store if the particle does not have
+  // contact candidate yet
+  if (candidates_container_it == contact_pair_candidates.end())
+    {
+      auto pair_it_bool = contact_pair_candidates.emplace(
+        particle_id,
+        ankerl::unordered_dense::map<
+          DEM::global_face_id,
+          std::tuple<Particles::ParticleIterator<dim>,
+                     Tensor<1, dim>,
+                     Point<dim>,
+                     DEM::global_face_id>>());
+
+      candidates_container_it = pair_it_bool.first;
+    }
+
+  // Store particle ids from the selected particle iterator
+  candidates_container_it->second.emplace(
+    boundary_cells_content.global_face_id,
+    std::make_tuple(particle_iterator,
+                    boundary_cells_content.normal_vector,
+                    boundary_cells_content.point_on_face,
+                    boundary_cells_content.boundary_id));
+}
+
+/**
+ * @brief Store the candidate particle-floating wall contact pairs.
+ *
+ * @param particle_iterator The particle iterator to start boundary id.
+ * @param floating_wall_id The floating wall id.
+ * @param contact_pair_candidates A map which will contain all the
+ * particle/wall pairs candidate.
+ */
+template <int dim>
+inline void
+store_candidates(
+  const typename Particles::ParticleHandler<
+    dim>::particle_iterator_range::iterator &particle_iterator,
+  const DEM::global_face_id                 &floating_wall_id,
+  typename DEM::dem_data_structures<dim>::particle_floating_wall_candidates
+    &contact_pair_candidates);
 
 #endif
diff --git a/include/dem/particle_wall_contact_force.h b/include/dem/particle_wall_contact_force.h
index fefddd6da6..550905246e 100644
--- a/include/dem/particle_wall_contact_force.h
+++ b/include/dem/particle_wall_contact_force.h
@@ -21,9 +21,9 @@
 #include <core/dem_properties.h>
 #include <core/serial_solid.h>
 
+#include <dem/contact_info.h>
 #include <dem/data_containers.h>
 #include <dem/dem_solver_parameters.h>
-#include <dem/particle_wall_contact_info.h>
 
 #include <boost/math/special_functions.hpp>
 #include <boost/range/adaptor/map.hpp>
diff --git a/include/dem/particle_wall_dmt_force.h b/include/dem/particle_wall_dmt_force.h
index 36a46fe04f..32eacfe309 100644
--- a/include/dem/particle_wall_dmt_force.h
+++ b/include/dem/particle_wall_dmt_force.h
@@ -20,7 +20,6 @@
 #include <core/dem_properties.h>
 
 #include <dem/dem_solver_parameters.h>
-#include <dem/particle_wall_contact_info.h>
 #include <dem/particle_wall_nonlinear_force.h>
 
 #include <deal.II/particles/particle.h>
diff --git a/include/dem/particle_wall_fine_search.h b/include/dem/particle_wall_fine_search.h
index c299fbc477..016635b7fb 100644
--- a/include/dem/particle_wall_fine_search.h
+++ b/include/dem/particle_wall_fine_search.h
@@ -17,98 +17,72 @@
 #ifndef lethe_particle_wall_fine_search_h
 #define lethe_particle_wall_fine_search_h
 
-#include <core/dem_properties.h>
-
 #include <dem/data_containers.h>
 #include <dem/dem_solver_parameters.h>
-#include <dem/particle_wall_contact_info.h>
-
-#include <deal.II/particles/particle.h>
-#include <deal.II/particles/particle_handler.h>
-#include <deal.II/particles/particle_iterator.h>
-
-#include <iostream>
-#include <vector>
-
-using namespace dealii;
 
 /**
- * @brief This class is used for fine particle-wall contact search. Fine search
- * is used to find all the particles which are physically in contact with
- * system boundaries and obtain all the required information for calculation
- * of the corresponding particle-wall contact force.
+ * @brief Iterate over the contact candidates from particle-wall broad search
+ * (particle_wall_contact_pair_candidates) to add new contact pairs to the
+ * particle_wall_pairs_in_contact container
+ *
+ * @param particle_wall_contact_pair_candidates The output of particle-wall broad search
+ * which shows contact pair candidates
+ * @param particle_wall_pairs_in_contact An unordered_map of maps which stores
+ * all the particle-wall pairs which are physically in contact, and the
+ * contact information in a struct. Note that the size of this unordered map
+ * is equal to the number of particles
  */
-
 template <int dim>
-class ParticleWallFineSearch
-{
-public:
-  ParticleWallFineSearch();
-
-  /**
-   * @brief Iterate over the contact candidates from particle-wall broad search
-   * (particle_wall_contact_pair_candidates) to add new contact pairs to the
-   * particle_wall_pairs_in_contact container
-   *
-   * @param particle_wall_contact_pair_candidates The output of particle-wall broad search
-   * which shows contact pair candidates
-   * @param particle_wall_pairs_in_contact An unordered_map of maps which stores
-   * all the particle-wall pairs which are physically in contact, and the
-   * contact information in a struct. Note that the size of this unordered map
-   * is equal to the number of particles
-   */
-
-  void
-  particle_wall_fine_search(
-    const typename DEM::dem_data_structures<dim>::particle_wall_candidates
-      &particle_wall_contact_pair_candidates,
-    typename DEM::dem_data_structures<dim>::particle_wall_in_contact
-      &particle_wall_pairs_in_contact);
-
-  /**
-   * @brief Iterate over the contact candidates from particle-floating wall broad
-   * search (particle_floating_wall_candidates) to add new contact pairs to the
-   * particle_floating_wall_pairs_in_contact container
-   *
-   * @param particle_floating_wall_pair_candidates The output of particle-floating wall
-   * broad search which shows contact pair candidates
-   * @param floating_wall_properties Properties of floating walls defined in the
-   * parameter handler
-   * @param simulation_time Simulation time
-   * @param particle_floating_wall_pairs_in_contact An unordered_map of maps which stores
-   * all the particle-floating wall pairs which are in contact, and
-   * the contact information in a struct. Note that the size of this unordered
-   * map is equal to the number of particles
-   */
-  void
-  particle_floating_wall_fine_search(
-    const typename DEM::dem_data_structures<
-      dim>::particle_floating_wall_candidates
-      &particle_floating_wall_contact_candidates,
-    const Parameters::Lagrangian::FloatingWalls<dim> &floating_wall_properties,
-    const double                                      simulation_time,
-    typename DEM::dem_data_structures<dim>::particle_wall_in_contact
-      &particle_floating_wall_in_contact);
+void
+particle_wall_fine_search(
+  const typename DEM::dem_data_structures<dim>::particle_wall_candidates
+    &particle_wall_contact_pair_candidates,
+  typename DEM::dem_data_structures<dim>::particle_wall_in_contact
+    &particle_wall_pairs_in_contact);
 
+/**
+ * @brief Iterate over the contact candidates from particle-floating wall broad
+ * search (particle_floating_wall_candidates) to add new contact pairs to the
+ * particle_floating_wall_pairs_in_contact container
+ *
+ * @param particle_floating_wall_pair_candidates The output of particle-floating wall
+ * broad search which shows contact pair candidates
+ * @param floating_wall_properties Properties of floating walls defined in the
+ * parameter handler
+ * @param simulation_time Simulation time
+ * @param particle_floating_wall_pairs_in_contact An unordered_map of maps which stores
+ * all the particle-floating wall pairs which are in contact, and
+ * the contact information in a struct. Note that the size of this unordered
+ * map is equal to the number of particles
+ */
+template <int dim>
+void
+particle_floating_wall_fine_search(
+  const typename DEM::dem_data_structures<
+    dim>::particle_floating_wall_candidates
+    &particle_floating_wall_contact_candidates,
+  const Parameters::Lagrangian::FloatingWalls<dim> &floating_wall_properties,
+  const double                                      simulation_time,
+  typename DEM::dem_data_structures<dim>::particle_wall_in_contact
+    &particle_floating_wall_in_contact);
 
-  /**
-   * @brief Iterate over the contact candidates from particle-floating mesh broad
-   * search (particle_floating_mesh_contact_candidates) to add new contact pairs
-   * to the particle_floating_mesh_in_contact container
-   *
-   * @param particle_floating_mesh_contact_candidates The output of particle-floating mesh
-   * broad search which shows contact pair candidates
-   * @param particle_floating_mesh_in_contact An map of maps which stores
-   * all the particle-floating mesh pairs which are in contact
-   */
-
-  void
-  particle_floating_mesh_fine_search(
-    const typename DEM::dem_data_structures<
-      dim>::particle_floating_mesh_candidates
-      &particle_floating_mesh_contact_candidates,
-    typename DEM::dem_data_structures<dim>::particle_floating_mesh_in_contact
-      &particle_floating_mesh_in_contact);
-};
+/**
+ * @brief Iterate over the contact candidates from particle-floating mesh broad
+ * search (particle_floating_mesh_contact_candidates) to add new contact pairs
+ * to the particle_floating_mesh_in_contact container
+ *
+ * @param particle_floating_mesh_contact_candidates The output of particle-floating mesh
+ * broad search which shows contact pair candidates
+ * @param particle_floating_mesh_in_contact An map of maps which stores
+ * all the particle-floating mesh pairs which are in contact
+ */
+template <int dim>
+void
+particle_floating_mesh_fine_search(
+  const typename DEM::dem_data_structures<
+    dim>::particle_floating_mesh_candidates
+    &particle_floating_mesh_contact_candidates,
+  typename DEM::dem_data_structures<dim>::particle_floating_mesh_in_contact
+    &particle_floating_mesh_in_contact);
 
 #endif
diff --git a/include/dem/particle_wall_jkr_force.h b/include/dem/particle_wall_jkr_force.h
index 20b6113154..1ef9c9c63f 100644
--- a/include/dem/particle_wall_jkr_force.h
+++ b/include/dem/particle_wall_jkr_force.h
@@ -21,7 +21,6 @@
 
 #include <dem/dem_solver_parameters.h>
 #include <dem/particle_wall_contact_force.h>
-#include <dem/particle_wall_contact_info.h>
 
 #include <deal.II/particles/particle.h>
 
diff --git a/include/dem/particle_wall_linear_force.h b/include/dem/particle_wall_linear_force.h
index 2ddcbe6269..1244eaed39 100644
--- a/include/dem/particle_wall_linear_force.h
+++ b/include/dem/particle_wall_linear_force.h
@@ -21,7 +21,6 @@
 
 #include <dem/dem_solver_parameters.h>
 #include <dem/particle_wall_contact_force.h>
-#include <dem/particle_wall_contact_info.h>
 
 #include <deal.II/particles/particle.h>
 
diff --git a/include/dem/particle_wall_nonlinear_force.h b/include/dem/particle_wall_nonlinear_force.h
index 5493367b4a..11909fd76b 100644
--- a/include/dem/particle_wall_nonlinear_force.h
+++ b/include/dem/particle_wall_nonlinear_force.h
@@ -21,7 +21,6 @@
 
 #include <dem/dem_solver_parameters.h>
 #include <dem/particle_wall_contact_force.h>
-#include <dem/particle_wall_contact_info.h>
 
 #include <deal.II/particles/particle.h>
 
diff --git a/include/dem/update_fine_search_candidates.h b/include/dem/update_fine_search_candidates.h
index c9a69fbc75..7d511297bd 100644
--- a/include/dem/update_fine_search_candidates.h
+++ b/include/dem/update_fine_search_candidates.h
@@ -19,8 +19,7 @@
 
 #include <dem/contact_type.h>
 #include <dem/data_containers.h>
-#include <dem/particle_particle_contact_info.h>
-#include <dem/particle_wall_contact_info.h>
+
 
 /**
  * @brief Manage removing repetitions and adding new contact pairs to the
diff --git a/include/dem/update_local_particle_containers.h b/include/dem/update_local_particle_containers.h
index 9d9ec5d54f..3407c09070 100644
--- a/include/dem/update_local_particle_containers.h
+++ b/include/dem/update_local_particle_containers.h
@@ -17,9 +17,9 @@
 #ifndef lethe_update_local_particle_containers_h
 #define lethe_update_local_particle_containers_h
 
+#include <dem/contact_info.h>
 #include <dem/contact_type.h>
 #include <dem/data_containers.h>
-#include <dem/particle_point_line_contact_info_struct.h>
 
 #include <deal.II/particles/particle_handler.h>
 
diff --git a/include/fem-dem/cfd_dem_coupling.h b/include/fem-dem/cfd_dem_coupling.h
index d0dd624dee..68dc370ed4 100644
--- a/include/fem-dem/cfd_dem_coupling.h
+++ b/include/fem-dem/cfd_dem_coupling.h
@@ -266,11 +266,10 @@ class CFDDEMSolver : public FluidDynamicsVANS<dim>
   PVDHandler grid_pvdhandler;
   PVDHandler particles_pvdhandler;
 
-  DEMSolverParameters<dim>      dem_parameters;
-  double                        dem_time_step;
-  const unsigned int            this_mpi_process;
-  const unsigned int            n_mpi_processes;
-  LagrangianPostProcessing<dim> dem_post_processing_object;
+  DEMSolverParameters<dim> dem_parameters;
+  double                   dem_time_step;
+  const unsigned int       this_mpi_process;
+  const unsigned int       n_mpi_processes;
 
   /// Post-processing variables to output total fluid volume and total particles
   /// volume
diff --git a/include/fem-dem/ib_particles_dem.h b/include/fem-dem/ib_particles_dem.h
index 03c04b80e2..9b29734bcc 100644
--- a/include/fem-dem/ib_particles_dem.h
+++ b/include/fem-dem/ib_particles_dem.h
@@ -15,23 +15,21 @@
  *
  */
 
+#ifndef lethe_ib_particles_dem_h
+#define lethe_ib_particles_dem_h
+
 #include <core/ib_particle.h>
 #include <core/ib_stencil.h>
 #include <core/lethe_grid_tools.h>
 
 #include <dem/particle_particle_contact_force.h>
-#include <dem/particle_particle_contact_info.h>
 #include <dem/particle_wall_contact_force.h>
-#include <dem/particle_wall_contact_info.h>
 
 #include <deal.II/base/tensor.h>
 
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_tools.h>
 
-#ifndef lethe_ib_particles_dem_h
-#  define lethe_ib_particles_dem_h
-
 using namespace dealii;
 
 /**
diff --git a/source/dem/CMakeLists.txt b/source/dem/CMakeLists.txt
index a83a37072b..5a0e9b33e4 100644
--- a/source/dem/CMakeLists.txt
+++ b/source/dem/CMakeLists.txt
@@ -50,6 +50,7 @@ add_library(lethe-dem
   # Headers
   ../../include/dem/adaptive_sparse_contacts.h
   ../../include/dem/boundary_cells_info_struct.h
+  ../../include/dem/contact_info.h
   ../../include/dem/contact_type.h
   ../../include/dem/data_containers.h
   ../../include/dem/dem.h
@@ -76,15 +77,12 @@ add_library(lethe-dem
   ../../include/dem/output_force_torque_calculation.h
   ../../include/dem/particle_particle_broad_search.h
   ../../include/dem/particle_particle_contact_force.h
-  ../../include/dem/particle_particle_contact_info.h
   ../../include/dem/particle_particle_fine_search.h
   ../../include/dem/particle_point_line_broad_search.h
   ../../include/dem/particle_point_line_contact_force.h
-  ../../include/dem/particle_point_line_contact_info_struct.h
   ../../include/dem/particle_point_line_fine_search.h
   ../../include/dem/particle_wall_broad_search.h
   ../../include/dem/particle_wall_contact_force.h
-  ../../include/dem/particle_wall_contact_info.h
   ../../include/dem/particle_wall_dmt_force.h
   ../../include/dem/particle_wall_fine_search.h
   ../../include/dem/particle_wall_jkr_force.h
diff --git a/source/dem/dem.cc b/source/dem/dem.cc
index 86c377c1f2..326db22513 100644
--- a/source/dem/dem.cc
+++ b/source/dem/dem.cc
@@ -55,8 +55,8 @@ DEMSolver<dim>::DEMSolver(DEMSolverParameters<dim> dem_parameters)
 template <int dim>
 void
 DEMSolver<dim>::setup_parameters()
-{ // Print simulation starting information
-
+{
+  // Print simulation starting information
   pcout << std::endl;
   std::stringstream ss;
   ss << "Running on " << n_mpi_processes << " rank(s)";
@@ -620,7 +620,7 @@ DEMSolver<dim>::finish_simulation()
   // Outputting force and torques over boundary
   if (parameters.forces_torques.calculate_force_torque)
     {
-      write_forces_torques_output_results<dim>(
+      write_forces_torques_output_results(
         parameters.forces_torques.force_torque_output_name,
         parameters.forces_torques.output_frequency,
         triangulation.get_boundary_ids(),
@@ -703,16 +703,15 @@ DEMSolver<dim>::post_process_results()
   if (parameters.post_processing.Lagrangian_post_processing &&
       simulation_control->is_output_iteration())
     {
-      post_processing_object.write_post_processing_results(
-        triangulation,
-        grid_pvdhandler,
-        background_dh,
-        particle_handler,
-        parameters,
-        simulation_control->get_current_time(),
-        simulation_control->get_step_number(),
-        mpi_communicator,
-        sparse_contacts_object);
+      write_post_processing_results<dim>(triangulation,
+                                         grid_pvdhandler,
+                                         background_dh,
+                                         particle_handler,
+                                         parameters,
+                                         simulation_control->get_current_time(),
+                                         simulation_control->get_step_number(),
+                                         mpi_communicator,
+                                         sparse_contacts_object);
     }
 }
 
@@ -823,7 +822,6 @@ DEMSolver<dim>::sort_particles_into_subdomains_and_cells()
   particle_handler.exchange_ghost_particles(true);
 }
 
-
 template <int dim>
 void
 DEMSolver<dim>::solve()
diff --git a/source/dem/dem_contact_manager.cc b/source/dem/dem_contact_manager.cc
index 1ed19303c6..bc14a9264d 100644
--- a/source/dem/dem_contact_manager.cc
+++ b/source/dem/dem_contact_manager.cc
@@ -1,25 +1,27 @@
 #include <dem/dem_action_manager.h>
 #include <dem/dem_contact_manager.h>
 
+using namespace DEM;
+
 template <int dim>
 void
 DEMContactManager<dim>::execute_cell_neighbors_search(
   const parallel::distributed::Triangulation<dim> &triangulation,
-  const typename DEM::dem_data_structures<dim>::periodic_boundaries_cells_info
+  const typename dem_data_structures<dim>::periodic_boundaries_cells_info
     periodic_boundaries_cells_information)
 {
   // Get action manager
   auto action_manager = DEMActionManager::get_action_manager();
 
   // Find cell neighbors
-  cell_neighbors_object.find_cell_neighbors(triangulation,
-                                            cells_local_neighbor_list,
-                                            cells_ghost_neighbor_list);
+  find_cell_neighbors<dim>(triangulation,
+                           cells_local_neighbor_list,
+                           cells_ghost_neighbor_list);
 
   // Find cell periodic neighbors
   if (action_manager->check_periodic_boundaries_enabled())
     {
-      cell_neighbors_object.find_cell_periodic_neighbors(
+      find_cell_periodic_neighbors<dim>(
         triangulation,
         periodic_boundaries_cells_information,
         cells_local_periodic_neighbor_list,
@@ -30,8 +32,7 @@ DEMContactManager<dim>::execute_cell_neighbors_search(
   // Get total (with repetition) neighbors list for floating mesh.
   if (action_manager->check_solid_objects_enabled())
     {
-      cell_neighbors_object.find_full_cell_neighbors(triangulation,
-                                                     total_neighbor_list);
+      find_full_cell_neighbors<dim>(triangulation, total_neighbor_list);
     }
 }
 
@@ -43,8 +44,8 @@ DEMContactManager<dim>::update_contacts()
   // search step with local_contact_pair_candidates
   update_fine_search_candidates<
     dim,
-    typename DEM::dem_data_structures<dim>::adjacent_particle_pairs,
-    typename DEM::dem_data_structures<dim>::particle_particle_candidates,
+    typename dem_data_structures<dim>::adjacent_particle_pairs,
+    typename dem_data_structures<dim>::particle_particle_candidates,
     ContactType::local_particle_particle>(local_adjacent_particles,
                                           local_contact_pair_candidates);
 
@@ -52,8 +53,8 @@ DEMContactManager<dim>::update_contacts()
   // search step with global_contact_pair_candidates
   update_fine_search_candidates<
     dim,
-    typename DEM::dem_data_structures<dim>::adjacent_particle_pairs,
-    typename DEM::dem_data_structures<dim>::particle_particle_candidates,
+    typename dem_data_structures<dim>::adjacent_particle_pairs,
+    typename dem_data_structures<dim>::particle_particle_candidates,
     ContactType::ghost_particle_particle>(ghost_adjacent_particles,
                                           ghost_contact_pair_candidates);
 
@@ -65,8 +66,8 @@ DEMContactManager<dim>::update_contacts()
       // local_contact_pair_periodic_candidates
       update_fine_search_candidates<
         dim,
-        typename DEM::dem_data_structures<dim>::adjacent_particle_pairs,
-        typename DEM::dem_data_structures<dim>::particle_particle_candidates,
+        typename dem_data_structures<dim>::adjacent_particle_pairs,
+        typename dem_data_structures<dim>::particle_particle_candidates,
         ContactType::local_periodic_particle_particle>(
         local_periodic_adjacent_particles,
         local_contact_pair_periodic_candidates);
@@ -76,8 +77,8 @@ DEMContactManager<dim>::update_contacts()
       // ghost_contact_pair_periodic_candidates
       update_fine_search_candidates<
         dim,
-        typename DEM::dem_data_structures<dim>::adjacent_particle_pairs,
-        typename DEM::dem_data_structures<dim>::particle_particle_candidates,
+        typename dem_data_structures<dim>::adjacent_particle_pairs,
+        typename dem_data_structures<dim>::particle_particle_candidates,
         ContactType::ghost_periodic_particle_particle>(
         ghost_periodic_adjacent_particles,
         ghost_contact_pair_periodic_candidates);
@@ -87,8 +88,8 @@ DEMContactManager<dim>::update_contacts()
       // ghost_local_contact_pair_periodic_candidates
       update_fine_search_candidates<
         dim,
-        typename DEM::dem_data_structures<dim>::adjacent_particle_pairs,
-        typename DEM::dem_data_structures<dim>::particle_particle_candidates,
+        typename dem_data_structures<dim>::adjacent_particle_pairs,
+        typename dem_data_structures<dim>::particle_particle_candidates,
         ContactType::ghost_local_periodic_particle_particle>(
         ghost_local_periodic_adjacent_particles,
         ghost_local_contact_pair_periodic_candidates);
@@ -98,8 +99,8 @@ DEMContactManager<dim>::update_contacts()
   // search step with particle_wall_contact_candidates
   update_fine_search_candidates<
     dim,
-    typename DEM::dem_data_structures<dim>::particle_wall_in_contact,
-    typename DEM::dem_data_structures<dim>::particle_wall_candidates,
+    typename dem_data_structures<dim>::particle_wall_in_contact,
+    typename dem_data_structures<dim>::particle_wall_candidates,
     ContactType::particle_wall>(particle_wall_in_contact,
                                 particle_wall_candidates);
 
@@ -107,8 +108,8 @@ DEMContactManager<dim>::update_contacts()
   // of fine search step with particle_floating_wall_contact_candidates
   update_fine_search_candidates<
     dim,
-    typename DEM::dem_data_structures<dim>::particle_wall_in_contact,
-    typename DEM::dem_data_structures<dim>::particle_floating_wall_candidates,
+    typename dem_data_structures<dim>::particle_wall_in_contact,
+    typename dem_data_structures<dim>::particle_floating_wall_candidates,
     ContactType::particle_floating_wall>(particle_floating_wall_in_contact,
                                          particle_floating_wall_candidates);
 
@@ -120,9 +121,9 @@ DEMContactManager<dim>::update_contacts()
     {
       update_fine_search_candidates<
         dim,
-        typename DEM::dem_data_structures<
+        typename dem_data_structures<
           dim>::particle_floating_wall_from_mesh_in_contact,
-        typename DEM::dem_data_structures<
+        typename dem_data_structures<
           dim>::particle_floating_wall_from_mesh_candidates,
         ContactType::particle_floating_mesh>(
         particle_floating_mesh_in_contact[solid_counter],
@@ -141,14 +142,14 @@ DEMContactManager<dim>::update_local_particles_in_cells(
   // Update contact containers for local particle-particle pairs in contact
   update_contact_container_iterators<
     dim,
-    typename DEM::dem_data_structures<dim>::adjacent_particle_pairs,
+    typename dem_data_structures<dim>::adjacent_particle_pairs,
     ContactType::local_particle_particle>(local_adjacent_particles,
                                           particle_container);
 
   // Update contact containers for ghost particle-particle pairs in contact
   update_contact_container_iterators<
     dim,
-    typename DEM::dem_data_structures<dim>::adjacent_particle_pairs,
+    typename dem_data_structures<dim>::adjacent_particle_pairs,
     ContactType::ghost_particle_particle>(ghost_adjacent_particles,
                                           particle_container);
 
@@ -159,7 +160,7 @@ DEMContactManager<dim>::update_local_particles_in_cells(
       // pairs in contact
       update_contact_container_iterators<
         dim,
-        typename DEM::dem_data_structures<dim>::adjacent_particle_pairs,
+        typename dem_data_structures<dim>::adjacent_particle_pairs,
         ContactType::local_periodic_particle_particle>(
         local_periodic_adjacent_particles, particle_container);
 
@@ -167,7 +168,7 @@ DEMContactManager<dim>::update_local_particles_in_cells(
       // pairs in contact
       update_contact_container_iterators<
         dim,
-        typename DEM::dem_data_structures<dim>::adjacent_particle_pairs,
+        typename dem_data_structures<dim>::adjacent_particle_pairs,
         ContactType::ghost_periodic_particle_particle>(
         ghost_periodic_adjacent_particles, particle_container);
 
@@ -175,7 +176,7 @@ DEMContactManager<dim>::update_local_particles_in_cells(
       // pairs in contact
       update_contact_container_iterators<
         dim,
-        typename DEM::dem_data_structures<dim>::adjacent_particle_pairs,
+        typename dem_data_structures<dim>::adjacent_particle_pairs,
         ContactType::ghost_local_periodic_particle_particle>(
         ghost_local_periodic_adjacent_particles, particle_container);
     }
@@ -183,13 +184,13 @@ DEMContactManager<dim>::update_local_particles_in_cells(
   // Update contact containers for particle-wall pairs in contact
   update_contact_container_iterators<
     dim,
-    typename DEM::dem_data_structures<dim>::particle_wall_in_contact,
+    typename dem_data_structures<dim>::particle_wall_in_contact,
     ContactType::particle_wall>(particle_wall_in_contact, particle_container);
 
   // Update contact containers for particle-floating wall pairs in contact
   update_contact_container_iterators<
     dim,
-    typename DEM::dem_data_structures<dim>::particle_wall_in_contact,
+    typename dem_data_structures<dim>::particle_wall_in_contact,
     ContactType::particle_floating_wall>(particle_floating_wall_in_contact,
                                          particle_container);
 
@@ -201,7 +202,7 @@ DEMContactManager<dim>::update_local_particles_in_cells(
     {
       update_contact_container_iterators<
         dim,
-        typename DEM::dem_data_structures<
+        typename dem_data_structures<
           dim>::particle_floating_wall_from_mesh_in_contact,
         ContactType::particle_floating_mesh>(
         particle_floating_mesh_in_contact[solid_counter], particle_container);
@@ -210,13 +211,13 @@ DEMContactManager<dim>::update_local_particles_in_cells(
   // Update contact containers for particle-line pairs in contact
   update_contact_container_iterators<
     dim,
-    typename DEM::dem_data_structures<dim>::particle_point_line_contact_info,
-    ContactType::particle_point>(particle_lines_in_contact, particle_container);
+    typename dem_data_structures<dim>::particle_line_in_contact,
+    ContactType::particle_line>(particle_lines_in_contact, particle_container);
 
   // Update contact containers for particle-point pairs in contact
   update_contact_container_iterators<
     dim,
-    typename DEM::dem_data_structures<dim>::particle_point_line_contact_info,
+    typename dem_data_structures<dim>::particle_point_in_contact,
     ContactType::particle_point>(particle_points_in_contact,
                                  particle_container);
 }
@@ -233,28 +234,24 @@ DEMContactManager<dim>::execute_particle_particle_broad_search(
   // The first broad search is the default one for sparse contacts
   if (action_manager->use_default_broad_search_functions())
     {
-      particle_particle_broad_search_object
-        .find_particle_particle_contact_pairs(particle_handler, *this);
+      find_particle_particle_contact_pairs<dim>(particle_handler, *this);
 
       if (action_manager->check_periodic_boundaries_enabled())
         {
-          particle_particle_broad_search_object
-            .find_particle_particle_periodic_contact_pairs(particle_handler,
-                                                           *this);
+          find_particle_particle_periodic_contact_pairs<dim>(particle_handler,
+                                                             *this);
         }
     }
   else
     {
-      particle_particle_broad_search_object
-        .find_particle_particle_contact_pairs(particle_handler,
-                                              *this,
-                                              sparse_contacts_object);
+      find_particle_particle_contact_pairs<dim>(particle_handler,
+                                                *this,
+                                                sparse_contacts_object);
 
       if (action_manager->check_periodic_boundaries_enabled())
         {
-          particle_particle_broad_search_object
-            .find_particle_particle_periodic_contact_pairs(
-              particle_handler, *this, sparse_contacts_object);
+          find_particle_particle_periodic_contact_pairs<dim>(
+            particle_handler, *this, sparse_contacts_object);
         }
     }
 }
@@ -275,7 +272,7 @@ DEMContactManager<dim>::execute_particle_wall_broad_search(
   if (action_manager->use_default_broad_search_functions())
     {
       // Particle-wall contact candidates
-      particle_wall_broad_search_object.find_particle_wall_contact_pairs(
+      find_particle_wall_contact_pairs<dim>(
         boundary_cell_object.get_boundary_cells_information(),
         particle_handler,
         particle_wall_candidates);
@@ -283,45 +280,41 @@ DEMContactManager<dim>::execute_particle_wall_broad_search(
       // Particle-floating wall contact pairs
       if (floating_walls.floating_walls_number > 0)
         {
-          particle_wall_broad_search_object
-            .find_particle_floating_wall_contact_pairs(
-              boundary_cell_object.get_boundary_cells_with_floating_walls(),
-              particle_handler,
-              floating_walls,
-              simulation_time,
-              particle_floating_wall_candidates);
+          find_particle_floating_wall_contact_pairs<dim>(
+            boundary_cell_object.get_boundary_cells_with_floating_walls(),
+            particle_handler,
+            floating_walls,
+            simulation_time,
+            particle_floating_wall_candidates);
         }
 
       // Particle-floating mesh broad search
       if (action_manager->check_solid_objects_enabled())
         {
-          particle_wall_broad_search_object
-            .particle_solid_surfaces_contact_search(
-              solid_surfaces_mesh_info,
-              particle_handler,
-              particle_floating_mesh_candidates,
-              total_neighbor_list);
+          particle_solid_surfaces_contact_search<dim>(
+            solid_surfaces_mesh_info,
+            particle_handler,
+            particle_floating_mesh_candidates,
+            total_neighbor_list);
         }
 
-      particle_point_candidates =
-        particle_point_line_broad_search_object
-          .find_particle_point_contact_pairs(
-            particle_handler,
-            boundary_cell_object.get_boundary_cells_with_points());
+      find_particle_point_contact_pairs<dim>(
+        particle_handler,
+        boundary_cell_object.get_boundary_cells_with_points(),
+        particle_point_candidates);
 
       if constexpr (dim == 3)
         {
-          particle_line_candidates =
-            particle_point_line_broad_search_object
-              .find_particle_line_contact_pairs(
-                particle_handler,
-                boundary_cell_object.get_boundary_cells_with_lines());
+          find_particle_line_contact_pairs<dim>(
+            particle_handler,
+            boundary_cell_object.get_boundary_cells_with_lines(),
+            particle_line_candidates);
         }
     }
   else
     {
       // Particle-wall contact candidates
-      particle_wall_broad_search_object.find_particle_wall_contact_pairs(
+      find_particle_wall_contact_pairs<dim>(
         boundary_cell_object.get_boundary_cells_information(),
         particle_handler,
         particle_wall_candidates,
@@ -330,43 +323,39 @@ DEMContactManager<dim>::execute_particle_wall_broad_search(
       // Particle-floating wall contact pairs
       if (floating_walls.floating_walls_number > 0)
         {
-          particle_wall_broad_search_object
-            .find_particle_floating_wall_contact_pairs(
-              boundary_cell_object.get_boundary_cells_with_floating_walls(),
-              particle_handler,
-              floating_walls,
-              simulation_time,
-              particle_floating_wall_candidates,
-              sparse_contacts_object);
+          find_particle_floating_wall_contact_pairs<dim>(
+            boundary_cell_object.get_boundary_cells_with_floating_walls(),
+            particle_handler,
+            floating_walls,
+            simulation_time,
+            particle_floating_wall_candidates,
+            sparse_contacts_object);
         }
 
       // Particle-floating mesh broad search
       if (action_manager->check_solid_objects_enabled())
         {
-          particle_wall_broad_search_object
-            .particle_solid_surfaces_contact_search(
-              solid_surfaces_mesh_info,
-              particle_handler,
-              particle_floating_mesh_candidates,
-              total_neighbor_list,
-              sparse_contacts_object);
-        }
-
-      particle_point_candidates =
-        particle_point_line_broad_search_object
-          .find_particle_point_contact_pairs(
+          particle_solid_surfaces_contact_search<dim>(
+            solid_surfaces_mesh_info,
             particle_handler,
-            boundary_cell_object.get_boundary_cells_with_points(),
+            particle_floating_mesh_candidates,
+            total_neighbor_list,
             sparse_contacts_object);
+        }
+
+      find_particle_point_contact_pairs<dim>(
+        particle_handler,
+        boundary_cell_object.get_boundary_cells_with_points(),
+        particle_point_candidates,
+        sparse_contacts_object);
 
       if constexpr (dim == 3)
         {
-          particle_line_candidates =
-            particle_point_line_broad_search_object
-              .find_particle_line_contact_pairs(
-                particle_handler,
-                boundary_cell_object.get_boundary_cells_with_lines(),
-                sparse_contacts_object);
+          find_particle_line_contact_pairs<dim>(
+            particle_handler,
+            boundary_cell_object.get_boundary_cells_with_lines(),
+            particle_line_candidates,
+            sparse_contacts_object);
         }
     }
 }
@@ -378,40 +367,36 @@ DEMContactManager<dim>::execute_particle_particle_fine_search(
   const Tensor<1, dim> periodic_offset)
 {
   // Fine search for local particle-particle
-  particle_particle_fine_search_object.particle_particle_fine_search(
-    particle_container,
-    local_adjacent_particles,
-    local_contact_pair_candidates,
-    neighborhood_threshold);
+  particle_particle_fine_search<dim>(particle_container,
+                                     local_adjacent_particles,
+                                     local_contact_pair_candidates,
+                                     neighborhood_threshold);
 
   // Fine search for ghost particle-particle
-  particle_particle_fine_search_object.particle_particle_fine_search(
-    particle_container,
-    ghost_adjacent_particles,
-    ghost_contact_pair_candidates,
-    neighborhood_threshold);
+  particle_particle_fine_search<dim>(particle_container,
+                                     ghost_adjacent_particles,
+                                     ghost_contact_pair_candidates,
+                                     neighborhood_threshold);
 
   if (DEMActionManager::get_action_manager()
         ->check_periodic_boundaries_enabled())
     {
       // Fine search for local-local periodic particle-particle
-      particle_particle_fine_search_object.particle_particle_fine_search(
-        particle_container,
-        local_periodic_adjacent_particles,
-        local_contact_pair_periodic_candidates,
-        neighborhood_threshold,
-        periodic_offset);
+      particle_particle_fine_search<dim>(particle_container,
+                                         local_periodic_adjacent_particles,
+                                         local_contact_pair_periodic_candidates,
+                                         neighborhood_threshold,
+                                         periodic_offset);
 
       // Fine search for local-ghost periodic particle-particle
-      particle_particle_fine_search_object.particle_particle_fine_search(
-        particle_container,
-        ghost_periodic_adjacent_particles,
-        ghost_contact_pair_periodic_candidates,
-        neighborhood_threshold,
-        periodic_offset);
+      particle_particle_fine_search<dim>(particle_container,
+                                         ghost_periodic_adjacent_particles,
+                                         ghost_contact_pair_periodic_candidates,
+                                         neighborhood_threshold,
+                                         periodic_offset);
 
       // Fine search for ghost-local periodic particle-particle
-      particle_particle_fine_search_object.particle_particle_fine_search(
+      particle_particle_fine_search<dim>(
         particle_container,
         ghost_local_periodic_adjacent_particles,
         ghost_local_contact_pair_periodic_candidates,
@@ -428,13 +413,13 @@ DEMContactManager<dim>::execute_particle_wall_fine_search(
   const double                                      neighborhood_threshold)
 {
   // Particle - wall fine search
-  particle_wall_fine_search_object.particle_wall_fine_search(
-    particle_wall_candidates, particle_wall_in_contact);
+  particle_wall_fine_search<dim>(particle_wall_candidates,
+                                 particle_wall_in_contact);
 
   // Particle - floating wall fine search
   if (floating_walls.floating_walls_number > 0)
     {
-      particle_wall_fine_search_object.particle_floating_wall_fine_search(
+      particle_floating_wall_fine_search<dim>(
         particle_floating_wall_candidates,
         floating_walls,
         simulation_time,
@@ -444,19 +429,19 @@ DEMContactManager<dim>::execute_particle_wall_fine_search(
   // Particle - floating mesh fine search
   if (DEMActionManager::get_action_manager()->check_solid_objects_enabled())
     {
-      particle_wall_fine_search_object.particle_floating_mesh_fine_search(
+      particle_floating_mesh_fine_search<dim>(
         particle_floating_mesh_candidates, particle_floating_mesh_in_contact);
     }
 
-  particle_points_in_contact =
-    particle_point_line_fine_search_object.particle_point_fine_search(
-      particle_point_candidates, neighborhood_threshold);
+  particle_point_fine_search<dim>(particle_point_candidates,
+                                  neighborhood_threshold,
+                                  particle_points_in_contact);
 
   if constexpr (dim == 3)
     {
-      particle_lines_in_contact =
-        particle_point_line_fine_search_object.particle_line_fine_search(
-          particle_line_candidates, neighborhood_threshold);
+      particle_line_fine_search<dim>(particle_line_candidates,
+                                     neighborhood_threshold,
+                                     particle_lines_in_contact);
     }
 }
 
diff --git a/source/dem/find_boundary_cells_information.cc b/source/dem/find_boundary_cells_information.cc
index 4ef588aee6..4786f43c25 100644
--- a/source/dem/find_boundary_cells_information.cc
+++ b/source/dem/find_boundary_cells_information.cc
@@ -323,10 +323,10 @@ BoundaryCellsInformation<dim>::find_particle_point_and_line_contact_cells(
                             {
                               all_cells_with_boundary_lines[cell->id()
                                                               .to_string()]
-                                .insert(
-                                  {face->line_index(l),
-                                   std::make_pair(face->line(l)->vertex(0),
-                                                  face->line(l)->vertex(1))});
+                                .emplace(
+                                  face->line_index(l),
+                                  std::make_pair(face->line(l)->vertex(0),
+                                                 face->line(l)->vertex(1)));
                             }
                         }
                     }
@@ -378,11 +378,24 @@ BoundaryCellsInformation<dim>::find_particle_point_and_line_contact_cells(
                     {
                       for (auto &map_iterator : cell_boundary_lines)
                         {
-                          boundary_cells_with_lines.insert(
-                            {cell_id_string,
-                             std::make_tuple(cell,
-                                             map_iterator.second.first,
-                                             map_iterator.second.second)});
+                          Point<3> vertex_one, vertex_two;
+
+                          if constexpr (dim == 3)
+                            {
+                              vertex_one = map_iterator.second.first;
+                              vertex_two = map_iterator.second.second;
+                            }
+
+                          if constexpr (dim == 2)
+                            {
+                              vertex_one =
+                                point_nd_to_3d(map_iterator.second.first);
+                              vertex_two =
+                                point_nd_to_3d(map_iterator.second.second);
+                            }
+                          boundary_cells_with_lines.emplace(
+                            cell_id_string,
+                            cell_line_info<dim>{cell, vertex_one, vertex_two});
 
                           // Add the cell to local_cells_with_boundary_lines to
                           // be used in
@@ -428,9 +441,17 @@ BoundaryCellsInformation<dim>::find_particle_point_and_line_contact_cells(
                     {
                       if (boundary_vertices.count(face->vertex_index(v)) > 0)
                         {
-                          boundary_cells_with_points.insert(
-                            {cell->id().to_string(),
-                             std::make_pair(cell, face->vertex(v))});
+                          Point<3> vertex;
+
+                          if constexpr (dim == 3)
+                            vertex = face->vertex(v);
+
+                          if constexpr (dim == 2)
+                            vertex = point_nd_to_3d(face->vertex(v));
+
+                          boundary_cells_with_points.emplace(
+                            cell->id().to_string(),
+                            cell_point_info<dim>{cell, vertex});
                         }
                     }
                 }
diff --git a/source/dem/find_cell_neighbors.cc b/source/dem/find_cell_neighbors.cc
index b646083e84..0fcff7656c 100644
--- a/source/dem/find_cell_neighbors.cc
+++ b/source/dem/find_cell_neighbors.cc
@@ -1,21 +1,14 @@
 #include <dem/find_cell_neighbors.h>
 
-using namespace dealii;
+using namespace DEM;
 
-// The constructor of this class is empty
-template <int dim>
-FindCellNeighbors<dim>::FindCellNeighbors()
-{}
-
-// This function finds the neighbor list (without repetition) of all the active
-// cells in the triangulation
 template <int dim>
 void
-FindCellNeighbors<dim>::find_cell_neighbors(
+find_cell_neighbors(
   const parallel::distributed::Triangulation<dim> &triangulation,
-  typename DEM::dem_data_structures<dim>::cells_neighbor_list
+  typename dem_data_structures<dim>::cells_neighbor_list
     &cells_local_neighbor_list,
-  typename DEM::dem_data_structures<dim>::cells_neighbor_list
+  typename dem_data_structures<dim>::cells_neighbor_list
     &cells_ghost_neighbor_list)
 {
   // The output vectors of the function are cells_local_neighbor_list and
@@ -23,14 +16,14 @@ FindCellNeighbors<dim>::find_cell_neighbors(
   // of all local cells; while the second contains all the ghost cells of all
   // local cells. They are two vectors with size of the number of active cells.
   // The first elements of all vectors are the main cells
-  typename DEM::dem_data_structures<dim>::cell_vector local_neighbor_vector;
-  typename DEM::dem_data_structures<dim>::cell_vector ghost_neighbor_vector;
+  typename dem_data_structures<dim>::cell_vector local_neighbor_vector;
+  typename dem_data_structures<dim>::cell_vector ghost_neighbor_vector;
 
   // This vector is used to avoid repetition of adjacent cells. For instance if
   // cell B is recognized as the neighbor of cell A, cell A will not be added to
   // the neighbor list of cell B again. This is done using the total_cell_list
   // vector
-  typename DEM::dem_data_structures<dim>::cell_set total_cell_list;
+  typename dem_data_structures<dim>::cell_set total_cell_list;
 
   // For each cell, the cell vertices are found and used to find adjacent cells.
   // The reason is to find the cells located on the corners of the main cell.
@@ -104,26 +97,24 @@ FindCellNeighbors<dim>::find_cell_neighbors(
 
 template <int dim>
 void
-FindCellNeighbors<dim>::find_cell_periodic_neighbors(
+find_cell_periodic_neighbors(
   const parallel::distributed::Triangulation<dim> &triangulation,
-  const typename DEM::dem_data_structures<dim>::periodic_boundaries_cells_info
+  const typename dem_data_structures<dim>::periodic_boundaries_cells_info
     &periodic_boundaries_cells_information,
-  typename DEM::dem_data_structures<dim>::cells_neighbor_list
+  typename dem_data_structures<dim>::cells_neighbor_list
     &cells_local_periodic_neighbor_list,
-  typename DEM::dem_data_structures<dim>::cells_neighbor_list
+  typename dem_data_structures<dim>::cells_neighbor_list
     &cells_ghost_periodic_neighbor_list,
-  typename DEM::dem_data_structures<dim>::cells_neighbor_list
+  typename dem_data_structures<dim>::cells_neighbor_list
     &cells_ghost_local_periodic_neighbor_list)
 {
-  typename DEM::dem_data_structures<dim>::cell_vector
-    local_periodic_neighbor_vector;
-  typename DEM::dem_data_structures<dim>::cell_vector
-    ghost_periodic_neighbor_vector;
-  typename DEM::dem_data_structures<dim>::cell_vector
+  typename dem_data_structures<dim>::cell_vector local_periodic_neighbor_vector;
+  typename dem_data_structures<dim>::cell_vector ghost_periodic_neighbor_vector;
+  typename dem_data_structures<dim>::cell_vector
     ghost_local_periodic_neighbor_vector;
 
-  typename DEM::dem_data_structures<dim>::cell_set total_cell_list;
-  typename DEM::dem_data_structures<dim>::cell_set total_ghost_cell_list;
+  typename dem_data_structures<dim>::cell_set total_cell_list;
+  typename dem_data_structures<dim>::cell_set total_ghost_cell_list;
 
   // For each cell, the cell vertices are found and used to find adjacent cells.
   // The reason is to find the cells located on the corners of the main cell.
@@ -154,15 +145,14 @@ FindCellNeighbors<dim>::find_cell_periodic_neighbors(
           total_cell_list.insert(cell);
 
           // Empty list of periodic cell neighbor
-          typename DEM::dem_data_structures<dim>::cell_vector
-            periodic_neighbor_list;
+          typename dem_data_structures<dim>::cell_vector periodic_neighbor_list;
 
           // Get the periodic neighbor of the cell
-          get_periodic_neighbor_list(cell,
-                                     coinciding_vertex_groups,
-                                     vertex_to_coinciding_vertex_group,
-                                     v_to_c,
-                                     periodic_neighbor_list);
+          get_periodic_neighbor_list<dim>(cell,
+                                          coinciding_vertex_groups,
+                                          vertex_to_coinciding_vertex_group,
+                                          v_to_c,
+                                          periodic_neighbor_list);
 
           for (const auto &periodic_neighbor : periodic_neighbor_list)
             {
@@ -229,15 +219,14 @@ FindCellNeighbors<dim>::find_cell_periodic_neighbors(
           total_ghost_cell_list.insert(cell);
 
           // Empty list of periodic cell neighbor
-          typename DEM::dem_data_structures<dim>::cell_vector
-            periodic_neighbor_list;
+          typename dem_data_structures<dim>::cell_vector periodic_neighbor_list;
 
           // Get the periodic neighbor of the cell
-          get_periodic_neighbor_list(cell,
-                                     coinciding_vertex_groups,
-                                     vertex_to_coinciding_vertex_group,
-                                     v_to_c,
-                                     periodic_neighbor_list);
+          get_periodic_neighbor_list<dim>(cell,
+                                          coinciding_vertex_groups,
+                                          vertex_to_coinciding_vertex_group,
+                                          v_to_c,
+                                          periodic_neighbor_list);
 
           for (const auto &periodic_neighbor : periodic_neighbor_list)
             {
@@ -274,9 +263,9 @@ FindCellNeighbors<dim>::find_cell_periodic_neighbors(
 // the main cell to search for possible collisions with the floating mesh.
 template <int dim>
 void
-FindCellNeighbors<dim>::find_full_cell_neighbors(
+find_full_cell_neighbors(
   const parallel::distributed::Triangulation<dim> &triangulation,
-  typename DEM::dem_data_structures<dim>::cells_total_neighbor_list
+  typename dem_data_structures<dim>::cells_total_neighbor_list
     &cells_total_neighbor_list)
 {
   auto v_to_c = GridTools::vertex_to_cell_map(triangulation);
@@ -316,16 +305,17 @@ FindCellNeighbors<dim>::find_full_cell_neighbors(
         }
     }
 }
+
 template <int dim>
 void
-FindCellNeighbors<dim>::get_periodic_neighbor_list(
+get_periodic_neighbor_list(
   const typename Triangulation<dim>::active_cell_iterator &cell,
   const std::map<unsigned int, std::vector<unsigned int>>
                                              &coinciding_vertex_groups,
   const std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group,
   const std::vector<std::set<typename Triangulation<dim>::active_cell_iterator>>
-                                                      &v_to_c,
-  typename DEM::dem_data_structures<dim>::cell_vector &periodic_neighbor_list)
+                                                 &v_to_c,
+  typename dem_data_structures<dim>::cell_vector &periodic_neighbor_list)
 {
   // Loop over vertices of the cell
   for (unsigned int vertex = 0; vertex < cell->n_vertices(); ++vertex)
@@ -360,5 +350,54 @@ FindCellNeighbors<dim>::get_periodic_neighbor_list(
     }
 }
 
-template class FindCellNeighbors<2>;
-template class FindCellNeighbors<3>;
+template void
+find_cell_neighbors<2>(
+  const parallel::distributed::Triangulation<2> &triangulation,
+  typename dem_data_structures<2>::cells_neighbor_list
+    &cells_local_neighbor_list,
+  typename dem_data_structures<2>::cells_neighbor_list
+    &cells_ghost_neighbor_list);
+
+template void
+find_cell_neighbors<3>(
+  const parallel::distributed::Triangulation<3> &triangulation,
+  typename dem_data_structures<3>::cells_neighbor_list
+    &cells_local_neighbor_list,
+  typename dem_data_structures<3>::cells_neighbor_list
+    &cells_ghost_neighbor_list);
+
+template void
+find_cell_periodic_neighbors<2>(
+  const parallel::distributed::Triangulation<2> &triangulation,
+  const typename dem_data_structures<2>::periodic_boundaries_cells_info
+    &periodic_boundaries_cells_information,
+  typename dem_data_structures<2>::cells_neighbor_list
+    &cells_local_periodic_neighbor_list,
+  typename dem_data_structures<2>::cells_neighbor_list
+    &cells_ghost_periodic_neighbor_list,
+  typename dem_data_structures<2>::cells_neighbor_list
+    &cells_ghost_local_periodic_neighbor_list);
+
+template void
+find_cell_periodic_neighbors<3>(
+  const parallel::distributed::Triangulation<3> &triangulation,
+  const typename dem_data_structures<3>::periodic_boundaries_cells_info
+    &periodic_boundaries_cells_information,
+  typename dem_data_structures<3>::cells_neighbor_list
+    &cells_local_periodic_neighbor_list,
+  typename dem_data_structures<3>::cells_neighbor_list
+    &cells_ghost_periodic_neighbor_list,
+  typename dem_data_structures<3>::cells_neighbor_list
+    &cells_ghost_local_periodic_neighbor_list);
+
+template void
+find_full_cell_neighbors<2>(
+  const parallel::distributed::Triangulation<2> &triangulation,
+  typename dem_data_structures<2>::cells_total_neighbor_list
+    &cells_total_neighbor_list);
+
+template void
+find_full_cell_neighbors<3>(
+  const parallel::distributed::Triangulation<3> &triangulation,
+  typename dem_data_structures<3>::cells_total_neighbor_list
+    &cells_total_neighbor_list);
diff --git a/source/dem/find_contact_detection_step.cc b/source/dem/find_contact_detection_step.cc
index ad7f110306..50f34b3977 100644
--- a/source/dem/find_contact_detection_step.cc
+++ b/source/dem/find_contact_detection_step.cc
@@ -104,11 +104,11 @@ find_floating_mesh_mapping_step(
 }
 
 template void
-find_floating_mesh_mapping_step(
+find_floating_mesh_mapping_step<2>(
   const double smallest_contact_search_criterion,
-  std::vector<std::shared_ptr<SerialSolid<2, 3>>> solids);
+  std::vector<std::shared_ptr<SerialSolid<1, 2>>> solids);
 
 template void
-find_floating_mesh_mapping_step(
+find_floating_mesh_mapping_step<3>(
   const double smallest_contact_search_criterion,
-  std::vector<std::shared_ptr<SerialSolid<1, 2>>> solids);
+  std::vector<std::shared_ptr<SerialSolid<2, 3>>> solids);
diff --git a/source/dem/grid_motion.cc b/source/dem/grid_motion.cc
index 3d6789946b..4d460f5828 100644
--- a/source/dem/grid_motion.cc
+++ b/source/dem/grid_motion.cc
@@ -3,7 +3,6 @@
 
 #include <deal.II/grid/grid_tools.h>
 
-#include <boost/math/special_functions.hpp>
 #include <boost/range/adaptor/map.hpp>
 
 using namespace dealii;
diff --git a/source/dem/input_parameter_inspection.cc b/source/dem/input_parameter_inspection.cc
index 3acdd3ffae..f4335cef2e 100644
--- a/source/dem/input_parameter_inspection.cc
+++ b/source/dem/input_parameter_inspection.cc
@@ -79,7 +79,7 @@ input_parameter_inspection(const DEMSolverParameters<dim> &dem_parameters,
   if (parameters.grid_motion.motion_type ==
       Parameters::Lagrangian::GridMotion<dim>::MotionType::rotational)
     {
-      if (dim == 2)
+      if constexpr (dim == 2)
         {
           if (parameters.grid_motion.grid_rotational_axis != 0 &&
               parameters.grid_motion.grid_rotational_axis != 1)
@@ -87,7 +87,8 @@ input_parameter_inspection(const DEMSolverParameters<dim> &dem_parameters,
               "Specified grid rotational axis is not valid, use 0 for rotation around x"
               " axis or 1 for rotation around y axis.");
         }
-      else if (dim == 3)
+
+      if constexpr (dim == 3)
         {
           if (parameters.grid_motion.grid_rotational_axis != 0 &&
               parameters.grid_motion.grid_rotational_axis != 1 &&
@@ -100,13 +101,13 @@ input_parameter_inspection(const DEMSolverParameters<dim> &dem_parameters,
 }
 
 template void
-input_parameter_inspection(const DEMSolverParameters<2> &dem_parameters,
-                           const ConditionalOStream     &pcout,
-                           const std::vector<std::shared_ptr<Distribution>>
-                             &size_distribution_object_container);
+input_parameter_inspection<2>(const DEMSolverParameters<2> &dem_parameters,
+                              const ConditionalOStream     &pcout,
+                              const std::vector<std::shared_ptr<Distribution>>
+                                &size_distribution_object_container);
 
 template void
-input_parameter_inspection(const DEMSolverParameters<3> &dem_parameters,
-                           const ConditionalOStream     &pcout,
-                           const std::vector<std::shared_ptr<Distribution>>
-                             &size_distribution_object_container);
+input_parameter_inspection<3>(const DEMSolverParameters<3> &dem_parameters,
+                              const ConditionalOStream     &pcout,
+                              const std::vector<std::shared_ptr<Distribution>>
+                                &size_distribution_object_container);
diff --git a/source/dem/lagrangian_post_processing.cc b/source/dem/lagrangian_post_processing.cc
index 9e10aafa45..3fe71ecb0b 100644
--- a/source/dem/lagrangian_post_processing.cc
+++ b/source/dem/lagrangian_post_processing.cc
@@ -12,21 +12,16 @@
 
 using namespace dealii;
 
-template <int dim>
-LagrangianPostProcessing<dim>::LagrangianPostProcessing()
-{}
-
 template <int dim>
 void
-LagrangianPostProcessing<dim>::calculate_average_particles_velocity(
+calculate_average_particles_velocity(
   const parallel::distributed::Triangulation<dim> &triangulation,
-  const Particles::ParticleHandler<dim>           &particle_handler)
+  const Particles::ParticleHandler<dim>           &particle_handler,
+  Vector<double>                                  &velocity_average_x,
+  Vector<double>                                  &velocity_average_y,
+  Vector<double>                                  &velocity_average_z,
+  Vector<double>                                  &velocity_average_magnitude)
 {
-  velocity_average_x.reinit(triangulation.n_active_cells());
-  velocity_average_y.reinit(triangulation.n_active_cells());
-  velocity_average_z.reinit(triangulation.n_active_cells());
-  velocity_average_magnitude.reinit(triangulation.n_active_cells());
-
   // Iterating through the active cells in the triangulation
   for (const auto &cell : triangulation.active_cell_iterators())
     {
@@ -41,30 +36,30 @@ LagrangianPostProcessing<dim>::calculate_average_particles_velocity(
             cell_velocity_average[1];
 
           if constexpr (dim == 3)
-            velocity_average_z[cell->active_cell_index()] =
-              cell_velocity_average[2];
+            {
+              velocity_average_z[cell->active_cell_index()] =
+                cell_velocity_average[2];
+              velocity_average_magnitude[cell->active_cell_index()] =
+                sqrt(pow(cell_velocity_average[0], 2) +
+                     pow(cell_velocity_average[1], 2) +
+                     pow(cell_velocity_average[2], 2));
+            }
 
           if constexpr (dim == 2)
             velocity_average_magnitude[cell->active_cell_index()] =
               sqrt(pow(cell_velocity_average[0], 2) +
                    pow(cell_velocity_average[1], 2));
-          if constexpr (dim == 3)
-            velocity_average_magnitude[cell->active_cell_index()] =
-              sqrt(pow(cell_velocity_average[0], 2) +
-                   pow(cell_velocity_average[1], 2) +
-                   pow(cell_velocity_average[2], 2));
         }
     }
 }
 
 template <int dim>
 void
-LagrangianPostProcessing<dim>::calculate_average_granular_temperature(
+calculate_average_granular_temperature(
   const parallel::distributed::Triangulation<dim> &triangulation,
-  const Particles::ParticleHandler<dim>           &particle_handler)
+  const Particles::ParticleHandler<dim>           &particle_handler,
+  Vector<double>                                  &granular_temperature_average)
 {
-  granular_temperature_average.reinit(triangulation.n_active_cells());
-
   // Iterating through the active cells in the triangulation
   for (const auto &cell : triangulation.active_cell_iterators())
     {
@@ -130,7 +125,7 @@ LagrangianPostProcessing<dim>::calculate_average_granular_temperature(
 
 template <int dim>
 Tensor<1, dim>
-LagrangianPostProcessing<dim>::calculate_cell_average_particles_velocity(
+calculate_cell_average_particles_velocity(
   const typename parallel::distributed::Triangulation<dim>::cell_iterator &cell,
   const Particles::ParticleHandler<dim> &particle_handler)
 {
@@ -173,7 +168,7 @@ LagrangianPostProcessing<dim>::calculate_cell_average_particles_velocity(
 
 template <int dim>
 void
-LagrangianPostProcessing<dim>::write_post_processing_results(
+write_post_processing_results(
   const parallel::distributed::Triangulation<dim> &triangulation,
   PVDHandler                                      &grid_pvdhandler,
   const DoFHandler<dim>                           &background_dh,
@@ -193,53 +188,49 @@ LagrangianPostProcessing<dim>::write_post_processing_results(
   DataOut<dim> data_out;
   data_out.attach_dof_handler(background_dh);
 
-  std::vector<std::string> average_solution_names;
-
   // Write particles' average velocity
-  calculate_average_particles_velocity(triangulation, particle_handler);
-
-  average_solution_names.push_back("average_velocity_x");
-  average_solution_names.push_back("average_velocity_y");
-  if constexpr (dim == 3)
-    average_solution_names.push_back("average_velocity_z");
-  average_solution_names.push_back("average_velocity_magnitude");
-
-  average_solution_names.push_back("average_velocity_magnitude");
+  Vector<double> velocity_average_x(triangulation.n_active_cells());
+  Vector<double> velocity_average_y(triangulation.n_active_cells());
+  Vector<double> velocity_average_z(triangulation.n_active_cells());
+  Vector<double> velocity_average_magnitude(triangulation.n_active_cells());
+  calculate_average_particles_velocity<dim>(triangulation,
+                                            particle_handler,
+                                            velocity_average_x,
+                                            velocity_average_y,
+                                            velocity_average_z,
+                                            velocity_average_magnitude);
 
   data_out.add_data_vector(velocity_average_x,
-                           average_solution_names[0],
+                           "average_velocity_x",
                            DataOut<dim>::type_cell_data);
   data_out.add_data_vector(velocity_average_y,
-                           average_solution_names[1],
+                           "average_velocity_y",
                            DataOut<dim>::type_cell_data);
+
   if constexpr (dim == 3)
     data_out.add_data_vector(velocity_average_z,
-                             average_solution_names[2],
-                             DataOut<dim>::type_cell_data);
-  if constexpr (dim == 2)
-    data_out.add_data_vector(velocity_average_magnitude,
-                             average_solution_names[2],
-                             DataOut<dim>::type_cell_data);
-  if constexpr (dim == 3)
-    data_out.add_data_vector(velocity_average_magnitude,
-                             average_solution_names[3],
+                             "average_velocity_z",
                              DataOut<dim>::type_cell_data);
 
+  data_out.add_data_vector(velocity_average_magnitude,
+                           "average_velocity_magnitude",
+                           DataOut<dim>::type_cell_data);
 
   // Write particles' granular temperature
-  calculate_average_granular_temperature(triangulation, particle_handler);
-  average_solution_names.push_back("granular_temperature");
+  Vector<double> granular_temperature_average(triangulation.n_active_cells());
+  calculate_average_granular_temperature<dim>(triangulation,
+                                              particle_handler,
+                                              granular_temperature_average);
 
   data_out.add_data_vector(granular_temperature_average,
-                           average_solution_names.back(),
+                           "granular_temperature",
                            DataOut<dim>::type_cell_data);
 
-  average_solution_names.push_back("mobility_status");
-
+  // Write mobility status of cells
   Vector<float> mobility_status(triangulation.n_active_cells());
   sparse_contacts_object.get_mobility_status_vector(mobility_status);
   data_out.add_data_vector(mobility_status,
-                           average_solution_names.back(),
+                           "mobility_status",
                            DataOut<dim>::type_cell_data);
 
   // Attach the solution data to data_out object
@@ -263,5 +254,26 @@ LagrangianPostProcessing<dim>::write_post_processing_results(
                          mpi_communicator);
 }
 
-template class LagrangianPostProcessing<2>;
-template class LagrangianPostProcessing<3>;
+template void
+write_post_processing_results<2>(
+  const parallel::distributed::Triangulation<2> &triangulation,
+  PVDHandler                                    &grid_pvdhandler,
+  const DoFHandler<2>                           &background_dh,
+  const Particles::ParticleHandler<2>           &particle_handler,
+  const DEMSolverParameters<2>                  &dem_parameters,
+  const double                                   current_time,
+  const unsigned int                             step_number,
+  const MPI_Comm                                &mpi_communicator,
+  AdaptiveSparseContacts<2>                     &sparse_contacts_object);
+
+template void
+write_post_processing_results<3>(
+  const parallel::distributed::Triangulation<3> &triangulation,
+  PVDHandler                                    &grid_pvdhandler,
+  const DoFHandler<3>                           &background_dh,
+  const Particles::ParticleHandler<3>           &particle_handler,
+  const DEMSolverParameters<3>                  &dem_parameters,
+  const double                                   current_time,
+  const unsigned int                             step_number,
+  const MPI_Comm                                &mpi_communicator,
+  AdaptiveSparseContacts<3>                     &sparse_contacts_object);
diff --git a/source/dem/output_force_torque_calculation.cc b/source/dem/output_force_torque_calculation.cc
index 47e8838463..3975112f7c 100644
--- a/source/dem/output_force_torque_calculation.cc
+++ b/source/dem/output_force_torque_calculation.cc
@@ -26,3 +26,57 @@ write_forces_torques_output_locally(
 
   table.write_text(std::cout);
 }
+
+void
+write_forces_torques_output_results(
+  const std::string                                filename,
+  const unsigned int                               output_frequency,
+  const std::vector<unsigned int>                  boundary_index,
+  const double                                     time_step,
+  DEM::dem_data_structures<3>::vector_on_boundary &forces_boundary_information,
+  DEM::dem_data_structures<3>::vector_on_boundary &torques_boundary_information)
+{
+  unsigned int this_mpi_process =
+    Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  unsigned int n_mpi_processes =
+    Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  for (unsigned int i = 0; i < boundary_index.size(); i++)
+    {
+      std::string filename_boundary_description =
+        std::to_string(boundary_index[i]);
+      std::string filename_output =
+        filename + "_boundary_" + filename_boundary_description;
+      if (this_mpi_process == i || (i + 1) > n_mpi_processes)
+        {
+          TableHandler table;
+
+          for (unsigned int j = 0; j < forces_boundary_information.size();
+               j += output_frequency)
+            {
+              table.add_value("Force x", forces_boundary_information[j][i][0]);
+              table.add_value("Force y", forces_boundary_information[j][i][1]);
+              table.add_value("Force z", forces_boundary_information[j][i][2]);
+              table.add_value("Torque x",
+                              torques_boundary_information[j][i][0]);
+              table.add_value("Torque y",
+                              torques_boundary_information[j][i][1]);
+              table.add_value("Torque z",
+                              torques_boundary_information[j][i][2]);
+              table.add_value("Current time", j * time_step);
+            }
+
+          table.set_precision("Force x", 9);
+          table.set_precision("Force y", 9);
+          table.set_precision("Force z", 9);
+          table.set_precision("Torque x", 9);
+          table.set_precision("Torque y", 9);
+          table.set_precision("Torque z", 9);
+          table.set_precision("Current time", 9);
+
+          // output
+          std::ofstream out_file(filename_output);
+          table.write_text(out_file);
+          out_file.close();
+        }
+    }
+}
diff --git a/source/dem/particle_particle_broad_search.cc b/source/dem/particle_particle_broad_search.cc
index f609798757..27ed2b68ee 100644
--- a/source/dem/particle_particle_broad_search.cc
+++ b/source/dem/particle_particle_broad_search.cc
@@ -3,13 +3,9 @@
 
 using namespace dealii;
 
-template <int dim>
-ParticleParticleBroadSearch<dim>::ParticleParticleBroadSearch()
-{}
-
 template <int dim>
 void
-ParticleParticleBroadSearch<dim>::find_particle_particle_contact_pairs(
+find_particle_particle_contact_pairs(
   dealii::Particles::ParticleHandler<dim> &particle_handler,
   DEMContactManager<dim>                  &container_manager)
 {
@@ -53,10 +49,10 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_contact_pairs(
                particle_in_main_cell != particles_in_main_cell.end();
                ++particle_in_main_cell)
             {
-              store_candidates(particle_in_main_cell->get_id(),
-                               std::next(particle_in_main_cell, 1),
-                               particles_in_main_cell,
-                               local_contact_pair_candidates);
+              store_candidates<dim>(particle_in_main_cell->get_id(),
+                                    std::next(particle_in_main_cell, 1),
+                                    particles_in_main_cell,
+                                    local_contact_pair_candidates);
             }
 
           // Going through neighbor cells of the main cell
@@ -75,10 +71,10 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_contact_pairs(
                    particle_in_main_cell != particles_in_main_cell.end();
                    ++particle_in_main_cell)
                 {
-                  store_candidates(particle_in_main_cell->get_id(),
-                                   particles_in_neighbor_cell.begin(),
-                                   particles_in_neighbor_cell,
-                                   local_contact_pair_candidates);
+                  store_candidates<dim>(particle_in_main_cell->get_id(),
+                                        particles_in_neighbor_cell.begin(),
+                                        particles_in_neighbor_cell,
+                                        local_contact_pair_candidates);
                 }
             }
         }
@@ -122,10 +118,10 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_contact_pairs(
                    particle_in_main_cell != particles_in_main_cell.end();
                    ++particle_in_main_cell)
                 {
-                  store_candidates(particle_in_main_cell->get_id(),
-                                   particles_in_neighbor_cell.begin(),
-                                   particles_in_neighbor_cell,
-                                   ghost_contact_pair_candidates);
+                  store_candidates<dim>(particle_in_main_cell->get_id(),
+                                        particles_in_neighbor_cell.begin(),
+                                        particles_in_neighbor_cell,
+                                        ghost_contact_pair_candidates);
                 }
             }
         }
@@ -134,7 +130,7 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_contact_pairs(
 
 template <int dim>
 void
-ParticleParticleBroadSearch<dim>::find_particle_particle_contact_pairs(
+find_particle_particle_contact_pairs(
   dealii::Particles::ParticleHandler<dim> &particle_handler,
   DEMContactManager<dim>                  &container_manager,
   const AdaptiveSparseContacts<dim>       &sparse_contacts_object)
@@ -189,10 +185,10 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_contact_pairs(
                particle_in_main_cell != particles_in_main_cell.end();
                ++particle_in_main_cell)
             {
-              store_candidates(particle_in_main_cell->get_id(),
-                               std::next(particle_in_main_cell, 1),
-                               particles_in_main_cell,
-                               local_contact_pair_candidates);
+              store_candidates<dim>(particle_in_main_cell->get_id(),
+                                    std::next(particle_in_main_cell, 1),
+                                    particles_in_main_cell,
+                                    local_contact_pair_candidates);
             }
         }
 
@@ -229,10 +225,10 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_contact_pairs(
                particle_in_main_cell != particles_in_main_cell.end();
                ++particle_in_main_cell)
             {
-              store_candidates(particle_in_main_cell->get_id(),
-                               particles_in_neighbor_cell.begin(),
-                               particles_in_neighbor_cell,
-                               local_contact_pair_candidates);
+              store_candidates<dim>(particle_in_main_cell->get_id(),
+                                    particles_in_neighbor_cell.begin(),
+                                    particles_in_neighbor_cell,
+                                    local_contact_pair_candidates);
             }
         }
     }
@@ -298,10 +294,10 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_contact_pairs(
                particle_in_main_cell != particles_in_main_cell.end();
                ++particle_in_main_cell)
             {
-              store_candidates(particle_in_main_cell->get_id(),
-                               particles_in_neighbor_cell.begin(),
-                               particles_in_neighbor_cell,
-                               ghost_contact_pair_candidates);
+              store_candidates<dim>(particle_in_main_cell->get_id(),
+                                    particles_in_neighbor_cell.begin(),
+                                    particles_in_neighbor_cell,
+                                    ghost_contact_pair_candidates);
             }
         }
     }
@@ -310,7 +306,7 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_contact_pairs(
 
 template <int dim>
 void
-ParticleParticleBroadSearch<dim>::find_particle_particle_periodic_contact_pairs(
+find_particle_particle_periodic_contact_pairs(
   dealii::Particles::ParticleHandler<dim> &particle_handler,
   DEMContactManager<dim>                  &container_manager)
 {
@@ -373,10 +369,11 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_periodic_contact_pairs(
                    particle_in_main_cell != particles_in_main_cell.end();
                    ++particle_in_main_cell)
                 {
-                  store_candidates(particle_in_main_cell->get_id(),
-                                   particles_in_periodic_neighbor_cell.begin(),
-                                   particles_in_periodic_neighbor_cell,
-                                   local_contact_pair_periodic_candidates);
+                  store_candidates<dim>(
+                    particle_in_main_cell->get_id(),
+                    particles_in_periodic_neighbor_cell.begin(),
+                    particles_in_periodic_neighbor_cell,
+                    local_contact_pair_periodic_candidates);
                 }
             }
         }
@@ -425,10 +422,11 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_periodic_contact_pairs(
                    particle_in_main_cell != particles_in_main_cell.end();
                    ++particle_in_main_cell)
                 {
-                  store_candidates(particle_in_main_cell->get_id(),
-                                   particles_in_periodic_neighbor_cell.begin(),
-                                   particles_in_periodic_neighbor_cell,
-                                   ghost_contact_pair_periodic_candidates);
+                  store_candidates<dim>(
+                    particle_in_main_cell->get_id(),
+                    particles_in_periodic_neighbor_cell.begin(),
+                    particles_in_periodic_neighbor_cell,
+                    ghost_contact_pair_periodic_candidates);
                 }
             }
         }
@@ -474,7 +472,7 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_periodic_contact_pairs(
                    particle_in_main_cell != particles_in_main_cell.end();
                    ++particle_in_main_cell)
                 {
-                  store_candidates(
+                  store_candidates<dim>(
                     particle_in_main_cell->get_id(),
                     particles_in_periodic_neighbor_cell.begin(),
                     particles_in_periodic_neighbor_cell,
@@ -487,7 +485,7 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_periodic_contact_pairs(
 
 template <int dim>
 void
-ParticleParticleBroadSearch<dim>::find_particle_particle_periodic_contact_pairs(
+find_particle_particle_periodic_contact_pairs(
   dealii::Particles::ParticleHandler<dim> &particle_handler,
   DEMContactManager<dim>                  &container_manager,
   const AdaptiveSparseContacts<dim>       &sparse_contacts_object)
@@ -571,10 +569,10 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_periodic_contact_pairs(
                particle_in_main_cell != particles_in_main_cell.end();
                ++particle_in_main_cell)
             {
-              store_candidates(particle_in_main_cell->get_id(),
-                               particles_in_periodic_neighbor_cell.begin(),
-                               particles_in_periodic_neighbor_cell,
-                               local_contact_pair_periodic_candidates);
+              store_candidates<dim>(particle_in_main_cell->get_id(),
+                                    particles_in_periodic_neighbor_cell.begin(),
+                                    particles_in_periodic_neighbor_cell,
+                                    local_contact_pair_periodic_candidates);
             }
         }
     }
@@ -642,10 +640,10 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_periodic_contact_pairs(
                particle_in_main_cell != particles_in_main_cell.end();
                ++particle_in_main_cell)
             {
-              store_candidates(particle_in_main_cell->get_id(),
-                               particles_in_periodic_neighbor_cell.begin(),
-                               particles_in_periodic_neighbor_cell,
-                               ghost_contact_pair_periodic_candidates);
+              store_candidates<dim>(particle_in_main_cell->get_id(),
+                                    particles_in_periodic_neighbor_cell.begin(),
+                                    particles_in_periodic_neighbor_cell,
+                                    ghost_contact_pair_periodic_candidates);
             }
         }
     }
@@ -710,14 +708,93 @@ ParticleParticleBroadSearch<dim>::find_particle_particle_periodic_contact_pairs(
                particle_in_main_cell != particles_in_main_cell.end();
                ++particle_in_main_cell)
             {
-              store_candidates(particle_in_main_cell->get_id(),
-                               particles_in_periodic_neighbor_cell.begin(),
-                               particles_in_periodic_neighbor_cell,
-                               ghost_local_contact_pair_periodic_candidates);
+              store_candidates<dim>(
+                particle_in_main_cell->get_id(),
+                particles_in_periodic_neighbor_cell.begin(),
+                particles_in_periodic_neighbor_cell,
+                ghost_local_contact_pair_periodic_candidates);
             }
         }
     }
 }
 
-template class ParticleParticleBroadSearch<2>;
-template class ParticleParticleBroadSearch<3>;
+template <int dim>
+void
+store_candidates(
+  const types::particle_index &main_particle_id,
+  const typename Particles::ParticleHandler<
+    dim>::particle_iterator_range::iterator &particle_begin,
+  const typename Particles::ParticleHandler<dim>::particle_iterator_range
+    &particles_to_evaluate,
+  typename DEM::dem_data_structures<dim>::particle_particle_candidates
+    &contact_pair_candidates)
+{
+  // Find the contact candidate container of the main particle
+  auto candidates_container_it = contact_pair_candidates.find(main_particle_id);
+
+  // Reserve arbitrary vector capacity and store if the particle does not have
+  // contact candidate yet
+  if (candidates_container_it == contact_pair_candidates.end())
+    {
+      std::vector<types::particle_index> candidates_container;
+      candidates_container.reserve(40);
+
+      // Insert the empty vector and get the iterator to the inserted element
+      // prior storing the particle ids
+      auto pair_it_bool =
+        contact_pair_candidates.emplace(main_particle_id, candidates_container);
+      candidates_container_it = pair_it_bool.first;
+    }
+
+  // Store particle ids from the selected particle iterator
+  for (auto particle_iterator = particle_begin;
+       particle_iterator != particles_to_evaluate.end();
+       ++particle_iterator)
+    {
+      candidates_container_it->second.emplace_back(particle_iterator->get_id());
+    }
+}
+
+template void
+find_particle_particle_contact_pairs<2>(
+  dealii::Particles::ParticleHandler<2> &particle_handler,
+  DEMContactManager<2>                  &container_manager);
+
+template void
+find_particle_particle_contact_pairs<3>(
+  dealii::Particles::ParticleHandler<3> &particle_handler,
+  DEMContactManager<3>                  &container_manager);
+
+template void
+find_particle_particle_contact_pairs<2>(
+  dealii::Particles::ParticleHandler<2> &particle_handler,
+  DEMContactManager<2>                  &container_manager,
+  const AdaptiveSparseContacts<2>       &sparse_contacts_object);
+
+template void
+find_particle_particle_contact_pairs<3>(
+  dealii::Particles::ParticleHandler<3> &particle_handler,
+  DEMContactManager<3>                  &container_manager,
+  const AdaptiveSparseContacts<3>       &sparse_contacts_object);
+
+template void
+find_particle_particle_periodic_contact_pairs<2>(
+  dealii::Particles::ParticleHandler<2> &particle_handler,
+  DEMContactManager<2>                  &container_manager);
+
+template void
+find_particle_particle_periodic_contact_pairs<3>(
+  dealii::Particles::ParticleHandler<3> &particle_handler,
+  DEMContactManager<3>                  &container_manager);
+
+template void
+find_particle_particle_periodic_contact_pairs<2>(
+  dealii::Particles::ParticleHandler<2> &particle_handler,
+  DEMContactManager<2>                  &container_manager,
+  const AdaptiveSparseContacts<2>       &sparse_contacts_object);
+
+template void
+find_particle_particle_periodic_contact_pairs<3>(
+  dealii::Particles::ParticleHandler<3> &particle_handler,
+  DEMContactManager<3>                  &container_manager,
+  const AdaptiveSparseContacts<3>       &sparse_contacts_object);
diff --git a/source/dem/particle_particle_fine_search.cc b/source/dem/particle_particle_fine_search.cc
index 3cb59b720c..97db12a26b 100644
--- a/source/dem/particle_particle_fine_search.cc
+++ b/source/dem/particle_particle_fine_search.cc
@@ -1,15 +1,18 @@
+#include <core/dem_properties.h>
+
+#include <dem/contact_info.h>
 #include <dem/dem_contact_manager.h>
 #include <dem/particle_particle_fine_search.h>
 
-using namespace dealii;
+#include <deal.II/particles/particle.h>
 
-template <int dim>
-ParticleParticleFineSearch<dim>::ParticleParticleFineSearch()
-{}
+#include <boost/range/adaptor/map.hpp>
+
+using namespace dealii;
 
 template <int dim>
 void
-ParticleParticleFineSearch<dim>::particle_particle_fine_search(
+particle_particle_fine_search(
   const typename DEM::dem_data_structures<dim>::particle_index_iterator_map
     &particle_container,
   typename DEM::dem_data_structures<dim>::adjacent_particle_pairs
@@ -92,12 +95,32 @@ ParticleParticleFineSearch<dim>::particle_particle_fine_search(
 
               particle_one_contact_list.emplace(
                 particle_two_id,
-                particle_particle_contact_info<dim>(particle_one,
-                                                    particle_two));
+                particle_particle_contact_info<dim>{particle_one,
+                                                    particle_two,
+                                                    Tensor<1, 3>()});
             }
         }
     }
 }
 
-template class ParticleParticleFineSearch<2>;
-template class ParticleParticleFineSearch<3>;
+template void
+particle_particle_fine_search<2>(
+  typename DEM::dem_data_structures<2>::particle_index_iterator_map const
+    &particle_container,
+  typename DEM::dem_data_structures<2>::adjacent_particle_pairs
+    &adjacent_particles,
+  const typename DEM::dem_data_structures<2>::particle_particle_candidates
+                    &contact_pair_candidates,
+  const double       neighborhood_threshold,
+  const Tensor<1, 2> periodic_offset = Tensor<1, 2>());
+
+template void
+particle_particle_fine_search<3>(
+  typename DEM::dem_data_structures<3>::particle_index_iterator_map const
+    &particle_container,
+  typename DEM::dem_data_structures<3>::adjacent_particle_pairs
+    &adjacent_particles,
+  const typename DEM::dem_data_structures<3>::particle_particle_candidates
+                    &contact_pair_candidates,
+  const double       neighborhood_threshold,
+  const Tensor<1, 3> periodic_offset = Tensor<1, 3>());
diff --git a/source/dem/particle_point_line_broad_search.cc b/source/dem/particle_point_line_broad_search.cc
index 6c79f59da5..aa639d875c 100644
--- a/source/dem/particle_point_line_broad_search.cc
+++ b/source/dem/particle_point_line_broad_search.cc
@@ -1,25 +1,20 @@
 #include <dem/particle_point_line_broad_search.h>
 
 using namespace dealii;
-
-// Constructor
-template <int dim>
-ParticlePointLineBroadSearch<dim>::ParticlePointLineBroadSearch()
-{}
+using namespace DEM;
 
 // This function finds all the particle-point contact candidates
 template <int dim>
-typename DEM::dem_data_structures<dim>::particle_point_candidates
-ParticlePointLineBroadSearch<dim>::find_particle_point_contact_pairs(
+void
+find_particle_point_contact_pairs(
   const Particles::ParticleHandler<dim> &particle_handler,
-  const std::unordered_map<
-    std::string,
-    std::pair<typename Triangulation<dim>::active_cell_iterator, Point<dim>>>
-    &boundary_cells_with_points)
+  const std::unordered_map<std::string, cell_point_info<dim>>
+    &boundary_cells_with_points,
+  typename DEM::dem_data_structures<dim>::particle_point_candidates
+    &particle_point_contact_candidates)
 {
-  std::unordered_map<types::particle_index,
-                     std::pair<Particles::ParticleIterator<dim>, Point<dim>>>
-    particle_point_contact_candidates;
+  // Clear the candidate map
+  particle_point_contact_candidates.clear();
 
   // Defining and resetting a local particle-point candidate counter. This is
   // used as a key to the output map
@@ -34,19 +29,13 @@ ParticlePointLineBroadSearch<dim>::find_particle_point_contact_pairs(
        map_iterator != boundary_cells_with_points.end();
        ++map_iterator)
     {
-      auto cells_with_boundary_points_information = &map_iterator->second;
-
-      // Getting the cell and boundary vertex as local variables
-      auto cell_with_boundary_point =
-        cells_with_boundary_points_information->first;
-
-      Point<dim> vertex_location =
-        cells_with_boundary_points_information->second;
+      const cell_point_info<dim> &cells_with_boundary_points_information =
+        map_iterator->second;
 
       // Finding particles located in the corresponding cell
       typename Particles::ParticleHandler<dim>::particle_iterator_range
-        particles_in_cell =
-          particle_handler.particles_in_cell(cell_with_boundary_point);
+        particles_in_cell = particle_handler.particles_in_cell(
+          cells_with_boundary_points_information.cell);
 
       for (typename Particles::ParticleHandler<dim>::particle_iterator_range::
              iterator particles_in_cell_iterator = particles_in_cell.begin();
@@ -56,28 +45,28 @@ ParticlePointLineBroadSearch<dim>::find_particle_point_contact_pairs(
           // Making the pair and adding it to the
           // particle_point_contact_candidates map. This map is the output of
           // this function
-          particle_point_contact_candidates.insert(
-            {contact_candidate_counter,
-             std::make_pair(particles_in_cell_iterator, vertex_location)});
+          particle_point_contact_candidates.emplace(
+            contact_candidate_counter,
+            particle_point_contact_info<dim>{
+              particles_in_cell_iterator,
+              cells_with_boundary_points_information.point});
           ++contact_candidate_counter;
         }
     }
-  return particle_point_contact_candidates;
 }
 
 template <int dim>
-typename DEM::dem_data_structures<dim>::particle_point_candidates
-ParticlePointLineBroadSearch<dim>::find_particle_point_contact_pairs(
+void
+find_particle_point_contact_pairs(
   const Particles::ParticleHandler<dim> &particle_handler,
-  const std::unordered_map<
-    std::string,
-    std::pair<typename Triangulation<dim>::active_cell_iterator, Point<dim>>>
-                                    &boundary_cells_with_points,
+  const std::unordered_map<std::string, cell_point_info<dim>>
+    &boundary_cells_with_points,
+  typename DEM::dem_data_structures<dim>::particle_point_candidates
+                                    &particle_point_contact_candidates,
   const AdaptiveSparseContacts<dim> &sparse_contacts_object)
 {
-  std::unordered_map<types::particle_index,
-                     std::pair<Particles::ParticleIterator<dim>, Point<dim>>>
-    particle_point_contact_candidates;
+  // Clear the candidate map
+  particle_point_contact_candidates.clear();
 
   // Defining and resetting a local particle-point candidate counter. This is
   // used as a key to the output map
@@ -92,25 +81,20 @@ ParticlePointLineBroadSearch<dim>::find_particle_point_contact_pairs(
        map_iterator != boundary_cells_with_points.end();
        ++map_iterator)
     {
-      auto cells_with_boundary_points_information = &map_iterator->second;
-
-      // Getting the cell and boundary vertex as local variables
-      auto cell_with_boundary_point =
-        cells_with_boundary_points_information->first;
+      const cell_point_info<dim> &cells_with_boundary_points_information =
+        map_iterator->second;
 
       // If main cell has status other than mobile, skip to next cell
       unsigned int main_cell_mobility_status =
-        sparse_contacts_object.check_cell_mobility(cell_with_boundary_point);
+        sparse_contacts_object.check_cell_mobility(
+          cells_with_boundary_points_information.cell);
       if (main_cell_mobility_status != AdaptiveSparseContacts<dim>::mobile)
         continue;
 
-      Point<dim> vertex_location =
-        cells_with_boundary_points_information->second;
-
       // Finding particles located in the corresponding cell
       typename Particles::ParticleHandler<dim>::particle_iterator_range
-        particles_in_cell =
-          particle_handler.particles_in_cell(cell_with_boundary_point);
+        particles_in_cell = particle_handler.particles_in_cell(
+          cells_with_boundary_points_information.cell);
 
       for (typename Particles::ParticleHandler<dim>::particle_iterator_range::
              iterator particles_in_cell_iterator = particles_in_cell.begin();
@@ -120,30 +104,28 @@ ParticlePointLineBroadSearch<dim>::find_particle_point_contact_pairs(
           // Making the pair and adding it to the
           // particle_point_contact_candidates map. This map is the output of
           // this function
-          particle_point_contact_candidates.insert(
-            {contact_candidate_counter,
-             std::make_pair(particles_in_cell_iterator, vertex_location)});
+          particle_point_contact_candidates.emplace(
+            contact_candidate_counter,
+            particle_point_contact_info<dim>{
+              particles_in_cell_iterator,
+              cells_with_boundary_points_information.point});
           ++contact_candidate_counter;
         }
     }
-  return particle_point_contact_candidates;
 }
 
 // This function finds all the particle-line contact candidates
 template <int dim>
-typename DEM::dem_data_structures<dim>::particle_line_candidates
-ParticlePointLineBroadSearch<dim>::find_particle_line_contact_pairs(
+void
+find_particle_line_contact_pairs(
   const Particles::ParticleHandler<dim> &particle_handler,
-  const std::unordered_map<
-    std::string,
-    std::tuple<typename Triangulation<dim>::active_cell_iterator,
-               Point<dim>,
-               Point<dim>>> &boundary_cells_with_lines)
+  const std::unordered_map<std::string, cell_line_info<dim>>
+    &boundary_cells_with_lines,
+  typename DEM::dem_data_structures<dim>::particle_line_candidates
+    &particle_line_contact_candidates)
 {
-  std::unordered_map<
-    types::particle_index,
-    std::tuple<Particles::ParticleIterator<dim>, Point<dim>, Point<dim>>>
-    particle_line_contact_candidates;
+  // Clear the candidates map
+  particle_line_contact_candidates.clear();
 
   // Defining and resetting a local particle-line candidate counter. This is
   // used as a key to the output map
@@ -158,58 +140,43 @@ ParticlePointLineBroadSearch<dim>::find_particle_line_contact_pairs(
        map_iterator != boundary_cells_with_lines.end();
        ++map_iterator)
     {
-      auto cells_with_boundary_lines_information = &map_iterator->second;
-
-      // Getting the cell and the locations of boundary vertices (beginning and
-      // ending points of the boundary line) as local variables
-      auto cell_with_boundary_line =
-        std::get<0>(*cells_with_boundary_lines_information);
-
-      Point<dim> first_vertex_location =
-        std::get<1>(*cells_with_boundary_lines_information);
-      Point<dim> second_vertex_location =
-        std::get<2>(*cells_with_boundary_lines_information);
+      const cell_line_info<dim> &cells_with_boundary_lines_info =
+        map_iterator->second;
 
       // Finding particles located in the corresponding cell
       typename Particles::ParticleHandler<dim>::particle_iterator_range
-        particles_in_cell =
-          particle_handler.particles_in_cell(cell_with_boundary_line);
+        particles_in_cell = particle_handler.particles_in_cell(
+          cells_with_boundary_lines_info.cell);
 
-      for (typename Particles::ParticleHandler<dim>::particle_iterator_range::
-             iterator particles_in_cell_iterator = particles_in_cell.begin();
+      for (auto particles_in_cell_iterator = particles_in_cell.begin();
            particles_in_cell_iterator != particles_in_cell.end();
            ++particles_in_cell_iterator)
         {
-          // Making the tuple (particle, beginning and ending points of the
-          // boundary line) and adding it to the
-          // particle_line_contact_candidates map. This map is the output of
-          // this function
-          particle_line_contact_candidates.insert(
-            {contact_candidate_counter,
-             std::make_tuple(particles_in_cell_iterator,
-                             first_vertex_location,
-                             second_vertex_location)});
+          // Adding the particle line contact info to the
+          // particle_line_contact_candidates map.
+          particle_line_contact_candidates.emplace(
+            contact_candidate_counter,
+            particle_line_contact_info<dim>{
+              particles_in_cell_iterator,
+              cells_with_boundary_lines_info.point_one,
+              cells_with_boundary_lines_info.point_two});
           ++contact_candidate_counter;
         }
     }
-  return particle_line_contact_candidates;
 }
 
 template <int dim>
-typename DEM::dem_data_structures<dim>::particle_line_candidates
-ParticlePointLineBroadSearch<dim>::find_particle_line_contact_pairs(
+void
+find_particle_line_contact_pairs(
   const Particles::ParticleHandler<dim> &particle_handler,
-  const std::unordered_map<
-    std::string,
-    std::tuple<typename Triangulation<dim>::active_cell_iterator,
-               Point<dim>,
-               Point<dim>>>         &boundary_cells_with_lines,
+  const std::unordered_map<std::string, cell_line_info<dim>>
+    &boundary_cells_with_lines,
+  typename DEM::dem_data_structures<dim>::particle_line_candidates
+                                    &particle_line_contact_candidates,
   const AdaptiveSparseContacts<dim> &sparse_contacts_object)
 {
-  std::unordered_map<
-    types::particle_index,
-    std::tuple<Particles::ParticleIterator<dim>, Point<dim>, Point<dim>>>
-    particle_line_contact_candidates;
+  // Clear the candidates map
+  particle_line_contact_candidates.clear();
 
   // Defining and resetting a local particle-line candidate counter. This is
   // used as a key to the output map
@@ -224,47 +191,102 @@ ParticlePointLineBroadSearch<dim>::find_particle_line_contact_pairs(
        map_iterator != boundary_cells_with_lines.end();
        ++map_iterator)
     {
-      auto cells_with_boundary_lines_information = &map_iterator->second;
-
-      // Getting the cell and the locations of boundary vertices (beginning and
-      // ending points of the boundary line) as local variables
-      auto cell_with_boundary_line =
-        std::get<0>(*cells_with_boundary_lines_information);
+      const cell_line_info<dim> &cells_with_boundary_lines_info =
+        map_iterator->second;
 
       // If main cell has status other than mobile, skip to next cell
       unsigned int main_cell_mobility_status =
-        sparse_contacts_object.check_cell_mobility(cell_with_boundary_line);
+        sparse_contacts_object.check_cell_mobility(
+          cells_with_boundary_lines_info.cell);
       if (main_cell_mobility_status != AdaptiveSparseContacts<dim>::mobile)
         continue;
 
-      Point<dim> first_vertex_location =
-        std::get<1>(*cells_with_boundary_lines_information);
-      Point<dim> second_vertex_location =
-        std::get<2>(*cells_with_boundary_lines_information);
-
       // Finding particles located in the corresponding cell
       typename Particles::ParticleHandler<dim>::particle_iterator_range
-        particles_in_cell =
-          particle_handler.particles_in_cell(cell_with_boundary_line);
+        particles_in_cell = particle_handler.particles_in_cell(
+          cells_with_boundary_lines_info.cell);
 
       for (auto particles_in_cell_iterator = particles_in_cell.begin();
            particles_in_cell_iterator != particles_in_cell.end();
            ++particles_in_cell_iterator)
         {
-          // Making the tuple (particle, beginning and ending points of the
-          // boundary line) and adding it to the
-          // particle_line_contact_candidates map. This map is the output of
-          // this function
-          particle_line_contact_candidates.insert(
-            {contact_candidate_counter,
-             std::make_tuple(particles_in_cell_iterator,
-                             first_vertex_location,
-                             second_vertex_location)});
+          // Adding the particle line contact info to the
+          // particle_line_contact_candidates map.
+          particle_line_contact_candidates.emplace(
+            contact_candidate_counter,
+            particle_line_contact_info<dim>{
+              particles_in_cell_iterator,
+              cells_with_boundary_lines_info.point_one,
+              cells_with_boundary_lines_info.point_two});
           ++contact_candidate_counter;
         }
     }
-  return particle_line_contact_candidates;
 }
 
-template class ParticlePointLineBroadSearch<2>;
-template class ParticlePointLineBroadSearch<3>;
+template void
+find_particle_line_contact_pairs<2>(
+  const Particles::ParticleHandler<2> &particle_handler,
+  const std::unordered_map<std::string, cell_line_info<2>>
+    &boundary_cells_with_lines,
+  typename DEM::dem_data_structures<2>::particle_line_candidates
+    &particle_line_contact_candidates);
+
+template void
+find_particle_line_contact_pairs<3>(
+  const Particles::ParticleHandler<3> &particle_handler,
+  const std::unordered_map<std::string, cell_line_info<3>>
+    &boundary_cells_with_lines,
+  typename DEM::dem_data_structures<3>::particle_line_candidates
+    &particle_line_contact_candidates);
+
+template void
+find_particle_line_contact_pairs<2>(
+  const Particles::ParticleHandler<2> &particle_handler,
+  const std::unordered_map<std::string, cell_line_info<2>>
+    &boundary_cells_with_lines,
+  typename DEM::dem_data_structures<2>::particle_line_candidates
+                                  &particle_line_contact_candidates,
+  const AdaptiveSparseContacts<2> &sparse_contacts_object);
+
+template void
+find_particle_line_contact_pairs<3>(
+  const Particles::ParticleHandler<3> &particle_handler,
+  const std::unordered_map<std::string, cell_line_info<3>>
+    &boundary_cells_with_lines,
+  typename DEM::dem_data_structures<3>::particle_line_candidates
+                                  &particle_line_contact_candidates,
+  const AdaptiveSparseContacts<3> &sparse_contacts_object);
+
+template void
+find_particle_point_contact_pairs<2>(
+  const Particles::ParticleHandler<2> &particle_handler,
+  const std::unordered_map<std::string, cell_point_info<2>>
+    &boundary_cells_with_points,
+  typename DEM::dem_data_structures<2>::particle_point_candidates
+    &particle_point_contact_candidates);
+
+template void
+find_particle_point_contact_pairs<3>(
+  const Particles::ParticleHandler<3> &particle_handler,
+  const std::unordered_map<std::string, cell_point_info<3>>
+    &boundary_cells_with_points,
+  typename DEM::dem_data_structures<3>::particle_point_candidates
+    &particle_point_contact_candidates);
+
+template void
+find_particle_point_contact_pairs<2>(
+  const Particles::ParticleHandler<2> &particle_handler,
+  const std::unordered_map<std::string, cell_point_info<2>>
+    &boundary_cells_with_points,
+  typename DEM::dem_data_structures<2>::particle_point_candidates
+                                  &particle_point_contact_candidates,
+  const AdaptiveSparseContacts<2> &sparse_contacts_object);
+
+template void
+find_particle_point_contact_pairs<3>(
+  const Particles::ParticleHandler<3> &particle_handler,
+  const std::unordered_map<std::string, cell_point_info<3>>
+    &boundary_cells_with_points,
+  typename DEM::dem_data_structures<3>::particle_point_candidates
+                                  &particle_point_contact_candidates,
+  const AdaptiveSparseContacts<3> &sparse_contacts_object);
diff --git a/source/dem/particle_point_line_contact_force.cc b/source/dem/particle_point_line_contact_force.cc
index 0c7d5e4173..18711eae4e 100644
--- a/source/dem/particle_point_line_contact_force.cc
+++ b/source/dem/particle_point_line_contact_force.cc
@@ -14,7 +14,7 @@ ParticlePointLineForce<dim>::ParticlePointLineForce()
 template <int dim>
 void
 ParticlePointLineForce<dim>::calculate_particle_point_contact_force(
-  const typename DEM::dem_data_structures<dim>::particle_point_line_contact_info
+  const typename DEM::dem_data_structures<dim>::particle_point_in_contact
     *particle_point_pairs_in_contact,
   const Parameters::Lagrangian::LagrangianPhysicalProperties
                             &physical_properties,
@@ -42,7 +42,7 @@ ParticlePointLineForce<dim>::calculate_particle_point_contact_force(
       if constexpr (dim == 2)
         particle_location_3d = point_nd_to_3d(particle->get_location());
 
-      const Point<3> point = contact_information->point_one;
+      const Point<3> point = contact_information->point;
       double         normal_overlap =
         ((particle_properties[DEM::PropertiesIndex::dp]) / 2) -
         point.distance(particle_location_3d);
@@ -138,7 +138,7 @@ ParticlePointLineForce<dim>::calculate_particle_point_contact_force(
 template <int dim>
 void
 ParticlePointLineForce<dim>::calculate_particle_line_contact_force(
-  const typename DEM::dem_data_structures<dim>::particle_point_line_contact_info
+  const typename DEM::dem_data_structures<dim>::particle_line_in_contact
     *particle_line_pairs_in_contact,
   const Parameters::Lagrangian::LagrangianPhysicalProperties
                             &physical_properties,
diff --git a/source/dem/particle_point_line_fine_search.cc b/source/dem/particle_point_line_fine_search.cc
index 4725e2124f..ed6a19b358 100644
--- a/source/dem/particle_point_line_fine_search.cc
+++ b/source/dem/particle_point_line_fine_search.cc
@@ -1,29 +1,26 @@
+#include <core/dem_properties.h>
 #include <core/tensors_and_points_dimension_manipulation.h>
 
+#include <dem/contact_info.h>
 #include <dem/particle_point_line_fine_search.h>
 
+#include <deal.II/particles/particle_handler.h>
 using namespace dealii;
 
-template <int dim>
-ParticlePointLineFineSearch<dim>::ParticlePointLineFineSearch()
-{}
-
 // In this function, the output of particle-point broad search is investigated
 // to see if the pairs are in contact or not. If they are in contact, the normal
 // overlap, normal vector and normal relative velocity of the contact are
 // calculated. The output of this function is used for calculation of the
 // contact force
 template <int dim>
-typename DEM::dem_data_structures<dim>::particle_point_line_contact_info
-ParticlePointLineFineSearch<dim>::particle_point_fine_search(
+void
+particle_point_fine_search(
   const typename DEM::dem_data_structures<dim>::particle_point_candidates
               &particle_point_contact_candidates,
-  const double neighborhood_threshold)
+  const double neighborhood_threshold,
+  typename DEM::dem_data_structures<dim>::particle_point_in_contact
+    &particle_point_pairs_in_contact)
 {
-  std::unordered_map<types::particle_index,
-                     particle_point_line_contact_info_struct<dim>>
-    particle_point_pairs_in_contact;
-
   // Iterating over contact candidates from broad search. If a particle-point
   // pair is in contact (distance > 0) it is inserted into the output of this
   // function (particle_point_pairs_in_contact)
@@ -35,48 +32,37 @@ ParticlePointLineFineSearch<dim>::particle_point_fine_search(
     {
       // Get the value of the map (pair candidate) from the
       // contact_pair_candidates_iterator
-      auto pair_candidates = &contact_pair_candidates_iterator->second;
-
-      // Get the particle, particle properties and boundary vertex location once
-      // to improve efficiency
-      auto     particle            = std::get<0>(*pair_candidates);
-      auto     particle_properties = particle->get_properties();
-      Point<3> vertex_location_3d;
-
-      if constexpr (dim == 3)
-        vertex_location_3d = std::get<1>(*pair_candidates);
-
-      if constexpr (dim == 2)
-        vertex_location_3d = point_nd_to_3d(std::get<1>(*pair_candidates));
-
-      Point<3> particle_location_3d;
-      if constexpr (dim == 3)
-        particle_location_3d = particle->get_location();
-
-      if constexpr (dim == 2)
-        particle_location_3d = point_nd_to_3d(particle->get_location());
+      const particle_point_contact_info<dim> &pair_candidates =
+        contact_pair_candidates_iterator->second;
+
+      // Get the particle, particle diameter and boundary vertex location once
+      Particles::ParticleIterator<dim> particle = pair_candidates.particle;
+      double                           particle_diameter =
+        particle->get_properties()[DEM::PropertiesIndex::dp];
+
+      Point<3> vertex_location   = pair_candidates.point;
+      Point<3> particle_location = [&] {
+        if constexpr (dim == 3)
+          return particle->get_location();
+        else
+          return point_nd_to_3d(particle->get_location());
+      }();
 
       // Calculation of the square_distance between the particle and boundary
       // vertex
       const double square_distance =
-        ((particle_properties[DEM::PropertiesIndex::dp]) / 2) -
-        vertex_location_3d.distance_square(particle_location_3d);
+        (particle_diameter / 2.0) -
+        vertex_location.distance_square(particle_location);
 
-      // If the distance is larger than neighberhood threshold, then the
+      // If the distance is larger than neighborhood threshold, then the
       // particle-point pair are in contact
       if (square_distance > neighborhood_threshold)
         {
-          // Creating a sample from the particle_point_line_contact_info_struct
-          // and adding contact info to the sample
-          particle_point_line_contact_info_struct<dim> contact_info;
-          contact_info.particle  = particle;
-          contact_info.point_one = vertex_location_3d;
-
-          particle_point_pairs_in_contact.insert(
-            {particle->get_id(), contact_info});
+          // Adding contact info to the sample
+          particle_point_pairs_in_contact.emplace(particle->get_id(),
+                                                  pair_candidates);
         }
     }
-  return particle_point_pairs_in_contact;
 }
 
 // In this function, the output of particle-line broad search is investigated
@@ -85,19 +71,17 @@ ParticlePointLineFineSearch<dim>::particle_point_fine_search(
 // calculated. The output of this function is used for calculation of the
 // contact force
 template <int dim>
-typename DEM::dem_data_structures<dim>::particle_point_line_contact_info
-ParticlePointLineFineSearch<dim>::particle_line_fine_search(
+void
+particle_line_fine_search(
   const typename DEM::dem_data_structures<dim>::particle_line_candidates
               &particle_line_contact_candidates,
-  const double neighborhood_threshold)
+  const double neighborhood_threshold,
+  typename DEM::dem_data_structures<dim>::particle_line_in_contact
+    &particle_line_pairs_in_contact)
 {
-  std::unordered_map<types::particle_index,
-                     particle_point_line_contact_info_struct<dim>>
-    particle_line_pairs_in_contact;
-
   // Iterating over contact candidates from broad search. If a particle-line
   // pair is in contact (distance > 0) it is inserted into the output of this
-  // function (particle_line_pairs_in_contact)
+  // function (particle_line_pairs_in_contact);
   for (auto contact_pair_candidates_iterator =
          particle_line_contact_candidates.begin();
        contact_pair_candidates_iterator !=
@@ -106,71 +90,51 @@ ParticlePointLineFineSearch<dim>::particle_line_fine_search(
     {
       // Get the value of the map (pair candidate) from the
       // contact_pair_candidates_iterator
-      auto pair_candidates = &contact_pair_candidates_iterator->second;
-
-      // Get the particle, particle properties and the locations of beginning
-      // and ending vertices of the boundary line once to improve efficiency
-      auto     particle            = std::get<0>(*pair_candidates);
-      auto     particle_properties = particle->get_properties();
-      Point<3> vertex_one_location_3d;
-      Point<3> vertex_two_location_3d;
-
-
-      if constexpr (dim == 3)
-        {
-          vertex_one_location_3d = std::get<1>(*pair_candidates);
-          vertex_two_location_3d = std::get<2>(*pair_candidates);
-        }
-
-      if constexpr (dim == 2)
-        {
-          vertex_one_location_3d =
-            point_nd_to_3d(std::get<1>(*pair_candidates));
-          vertex_two_location_3d =
-            point_nd_to_3d(std::get<2>(*pair_candidates));
-        }
-
-      Point<3> particle_location_3d;
-      if constexpr (dim == 3)
-        particle_location_3d = particle->get_location();
-
-      if constexpr (dim == 2)
-        particle_location_3d = point_nd_to_3d(particle->get_location());
+      const particle_line_contact_info<dim> &pair_candidates =
+        contact_pair_candidates_iterator->second;
+
+      // Get the particle, particle diameter and the locations of beginning
+      // and ending vertices of the boundary line
+      Particles::ParticleIterator<dim> particle = pair_candidates.particle;
+      double                           particle_diameter =
+        particle->get_properties()[DEM::PropertiesIndex::dp];
+
+      Point<3> vertex_one_location = pair_candidates.point_one;
+      Point<3> vertex_two_location = pair_candidates.point_two;
+      Point<3> particle_location   = [&] {
+        if constexpr (dim == 3)
+          return particle->get_location();
+        else
+          return point_nd_to_3d(particle->get_location());
+      }();
 
       // For finding the particle-line distance, the projection of the particle
       // on the line should be obtained
-      Point<3> projection = find_projection_point(particle_location_3d,
-                                                  vertex_one_location_3d,
-                                                  vertex_two_location_3d);
+      Point<3> projection = find_projection_point(particle_location,
+                                                  vertex_one_location,
+                                                  vertex_two_location);
 
       // Calculation of the distance between the particle and boundary line
       const double square_distance =
-        ((particle_properties[DEM::PropertiesIndex::dp]) / 2) -
-        projection.distance_square(particle_location_3d);
+        (particle_diameter / 2.0) -
+        projection.distance_square(particle_location);
 
       // If the distance is positive, then the particle-line pair are in
       // contact
       if (square_distance > neighborhood_threshold)
         {
-          // Creating a sample from the particle_point_line_contact_info_struct
+          // Creating a sample from the particle_point_line_contact_info
           // and adding contact info to the sample
-          particle_point_line_contact_info_struct<dim> contact_info;
-          contact_info.particle  = particle;
-          contact_info.point_one = vertex_one_location_3d;
-          contact_info.point_two = vertex_two_location_3d;
-
-          particle_line_pairs_in_contact.insert(
-            {particle->get_id(), contact_info});
+          particle_line_pairs_in_contact.emplace(particle->get_id(),
+                                                 pair_candidates);
         }
     }
-  return particle_line_pairs_in_contact;
 }
 
-template <int dim>
 Point<3>
-ParticlePointLineFineSearch<dim>::find_projection_point(const Point<3> &point_p,
-                                                        const Point<3> &point_a,
-                                                        const Point<3> &point_b)
+find_projection_point(const Point<3> &point_p,
+                      const Point<3> &point_a,
+                      const Point<3> &point_b)
 {
   Tensor<1, 3> vector_ab = point_b - point_a;
   Tensor<1, 3> vector_ap = point_p - point_a;
@@ -181,5 +145,34 @@ ParticlePointLineFineSearch<dim>::find_projection_point(const Point<3> &point_p,
   return projection;
 }
 
-template class ParticlePointLineFineSearch<2>;
-template class ParticlePointLineFineSearch<3>;
+template void
+particle_point_fine_search<2>(
+  const typename DEM::dem_data_structures<2>::particle_point_candidates
+              &particle_point_contact_candidates,
+  const double neighborhood_threshold,
+  typename DEM::dem_data_structures<2>::particle_point_in_contact
+    &particle_point_pairs_in_contact);
+
+template void
+particle_point_fine_search<3>(
+  const typename DEM::dem_data_structures<3>::particle_point_candidates
+              &particle_point_contact_candidates,
+  const double neighborhood_threshold,
+  typename DEM::dem_data_structures<3>::particle_point_in_contact
+    &particle_point_pairs_in_contact);
+
+template void
+particle_line_fine_search<2>(
+  const typename DEM::dem_data_structures<2>::particle_line_candidates
+              &particle_line_contact_candidates,
+  const double neighborhood_threshold,
+  typename DEM::dem_data_structures<2>::particle_line_in_contact
+    &particle_line_pairs_in_contact);
+
+template void
+particle_line_fine_search<3>(
+  const typename DEM::dem_data_structures<3>::particle_line_candidates
+              &particle_line_contact_candidates,
+  const double neighborhood_threshold,
+  typename DEM::dem_data_structures<3>::particle_line_in_contact
+    &particle_line_pairs_in_contact);
diff --git a/source/dem/particle_wall_broad_search.cc b/source/dem/particle_wall_broad_search.cc
index cfd628d064..112554bb2b 100644
--- a/source/dem/particle_wall_broad_search.cc
+++ b/source/dem/particle_wall_broad_search.cc
@@ -2,13 +2,9 @@
 
 using namespace dealii;
 
-template <int dim>
-ParticleWallBroadSearch<dim>::ParticleWallBroadSearch()
-{}
-
 template <int dim>
 void
-ParticleWallBroadSearch<dim>::find_particle_wall_contact_pairs(
+find_particle_wall_contact_pairs(
   const std::map<int, boundary_cells_info_struct<dim>>
                                         &boundary_cells_information,
   const Particles::ParticleHandler<dim> &particle_handler,
@@ -58,7 +54,7 @@ ParticleWallBroadSearch<dim>::find_particle_wall_contact_pairs(
 
 template <int dim>
 void
-ParticleWallBroadSearch<dim>::find_particle_floating_wall_contact_pairs(
+find_particle_floating_wall_contact_pairs(
   const std::unordered_map<
     DEM::global_face_id,
     std::set<typename Triangulation<dim>::active_cell_iterator>>
@@ -115,9 +111,9 @@ ParticleWallBroadSearch<dim>::find_particle_floating_wall_contact_pairs(
                        particle_in_cell_iterator != particles_in_cell.end();
                        ++particle_in_cell_iterator)
                     {
-                      store_candidates(particle_in_cell_iterator,
-                                       floating_wall_id,
-                                       particle_floating_wall_candidates);
+                      store_candidates<dim>(particle_in_cell_iterator,
+                                            floating_wall_id,
+                                            particle_floating_wall_candidates);
                     }
                 }
             }
@@ -127,7 +123,7 @@ ParticleWallBroadSearch<dim>::find_particle_floating_wall_contact_pairs(
 
 template <int dim>
 void
-ParticleWallBroadSearch<dim>::particle_solid_surfaces_contact_search(
+particle_solid_surfaces_contact_search(
   const typename DEM::dem_data_structures<dim>::solid_surfaces_mesh_information
                                         &solid_surfaces_mesh_information,
   const Particles::ParticleHandler<dim> &particle_handler,
@@ -213,7 +209,7 @@ ParticleWallBroadSearch<dim>::particle_solid_surfaces_contact_search(
 
 template <int dim>
 void
-ParticleWallBroadSearch<dim>::find_particle_wall_contact_pairs(
+find_particle_wall_contact_pairs(
   const std::map<int, boundary_cells_info_struct<dim>>
                                         &boundary_cells_information,
   const Particles::ParticleHandler<dim> &particle_handler,
@@ -265,7 +261,7 @@ ParticleWallBroadSearch<dim>::find_particle_wall_contact_pairs(
 
 template <int dim>
 void
-ParticleWallBroadSearch<dim>::find_particle_floating_wall_contact_pairs(
+find_particle_floating_wall_contact_pairs(
   const std::unordered_map<
     types::global_dof_index,
     std::set<typename Triangulation<dim>::active_cell_iterator>>
@@ -326,9 +322,9 @@ ParticleWallBroadSearch<dim>::find_particle_floating_wall_contact_pairs(
                    particle_in_cell_iterator != particles_in_cell.end();
                    ++particle_in_cell_iterator)
                 {
-                  store_candidates(particle_in_cell_iterator,
-                                   floating_wall_id,
-                                   particle_floating_wall_candidates);
+                  store_candidates<dim>(particle_in_cell_iterator,
+                                        floating_wall_id,
+                                        particle_floating_wall_candidates);
                 }
             }
         }
@@ -337,7 +333,7 @@ ParticleWallBroadSearch<dim>::find_particle_floating_wall_contact_pairs(
 
 template <int dim>
 void
-ParticleWallBroadSearch<dim>::particle_solid_surfaces_contact_search(
+particle_solid_surfaces_contact_search(
   const typename DEM::dem_data_structures<dim>::solid_surfaces_mesh_information
                                         &solid_surfaces_mesh_information,
   const Particles::ParticleHandler<dim> &particle_handler,
@@ -419,6 +415,159 @@ ParticleWallBroadSearch<dim>::particle_solid_surfaces_contact_search(
     }
 }
 
+template <int dim>
+void
+store_candidates(
+  const typename Particles::ParticleHandler<
+    dim>::particle_iterator_range::iterator &particle_iterator,
+  const DEM::global_face_id                 &floating_wall_id,
+  typename DEM::dem_data_structures<dim>::particle_floating_wall_candidates
+    &contact_pair_candidates)
+{
+  // Find the contact candidate container of the particle
+  const types::particle_index particle_id = particle_iterator->get_id();
+  auto candidates_container_it = contact_pair_candidates.find(particle_id);
+
+  // Reserve arbitrary vector capacity and store if the particle does not have
+  // contact candidate yet
+  if (candidates_container_it == contact_pair_candidates.end())
+    {
+      auto pair_it_bool = contact_pair_candidates.emplace(
+        particle_id,
+        std::unordered_map<DEM::global_face_id,
+                           Particles::ParticleIterator<dim>>());
+
+      candidates_container_it = pair_it_bool.first;
+    }
+
+  // Store particle ids from the selected particle iterator
+  candidates_container_it->second.emplace(floating_wall_id, particle_iterator);
+}
+
 
-template class ParticleWallBroadSearch<2>;
-template class ParticleWallBroadSearch<3>;
+// Explicit Instantiation
+template void
+find_particle_wall_contact_pairs<2>(
+  const std::map<int, boundary_cells_info_struct<2>>
+                                      &boundary_cells_information,
+  const Particles::ParticleHandler<2> &particle_handler,
+  DEM::dem_data_structures<2>::particle_wall_candidates
+    &particle_wall_contact_candidates);
+
+template void
+find_particle_wall_contact_pairs<3>(
+  const std::map<int, boundary_cells_info_struct<3>>
+                                      &boundary_cells_information,
+  const Particles::ParticleHandler<3> &particle_handler,
+  DEM::dem_data_structures<3>::particle_wall_candidates
+    &particle_wall_contact_candidates);
+
+template void
+find_particle_wall_contact_pairs<2>(
+  const std::map<int, boundary_cells_info_struct<2>>
+                                      &boundary_cells_information,
+  const Particles::ParticleHandler<2> &particle_handler,
+  DEM::dem_data_structures<2>::particle_wall_candidates
+                                  &particle_wall_contact_candidates,
+  const AdaptiveSparseContacts<2> &sparse_contacts_object);
+
+template void
+find_particle_wall_contact_pairs<3>(
+  const std::map<int, boundary_cells_info_struct<3>>
+                                      &boundary_cells_information,
+  const Particles::ParticleHandler<3> &particle_handler,
+  DEM::dem_data_structures<3>::particle_wall_candidates
+                                  &particle_wall_contact_candidates,
+  const AdaptiveSparseContacts<3> &sparse_contacts_object);
+
+template void
+find_particle_floating_wall_contact_pairs<2>(
+  const std::unordered_map<
+    types::global_dof_index,
+    std::set<typename Triangulation<2>::active_cell_iterator>>
+                                      &boundary_cells_for_floating_walls,
+  const Particles::ParticleHandler<2> &particle_handler,
+  const Parameters::Lagrangian::FloatingWalls<2> &floating_wall_properties,
+  const double                                    simulation_time,
+  DEM::dem_data_structures<2>::particle_floating_wall_candidates
+    &particle_floating_wall_candidates);
+
+template void
+find_particle_floating_wall_contact_pairs<3>(
+  const std::unordered_map<
+    types::global_dof_index,
+    std::set<typename Triangulation<3>::active_cell_iterator>>
+                                      &boundary_cells_for_floating_walls,
+  const Particles::ParticleHandler<3> &particle_handler,
+  const Parameters::Lagrangian::FloatingWalls<3> &floating_wall_properties,
+  const double                                    simulation_time,
+  DEM::dem_data_structures<3>::particle_floating_wall_candidates
+    &particle_floating_wall_candidates);
+
+template void
+find_particle_floating_wall_contact_pairs<2>(
+  const std::unordered_map<
+    types::global_dof_index,
+    std::set<typename Triangulation<2>::active_cell_iterator>>
+                                      &boundary_cells_for_floating_walls,
+  const Particles::ParticleHandler<2> &particle_handler,
+  const Parameters::Lagrangian::FloatingWalls<2> &floating_wall_properties,
+  const double                                    simulation_time,
+  DEM::dem_data_structures<2>::particle_floating_wall_candidates
+                                  &particle_floating_wall_candidates,
+  const AdaptiveSparseContacts<2> &sparse_contacts_object);
+
+template void
+find_particle_floating_wall_contact_pairs<3>(
+  const std::unordered_map<
+    types::global_dof_index,
+    std::set<typename Triangulation<3>::active_cell_iterator>>
+                                      &boundary_cells_for_floating_walls,
+  const Particles::ParticleHandler<3> &particle_handler,
+  const Parameters::Lagrangian::FloatingWalls<3> &floating_wall_properties,
+  const double                                    simulation_time,
+  DEM::dem_data_structures<3>::particle_floating_wall_candidates
+                                  &particle_floating_wall_candidates,
+  const AdaptiveSparseContacts<3> &sparse_contacts_object);
+
+template void
+particle_solid_surfaces_contact_search<2>(
+  const DEM::dem_data_structures<2>::solid_surfaces_mesh_information
+                                      &solid_surfaces_mesh_information,
+  const Particles::ParticleHandler<2> &particle_handler,
+  DEM::dem_data_structures<2>::particle_floating_mesh_candidates
+    &particle_floating_mesh_contact_candidates,
+  DEM::dem_data_structures<2>::cells_total_neighbor_list
+    &cells_total_neighbor_list);
+
+template void
+particle_solid_surfaces_contact_search<3>(
+  const DEM::dem_data_structures<3>::solid_surfaces_mesh_information
+                                      &solid_surfaces_mesh_information,
+  const Particles::ParticleHandler<3> &particle_handler,
+  DEM::dem_data_structures<3>::particle_floating_mesh_candidates
+    &particle_floating_mesh_contact_candidates,
+  DEM::dem_data_structures<3>::cells_total_neighbor_list
+    &cells_total_neighbor_list);
+
+template void
+particle_solid_surfaces_contact_search<2>(
+  const DEM::dem_data_structures<2>::solid_surfaces_mesh_information
+                                      &solid_surfaces_mesh_information,
+  const Particles::ParticleHandler<2> &particle_handler,
+  DEM::dem_data_structures<2>::particle_floating_mesh_candidates
+    &particle_floating_mesh_contact_candidates,
+  DEM::dem_data_structures<2>::cells_total_neighbor_list
+                                  &cells_total_neighbor_list,
+  const AdaptiveSparseContacts<2> &sparse_contacts_object);
+
+template void
+particle_solid_surfaces_contact_search<3>(
+  const DEM::dem_data_structures<3>::solid_surfaces_mesh_information
+                                      &solid_surfaces_mesh_information,
+  const Particles::ParticleHandler<3> &particle_handler,
+  DEM::dem_data_structures<3>::particle_floating_mesh_candidates
+    &particle_floating_mesh_contact_candidates,
+  DEM::dem_data_structures<3>::cells_total_neighbor_list
+                                  &cells_total_neighbor_list,
+  const AdaptiveSparseContacts<3> &sparse_contacts_object);
diff --git a/source/dem/particle_wall_fine_search.cc b/source/dem/particle_wall_fine_search.cc
index 4a01a5e0db..38d153007a 100644
--- a/source/dem/particle_wall_fine_search.cc
+++ b/source/dem/particle_wall_fine_search.cc
@@ -1,22 +1,23 @@
 #include <core/tensors_and_points_dimension_manipulation.h>
 
+#include <dem/contact_info.h>
 #include <dem/particle_wall_fine_search.h>
 
+#include <deal.II/particles/particle_handler.h>
+
 #include <boost/math/special_functions/sign.hpp>
 
+#include <vector>
 
 using namespace dealii;
-
-template <int dim>
-ParticleWallFineSearch<dim>::ParticleWallFineSearch()
-{}
+using namespace DEM;
 
 template <int dim>
 void
-ParticleWallFineSearch<dim>::particle_wall_fine_search(
-  const typename DEM::dem_data_structures<dim>::particle_wall_candidates
+particle_wall_fine_search(
+  const typename dem_data_structures<dim>::particle_wall_candidates
     &particle_wall_contact_pair_candidates,
-  typename DEM::dem_data_structures<dim>::particle_wall_in_contact
+  typename dem_data_structures<dim>::particle_wall_in_contact
     &particle_wall_pairs_in_contact)
 {
   // Iterating over contact candidates from broad search and adding the pairs to
@@ -66,12 +67,12 @@ ParticleWallFineSearch<dim>::particle_wall_fine_search(
 
 template <int dim>
 void
-ParticleWallFineSearch<dim>::particle_floating_wall_fine_search(
-  const typename DEM::dem_data_structures<
-    dim>::particle_floating_wall_candidates &particle_floating_wall_candidates,
+particle_floating_wall_fine_search(
+  const typename dem_data_structures<dim>::particle_floating_wall_candidates
+    &particle_floating_wall_candidates,
   const Parameters::Lagrangian::FloatingWalls<dim> &floating_wall_properties,
   const double                                      simulation_time,
-  typename DEM::dem_data_structures<dim>::particle_wall_in_contact
+  typename dem_data_structures<dim>::particle_wall_in_contact
     &particle_floating_wall_pairs_in_contact)
 {
   // Reading floating wall properties
@@ -160,11 +161,10 @@ ParticleWallFineSearch<dim>::particle_floating_wall_fine_search(
 
 template <int dim>
 void
-ParticleWallFineSearch<dim>::particle_floating_mesh_fine_search(
-  const typename DEM::dem_data_structures<
-    dim>::particle_floating_mesh_candidates
+particle_floating_mesh_fine_search(
+  const typename dem_data_structures<dim>::particle_floating_mesh_candidates
     &particle_floating_mesh_contact_candidates,
-  typename DEM::dem_data_structures<dim>::particle_floating_mesh_in_contact
+  typename dem_data_structures<dim>::particle_floating_mesh_in_contact
     &particle_floating_mesh_in_contact)
 {
   for (unsigned int solid_counter = 0;
@@ -194,5 +194,48 @@ ParticleWallFineSearch<dim>::particle_floating_mesh_fine_search(
     }
 }
 
-template class ParticleWallFineSearch<2>;
-template class ParticleWallFineSearch<3>;
+template void
+particle_wall_fine_search<2>(
+  const typename dem_data_structures<2>::particle_wall_candidates
+    &particle_wall_contact_pair_candidates,
+  typename dem_data_structures<2>::particle_wall_in_contact
+    &particle_wall_pairs_in_contact);
+
+template void
+particle_wall_fine_search<3>(
+  const typename dem_data_structures<3>::particle_wall_candidates
+    &particle_wall_contact_pair_candidates,
+  typename dem_data_structures<3>::particle_wall_in_contact
+    &particle_wall_pairs_in_contact);
+
+template void
+particle_floating_wall_fine_search<2>(
+  const typename dem_data_structures<2>::particle_floating_wall_candidates
+    &particle_floating_wall_candidates,
+  const Parameters::Lagrangian::FloatingWalls<2> &floating_wall_properties,
+  const double                                    simulation_time,
+  typename dem_data_structures<2>::particle_wall_in_contact
+    &particle_floating_wall_pairs_in_contact);
+
+template void
+particle_floating_wall_fine_search<3>(
+  const typename dem_data_structures<3>::particle_floating_wall_candidates
+    &particle_floating_wall_candidates,
+  const Parameters::Lagrangian::FloatingWalls<3> &floating_wall_properties,
+  const double                                    simulation_time,
+  typename dem_data_structures<3>::particle_wall_in_contact
+    &particle_floating_wall_pairs_in_contact);
+
+template void
+particle_floating_mesh_fine_search<2>(
+  const typename dem_data_structures<2>::particle_floating_mesh_candidates
+    &particle_floating_mesh_contact_candidates,
+  typename dem_data_structures<2>::particle_floating_mesh_in_contact
+    &particle_floating_mesh_in_contact);
+
+template void
+particle_floating_mesh_fine_search<3>(
+  const typename dem_data_structures<3>::particle_floating_mesh_candidates
+    &particle_floating_mesh_contact_candidates,
+  typename dem_data_structures<3>::particle_floating_mesh_in_contact
+    &particle_floating_mesh_in_contact);
diff --git a/source/dem/update_fine_search_candidates.cc b/source/dem/update_fine_search_candidates.cc
index e31eb99ae4..a7687d760b 100644
--- a/source/dem/update_fine_search_candidates.cc
+++ b/source/dem/update_fine_search_candidates.cc
@@ -1,4 +1,4 @@
-#include <dem/contact_type.h>
+#include <dem/contact_info.h>
 #include <dem/update_fine_search_candidates.h>
 
 using namespace dealii;
diff --git a/source/dem/update_local_particle_containers.cc b/source/dem/update_local_particle_containers.cc
index 518b3f25a0..273504cd98 100644
--- a/source/dem/update_local_particle_containers.cc
+++ b/source/dem/update_local_particle_containers.cc
@@ -367,9 +367,9 @@ update_contact_container_iterators<
 template void
 update_contact_container_iterators<
   2,
-  typename DEM::dem_data_structures<2>::particle_point_line_contact_info,
+  typename DEM::dem_data_structures<2>::particle_point_in_contact,
   ContactType::particle_point>(
-  typename DEM::dem_data_structures<2>::particle_point_line_contact_info
+  typename DEM::dem_data_structures<2>::particle_point_in_contact
     &pairs_in_contact,
   const typename DEM::dem_data_structures<2>::particle_index_iterator_map
     &particle_container);
@@ -377,9 +377,9 @@ update_contact_container_iterators<
 template void
 update_contact_container_iterators<
   3,
-  typename DEM::dem_data_structures<3>::particle_point_line_contact_info,
+  typename DEM::dem_data_structures<3>::particle_point_in_contact,
   ContactType::particle_point>(
-  typename DEM::dem_data_structures<3>::particle_point_line_contact_info
+  typename DEM::dem_data_structures<3>::particle_point_in_contact
     &pairs_in_contact,
   const typename DEM::dem_data_structures<3>::particle_index_iterator_map
     &particle_container);
@@ -388,9 +388,9 @@ update_contact_container_iterators<
 template void
 update_contact_container_iterators<
   2,
-  typename DEM::dem_data_structures<2>::particle_point_line_contact_info,
+  typename DEM::dem_data_structures<2>::particle_line_in_contact,
   ContactType::particle_line>(
-  typename DEM::dem_data_structures<2>::particle_point_line_contact_info
+  typename DEM::dem_data_structures<2>::particle_line_in_contact
     &pairs_in_contact,
   const typename DEM::dem_data_structures<2>::particle_index_iterator_map
     &particle_container);
@@ -398,9 +398,9 @@ update_contact_container_iterators<
 template void
 update_contact_container_iterators<
   3,
-  typename DEM::dem_data_structures<3>::particle_point_line_contact_info,
+  typename DEM::dem_data_structures<3>::particle_line_in_contact,
   ContactType::particle_line>(
-  typename DEM::dem_data_structures<3>::particle_point_line_contact_info
+  typename DEM::dem_data_structures<3>::particle_line_in_contact
     &pairs_in_contact,
   const typename DEM::dem_data_structures<3>::particle_index_iterator_map
     &particle_container);
diff --git a/source/fem-dem/cfd_dem_coupling.cc b/source/fem-dem/cfd_dem_coupling.cc
index 5aa3b4ddbd..b02b39998b 100644
--- a/source/fem-dem/cfd_dem_coupling.cc
+++ b/source/fem-dem/cfd_dem_coupling.cc
@@ -947,8 +947,7 @@ CFDDEMSolver<dim>::report_particle_statistics()
         dynamic_cast<parallel::distributed::Triangulation<dim> *>(
           &*this->triangulation);
 
-
-      dem_post_processing_object.write_post_processing_results(
+      write_post_processing_results<dim>(
         *parallel_triangulation,
         grid_pvdhandler,
         this->dof_handler,
diff --git a/tests/dem/find_cell_neighbors.cc b/tests/dem/find_cell_neighbors.cc
index abd5c3454c..10c708dc7f 100644
--- a/tests/dem/find_cell_neighbors.cc
+++ b/tests/dem/find_cell_neighbors.cc
@@ -52,10 +52,9 @@ test()
   std::vector<std::vector<typename Triangulation<dim>::active_cell_iterator>>
     cells_ghost_neighbor_list;
 
-  FindCellNeighbors<dim> cell_neighbor_object;
-  cell_neighbor_object.find_cell_neighbors(triangulation,
-                                           cells_local_neighbor_list,
-                                           cells_ghost_neighbor_list);
+  find_cell_neighbors<dim>(triangulation,
+                           cells_local_neighbor_list,
+                           cells_ghost_neighbor_list);
 
   // Output
   int i = 0;
diff --git a/tests/dem/find_contact_pairs.cc b/tests/dem/find_contact_pairs.cc
index 0b74e954f7..a16095d69c 100644
--- a/tests/dem/find_contact_pairs.cc
+++ b/tests/dem/find_contact_pairs.cc
@@ -58,16 +58,14 @@ test()
 
   MappingQ1<dim> mapping;
 
-  DEMContactManager<dim>          container_manager;
+  DEMContactManager<dim>          contact_manager;
   Particles::ParticleHandler<dim> particle_handler(triangulation, mapping);
 
   // Finding cell neighbors list, it is required for finding the broad search
-  // pairs in the container_manager
-  FindCellNeighbors<dim> cell_neighbor_object;
-  cell_neighbor_object.find_cell_neighbors(
-    triangulation,
-    container_manager.cells_local_neighbor_list,
-    container_manager.cells_ghost_neighbor_list);
+  // pairs in the contact_manager
+  find_cell_neighbors<dim>(triangulation,
+                           contact_manager.cells_local_neighbor_list,
+                           contact_manager.cells_ghost_neighbor_list);
 
   // inserting three particles at x = -0.4 , x = 0.4 and x = 0.8
   // which means they are inserted in three adjacent cells in x direction
@@ -107,13 +105,13 @@ test()
   AdaptiveSparseContacts<dim> dummy_adaptive_sparse_contacts;
 
   // Calling broad search function
-  container_manager.execute_particle_particle_broad_search(
+  contact_manager.execute_particle_particle_broad_search(
     particle_handler, dummy_adaptive_sparse_contacts);
 
   // Output
   for (auto pairs_iterator =
-         container_manager.local_contact_pair_candidates.begin();
-       pairs_iterator != container_manager.local_contact_pair_candidates.end();
+         contact_manager.local_contact_pair_candidates.begin();
+       pairs_iterator != contact_manager.local_contact_pair_candidates.end();
        ++pairs_iterator)
     {
       unsigned int first_particle_id = pairs_iterator->first;
@@ -127,8 +125,8 @@ test()
     }
 
   for (auto pairs_iterator =
-         container_manager.local_contact_pair_candidates.begin();
-       pairs_iterator != container_manager.local_contact_pair_candidates.end();
+         contact_manager.local_contact_pair_candidates.begin();
+       pairs_iterator != contact_manager.local_contact_pair_candidates.end();
        ++pairs_iterator)
     {
       unsigned int first_particle_id = pairs_iterator->first;
diff --git a/tests/dem/find_full_cell_neighbors.cc b/tests/dem/find_full_cell_neighbors.cc
index 09c3c88f85..d21336dd96 100644
--- a/tests/dem/find_full_cell_neighbors.cc
+++ b/tests/dem/find_full_cell_neighbors.cc
@@ -51,9 +51,7 @@ test()
   typename DEM::dem_data_structures<dim>::cells_total_neighbor_list
     cells_total_neighbor_list;
 
-  FindCellNeighbors<dim> cell_neighbor_object;
-  cell_neighbor_object.find_full_cell_neighbors(triangulation,
-                                                cells_total_neighbor_list);
+  find_full_cell_neighbors<dim>(triangulation, cells_total_neighbor_list);
 
   // Output
   int i = 0;
diff --git a/tests/dem/normal_force.cc b/tests/dem/normal_force.cc
index ead7bc4190..369f08c5d6 100644
--- a/tests/dem/normal_force.cc
+++ b/tests/dem/normal_force.cc
@@ -152,16 +152,14 @@ test()
                        Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
 
   // P-W broad search
-  ParticleWallBroadSearch<dim> particle_wall_broad_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_candidates
     particle_wall_contact_list;
-  particle_wall_broad_search_object.find_particle_wall_contact_pairs(
+  find_particle_wall_contact_pairs<dim>(
     boundary_cells_object.get_boundary_cells_information(),
     particle_handler,
     particle_wall_contact_list);
 
   // Particle-Wall fine search
-  ParticleWallFineSearch<dim> particle_wall_fine_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_in_contact
                                   particle_wall_contact_information;
   ParticleWallNonLinearForce<dim> particle_wall_force_object(dem_parameters);
@@ -191,8 +189,8 @@ test()
       else
         {
           // If particle and wall are in contact
-          particle_wall_fine_search_object.particle_wall_fine_search(
-            particle_wall_contact_list, particle_wall_contact_information);
+          particle_wall_fine_search<dim>(particle_wall_contact_list,
+                                         particle_wall_contact_information);
           auto particle_wall_pairs_in_contact_iterator =
             &particle_wall_contact_information.begin()->second;
           auto particle_wall_contact_information_iterator =
diff --git a/tests/dem/particle_particle_contact_force_linear.cc b/tests/dem/particle_particle_contact_force_linear.cc
index e9255299da..3f7373b762 100644
--- a/tests/dem/particle_particle_contact_force_linear.cc
+++ b/tests/dem/particle_particle_contact_force_linear.cc
@@ -36,11 +36,7 @@
 #include <core/dem_properties.h>
 
 #include <dem/dem_contact_manager.h>
-#include <dem/dem_solver_parameters.h>
-#include <dem/find_cell_neighbors.h>
-#include <dem/particle_particle_broad_search.h>
 #include <dem/particle_particle_contact_force.h>
-#include <dem/particle_particle_fine_search.h>
 
 
 // Tests (with common definitions)
@@ -92,14 +88,12 @@ test()
 
   // Creating containers manager for finding cell neighbor and also broad and
   // fine particle-particle search objects
-  DEMContactManager<dim> container_manager;
+  DEMContactManager<dim> contact_manager;
 
   // Finding cell neighbors
-  FindCellNeighbors<dim> cell_neighbor_object;
-  cell_neighbor_object.find_cell_neighbors(
-    triangulation,
-    container_manager.cells_local_neighbor_list,
-    container_manager.cells_ghost_neighbor_list);
+  find_cell_neighbors<dim>(triangulation,
+                           contact_manager.cells_local_neighbor_list,
+                           contact_manager.cells_ghost_neighbor_list);
 
   // Inserting two particles in contact
   Point<3> position1 = {0.4, 0, 0};
@@ -155,18 +149,17 @@ test()
        particle_iterator != particle_handler.end();
        ++particle_iterator)
     {
-      container_manager.particle_container[particle_iterator->get_id()] =
+      contact_manager.particle_container[particle_iterator->get_id()] =
         particle_iterator;
     }
 
   // Dummy Adaptive sparse contacts object and particle-particle broad search
   AdaptiveSparseContacts<dim> dummy_adaptive_sparse_contacts;
-  container_manager.execute_particle_particle_broad_search(
+  contact_manager.execute_particle_particle_broad_search(
     particle_handler, dummy_adaptive_sparse_contacts);
 
   // Calling fine search
-  container_manager.execute_particle_particle_fine_search(
-    neighborhood_threshold);
+  contact_manager.execute_particle_particle_fine_search(neighborhood_threshold);
 
   // Calling linear force
   ParticleParticleContactForce<
@@ -174,8 +167,10 @@ test()
     Parameters::Lagrangian::ParticleParticleContactForceModel::linear,
     Parameters::Lagrangian::RollingResistanceMethod::constant_resistance>
     linear_force_object(dem_parameters);
-  linear_force_object.calculate_particle_particle_contact_force(
-    container_manager, dt, torque, force);
+  linear_force_object.calculate_particle_particle_contact_force(contact_manager,
+                                                                dt,
+                                                                torque,
+                                                                force);
 
   // Output
   auto particle = particle_handler.begin();
diff --git a/tests/dem/particle_particle_contact_force_nonlinear.cc b/tests/dem/particle_particle_contact_force_nonlinear.cc
index 2e3b22a6b9..17aa7322ab 100644
--- a/tests/dem/particle_particle_contact_force_nonlinear.cc
+++ b/tests/dem/particle_particle_contact_force_nonlinear.cc
@@ -36,11 +36,7 @@
 #include <core/dem_properties.h>
 
 #include <dem/dem_contact_manager.h>
-#include <dem/dem_solver_parameters.h>
-#include <dem/find_cell_neighbors.h>
-#include <dem/particle_particle_broad_search.h>
 #include <dem/particle_particle_contact_force.h>
-#include <dem/particle_particle_fine_search.h>
 
 // Tests (with common definitions)
 #include <../tests/tests.h>
@@ -91,14 +87,12 @@ test()
 
   // Creating containers manager for finding cell neighbor and also broad and
   // fine particle-particle search objects
-  DEMContactManager<dim> container_manager;
+  DEMContactManager<dim> contact_manager;
 
   // Finding cell neighbors
-  FindCellNeighbors<dim> cell_neighbor_object;
-  cell_neighbor_object.find_cell_neighbors(
-    triangulation,
-    container_manager.cells_local_neighbor_list,
-    container_manager.cells_ghost_neighbor_list);
+  find_cell_neighbors<dim>(triangulation,
+                           contact_manager.cells_local_neighbor_list,
+                           contact_manager.cells_ghost_neighbor_list);
 
 
   // Inserting two particles in contact
@@ -154,18 +148,17 @@ test()
        particle_iterator != particle_handler.end();
        ++particle_iterator)
     {
-      container_manager.particle_container[particle_iterator->get_id()] =
+      contact_manager.particle_container[particle_iterator->get_id()] =
         particle_iterator;
     }
 
   // Dummy Adaptive sparse contacts object and particle-particle broad search
   AdaptiveSparseContacts<dim> dummy_adaptive_sparse_contacts;
-  container_manager.execute_particle_particle_broad_search(
+  contact_manager.execute_particle_particle_broad_search(
     particle_handler, dummy_adaptive_sparse_contacts);
 
   // Calling fine search
-  container_manager.execute_particle_particle_fine_search(
-    neighborhood_threshold);
+  contact_manager.execute_particle_particle_fine_search(neighborhood_threshold);
 
   // Calling linear force
   ParticleParticleContactForce<
@@ -175,7 +168,7 @@ test()
     Parameters::Lagrangian::RollingResistanceMethod::constant_resistance>
     nonlinear_force_object(dem_parameters);
   nonlinear_force_object.calculate_particle_particle_contact_force(
-    container_manager, dt, torque, force);
+    contact_manager, dt, torque, force);
 
   // Output
   auto particle = particle_handler.begin();
diff --git a/tests/dem/particle_particle_contact_on_two_processors.cc b/tests/dem/particle_particle_contact_on_two_processors.cc
index 896f5aaeb8..33b34c08d3 100644
--- a/tests/dem/particle_particle_contact_on_two_processors.cc
+++ b/tests/dem/particle_particle_contact_on_two_processors.cc
@@ -117,23 +117,15 @@ test()
 
   typename dem_data_structures<2>::particle_index_iterator_map
     local_particle_container;
-  typename dem_data_structures<2>::adjacent_particle_pairs
-    cleared_local_adjacent_particles;
-  typename dem_data_structures<2>::adjacent_particle_pairs
-    cleared_ghost_adjacent_particles;
 
-  DEMContactManager<dim> container_manager;
+  DEMContactManager<dim> contact_manager;
 
   // Finding cell neighbors
-  FindCellNeighbors<dim> cell_neighbor_object;
-  cell_neighbor_object.find_cell_neighbors(
-    triangulation,
-    container_manager.cells_local_neighbor_list,
-    container_manager.cells_ghost_neighbor_list);
-
-  // Creating broad search, fine search and particle-particle force objects
-  ParticleParticleBroadSearch<dim> broad_search_object;
-  ParticleParticleFineSearch<dim>  fine_search_object;
+  find_cell_neighbors<dim>(triangulation,
+                           contact_manager.cells_local_neighbor_list,
+                           contact_manager.cells_ghost_neighbor_list);
+
+  // Creating particle-particle force objects
   ParticleParticleContactForce<
     dim,
     Parameters::Lagrangian::ParticleParticleContactForceModel::
@@ -210,27 +202,27 @@ test()
 
       particle_handler.exchange_ghost_particles();
 
-      container_manager.update_local_particles_in_cells(particle_handler);
+      contact_manager.update_local_particles_in_cells(particle_handler);
 
       // Dummy Adaptive sparse contacts object and particle-particle broad
       // search
       AdaptiveSparseContacts<dim> dummy_adaptive_sparse_contacts;
-      container_manager.execute_particle_particle_broad_search(
+      contact_manager.execute_particle_particle_broad_search(
         particle_handler, dummy_adaptive_sparse_contacts);
 
       // Calling fine search
-      container_manager.execute_particle_particle_fine_search(
+      contact_manager.execute_particle_particle_fine_search(
         neighborhood_threshold);
 
       // Integration
       // Calling non-linear force
       nonlinear_force_object.calculate_particle_particle_contact_force(
-        container_manager, dt, torque, force);
+        contact_manager, dt, torque, force);
 
       // Integration
       integrator_object.integrate(particle_handler, g, dt, torque, force, MOI);
 
-      container_manager.update_contacts();
+      contact_manager.update_contacts();
 
       if (iteration % output_frequency == 0)
         {
diff --git a/tests/dem/particle_particle_fine_search.cc b/tests/dem/particle_particle_fine_search.cc
index 2a46f07823..214220d307 100644
--- a/tests/dem/particle_particle_fine_search.cc
+++ b/tests/dem/particle_particle_fine_search.cc
@@ -33,10 +33,10 @@
 #include <deal.II/particles/particle_iterator.h>
 
 // Lethe
+#include <dem/contact_info.h>
 #include <dem/dem_contact_manager.h>
 #include <dem/find_cell_neighbors.h>
 #include <dem/particle_particle_broad_search.h>
-#include <dem/particle_particle_contact_info.h>
 #include <dem/particle_particle_fine_search.h>
 
 // Tests (with common definitions)
@@ -63,14 +63,12 @@ test()
   Particles::ParticleHandler<dim> particle_handler(
     triangulation, mapping, DEM::get_number_properties());
 
-  DEMContactManager<dim> container_manager;
+  DEMContactManager<dim> contact_manager;
 
   // Finding cell neighbors
-  FindCellNeighbors<dim> cell_neighbor_object;
-  cell_neighbor_object.find_cell_neighbors(
-    triangulation,
-    container_manager.cells_local_neighbor_list,
-    container_manager.cells_ghost_neighbor_list);
+  find_cell_neighbors<dim>(triangulation,
+                           contact_manager.cells_local_neighbor_list,
+                           contact_manager.cells_ghost_neighbor_list);
 
   // Inserting two particles in contact
   Point<3> position1 = {0.4, 0, 0};
@@ -119,24 +117,23 @@ test()
        particle_iterator != particle_handler.end();
        ++particle_iterator)
     {
-      container_manager.particle_container[particle_iterator->get_id()] =
+      contact_manager.particle_container[particle_iterator->get_id()] =
         particle_iterator;
     }
 
   // Dummy Adaptive sparse contacts object and particle-particle broad search
   AdaptiveSparseContacts<dim> dummy_adaptive_sparse_contacts;
-  container_manager.execute_particle_particle_broad_search(
+  contact_manager.execute_particle_particle_broad_search(
     particle_handler, dummy_adaptive_sparse_contacts);
 
   // Calling fine search
-  container_manager.execute_particle_particle_fine_search(
-    neighborhood_threshold);
+  contact_manager.execute_particle_particle_fine_search(neighborhood_threshold);
 
   // Output
   for (auto adjacent_particles_iterator =
-         container_manager.local_adjacent_particles.begin();
+         contact_manager.local_adjacent_particles.begin();
        adjacent_particles_iterator !=
-       container_manager.local_adjacent_particles.end();
+       contact_manager.local_adjacent_particles.end();
        ++adjacent_particles_iterator)
     {
       auto information_iterator =
diff --git a/tests/dem/particle_point_contact.cc b/tests/dem/particle_point_contact.cc
index 1eeb01953a..82dd99a057 100644
--- a/tests/dem/particle_point_contact.cc
+++ b/tests/dem/particle_point_contact.cc
@@ -36,6 +36,8 @@
 // Lethe
 #include <core/dem_properties.h>
 
+#include <dem/contact_info.h>
+#include <dem/contact_type.h>
 #include <dem/dem_solver_parameters.h>
 #include <dem/find_boundary_cells_information.h>
 #include <dem/particle_point_line_broad_search.h>
@@ -123,14 +125,11 @@ test()
                        Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
 
   // Particle-point broad search
-  std::unordered_map<unsigned int,
-                     std::pair<Particles::ParticleIterator<dim>, Point<dim>>>
-                                    contact_candidates;
-  ParticlePointLineBroadSearch<dim> broad_search_object;
+  typename DEM::dem_data_structures<dim>::particle_point_candidates
+    contact_candidates;
 
   // Particle-point fine search
-  ParticlePointLineFineSearch<dim> fine_search_object;
-  std::unordered_map<unsigned int, particle_point_line_contact_info_struct<dim>>
+  typename DEM::dem_data_structures<dim>::particle_point_in_contact
     contact_information;
 
   ParticlePointLineForce<dim>   force_object;
@@ -156,14 +155,14 @@ test()
       force[particle->get_id()][1] = 0;
       force[particle->get_id()][2] = 0;
 
-      contact_candidates =
-        broad_search_object.find_particle_point_contact_pairs(
-          particle_handler,
-          boundary_cells_object.get_boundary_cells_with_points());
+      find_particle_point_contact_pairs<dim>(
+        particle_handler,
+        boundary_cells_object.get_boundary_cells_with_points(),
+        contact_candidates);
 
-      contact_information =
-        fine_search_object.particle_point_fine_search(contact_candidates,
-                                                      neighborhood_threshold);
+      particle_point_fine_search<dim>(contact_candidates,
+                                      neighborhood_threshold,
+                                      contact_information);
 
       force_object.calculate_particle_point_contact_force(
         &contact_information,
diff --git a/tests/dem/particle_wall_contact_force_linear.cc b/tests/dem/particle_wall_contact_force_linear.cc
index 4004eca209..da3603d9e3 100644
--- a/tests/dem/particle_wall_contact_force_linear.cc
+++ b/tests/dem/particle_wall_contact_force_linear.cc
@@ -148,27 +148,24 @@ test()
     ConditionalOStream(std::cout,
                        Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
 
-  // Calling broad search
-  ParticleWallBroadSearch<dim> broad_search_object;
+
   // P-W broad search
-  ParticleWallBroadSearch<dim> particle_wall_broad_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_candidates
     particle_wall_contact_list;
-  particle_wall_broad_search_object.find_particle_wall_contact_pairs(
+  find_particle_wall_contact_pairs<dim>(
     boundary_cells_object.get_boundary_cells_information(),
     particle_handler,
     particle_wall_contact_list);
-  broad_search_object.find_particle_wall_contact_pairs(
+  find_particle_wall_contact_pairs<dim>(
     boundary_cells_object.get_boundary_cells_information(),
     particle_handler,
     particle_wall_contact_list);
 
   // Calling fine search
-  ParticleWallFineSearch<dim> fine_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_in_contact
     particle_wall_contact_information;
-  fine_search_object.particle_wall_fine_search(
-    particle_wall_contact_list, particle_wall_contact_information);
+  particle_wall_fine_search<dim>(particle_wall_contact_list,
+                                 particle_wall_contact_information);
 
   // Calling linear force
   ParticleWallLinearForce<dim> force_object(dem_parameters);
diff --git a/tests/dem/particle_wall_contact_force_nonlinear.cc b/tests/dem/particle_wall_contact_force_nonlinear.cc
index 4ab8674166..02dad5b8f7 100644
--- a/tests/dem/particle_wall_contact_force_nonlinear.cc
+++ b/tests/dem/particle_wall_contact_force_nonlinear.cc
@@ -149,27 +149,23 @@ test()
     ConditionalOStream(std::cout,
                        Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
 
-  // Calling broad search
-  ParticleWallBroadSearch<dim> broad_search_object;
   // P-W broad search
-  ParticleWallBroadSearch<dim> particle_wall_broad_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_candidates
     particle_wall_contact_list;
-  particle_wall_broad_search_object.find_particle_wall_contact_pairs(
+  find_particle_wall_contact_pairs<dim>(
     boundary_cells_object.get_boundary_cells_information(),
     particle_handler,
     particle_wall_contact_list);
-  broad_search_object.find_particle_wall_contact_pairs(
+  find_particle_wall_contact_pairs<dim>(
     boundary_cells_object.get_boundary_cells_information(),
     particle_handler,
     particle_wall_contact_list);
 
   // Calling fine search
-  ParticleWallFineSearch<dim> fine_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_in_contact
     particle_wall_contact_information;
-  fine_search_object.particle_wall_fine_search(
-    particle_wall_contact_list, particle_wall_contact_information);
+  particle_wall_fine_search<dim>(particle_wall_contact_list,
+                                 particle_wall_contact_information);
 
   // Calling non-linear force
   ParticleWallNonLinearForce<dim> force_object(dem_parameters);
diff --git a/tests/dem/particle_wall_contact_pairs.cc b/tests/dem/particle_wall_contact_pairs.cc
index 487d83fa4f..a93460466f 100644
--- a/tests/dem/particle_wall_contact_pairs.cc
+++ b/tests/dem/particle_wall_contact_pairs.cc
@@ -99,17 +99,14 @@ test()
     ConditionalOStream(std::cout,
                        Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
 
-  // Calling particle-wall broad search
-  ParticleWallBroadSearch<dim> broad_search_object;
   // P-W broad search
-  ParticleWallBroadSearch<dim> particle_wall_broad_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_candidates
     particle_wall_contact_list;
-  particle_wall_broad_search_object.find_particle_wall_contact_pairs(
+  find_particle_wall_contact_pairs<dim>(
     boundary_cells_object.get_boundary_cells_information(),
     particle_handler,
     particle_wall_contact_list);
-  broad_search_object.find_particle_wall_contact_pairs(
+  find_particle_wall_contact_pairs<dim>(
     boundary_cells_object.get_boundary_cells_information(),
     particle_handler,
     particle_wall_contact_list);
diff --git a/tests/dem/particle_wall_fine_search.cc b/tests/dem/particle_wall_fine_search.cc
index afc2f65b34..478bef66b5 100644
--- a/tests/dem/particle_wall_fine_search.cc
+++ b/tests/dem/particle_wall_fine_search.cc
@@ -36,7 +36,6 @@
 // Lethe
 #include <dem/find_boundary_cells_information.h>
 #include <dem/particle_wall_broad_search.h>
-#include <dem/particle_wall_contact_info.h>
 #include <dem/particle_wall_fine_search.h>
 
 // Tests (with common definitions)
@@ -95,20 +94,18 @@ test()
                        Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
 
   // Calling particle-wall broad search
-  ParticleWallBroadSearch<dim> broad_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_candidates
     particle_wall_contact_list;
-  broad_search_object.find_particle_wall_contact_pairs(
+  find_particle_wall_contact_pairs<dim>(
     boundary_cells_object.get_boundary_cells_information(),
     particle_handler,
     particle_wall_contact_list);
 
   // Calling particle-wall fine search
-  ParticleWallFineSearch<dim> fine_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_in_contact
     particle_wall_contact_information;
-  fine_search_object.particle_wall_fine_search(
-    particle_wall_contact_list, particle_wall_contact_information);
+  particle_wall_fine_search<dim>(particle_wall_contact_list,
+                                 particle_wall_contact_information);
 
   // Output
   for (auto particle_wall_contact_information_iterator =
diff --git a/tests/dem/post_collision_velocity.cc b/tests/dem/post_collision_velocity.cc
index 3c222b91ee..023e7057ef 100644
--- a/tests/dem/post_collision_velocity.cc
+++ b/tests/dem/post_collision_velocity.cc
@@ -154,16 +154,14 @@ test(double coefficient_of_restitution)
                        Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
 
   // P-W broad search
-  ParticleWallBroadSearch<dim> particle_wall_broad_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_candidates
     particle_wall_contact_list;
-  particle_wall_broad_search_object.find_particle_wall_contact_pairs(
+  find_particle_wall_contact_pairs<dim>(
     boundary_cells_object.get_boundary_cells_information(),
     particle_handler,
     particle_wall_contact_list);
 
   // P-W fine search
-  ParticleWallFineSearch<dim> particle_wall_fine_search_object;
   typename DEM::dem_data_structures<dim>::particle_wall_in_contact
                                   particle_wall_contact_information;
   ParticleWallNonLinearForce<dim> particle_wall_force_object(dem_parameters);
@@ -174,8 +172,8 @@ test(double coefficient_of_restitution)
   for (double time = 0; time < 0.0001; time += dt)
     {
       // If particle and wall are in contact
-      particle_wall_fine_search_object.particle_wall_fine_search(
-        particle_wall_contact_list, particle_wall_contact_information);
+      particle_wall_fine_search<dim>(particle_wall_contact_list,
+                                     particle_wall_contact_information);
 
       particle_wall_force_object.calculate_particle_wall_contact_force(
         particle_wall_contact_information, dt, torque, force);
diff --git a/tests/dem/two_particles_multiple_contacts_parallel.cc b/tests/dem/two_particles_multiple_contacts_parallel.cc
index e7a6035a3e..1ab442137a 100644
--- a/tests/dem/two_particles_multiple_contacts_parallel.cc
+++ b/tests/dem/two_particles_multiple_contacts_parallel.cc
@@ -93,23 +93,15 @@ test()
 
   typename dem_data_structures<2>::particle_index_iterator_map
     local_particle_container;
-  typename dem_data_structures<2>::adjacent_particle_pairs
-    cleared_local_adjacent_particles;
-  typename dem_data_structures<2>::adjacent_particle_pairs
-    cleared_ghost_adjacent_particles;
 
-  DEMContactManager<dim> container_manager;
+  DEMContactManager<dim> contact_manager;
 
   // Finding cell neighbors
-  FindCellNeighbors<dim> cell_neighbor_object;
-  cell_neighbor_object.find_cell_neighbors(
-    triangulation,
-    container_manager.cells_local_neighbor_list,
-    container_manager.cells_ghost_neighbor_list);
-
-  // Creating broad search, fine search and particle-particle force objects
-  ParticleParticleBroadSearch<dim> broad_search_object;
-  ParticleParticleFineSearch<dim>  fine_search_object;
+  find_cell_neighbors<dim>(triangulation,
+                           contact_manager.cells_local_neighbor_list,
+                           contact_manager.cells_ghost_neighbor_list);
+
+  // Particle-particle force objects
   ParticleParticleContactForce<
     dim,
     Parameters::Lagrangian::ParticleParticleContactForceModel::
@@ -186,22 +178,22 @@ test()
 
       particle_handler.exchange_ghost_particles();
 
-      container_manager.update_local_particles_in_cells(particle_handler);
+      contact_manager.update_local_particles_in_cells(particle_handler);
 
       // Dummy Adaptive sparse contacts object and particle-particle broad
       // search
       AdaptiveSparseContacts<dim> dummy_adaptive_sparse_contacts;
-      container_manager.execute_particle_particle_broad_search(
+      contact_manager.execute_particle_particle_broad_search(
         particle_handler, dummy_adaptive_sparse_contacts);
 
       // Calling fine search
-      container_manager.execute_particle_particle_fine_search(
+      contact_manager.execute_particle_particle_fine_search(
         neighborhood_threshold);
 
       // Integration
       // Calling non-linear force
       nonlinear_force_object.calculate_particle_particle_contact_force(
-        container_manager, dt, torque, force);
+        contact_manager, dt, torque, force);
 
       // Store force before integration for proc 1
       // TODO - Improve this in the future, this is not clean.
@@ -211,7 +203,7 @@ test()
       // Integration
       integrator_object.integrate(particle_handler, g, dt, torque, force, MOI);
 
-      container_manager.update_contacts();
+      contact_manager.update_contacts();
 
       if (iteration % output_frequency == 0)
         {