Finished adding core functions into new class

include/fem-dem/cfd_dem_simulation_parameters.h

@@ -39,7 +39,7 @@ class CFDDEMSimulationParameters
   DEMSolverParameters<dim>  dem_parameters;
 
   std::shared_ptr<Parameters::VoidFractionParameters<dim>> void_fraction;
-  Parameters::CFDDEM                             cfd_dem;
+  Parameters::CFDDEM                                       cfd_dem;
 
   void
   declare(ParameterHandler             &prm,

include/fem-dem/fluid_dynamics_vans.h

@@ -15,8 +15,8 @@
  *
  */
 
-#ifndef lethe_gls_vans_h
-#  define lethe_gls_vans_h
+#ifndef lethe_fluid_dynamics_vans_h
+#define lethe_fluid_dynamics_vans_h
 
 #include <core/bdf.h>
 #include <core/dem_properties.h>
@@ -43,8 +43,6 @@
 #include <boost/archive/text_iarchive.hpp>
 #include <boost/archive/text_oarchive.hpp>
 
-#ifndef lethe_fluid_dynamics_vans_h
-#define lethe_fluid_dynamics_vans_h
 
 using namespace dealii;
 

include/fem-dem/void_fraction.h

@@ -1,23 +1,23 @@
 /* ---------------------------------------------------------------------
-*
-* 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.
-*
-* ---------------------------------------------------------------------
-*
-*/
+ *
+ * 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_void_fraction_h
-#  define lethe_void_fraction_h
+#define lethe_void_fraction_h
 
 #include <core/vector.h>
 
@@ -30,42 +30,129 @@
 #include <deal.II/dofs/dof_handler.h>
 
 #include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/mapping_fe.h>
 
-#include <deal.II/lac/trilinos_precondition.h>
 #include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/trilinos_precondition.h>
 #include <deal.II/lac/trilinos_sparse_matrix.h>
 
+#include <deal.II/particles/particle_handler.h>
+
 using namespace dealii;
 
 /**
- * @brief Base for all of the void fraction calculators. The void fraction is used as an auxiliary fields for the solvers that solve the Volume-Averaged Navier-Stokes equations.
+ * @brief Calculates the area of intersection between a circular (2D) particle and a circle
+ *
+ * @param r_particle Radius of the particle
+ *
+ * @param r_circle Radius of the circle
+ *
+ * @param neighbor_distance Distance between the particle and the circle
+ */
+inline double
+particle_circle_intersection_2d(double r_particle,
+                                double r_circle,
+                                double neighbor_distance)
+{
+  return pow(r_particle, 2) * Utilities::fixed_power<-1, double>(
+                                cos((pow(neighbor_distance, 2) +
+                                     pow(r_particle, 2) - pow(r_circle, 2)) /
+                                    (2 * neighbor_distance * r_particle))) +
+         Utilities::fixed_power<2, double>(r_circle) *
+           Utilities::fixed_power<-1, double>(
+             cos((pow(neighbor_distance, 2) - pow(r_particle, 2) +
+                  pow(r_circle, 2)) /
+                 (2 * neighbor_distance * r_circle))) -
+         0.5 * sqrt((-neighbor_distance + r_particle + r_circle) *
+                    (neighbor_distance + r_particle - r_circle) *
+                    (neighbor_distance - r_particle + r_circle) *
+                    (neighbor_distance + r_particle + r_circle));
+}
+
+/**
+ * @brief Calculates the volume of intersection between a spherical (3D) particle and a sphere
+ *
+ * @param r_particle Radius of the particle
+ *
+ * @param r_sphere Radius of the sphere
+ *
+ * @param neighbor_distance Distance between the particle and the sphere
+ */
+
+inline double
+particle_sphere_intersection_3d(double r_particle,
+                                double r_sphere,
+                                double neighbor_distance)
+{
+  return M_PI *
+         Utilities::fixed_power<2, double>(r_sphere + r_particle -
+                                           neighbor_distance) *
+         (Utilities::fixed_power<2, double>(neighbor_distance) +
+          (2 * neighbor_distance * r_particle) -
+          (3 * Utilities::fixed_power<2, double>(r_particle)) +
+          (2 * neighbor_distance * r_sphere) + (6 * r_sphere * r_particle) -
+          (3 * Utilities::fixed_power<2, double>(r_sphere))) /
+         (12 * neighbor_distance);
+}
+
+/**
+ * @brief Void fraction calculator. This class manages the calculation of the
+ * void fraction which is used as an auxiliary field for the solvers that solve
+ * the Volume-Averaged Navier-Stokes equations. The present architecture
+ * of the class support all of the different void fraction calculation methods
+ * within a single class instead of building a class hierarchy. This is
+ * because there are only a few void fraction calculation strategies.
  *
  * @tparam dim An integer that denotes the number of spatial dimensions.
  */
 template <int dim>
 class VoidFractionBase
 {
-  VoidFractionBase(  parallel::DistributedTriangulationBase<dim> * triangulation,
-                   const Parameters::VoidFractionParameters<dim> & input_parameters, const unsigned int fe_degree):
-    dof_handler(*triangulation),
-    parameters(input_parameters),
-    fe(fe_degree)
-  {;}
+  VoidFractionBase(
+    parallel::DistributedTriangulationBase<dim>   *triangulation,
+    const Parameters::VoidFractionParameters<dim> &input_parameters,
+    const Parameters::LinearSolver                &linear_solver_parameters,
+    const Particles::ParticleHandler<dim>         *particle_handler,
+    const unsigned int                             fe_degree,
+    const bool                                     simplex)
+    : dof_handler(*triangulation)
+    , void_fraction_parameters(input_parameters)
+    , particle_handler(particle_handler)
+  {
+    if (simplex)
+      {
+        fe         = std::make_shared<FE_SimplexP<dim>>(fe_degree);
+        mapping    = std::make_shared<MappingFE<dim>>(*fe);
+        quadrature = std::make_shared<QGaussSimplex<dim>>(fe->degree + 1);
+      }
+    else
+      {
+        // Usual case, for quad/hex meshes
+        fe         = std::make_shared<FE_Q<dim>>(fe_degree);
+        mapping    = std::make_shared<MappingQ<dim>>(fe->degree);
+        quadrature = std::make_shared<QGauss<dim>>(fe->degree + 1);
+      }
+  }
 
   /**
-  * @brief Setup the degrees of freedom and establishes the constraints of the void fraction systems.
+   * @brief Setup the degrees of freedom and establishes the constraints of the void fraction systems.
    *
    */
-  void setup_dofs();
+  void
+  setup_dofs();
+
 
   /**
-  * @brief Assembles the diagonal of the mass matrix to impose constraints on the system
-  *
-  *  @param diagonal_mass_matrix The matrix for which the diagonal entries will be filled with the diagonal of the mass matrix
+   * @brief Assembles the diagonal of the mass matrix to impose constraints on the system
+   *
+   *  @param diagonal_mass_matrix The matrix for which the diagonal entries will be filled with the diagonal of the mass matrix
    *
    *  @todo Establish if it is really necessary to keep this as a matrix and not as a vector since a diagonal is nothing more than a vector.
-  */
-  void assemble_mass_matrix_diagonal(TrilinosWrappers::SparseMatrix &diagonal_mass_matrix);
+   */
+  void
+  assemble_mass_matrix_diagonal(
+    TrilinosWrappers::SparseMatrix &diagonal_mass_matrix);
 
   /**
    * @brief Calculates the void fraction
@@ -78,17 +165,65 @@ class VoidFractionBase
 
 
 
-protected:
+private:
+  /**
+   * @brief Calculates the void fraction using a function. This is a straightforward usage of VectorTools.
+   *
+   * @param time current time for which the void fraction is to be calculated.
+   *
+   */
+  void
+  calculate_void_fraction_function(const double time);
+
+  /**
+   * @brief Calculates the void fraction using the particle centered method.
+   *
+   */
+  void
+  calculate_void_fraction_particle_centered_method();
+
+  /**
+   * @brief Calculates the void fraction using the satellite point method.
+   *
+   */
+  void
+  calculate_void_fraction_satellite_point_method();
+
+  /**
+   * @brief Calculates the void fraction using the Quadrature-Centered Method (QCM).
+   *
+   */
+  void
+  calculate_void_fraction_quadrature_centered_method();
+
+  /**
+   * @brief Solve the linear system resulting from the assemblies
+   *
+   */
+  void
+  solve_linear_system();
+
+  /// Finite element for the void fraction
+  std::shared_ptr<FiniteElement<dim>> fe;
+
+  /// Mapping for the void fraction
+  std::shared_ptr<Mapping<dim>> mapping;
+
+  /// Quadrature for the void fraction
+  std::shared_ptr<Quadrature<dim>> quadrature;
 
   /// Parameters for the calculation of the void fraction
-  const Parameters::VoidFractionParameters<dim> parameters;
+  const Parameters::VoidFractionParameters<dim> void_fraction_parameters;
+
+  /// Linear solvers for the calculation of the void fraction
+  const Parameters::LinearSolver linear_solver_parameters;
+
+  /// Particle handler used when the void fraction depends on particles
+  const Particles::ParticleHandler<dim> *particle_handler;
 
   /// DoFHandler that manages the void fraction
   DoFHandler<dim> dof_handler;
 
-  /// FE for the void fraction. Currently only FE_Q (Lagrange polynomials) are supported.
-  FE_Q<dim>       fe;
-
   /// Index set for the locally owned degree of freedoms
   IndexSet locally_owned_dofs;
 
@@ -108,37 +243,61 @@ class VoidFractionBase
   TrilinosWrappers::SparseMatrix complete_system_matrix_void_fraction;
   GlobalVectorType               complete_system_rhs_void_fraction;
 
-  /// Mass matrix used to constraint the value of the void fraction to be bounded
+  /// Mass matrix used to constraint the value of the void fraction to be
+  /// bounded
   TrilinosWrappers::SparseMatrix mass_matrix;
 
-  /// Mass matrix diagonal used to constraint the value of the void fraction to be bounded
-  GlobalVectorType               diagonal_of_mass_matrix;
+  /// Mass matrix diagonal used to constraint the value of the void fraction to
+  /// be bounded
+  GlobalVectorType diagonal_of_mass_matrix;
 
   /// ??? BB Not fully sure of this yet
-  IndexSet                       active_set;
+  IndexSet active_set;
 
-  /// Preconditioner used for the solution of the smoothed L2 projection of the void fraction
+  /// Preconditioner used for the solution of the smoothed L2 projection of the
+  /// void fraction
   std::shared_ptr<TrilinosWrappers::PreconditionILU> ilu_preconditioner;
 
-  /// Constraints used for the boundary conditions of the void fraction. Currently, this is only used to establish periodic void fractions.
-  AffineConstraints<double>                          void_fraction_constraints;
+  /// Constraints used for the boundary conditions of the void fraction.
+  /// Currently, this is only used to establish periodic void fractions.
+  AffineConstraints<double> void_fraction_constraints;
 
-  /// System matrix used to assembled the smoothed L2 projection of the void fraction
+  /// System matrix used to assembled the smoothed L2 projection of the void
+  /// fraction
   TrilinosWrappers::SparseMatrix system_matrix_void_fraction;
 
-  /// Right-hand side used to assemble the smoothed L2 projection of the void fraction
-  GlobalVectorType               system_rhs_void_fraction;
-};
+  /// Right-hand side used to assemble the smoothed L2 projection of the void
+  /// fraction
+  GlobalVectorType system_rhs_void_fraction;
 
-/**
- * @brief Calculate the void fraction using an analytical function (time and space). Although the void fraction is known analytically, it is still calculated on the grid using a DofHandler to ensure consistency with the underlying finite element scheme.
- *
- * @tparam dim An integer that denotes the number of spatial dimensions.
- */
-template <int dim>
-class VoidFractionFunction : public VoidFractionBase<dim>
-{
+  /// Vertices to cell map, this is used in the QCM and the satellite point
+  /// method to access the neighboring cells of a cell.
+  std::map<unsigned int,
+           std::set<typename DoFHandler<dim>::active_cell_iterator>>
+    vertices_to_cell;
 
+
+  /**
+   * Member related to boundary conditions. At the moment, only a single
+   * boundary condition is supported.
+   */
+
+  /// Boolean to indicate if the mesh has periodic boundaries
+  bool has_periodic_boundaries;
+
+  /// Offset for the periodic boundary condition
+  Tensor<1, dim> periodic_offset;
+
+  /// Direction associated of the periodic boundary condition
+  unsigned int periodic_direction;
+
+  /// Vertices to periodic cell map, this is used in the QCM and the satellite
+  /// point method to access the neighboring cells of a cell.
+  std::map<unsigned int,
+           std::set<typename DoFHandler<dim>::active_cell_iterator>>
+    vertices_to_periodic_cell;
 };
 
+
+
 #endif

source/fem-dem/fluid_dynamics_vans.cc

@@ -247,6 +247,7 @@ FluidDynamicsVANS<dim>::calculate_void_fraction(const double time,
   announce_string(this->pcout, "Void Fraction");
 
   TimerOutput::Scope t(this->computing_timer, "Calculate void fraction");
+
   if (this->cfd_dem_simulation_parameters.void_fraction->mode ==
       Parameters::VoidFractionMode::function)
     {