dem_structures_coupling_utilities.h

Go to the documentation of this file.

/*
 * Author: Miguel Angel Celigueta
 *
 *  maceli@cimne.upc.edu
 */
 
#ifndef KRATOS_STRUCTURES_DEM_COUPLING_UTILITIES_H
#define KRATOS_STRUCTURES_DEM_COUPLING_UTILITIES_H
// /* External includes */
 
// System includes
 
// Project includes
#include "includes/variables.h"
 
/* System includes */
#include <limits>
#include <iostream>
#include <iomanip>
 
/* External includes */
#ifdef _OPENMP
#include <omp.h>
#endif
 
/* Project includes */
#include "includes/define.h"
#include "includes/model_part.h"
#include "custom_conditions/RigidFace.h"
#include "custom_conditions/RigidEdge.h"
#include "DEM_application_variables.h"
#include "dem_structures_coupling_application_variables.h"
#include "custom_elements/spheric_continuum_particle.h"
 
namespace Kratos
{
class DemStructuresCouplingUtilities
{
public:
typedef ModelPart::NodesContainerType::ContainerType::iterator NodesIteratorType;
 
KRATOS_CLASS_POINTER_DEFINITION(DemStructuresCouplingUtilities);
 
DemStructuresCouplingUtilities(){}
 
virtual ~DemStructuresCouplingUtilities(){}
 
//***************************************************************************************************************
//***************************************************************************************************************
 
void TransferStructuresSkinToDem(ModelPart& r_source_model_part, ModelPart& r_destination_model_part, Properties::Pointer props) {
 
    // std::string error = CheckProvidedProperties(props);
    const int dimension = r_source_model_part.GetProcessInfo()[DOMAIN_SIZE];
 
    // if (error != "all_ok") KRATOS_ERROR << "The Dem Walls ModelPart has no valid Properties. Missing " << error << " . Exiting." << std::endl;
 
    r_destination_model_part.Conditions().Sort();
    int id = 1;
 
    if (r_destination_model_part.Conditions().size()) id = (r_destination_model_part.ConditionsEnd()-1)->Id() + 1;
 
    ModelPart::ConditionsContainerType& source_conditions = r_source_model_part.Conditions();
 
    // Adding conditions
    for (unsigned int i = 0; i < source_conditions.size(); i++) {
        ModelPart::ConditionsContainerType::iterator it = r_source_model_part.ConditionsBegin() + i;
        Geometry< Node >::Pointer p_geometry =  it->pGetGeometry();
        Condition::Pointer cond;
        if (dimension == 2) {
            cond = Condition::Pointer(new RigidEdge2D(id, p_geometry, props));
        } else {
            cond = Condition::Pointer(new RigidFace3D(id, p_geometry, props));
        }
 
        cond->Set(DEMFlags::STICKY, true);
        r_destination_model_part.AddCondition(cond); //TODO: add all of them in a single sentence! AddConditions. Use a temporary PointerVector as a list (not std::vector!).
        id++;
    }
 
    // Adding nodes
    r_destination_model_part.AddNodes(r_source_model_part.NodesBegin(), r_source_model_part.NodesEnd());
}
 
std::string CheckProvidedProperties(Properties::Pointer props) {
    std::vector<const Variable<double>* > list_of_variables_double_to_check = {&STATIC_FRICTION, &DYNAMIC_FRICTION, &FRICTION_DECAY, &WALL_COHESION, &SEVERITY_OF_WEAR, &IMPACT_WEAR_SEVERITY, &BRINELL_HARDNESS, &YOUNG_MODULUS, &POISSON_RATIO};
    std::vector<const Variable<bool>* > list_of_variables_bool_to_check = {&COMPUTE_WEAR};
    for (int i=0; i<(int)list_of_variables_double_to_check.size(); i++) {
        if(!props->Has(*list_of_variables_double_to_check[i])) return list_of_variables_double_to_check[i]->Name();
    }
    for (int i=0; i<(int)list_of_variables_bool_to_check.size(); i++) {
        if(!props->Has(*list_of_variables_bool_to_check[i])) return list_of_variables_bool_to_check[i]->Name();
    }
    return "all_ok";
}
 
void SmoothLoadTrasferredToFem(ModelPart& r_model_part, const double portion_of_the_force_which_is_new) {
    #pragma omp parallel for
    for (int i=0; i<(int)r_model_part.Nodes().size(); i++) {
        auto node_it = r_model_part.NodesBegin() + i;
        array_1d<double, 3> averaged_force;
        array_1d<double, 3>& node_dem_load = node_it->FastGetSolutionStepValue(DEM_SURFACE_LOAD);
        noalias(averaged_force) = portion_of_the_force_which_is_new * node_dem_load + (1.0 - portion_of_the_force_which_is_new) * node_it->FastGetSolutionStepValue(DEM_SURFACE_LOAD, 1);
        noalias(node_dem_load) = averaged_force;
    }
}
 
double ComputeSandProduction(ModelPart& dem_model_part, ModelPart& outer_walls_model_part, const double time) {
 
    const std::string sand_prod_filename = "sand_production_graph.txt";
    static std::ofstream ofs_sand_prod_file;
    static bool first_time_entered = true;
    if (first_time_entered) {
        ofs_sand_prod_file.open(sand_prod_filename, std::ofstream::out | std::ofstream::trunc);
        first_time_entered = false;
    }
 
    ModelPart::ElementsContainerType& pElements = dem_model_part.GetCommunicator().LocalMesh().Elements();
    double current_total_mass_in_grams = 0.0;
 
    for (unsigned int k = 0; k < pElements.size(); k++) {
 
        ModelPart::ElementsContainerType::iterator it = pElements.ptr_begin() + k;
        Element* raw_p_element = &(*it);
        SphericParticle* p_sphere = dynamic_cast<SphericParticle*>(raw_p_element);
        if (p_sphere->Is(ISOLATED)) continue;
        const double particle_density = p_sphere->GetDensity();
        const double particle_volume = p_sphere->CalculateVolume();
        current_total_mass_in_grams += particle_volume * particle_density * 1.0e3;
    }
    static const double initial_total_mass_in_grams = current_total_mass_in_grams;
    const double cumulative_sand_mass_in_grams = initial_total_mass_in_grams - current_total_mass_in_grams;
 
    //ModelPart::ConditionsContainerType::iterator condition_begin = outer_walls_model_part.ConditionsBegin();
    //const double face_pressure_in_psi = condition_begin->GetValue(POSITIVE_FACE_PRESSURE) * 0.000145;
    ProcessInfo& r_process_info = dem_model_part.GetProcessInfo();
    const double Pascals_to_psi_factor = 0.000145;
    const double face_pressure_in_psi = fabs(r_process_info[TARGET_STRESS_Z]) * Pascals_to_psi_factor;
 
    static std::ofstream sand_prod_file("sand_production_graph.txt", std::ios_base::out | std::ios_base::app);
    sand_prod_file << time << " " << face_pressure_in_psi << " " << cumulative_sand_mass_in_grams << '\n';
    sand_prod_file.flush();
 
    return cumulative_sand_mass_in_grams;
}
 
void MarkBrokenSpheres(ModelPart& dem_model_part) {
 
    ModelPart::ElementsContainerType& pElements = dem_model_part.GetCommunicator().LocalMesh().Elements();
 
    for (unsigned int k = 0; k < pElements.size(); k++) {
 
        ModelPart::ElementsContainerType::iterator it = pElements.ptr_begin() + k;
        Element* raw_p_element = &(*it);
        SphericContinuumParticle* p_sphere = dynamic_cast<SphericContinuumParticle*>(raw_p_element);
        if (p_sphere->Is(ISOLATED)) continue;
        bool go_to_next_particle = false;
 
        for (unsigned int i = 0; i < p_sphere->mContinuumInitialNeighborsSize; i++) {
 
            if (!p_sphere->mIniNeighbourFailureId[i]) {
                go_to_next_particle = true;
                break;
            }
        }
        if (go_to_next_particle) continue;
        else p_sphere->Set(ISOLATED, true);
    }
}
 
double ComputeSandProductionWithDepthFirstSearchNonRecursiveImplementation(ModelPart& dem_model_part, ModelPart& outer_walls_model_part, const double time) {
 
    const std::string sand_prod_filename = "sand_production_graph_with_chunks_non_recursive.txt";
    static std::ofstream ofs_sand_prod_file;
    const std::string granulometry_distr_filename = "granulometry_distribution.txt";
    static std::ofstream ofs_granulometry_distr_file;
    static bool first_time_entered = true;
    if (first_time_entered) {
        ofs_sand_prod_file.open(sand_prod_filename, std::ofstream::out | std::ofstream::trunc);
        ofs_granulometry_distr_file.open(granulometry_distr_filename, std::ofstream::out | std::ofstream::trunc);
        first_time_entered = false;
    }
 
    ModelPart::ElementsContainerType& pElements = dem_model_part.GetCommunicator().LocalMesh().Elements();
 
    std::vector<double> chunks_masses;
 
    for (unsigned int k = 0; k < pElements.size(); k++) {
        ModelPart::ElementsContainerType::iterator it = pElements.ptr_begin() + k;
        it->Set(VISITED, false);
    }
 
    std::vector<SphericContinuumParticle*> stack_of_particles_to_check;
 
    for (unsigned int k = 0; k < pElements.size(); k++) {
        ModelPart::ElementsContainerType::iterator it = pElements.ptr_begin() + k;
        Element* raw_p_element = &(*it);
        SphericContinuumParticle* p_sphere = dynamic_cast<SphericContinuumParticle*>(raw_p_element);
        double this_chunk_mass = 0.0;
        stack_of_particles_to_check.push_back(p_sphere);
        while (stack_of_particles_to_check.size()) {
            SphericContinuumParticle* current_particle = stack_of_particles_to_check.back();
            stack_of_particles_to_check.pop_back();
            if (current_particle->Is(VISITED)) continue;
            const double particle_density = current_particle->GetDensity();
            const double particle_volume = current_particle->CalculateVolume();
            this_chunk_mass += particle_volume * particle_density * 1.0e3;
 
            current_particle->Set(VISITED, true);
 
            for (size_t i = 0; i < current_particle->mContinuumInitialNeighborsSize; i++) {
                SphericParticle* p_neighbour_sphere = current_particle->mNeighbourElements[i];
                if (p_neighbour_sphere == NULL) continue;
 
                if (p_neighbour_sphere->Is(VISITED)) continue; //not necessary, but saves increasing and decreasing stack_of_particles_to_check's size
                if (current_particle->mIniNeighbourFailureId[i]) continue;
 
                auto existing_element_it = dem_model_part.GetMesh(0).Elements().find(p_neighbour_sphere->Id());
                if (existing_element_it == dem_model_part.GetMesh(0).ElementsEnd()) continue;
 
                SphericContinuumParticle* p_neigh_cont_sphere = dynamic_cast<SphericContinuumParticle*>(p_neighbour_sphere);
                stack_of_particles_to_check.push_back(p_neigh_cont_sphere);
            }
        }
        if (this_chunk_mass) chunks_masses.push_back(this_chunk_mass);
    }
 
    const double max_mass_of_a_single_chunck = *std::max_element(chunks_masses.begin(), chunks_masses.end());
    const double current_total_mass_in_grams = max_mass_of_a_single_chunck;
    static const double initial_total_mass_in_grams = current_total_mass_in_grams;
    const double cumulative_sand_mass_in_grams = initial_total_mass_in_grams - current_total_mass_in_grams;
 
    ProcessInfo& r_process_info = dem_model_part.GetProcessInfo();
    const double Pascals_to_psi_factor = 0.000145;
    const double face_pressure_in_psi = fabs(r_process_info[TARGET_STRESS_Z]) * Pascals_to_psi_factor;
 
    ofs_sand_prod_file << time << " " << face_pressure_in_psi << " " << cumulative_sand_mass_in_grams << '\n';
    ofs_sand_prod_file.flush();
 
    unsigned int number_of_time_steps_between_granulometry_prints = 1e9; //TODO: For the time being, we are not interested in printing this
    static unsigned int printing_counter = 0;
    if (printing_counter == number_of_time_steps_between_granulometry_prints) {
        ofs_granulometry_distr_file << time;
        for (unsigned int k = 0; k < chunks_masses.size(); k++) {
            ofs_granulometry_distr_file << " " << chunks_masses[k];
        }
        ofs_granulometry_distr_file << '\n';
        printing_counter = 0;
    }
    printing_counter++;
    ofs_granulometry_distr_file.flush();
 
    return cumulative_sand_mass_in_grams;
}
 
void ComputeSandProductionWithDepthFirstSearch(ModelPart& dem_model_part, ModelPart& outer_walls_model_part, const double time) {
 
    const std::string filename = "sand_production_graph_with_chunks.txt";
    std::ifstream ifile(filename.c_str());
    static bool first_time_entered = true;
    if ((bool) ifile && first_time_entered) {
        std::remove(filename.c_str());
        first_time_entered = false;
    }
 
    ModelPart::ElementsContainerType& pElements = dem_model_part.GetCommunicator().LocalMesh().Elements();
 
    std::vector<double> chunks_masses;
 
    for (unsigned int k = 0; k < pElements.size(); k++) {
        ModelPart::ElementsContainerType::iterator it = pElements.ptr_begin() + k;
        it->Set(VISITED, false);
    }
 
    for (unsigned int k = 0; k < pElements.size(); k++) {
        ModelPart::ElementsContainerType::iterator it = pElements.ptr_begin() + k;
        Element* raw_p_element = &(*it);
        SphericContinuumParticle* p_sphere = dynamic_cast<SphericContinuumParticle*>(raw_p_element);
        double this_chunk_mass = 0.0;
        if( it->IsNot(VISITED) ) {
            DepthFirstSearchVisit(p_sphere, this_chunk_mass);
            chunks_masses.push_back(this_chunk_mass);
        }
    }
 
    const double max_mass_of_a_single_chunck = *std::max_element(chunks_masses.begin(), chunks_masses.end());
    const double current_total_mass_in_grams = max_mass_of_a_single_chunck;
    static const double initial_total_mass_in_grams = current_total_mass_in_grams;
    const double cumulative_sand_mass_in_grams = initial_total_mass_in_grams - current_total_mass_in_grams;
 
    ModelPart::ConditionsContainerType::iterator condition_begin = outer_walls_model_part.ConditionsBegin();
    const double Pascals_to_psi_factor = 0.000145;
    const double face_pressure_in_psi = condition_begin->GetValue(POSITIVE_FACE_PRESSURE) * Pascals_to_psi_factor;
 
    static std::ofstream sand_prod_file(filename, std::ios_base::out | std::ios_base::app);
    sand_prod_file << time << " " << face_pressure_in_psi << " " << cumulative_sand_mass_in_grams << '\n';
    sand_prod_file.flush();
}
 
void DepthFirstSearchVisit(SphericContinuumParticle* p_sphere, double& this_chunk_mass) {
    p_sphere->Set(VISITED, true);
    const double particle_radius = p_sphere->GetRadius();
    const double particle_density = p_sphere->GetDensity();
    this_chunk_mass += (4.0/3.0) * Globals::Pi * particle_density * particle_radius * particle_radius * particle_radius * 1000.0;
    for (size_t i=0; i<p_sphere->mContinuumInitialNeighborsSize; i++) {
        SphericParticle* p_neighbour_sphere = p_sphere->mNeighbourElements[i];
        if (p_neighbour_sphere==NULL) continue;
        if (p_sphere->mIniNeighbourFailureId[i]) continue;
        if (p_neighbour_sphere->IsNot(VISITED)) {
            SphericContinuumParticle* p_neigh_cont_sphere = dynamic_cast<SphericContinuumParticle*>(p_neighbour_sphere);
            DepthFirstSearchVisit(p_neigh_cont_sphere, this_chunk_mass);
        }
    }
}
 
void ComputeTriaxialSandProduction(ModelPart& dem_model_part, ModelPart& outer_walls_model_part_1, ModelPart& outer_walls_model_part_2, const double time) {
 
    const std::string filename = "sand_production_graph.txt";
    std::ifstream ifile(filename.c_str());
    static bool first_time_entered = true;
    if ((bool) ifile && first_time_entered) {
        std::remove(filename.c_str());
        first_time_entered = false;
    }
 
    ModelPart::ElementsContainerType& pElements = dem_model_part.GetCommunicator().LocalMesh().Elements();
    double current_total_mass_in_grams = 0.0;
 
    for (unsigned int k = 0; k < pElements.size(); k++) {
 
        ModelPart::ElementsContainerType::iterator it = pElements.ptr_begin() + k;
        Element* raw_p_element = &(*it);
        SphericParticle* p_sphere = dynamic_cast<SphericParticle*>(raw_p_element);
        if (p_sphere->Is(ISOLATED)) continue;
        const double particle_radius = p_sphere->GetRadius();
        const double particle_density = p_sphere->GetDensity();
        current_total_mass_in_grams += (4.0/3.0) * Globals::Pi * particle_density * particle_radius * particle_radius * particle_radius * 1000.0;
    }
    static const double initial_total_mass_in_grams = current_total_mass_in_grams;
    const double cumulative_sand_mass_in_grams = initial_total_mass_in_grams - current_total_mass_in_grams;
 
    ModelPart::ConditionsContainerType::iterator condition_begin_1 = outer_walls_model_part_1.ConditionsBegin();
    ModelPart::ConditionsContainerType::iterator condition_begin_2 = outer_walls_model_part_2.ConditionsBegin();
 
    const double Pascals_to_psi_factor = 0.000145;
    const double face_pressure_in_psi = (condition_begin_1->GetValue(POSITIVE_FACE_PRESSURE) +
                                         condition_begin_2->GetValue(POSITIVE_FACE_PRESSURE) +
                                         3.45e6) * Pascals_to_psi_factor * 0.33333333333333; // 3.45e6 is the sigma_z constant pressure
 
    static std::ofstream sand_prod_file(filename, std::ios_base::out | std::ios_base::app);
    sand_prod_file << time << " " << face_pressure_in_psi << " " << cumulative_sand_mass_in_grams << '\n';
    sand_prod_file.flush();
}
 
//***************************************************************************************************************
//***************************************************************************************************************
 
 
virtual std::string Info() const
{
    return "";
}
 
 
virtual void PrintInfo(std::ostream& rOStream) const
{
}
 
 
virtual void PrintData(std::ostream& rOStream) const
{
}
 
protected:
 
private:
 
DemStructuresCouplingUtilities & operator=(DemStructuresCouplingUtilities const& rOther);
 
 
 
}; // Class DemStructuresCouplingUtilities
 
}  // namespace Python.
 
#endif // KRATOS_STRUCTURES_DEM_COUPLING_UTILITIES_H