Merge pull request #12539 from KratosMultiphysics/dem/improvement_part1

[DEMApplication] Added new features and fixed a bug related to the periodic boundary [improvement part 1]

applications/DEMApplication/CMakeLists.txt

@@ -92,6 +92,7 @@ set(KRATOS_DEM_APPLICATION_CORE
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_smooth_joint_CL.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_constitutive/DEM_parallel_bond_for_membrane_CL.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/inlet.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/fast_filling_creator.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/force_based_inlet.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/create_and_destroy.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/properties_proxies.cpp

applications/DEMApplication/custom_python/add_custom_utilities_to_python.cpp

@@ -25,6 +25,7 @@
 #include "custom_utilities/dem_fem_search.h"
 #include "custom_utilities/dem_fem_utilities.h"
 #include "custom_utilities/inlet.h"
+#include "custom_utilities/fast_filling_creator.h"
 #include "custom_utilities/force_based_inlet.h"
 #include "custom_utilities/reorder_consecutive_from_given_ids_model_part_io.h"
 #include "custom_utilities/AuxiliaryUtilities.h"
@@ -204,6 +205,13 @@ void AddCustomUtilitiesToPython(pybind11::module& m) {
         .def("GetMaxRadius", &DEM_Inlet::GetMaxRadius)
         ;
 
+    py::class_<Fast_Filling_Creator, Fast_Filling_Creator::Pointer>(m, "Fast_Filling_Creator")
+        .def(py::init<const int>())
+        .def(py::init<Parameters&, const int>())
+        .def("GetRandomParticleRadius", &Fast_Filling_Creator::GetRandomParticleRadius, py::arg("creator_destructor"))
+        .def("CheckHasIndentationOrNot", &Fast_Filling_Creator::CheckHasIndentationOrNot)
+        ;
+
     py::class_<DEM_Force_Based_Inlet, DEM_Force_Based_Inlet::Pointer, DEM_Inlet>(m, "DEM_Force_Based_Inlet")
         .def(py::init<ModelPart&, array_1d<double, 3>, const int>())
         .def(py::init<ModelPart&, array_1d<double, 3>>())

applications/DEMApplication/custom_strategies/strategies/explicit_solver_continuum.cpp

@@ -654,9 +654,19 @@ namespace Kratos {
     void ContinuumExplicitSolverStrategy::SetSearchRadiiOnAllParticles(ModelPart& r_model_part, const double added_search_distance, const double amplification) {
         KRATOS_TRY
         const int number_of_elements = r_model_part.GetCommunicator().LocalMesh().NumberOfElements();
-        #pragma omp parallel for
-        for (int i = 0; i < number_of_elements; i++) {
-            mListOfSphericContinuumParticles[i]->SetSearchRadius(amplification * mListOfSphericContinuumParticles[i]->mLocalRadiusAmplificationFactor * (added_search_distance + mListOfSphericContinuumParticles[i]->GetRadius()));
+        if (GetDeltaOption() == 3){
+            // In this case, the parameter "added_search_distance" is actually a multiplier for getting the added_search_distance
+            const double search_radius_multiplier = added_search_distance;
+            #pragma omp parallel for
+            for (int i = 0; i < number_of_elements; i++) {
+                mListOfSphericContinuumParticles[i]->SetSearchRadius(amplification * mListOfSphericContinuumParticles[i]->mLocalRadiusAmplificationFactor * ((1 + search_radius_multiplier) * mListOfSphericContinuumParticles[i]->GetRadius()));
+            }
+        }
+        else{
+            #pragma omp parallel for
+            for (int i = 0; i < number_of_elements; i++) {
+                mListOfSphericContinuumParticles[i]->SetSearchRadius(amplification * mListOfSphericContinuumParticles[i]->mLocalRadiusAmplificationFactor * (added_search_distance + mListOfSphericContinuumParticles[i]->GetRadius()));
+            }
         }
         KRATOS_CATCH("")
     }

applications/DEMApplication/custom_strategies/strategies/explicit_solver_strategy.cpp

@@ -1291,9 +1291,18 @@ namespace Kratos {
         KRATOS_TRY
 
         int number_of_elements = r_model_part.GetCommunicator().LocalMesh().ElementsArray().end() - r_model_part.GetCommunicator().LocalMesh().ElementsArray().begin();
-        IndexPartition<unsigned int>(number_of_elements).for_each([&](unsigned int i) {
-            mListOfSphericParticles[i]->SetSearchRadius(amplification * (added_search_distance + mListOfSphericParticles[i]->GetRadius()));
-        });
+        if (GetDeltaOption() == 3){
+            // In this case, the parameter "added_search_distance" is actually a multiplier for getting the added_search_distance
+            const double search_radius_multiplier = added_search_distance;
+            IndexPartition<unsigned int>(number_of_elements).for_each([&](unsigned int i) {
+                mListOfSphericParticles[i]->SetSearchRadius(amplification * ((1 + search_radius_multiplier) *  mListOfSphericParticles[i]->GetRadius()));
+            });
+        }
+        else{
+            IndexPartition<unsigned int>(number_of_elements).for_each([&](unsigned int i) {
+                mListOfSphericParticles[i]->SetSearchRadius(amplification * (added_search_distance + mListOfSphericParticles[i]->GetRadius()));
+            });
+        }
 
         KRATOS_CATCH("")
     }

applications/DEMApplication/custom_utilities/create_and_destroy.cpp

@@ -332,6 +332,27 @@ namespace Kratos {
         KRATOS_CATCH("")
     }
 
+    double ParticleCreatorDestructor::SelectRadius(Parameters r_sub_model_part_with_parameters,
+                                                   std::map<std::string, std::unique_ptr<RandomVariable>>& r_random_variables_map){
+
+        KRATOS_TRY
+
+        double radius = r_sub_model_part_with_parameters["RADIUS"].GetDouble();
+        const double& max_radius = r_sub_model_part_with_parameters["MAXIMUM_RADIUS"].GetDouble();
+        const std::string& distribution_type = r_sub_model_part_with_parameters["PROBABILITY_DISTRIBUTION"].GetString();
+
+        const double& std_deviation = r_sub_model_part_with_parameters["STANDARD_DEVIATION"].GetDouble();
+        const double& min_radius = r_sub_model_part_with_parameters["MINIMUM_RADIUS"].GetDouble();
+
+        if (distribution_type == "normal") radius = rand_normal(radius, std_deviation, max_radius, min_radius);
+        else if (distribution_type == "lognormal") radius = rand_lognormal(radius, std_deviation, max_radius, min_radius);
+        else if (distribution_type == "piecewise_linear" || distribution_type == "discrete") radius = r_random_variables_map[r_sub_model_part_with_parameters["NAME"].GetString()]->Sample();
+        else KRATOS_ERROR << "Unknown probability distribution in submodelpart " << r_sub_model_part_with_parameters["NAME"].GetString() << std::endl;
+
+        return radius;
+        KRATOS_CATCH("")
+    }
+
     SphericParticle* ParticleCreatorDestructor::ElementCreatorWithPhysicalParameters(ModelPart& r_modelpart,
                                                                         int r_Elem_Id,
                                                                         Node ::Pointer reference_node,

applications/DEMApplication/custom_utilities/create_and_destroy.h

@@ -68,6 +68,8 @@ friend class ExplicitSolverStrategy;
     virtual double SelectRadius(bool initial,
                                 ModelPart& r_sub_model_part_with_parameters,
                                 std::map<std::string, std::unique_ptr<RandomVariable>>& r_random_variables_map);
+    virtual double SelectRadius(Parameters r_sub_model_part_with_parameters,
+                                std::map<std::string, std::unique_ptr<RandomVariable>>& r_random_variables_map);
 
     void NodeCreatorWithPhysicalParameters(ModelPart& r_modelpart,
                                            Node ::Pointer& pnew_node,

applications/DEMApplication/custom_utilities/fast_filling_creator.cpp

@@ -0,0 +1,83 @@
+/////////////////////////////////////////////////
+// Main author: Chengshun Shang (CIMNE)
+// Email: [email protected]
+// Date: Sep 2023
+/////////////////////////////////////////////////
+
+#include <string>
+#include <iostream>
+#include <random>
+#include <functional>
+
+#include "fast_filling_creator.h"
+#include "create_and_destroy.h"
+#include "custom_elements/spheric_continuum_particle.h"
+#include "custom_elements/cluster3D.h"
+#include "custom_constitutive/DEM_discontinuum_constitutive_law.h"
+#include "custom_constitutive/DEM_continuum_constitutive_law.h"
+#include "dem_fem_utilities.h"
+#include "GeometryFunctions.h"
+
+
+namespace Kratos {
+
+    /// Constructor
+
+    Fast_Filling_Creator::Fast_Filling_Creator(const int seed):
+                        Fast_Filling_Creator(Parameters(R"({})"), seed){}
+
+    Fast_Filling_Creator::Fast_Filling_Creator(const Parameters& r_fast_filling_creator_settings, const int seed):
+     mFastFillingCreatorSettings(Parameters(r_fast_filling_creator_settings))
+        {
+        std::mt19937 gen(seed);
+        mGenerator = gen;
+
+        mTotalNumberOfParticlesInjected = 0;
+        mTotalMassInjected = 0.0;
+    }
+
+    double Fast_Filling_Creator::GetRandomParticleRadius(ParticleCreatorDestructor& creator) {
+
+        KRATOS_TRY
+
+        if (mFastFillingCreatorSettings["PROBABILITY_DISTRIBUTION"].GetString() == "piecewise_linear" || mFastFillingCreatorSettings["PROBABILITY_DISTRIBUTION"].GetString() == "discrete"){
+
+            const Parameters& rv_settings = mFastFillingCreatorSettings["random_variable_settings"];
+            int seed = rv_settings["seed"].GetInt();
+            if (!rv_settings["do_use_seed"].GetBool()){
+                seed = std::random_device{}();
+            }
+            if (mFastFillingCreatorSettings["PROBABILITY_DISTRIBUTION"].GetString() == "piecewise_linear"){
+
+                mInletsRandomVariables[mFastFillingCreatorSettings["NAME"].GetString()] = std::unique_ptr<PiecewiseLinearRandomVariable>(new PiecewiseLinearRandomVariable(rv_settings, seed));
+            }
+
+            else if (mFastFillingCreatorSettings["PROBABILITY_DISTRIBUTION"].GetString() == "discrete"){
+                mInletsRandomVariables[mFastFillingCreatorSettings["NAME"].GetString()] = std::unique_ptr<DiscreteRandomVariable>(new DiscreteRandomVariable(rv_settings, seed));
+            }
+
+            else {
+                KRATOS_ERROR << "Unknown Fast Filling Creator random variable: " << mFastFillingCreatorSettings["PROBABILITY_DISTRIBUTION"].GetString() << ".";
+            }
+        }
+
+        double radius = creator.SelectRadius(mFastFillingCreatorSettings, mInletsRandomVariables);
+
+        return radius;
+
+        KRATOS_CATCH("")
+    }
+
+    bool Fast_Filling_Creator::CheckHasIndentationOrNot(const double x_1, const double y_1, const double z_1, const double r_1,
+                                              const double x_2, const double y_2, const double z_2, const double r_2){
+
+        KRATOS_TRY
+
+        double distance = std::sqrt(std::pow(x_1 - x_2, 2) + std::pow(y_1 - y_2, 2) + std::pow(z_1 - z_2, 2));
+
+        return distance <= (r_1 + r_2); 
+
+        KRATOS_CATCH("")                                    
+    }
+
+} // namespace Kratos

applications/DEMApplication/custom_utilities/fast_filling_creator.h

@@ -0,0 +1,80 @@
+/////////////////////////////////////////////////
+// Main author: Chengshun Shang (CIMNE)
+// Email: [email protected]
+// Date: Sep 2023
+/////////////////////////////////////////////////
+
+#if !defined(FAST_FILLING_PARTICLES_H)
+#define FAST_FILLING_PARTICLES_H
+
+// System includes
+#include <string>
+#include <iostream>
+#include <random>
+
+// Project includes
+#include "includes/define.h"
+#include "includes/variables.h"
+#include "includes/node.h"
+#include "includes/element.h"
+#include "geometries/geometry.h"
+#include "includes/properties.h"
+#include "includes/process_info.h"
+#include "includes/indexed_object.h"
+#include "containers/global_pointers_vector.h"
+#include "includes/constitutive_law.h"
+#include "includes/condition.h"
+#include "custom_elements/discrete_element.h"
+#include "custom_utilities/AuxiliaryFunctions.h"
+#include "custom_utilities/random_variable.h"
+#include "custom_utilities/piecewise_linear_random_variable.h"
+#include "custom_utilities/discrete_random_variable.h"
+#include "custom_utilities/properties_proxies.h"
+#include "custom_elements/spheric_particle.h"
+
+namespace Kratos {
+
+    class ParticleCreatorDestructor;
+
+    class  KRATOS_API(DEM_APPLICATION) Fast_Filling_Creator
+    {
+    public:
+
+        typedef GlobalPointersVector<Element >::iterator ParticleWeakIteratorType;
+        typedef GlobalPointersVector<Element> ParticleWeakVectorType;
+        typedef ModelPart::ElementsContainerType ElementsArrayType;
+
+        KRATOS_CLASS_POINTER_DEFINITION(Fast_Filling_Creator);
+
+        /// Constructor:
+        Fast_Filling_Creator(const int seed=42);
+        Fast_Filling_Creator(const Parameters& r_fast_filling_creator_settings, const int seed=42);
+
+        /// Destructor.
+        virtual ~Fast_Filling_Creator(){}
+
+        Fast_Filling_Creator(const Fast_Filling_Creator&) = delete;
+        Fast_Filling_Creator& operator=(const Fast_Filling_Creator&) = delete;
+
+        virtual double GetRandomParticleRadius(ParticleCreatorDestructor& creator);
+        virtual bool CheckHasIndentationOrNot(const double x_1, const double y_1, const double z_1, const double r_1,
+                                              const double x_2, const double y_2, const double z_2, const double r_2);
+
+    protected:
+
+    private:
+
+        int  mTotalNumberOfParticlesInjected;
+        double mTotalMassInjected;
+        //int mSeed;
+        std::vector<double> mMassInjected;
+        std::mt19937 mGenerator;
+
+        std::map<std::string, std::unique_ptr<RandomVariable>> mInletsRandomVariables;
+        std::map<std::string, Parameters> mInletsRandomSettings;
+        Parameters mFastFillingCreatorSettings;
+    };
+}// namespace Kratos.
+
+#endif // FAST_FILLING_PARTICLES_H defined
+