diff --git a/include/cell.h b/include/cell.h index 0d1f1a39f..ae52573fa 100644 --- a/include/cell.h +++ b/include/cell.h @@ -83,12 +83,6 @@ class Cell { //! Return the status of a cell: active (if a particle is present) bool status() const { return particles_.size(); } - //! Assign solving status - void assign_solving_status(bool status) { solving_status_ = status; } - - //! Return solving status - bool solving_status() const { return solving_status_; } - //! Return particles_ std::vector particles() const { return particles_; } @@ -223,28 +217,50 @@ class Cell { //! Return previous mpi rank unsigned previous_mpirank() const; + /** + * \defgroup MultiPhase Functions dealing with multi-phase MPM + */ + /**@{*/ + + //! Assign solving status + //! \ingroup MultiPhase + //! \param[in] status Cell solving status for parallel free-surface detection + void assign_solving_status(bool status) { solving_status_ = status; } + + //! Return solving status of a cell: active (if a particle is present in all + //! MPI rank) + //! \ingroup MultiPhase + bool solving_status() const { return solving_status_; } + //! Assign free surface + //! \ingroup MultiPhase //! \param[in] free_surface boolean indicating free surface cell void assign_free_surface(bool free_surface) { free_surface_ = free_surface; }; //! Return free surface bool + //! \ingroup MultiPhase //! \retval free_surface_ indicating free surface cell bool free_surface() const { return free_surface_; }; //! Assign volume traction to node + //! \ingroup MultiPhase //! \param[in] volume_fraction cell volume fraction void assign_volume_fraction(double volume_fraction) { volume_fraction_ = volume_fraction; }; //! Return cell volume fraction + //! \ingroup MultiPhase //! \retval volume_fraction_ cell volume fraction double volume_fraction() const { return volume_fraction_; }; //! Map cell volume to the nodes + //! \ingroup MultiPhase //! \param[in] phase to map volume void map_cell_volume_to_nodes(unsigned phase); + /**@}*/ + private: //! Approximately check if a point is in a cell //! \param[in] point Coordinates of point diff --git a/include/cell.tcc b/include/cell.tcc index 43b9da653..2bbbc5731 100644 --- a/include/cell.tcc +++ b/include/cell.tcc @@ -850,7 +850,7 @@ template void mpm::Cell::map_cell_volume_to_nodes(unsigned phase) { if (this->status()) { // Check if cell volume is set - if (volume_ == std::numeric_limits::lowest()) + if (volume_ <= std::numeric_limits::lowest()) this->compute_volume(); for (unsigned i = 0; i < nodes_.size(); ++i) { diff --git a/include/mesh.h b/include/mesh.h index 00869f585..3a31fb75b 100644 --- a/include/mesh.h +++ b/include/mesh.h @@ -481,7 +481,19 @@ class Mesh { //! Inject particles void inject_particles(double current_time); + // Create the nodal properties' map + void create_nodal_properties(); + + // Initialise the nodal properties' map + void initialise_nodal_properties(); + + /** + * \defgroup MultiPhase Functions dealing with multi-phase MPM + */ + /**@{*/ + //! Compute free surface + //! \ingroup MultiPhase //! \param[in] method Type of method to use //! \param[in] volume_tolerance for volume_fraction approach //! \retval status Status of compute_free_surface @@ -490,6 +502,7 @@ class Mesh { double volume_tolerance = std::numeric_limits::epsilon()); //! Compute free surface by density method + //! \ingroup MultiPhase //! \details Using simple approach of volume fraction approach as (Kularathna //! & Soga, 2017) and density ratio comparison (Hamad, 2015). This method is //! fast, but less accurate. @@ -499,6 +512,7 @@ class Mesh { double volume_tolerance = std::numeric_limits::epsilon()); //! Compute free surface by geometry method + //! \ingroup MultiPhase //! \details Using a more expensive approach using neighbouring particles and //! current geometry. This method combine multiple checks in order to simplify //! and fasten the process: (1) Volume fraction approach as (Kularathna & Soga @@ -510,30 +524,30 @@ class Mesh { double volume_tolerance = std::numeric_limits::epsilon()); //! Get free surface node set + //! \ingroup MultiPhase //! \retval id_set Set of free surface node ids std::set free_surface_nodes(); //! Get free surface cell set + //! \ingroup MultiPhase //! \retval id_set Set of free surface cell ids std::set free_surface_cells(); //! Get free surface particle set + //! \ingroup MultiPhase //! \retval status Status of compute_free_surface //! \retval id_set Set of free surface particle ids std::set free_surface_particles(); - // Create the nodal properties' map - void create_nodal_properties(); - - // Initialise the nodal properties' map - void initialise_nodal_properties(); - //! Assign particles pore pressures + //! \ingroup MultiPhase //! \param[in] particle_pore_pressure Initial pore pressure of particle bool assign_particles_pore_pressures( const std::vector>& particle_pore_pressures); + /**@}*/ + private: // Read particles from file //! \param[in] pset_id Set ID of the particles @@ -606,5 +620,6 @@ class Mesh { } // namespace mpm #include "mesh.tcc" +#include "mesh_multiphase.tcc" #endif // MPM_MESH_H_ diff --git a/include/mesh.tcc b/include/mesh.tcc index 0d024d7e7..7411235c5 100644 --- a/include/mesh.tcc +++ b/include/mesh.tcc @@ -2092,504 +2092,9 @@ void mpm::Mesh::create_nodal_properties() { } } -//! Compute free surface cells, nodes, and particles -template -bool mpm::Mesh::compute_free_surface(const std::string& method, - double volume_tolerance) { - if (method == "density") { - return this->compute_free_surface_by_density(volume_tolerance); - } else if (method == "geometry") { - return this->compute_free_surface_by_geometry(volume_tolerance); - } else { - console_->info( - "The selected free-surface computation method: {}\n is not " - "available. " - "Using density approach as default method.", - method); - return this->compute_free_surface_by_density(volume_tolerance); - } -} - -//! Compute free surface cells, nodes, and particles by density and geometry -template -bool mpm::Mesh::compute_free_surface_by_geometry( - double volume_tolerance) { - bool status = true; - try { - // Reset free surface cell and particles - this->iterate_over_cells(std::bind(&mpm::Cell::assign_free_surface, - std::placeholders::_1, false)); - - VectorDim temp_normal; - temp_normal.setZero(); - this->iterate_over_particles_predicate( - std::bind(&mpm::ParticleBase::assign_normal, - std::placeholders::_1, temp_normal), - std::bind(&mpm::ParticleBase::free_surface, - std::placeholders::_1)); - - this->iterate_over_particles( - std::bind(&mpm::ParticleBase::assign_free_surface, - std::placeholders::_1, false)); - - // Reset volume fraction - this->iterate_over_cells(std::bind(&mpm::Cell::assign_volume_fraction, - std::placeholders::_1, 0.0)); - - // Compute and assign volume fraction to each cell - for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); - ++citr) { - if ((*citr)->status()) { - // Compute volume fraction - double cell_volume_fraction = 0.0; - for (const auto p_id : (*citr)->particles()) - cell_volume_fraction += map_particles_[p_id]->volume(); - - cell_volume_fraction = cell_volume_fraction / (*citr)->volume(); - (*citr)->assign_volume_fraction(cell_volume_fraction); - } - } - - // First, we detect the cell with possible free surfaces - // Compute boundary cells and nodes based on geometry - std::set free_surface_candidate_cells; - for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); - ++citr) { - // Cell contains particles - if ((*citr)->status()) { - bool candidate_cell = false; - const auto& node_id = (*citr)->nodes_id(); - if ((*citr)->volume_fraction() < volume_tolerance) { - candidate_cell = true; - for (const auto id : node_id) { - map_nodes_[id]->assign_free_surface(true); - } - } else { - // Loop over neighbouring cells - for (const auto neighbour_cell_id : (*citr)->neighbours()) { - if (!map_cells_[neighbour_cell_id]->status()) { - candidate_cell = true; - const auto& n_node_id = map_cells_[neighbour_cell_id]->nodes_id(); - - // Detect common node id - std::set common_node_id; - std::set_intersection( - node_id.begin(), node_id.end(), n_node_id.begin(), - n_node_id.end(), - std::inserter(common_node_id, common_node_id.begin())); - - // Assign free surface nodes - if (!common_node_id.empty()) { - for (const auto common_id : common_node_id) { - map_nodes_[common_id]->assign_free_surface(true); - } - } - } - } - } - - // Assign free surface cell - if (candidate_cell) { - (*citr)->assign_free_surface(true); - free_surface_candidate_cells.insert((*citr)->id()); - } - } - } - - // Compute particle neighbours for particles at candidate cells - std::vector free_surface_candidate_particles_first; - for (const auto cell_id : free_surface_candidate_cells) { - this->find_particle_neighbours(map_cells_[cell_id]); - const auto& particle_ids = map_cells_[cell_id]->particles(); - free_surface_candidate_particles_first.insert( - free_surface_candidate_particles_first.end(), particle_ids.begin(), - particle_ids.end()); - } - - // Compute boundary particles based on density function - // Lump cell volume to nodes - this->iterate_over_cells(std::bind( - &mpm::Cell::map_cell_volume_to_nodes, std::placeholders::_1, 0)); - - // Compute nodal value of mass density - this->iterate_over_nodes_predicate( - std::bind(&mpm::NodeBase::compute_density, std::placeholders::_1), - std::bind(&mpm::NodeBase::status, std::placeholders::_1)); - - std::set free_surface_candidate_particles_second; - for (const auto p_id : free_surface_candidate_particles_first) { - const auto& particle = map_particles_[p_id]; - status = particle->compute_free_surface_by_density(); - if (status) free_surface_candidate_particles_second.insert(p_id); - } - - // Find free surface particles through geometry - for (const auto p_id : free_surface_candidate_particles_second) { - // Initialize renormalization matrix - Eigen::Matrix renormalization_matrix_inv; - renormalization_matrix_inv.setZero(); - - // Loop over neighbours - const auto& particle = map_particles_[p_id]; - const auto& p_coord = particle->coordinates(); - const auto& neighbour_particles = particle->neighbours(); - const double smoothing_length = 1.33 * particle->diameter(); - for (const auto neighbour_particle_id : neighbour_particles) { - const auto& n_coord = - map_particles_[neighbour_particle_id]->coordinates(); - const VectorDim rel_coord = n_coord - p_coord; - - // Compute kernel gradient - const VectorDim kernel_gradient = - mpm::RadialBasisFunction::gradient(smoothing_length, - -rel_coord, "gaussian"); - - // Inverse of renormalization matrix B - renormalization_matrix_inv += - (particle->mass() / particle->mass_density()) * kernel_gradient * - rel_coord.transpose(); - } - - // Compute lambda: minimum eigenvalue of B_inverse - Eigen::SelfAdjointEigenSolver es( - renormalization_matrix_inv); - double lambda = es.eigenvalues().minCoeff(); - - // Categorize particle based on lambda - bool free_surface = false; - bool secondary_check = false; - bool interior = false; - if (lambda <= 0.2) - free_surface = true; - else if (lambda > 0.2 && lambda <= 0.75) - secondary_check = true; - else - interior = true; - - // Compute numerical normal vector - VectorDim normal; - normal.setZero(); - if (!interior) { - VectorDim temporary_vec; - temporary_vec.setZero(); - for (const auto neighbour_particle_id : neighbour_particles) { - const auto& n_coord = - map_particles_[neighbour_particle_id]->coordinates(); - const VectorDim rel_coord = n_coord - p_coord; - - // Compute kernel gradient - const VectorDim kernel_gradient = - mpm::RadialBasisFunction::gradient(smoothing_length, - -rel_coord, "gaussian"); - - // Sum of kernel by volume - temporary_vec += - (particle->mass() / particle->mass_density()) * kernel_gradient; - } - normal = -renormalization_matrix_inv.inverse() * temporary_vec; - if (normal.norm() > std::numeric_limits::epsilon()) - normal.normalize(); - else - normal.setZero(); - } - - // If secondary check is needed - if (secondary_check) { - // Construct scanning region - // TODO: spacing distance should be a function of porosity - const double spacing_distance = smoothing_length; - VectorDim t_coord = p_coord + spacing_distance * normal; - - // Check all neighbours - for (const auto neighbour_particle_id : neighbour_particles) { - const auto& n_coord = - map_particles_[neighbour_particle_id]->coordinates(); - const VectorDim rel_coord_np = n_coord - p_coord; - const double distance_np = rel_coord_np.norm(); - const VectorDim rel_coord_nt = n_coord - t_coord; - const double distance_nt = rel_coord_nt.norm(); - - free_surface = true; - if (distance_np < std::sqrt(2) * spacing_distance) { - if (std::acos(normal.dot(rel_coord_np) / distance_np) < M_PI / 4) { - free_surface = false; - break; - } - } else { - if (distance_nt < spacing_distance) { - free_surface = false; - break; - } - } - } - } - - // Assign normal only to validated free surface - if (free_surface) { - particle->assign_free_surface(true); - particle->assign_normal(normal); - } - } - } catch (std::exception& exception) { - console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); - } - return status; -} - -//! Compute free surface cells, nodes, and particles by density -template -bool mpm::Mesh::compute_free_surface_by_density(double volume_tolerance) { - bool status = true; - try { - - // Get number of MPI ranks - int mpi_rank = 0; - int mpi_size = 1; - -#ifdef USE_MPI - MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); - MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); - - // Reset solving status - this->iterate_over_cells(std::bind(&mpm::Cell::assign_solving_status, - std::placeholders::_1, false)); - - // Initialize pointer of booleans for send and receive - bool* send_cell_solving_status = new bool[ncells()]; - memset(send_cell_solving_status, 0, ncells() * sizeof(bool)); - bool* receive_cell_solving_status = new bool[ncells()]; - - for (auto citr = cells_.cbegin(); citr != cells_.cend(); ++citr) - if ((*citr)->status()) - // Assign solving status for MPI solver - send_cell_solving_status[(*citr)->id()] = true; - - MPI_Allreduce(send_cell_solving_status, receive_cell_solving_status, - ncells(), MPI_CXX_BOOL, MPI_LOR, MPI_COMM_WORLD); - - for (auto citr = cells_.cbegin(); citr != cells_.cend(); ++citr) { - if (receive_cell_solving_status[(*citr)->id()]) { - // Assign solving status for MPI solver - (*citr)->assign_solving_status(true); - } - } - - delete[] send_cell_solving_status; - delete[] receive_cell_solving_status; -#endif - - // Reset free surface cell - this->iterate_over_cells(std::bind(&mpm::Cell::assign_free_surface, - std::placeholders::_1, false)); - - // Reset volume fraction - this->iterate_over_cells(std::bind(&mpm::Cell::assign_volume_fraction, - std::placeholders::_1, 0.0)); - - // Reset free surface particle - this->iterate_over_particles( - std::bind(&mpm::ParticleBase::assign_free_surface, - std::placeholders::_1, false)); - - // Compute and assign volume fraction to each cell - for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); - ++citr) { - if ((*citr)->status()) { - // Compute volume fraction - double cell_volume_fraction = 0.0; - for (const auto p_id : (*citr)->particles()) - cell_volume_fraction += map_particles_[p_id]->volume(); - - cell_volume_fraction = cell_volume_fraction / (*citr)->volume(); - (*citr)->assign_volume_fraction(cell_volume_fraction); - } - } - -#ifdef USE_MPI - // Initialize vector of double for send and receive - std::vector send_cell_vol_fraction; - send_cell_vol_fraction.resize(ncells()); - std::vector receive_cell_vol_fraction; - receive_cell_vol_fraction.resize(ncells()); - - for (auto citr = cells_.cbegin(); citr != cells_.cend(); ++citr) - if ((*citr)->status()) - // Assign volume_fraction for MPI solver - send_cell_vol_fraction[(*citr)->id()] = (*citr)->volume_fraction(); - - MPI_Allreduce(send_cell_vol_fraction.data(), - receive_cell_vol_fraction.data(), ncells(), MPI_DOUBLE, - MPI_SUM, MPI_COMM_WORLD); - - for (auto citr = cells_.cbegin(); citr != cells_.cend(); ++citr) { - // Assign volume_fraction for MPI solver - (*citr)->assign_volume_fraction(receive_cell_vol_fraction[(*citr)->id()]); - } -#endif - - // Compute boundary cells and nodes based on geometry - for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); - ++citr) { - - if ((*citr)->status()) { - bool cell_at_interface = false; - const auto& node_id = (*citr)->nodes_id(); - bool internal = true; - - //! Check internal cell - for (const auto neighbour_cell_id : (*citr)->neighbours()) { -#if USE_MPI - if (!map_cells_[neighbour_cell_id]->solving_status()) { - internal = false; - break; - } -#else - if (!map_cells_[neighbour_cell_id]->status()) { - internal = false; - break; - } -#endif - } - - //! Check volume fraction only for boundary cell - if (!internal) { - if ((*citr)->volume_fraction() < volume_tolerance) { - cell_at_interface = true; - for (const auto id : node_id) { - map_nodes_[id]->assign_free_surface(cell_at_interface); - } - } else { - for (const auto neighbour_cell_id : (*citr)->neighbours()) { - if (map_cells_[neighbour_cell_id]->volume_fraction() < - volume_tolerance) { - cell_at_interface = true; - const auto& n_node_id = - map_cells_[neighbour_cell_id]->nodes_id(); - - // Detect common node id - std::set common_node_id; - std::set_intersection( - node_id.begin(), node_id.end(), n_node_id.begin(), - n_node_id.end(), - std::inserter(common_node_id, common_node_id.begin())); - - // Assign free surface nodes - if (!common_node_id.empty()) { - for (const auto common_id : common_node_id) { - map_nodes_[common_id]->assign_free_surface(true); - } - } - } - } - } - - // Assign free surface cell - if (cell_at_interface) { - (*citr)->assign_free_surface(cell_at_interface); - } - } - } - } - - // Compute boundary particles based on density function - // Lump cell volume to nodes - this->iterate_over_cells( - std::bind(&mpm::Cell::map_cell_volume_to_nodes, - std::placeholders::_1, mpm::ParticlePhase::SinglePhase)); - -#ifdef USE_MPI - // Run if there is more than a single MPI task - if (mpi_size > 1) { - // MPI all reduce nodal volume - this->template nodal_halo_exchange( - std::bind(&mpm::NodeBase::volume, std::placeholders::_1, - mpm::ParticlePhase::SinglePhase), - std::bind(&mpm::NodeBase::update_volume, std::placeholders::_1, - false, mpm::ParticlePhase::SinglePhase, - std::placeholders::_2)); - } -#endif - - // Compute nodal value of mass density - this->iterate_over_nodes_predicate( - std::bind(&mpm::NodeBase::compute_density, std::placeholders::_1), - std::bind(&mpm::NodeBase::status, std::placeholders::_1)); - - // Evaluate free surface particles - for (auto pitr = this->particles_.cbegin(); pitr != this->particles_.cend(); - ++pitr) { - bool status = (*pitr)->compute_free_surface_by_density(); - if (status) { - (*pitr)->assign_free_surface(status); - } - } - } catch (std::exception& exception) { - console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); - } - return status; -} - -//! Get free surface node set -template -std::set mpm::Mesh::free_surface_nodes() { - std::set id_set; - for (auto nitr = this->nodes_.cbegin(); nitr != this->nodes_.cend(); ++nitr) - if ((*nitr)->free_surface()) id_set.insert((*nitr)->id()); - return id_set; -} - -//! Get free surface cell set -template -std::set mpm::Mesh::free_surface_cells() { - std::set id_set; - for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); ++citr) - if ((*citr)->free_surface()) id_set.insert((*citr)->id()); - return id_set; -} - -//! Get free surface particle set -template -std::set mpm::Mesh::free_surface_particles() { - std::set id_set; - for (auto pitr = this->particles_.cbegin(); pitr != this->particles_.cend(); - ++pitr) - if ((*pitr)->free_surface()) id_set.insert((*pitr)->id()); - return id_set; -} - // Initialise the nodal properties' map template void mpm::Mesh::initialise_nodal_properties() { // Call initialise_properties function from the nodal properties nodal_properties_->initialise_nodal_properties(); -} - -//! Assign particle pore pressures -template -bool mpm::Mesh::assign_particles_pore_pressures( - const std::vector>& - particle_pore_pressures) { - bool status = true; - - try { - if (!particles_.size()) - throw std::runtime_error( - "No particles have been assigned in mesh, cannot assign pore " - "pressures"); - - for (const auto& particle_pore_pressure : particle_pore_pressures) { - // Particle id - mpm::Index pid = std::get<0>(particle_pore_pressure); - // Pore pressure - double pore_pressure = std::get<1>(particle_pore_pressure); - - if (map_particles_.find(pid) != map_particles_.end()) - map_particles_[pid]->assign_pressure(pore_pressure, - mpm::ParticlePhase::Liquid); - } - } catch (std::exception& exception) { - console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); - status = false; - } - return status; } \ No newline at end of file diff --git a/include/mesh_multiphase.tcc b/include/mesh_multiphase.tcc new file mode 100644 index 000000000..4a5b0a830 --- /dev/null +++ b/include/mesh_multiphase.tcc @@ -0,0 +1,494 @@ +//! Compute free surface cells, nodes, and particles +template +bool mpm::Mesh::compute_free_surface(const std::string& method, + double volume_tolerance) { + if (method == "density") { + return this->compute_free_surface_by_density(volume_tolerance); + } else if (method == "geometry") { + return this->compute_free_surface_by_geometry(volume_tolerance); + } else { + console_->info( + "The selected free-surface computation method: {}\n is not " + "available. " + "Using density approach as default method.", + method); + return this->compute_free_surface_by_density(volume_tolerance); + } +} + +//! Compute free surface cells, nodes, and particles by density and geometry +template +bool mpm::Mesh::compute_free_surface_by_geometry( + double volume_tolerance) { + bool status = true; + try { + // Reset free surface cell and particles + this->iterate_over_cells(std::bind(&mpm::Cell::assign_free_surface, + std::placeholders::_1, false)); + + VectorDim temp_normal; + temp_normal.setZero(); + this->iterate_over_particles_predicate( + std::bind(&mpm::ParticleBase::assign_normal, + std::placeholders::_1, temp_normal), + std::bind(&mpm::ParticleBase::free_surface, + std::placeholders::_1)); + + this->iterate_over_particles( + std::bind(&mpm::ParticleBase::assign_free_surface, + std::placeholders::_1, false)); + + // Reset volume fraction + this->iterate_over_cells(std::bind(&mpm::Cell::assign_volume_fraction, + std::placeholders::_1, 0.0)); + + // Compute and assign volume fraction to each cell + for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); + ++citr) { + if ((*citr)->status()) { + // Compute volume fraction + double cell_volume_fraction = 0.0; + for (const auto p_id : (*citr)->particles()) + cell_volume_fraction += map_particles_[p_id]->volume(); + + cell_volume_fraction = cell_volume_fraction / (*citr)->volume(); + (*citr)->assign_volume_fraction(cell_volume_fraction); + } + } + + // First, we detect the cell with possible free surfaces + // Compute boundary cells and nodes based on geometry + std::set free_surface_candidate_cells; + for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); + ++citr) { + // Cell contains particles + if ((*citr)->status()) { + bool candidate_cell = false; + const auto& node_id = (*citr)->nodes_id(); + if ((*citr)->volume_fraction() < volume_tolerance) { + candidate_cell = true; + for (const auto id : node_id) { + map_nodes_[id]->assign_free_surface(true); + } + } else { + // Loop over neighbouring cells + for (const auto neighbour_cell_id : (*citr)->neighbours()) { + if (!map_cells_[neighbour_cell_id]->status()) { + candidate_cell = true; + const auto& n_node_id = map_cells_[neighbour_cell_id]->nodes_id(); + + // Detect common node id + std::set common_node_id; + std::set_intersection( + node_id.begin(), node_id.end(), n_node_id.begin(), + n_node_id.end(), + std::inserter(common_node_id, common_node_id.begin())); + + // Assign free surface nodes + if (!common_node_id.empty()) { + for (const auto common_id : common_node_id) { + map_nodes_[common_id]->assign_free_surface(true); + } + } + } + } + } + + // Assign free surface cell + if (candidate_cell) { + (*citr)->assign_free_surface(true); + free_surface_candidate_cells.insert((*citr)->id()); + } + } + } + + // Compute particle neighbours for particles at candidate cells + std::vector free_surface_candidate_particles_first; + for (const auto cell_id : free_surface_candidate_cells) { + this->find_particle_neighbours(map_cells_[cell_id]); + const auto& particle_ids = map_cells_[cell_id]->particles(); + free_surface_candidate_particles_first.insert( + free_surface_candidate_particles_first.end(), particle_ids.begin(), + particle_ids.end()); + } + + // Compute boundary particles based on density function + // Lump cell volume to nodes + this->iterate_over_cells(std::bind( + &mpm::Cell::map_cell_volume_to_nodes, std::placeholders::_1, 0)); + + // Compute nodal value of mass density + this->iterate_over_nodes_predicate( + std::bind(&mpm::NodeBase::compute_density, std::placeholders::_1), + std::bind(&mpm::NodeBase::status, std::placeholders::_1)); + + std::set free_surface_candidate_particles_second; + for (const auto p_id : free_surface_candidate_particles_first) { + const auto& particle = map_particles_[p_id]; + status = particle->compute_free_surface_by_density(); + if (status) free_surface_candidate_particles_second.insert(p_id); + } + + // Find free surface particles through geometry + for (const auto p_id : free_surface_candidate_particles_second) { + // Initialize renormalization matrix + Eigen::Matrix renormalization_matrix_inv; + renormalization_matrix_inv.setZero(); + + // Loop over neighbours + const auto& particle = map_particles_[p_id]; + const auto& p_coord = particle->coordinates(); + const auto& neighbour_particles = particle->neighbours(); + const double smoothing_length = 1.33 * particle->diameter(); + for (const auto neighbour_particle_id : neighbour_particles) { + const auto& n_coord = + map_particles_[neighbour_particle_id]->coordinates(); + const VectorDim rel_coord = n_coord - p_coord; + + // Compute kernel gradient + const VectorDim kernel_gradient = + mpm::RadialBasisFunction::gradient(smoothing_length, + -rel_coord, "gaussian"); + + // Inverse of renormalization matrix B + renormalization_matrix_inv += + (particle->mass() / particle->mass_density()) * kernel_gradient * + rel_coord.transpose(); + } + + // Compute lambda: minimum eigenvalue of B_inverse + Eigen::SelfAdjointEigenSolver es( + renormalization_matrix_inv); + double lambda = es.eigenvalues().minCoeff(); + + // Categorize particle based on lambda + bool free_surface = false; + bool secondary_check = false; + bool interior = false; + if (lambda <= 0.2) + free_surface = true; + else if (lambda > 0.2 && lambda <= 0.75) + secondary_check = true; + else + interior = true; + + // Compute numerical normal vector + VectorDim normal; + normal.setZero(); + if (!interior) { + VectorDim temporary_vec; + temporary_vec.setZero(); + for (const auto neighbour_particle_id : neighbour_particles) { + const auto& n_coord = + map_particles_[neighbour_particle_id]->coordinates(); + const VectorDim rel_coord = n_coord - p_coord; + + // Compute kernel gradient + const VectorDim kernel_gradient = + mpm::RadialBasisFunction::gradient(smoothing_length, + -rel_coord, "gaussian"); + + // Sum of kernel by volume + temporary_vec += + (particle->mass() / particle->mass_density()) * kernel_gradient; + } + normal = -renormalization_matrix_inv.inverse() * temporary_vec; + if (normal.norm() > std::numeric_limits::epsilon()) + normal.normalize(); + else + normal.setZero(); + } + + // If secondary check is needed + if (secondary_check) { + // Construct scanning region + // TODO: spacing distance should be a function of porosity + const double spacing_distance = smoothing_length; + VectorDim t_coord = p_coord + spacing_distance * normal; + + // Check all neighbours + for (const auto neighbour_particle_id : neighbour_particles) { + const auto& n_coord = + map_particles_[neighbour_particle_id]->coordinates(); + const VectorDim rel_coord_np = n_coord - p_coord; + const double distance_np = rel_coord_np.norm(); + const VectorDim rel_coord_nt = n_coord - t_coord; + const double distance_nt = rel_coord_nt.norm(); + + free_surface = true; + if (distance_np < std::sqrt(2) * spacing_distance) { + if (std::acos(normal.dot(rel_coord_np) / distance_np) < M_PI / 4) { + free_surface = false; + break; + } + } else { + if (distance_nt < spacing_distance) { + free_surface = false; + break; + } + } + } + } + + // Assign normal only to validated free surface + if (free_surface) { + particle->assign_free_surface(true); + particle->assign_normal(normal); + } + } + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + } + return status; +} + +//! Compute free surface cells, nodes, and particles by density +template +bool mpm::Mesh::compute_free_surface_by_density(double volume_tolerance) { + bool status = true; + try { + + // Get number of MPI ranks + int mpi_rank = 0; + int mpi_size = 1; + +#ifdef USE_MPI + MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); + + // Reset solving status + this->iterate_over_cells(std::bind(&mpm::Cell::assign_solving_status, + std::placeholders::_1, false)); + + // Initialize pointer of booleans for send and receive + bool* send_cell_solving_status = new bool[ncells()]; + memset(send_cell_solving_status, 0, ncells() * sizeof(bool)); + bool* receive_cell_solving_status = new bool[ncells()]; + + for (auto citr = cells_.cbegin(); citr != cells_.cend(); ++citr) + if ((*citr)->status()) + // Assign solving status for MPI solver + send_cell_solving_status[(*citr)->id()] = true; + + MPI_Allreduce(send_cell_solving_status, receive_cell_solving_status, + ncells(), MPI_CXX_BOOL, MPI_LOR, MPI_COMM_WORLD); + + for (auto citr = cells_.cbegin(); citr != cells_.cend(); ++citr) { + if (receive_cell_solving_status[(*citr)->id()]) { + // Assign solving status for MPI solver + (*citr)->assign_solving_status(true); + } + } + + delete[] send_cell_solving_status; + delete[] receive_cell_solving_status; +#endif + + // Reset free surface cell + this->iterate_over_cells(std::bind(&mpm::Cell::assign_free_surface, + std::placeholders::_1, false)); + + // Reset volume fraction + this->iterate_over_cells(std::bind(&mpm::Cell::assign_volume_fraction, + std::placeholders::_1, 0.0)); + + // Reset free surface particle + this->iterate_over_particles( + std::bind(&mpm::ParticleBase::assign_free_surface, + std::placeholders::_1, false)); + + // Compute and assign volume fraction to each cell + for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); + ++citr) { + if ((*citr)->status()) { + // Compute volume fraction + double cell_volume_fraction = 0.0; + for (const auto p_id : (*citr)->particles()) + cell_volume_fraction += map_particles_[p_id]->volume(); + + cell_volume_fraction = cell_volume_fraction / (*citr)->volume(); + (*citr)->assign_volume_fraction(cell_volume_fraction); + } + } + +#ifdef USE_MPI + // Initialize vector of double for send and receive + std::vector send_cell_vol_fraction; + send_cell_vol_fraction.resize(ncells()); + std::vector receive_cell_vol_fraction; + receive_cell_vol_fraction.resize(ncells()); + + for (auto citr = cells_.cbegin(); citr != cells_.cend(); ++citr) + if ((*citr)->status()) + // Assign volume_fraction for MPI solver + send_cell_vol_fraction[(*citr)->id()] = (*citr)->volume_fraction(); + + MPI_Allreduce(send_cell_vol_fraction.data(), + receive_cell_vol_fraction.data(), ncells(), MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); + + for (auto citr = cells_.cbegin(); citr != cells_.cend(); ++citr) { + // Assign volume_fraction for MPI solver + (*citr)->assign_volume_fraction(receive_cell_vol_fraction[(*citr)->id()]); + } +#endif + + // Compute boundary cells and nodes based on geometry + for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); + ++citr) { + + if ((*citr)->status()) { + bool cell_at_interface = false; + const auto& node_id = (*citr)->nodes_id(); + bool internal = true; + + //! Check internal cell + for (const auto neighbour_cell_id : (*citr)->neighbours()) { +#if USE_MPI + if (!map_cells_[neighbour_cell_id]->solving_status()) { + internal = false; + break; + } +#else + if (!map_cells_[neighbour_cell_id]->status()) { + internal = false; + break; + } +#endif + } + + //! Check volume fraction only for boundary cell + if (!internal) { + if ((*citr)->volume_fraction() < volume_tolerance) { + cell_at_interface = true; + for (const auto id : node_id) { + map_nodes_[id]->assign_free_surface(cell_at_interface); + } + } else { + for (const auto neighbour_cell_id : (*citr)->neighbours()) { + if (map_cells_[neighbour_cell_id]->volume_fraction() < + volume_tolerance) { + cell_at_interface = true; + const auto& n_node_id = + map_cells_[neighbour_cell_id]->nodes_id(); + + // Detect common node id + std::set common_node_id; + std::set_intersection( + node_id.begin(), node_id.end(), n_node_id.begin(), + n_node_id.end(), + std::inserter(common_node_id, common_node_id.begin())); + + // Assign free surface nodes + if (!common_node_id.empty()) { + for (const auto common_id : common_node_id) { + map_nodes_[common_id]->assign_free_surface(true); + } + } + } + } + } + + // Assign free surface cell + if (cell_at_interface) { + (*citr)->assign_free_surface(cell_at_interface); + } + } + } + } + + // Compute boundary particles based on density function + // Lump cell volume to nodes + this->iterate_over_cells( + std::bind(&mpm::Cell::map_cell_volume_to_nodes, + std::placeholders::_1, mpm::ParticlePhase::SinglePhase)); + +#ifdef USE_MPI + // Run if there is more than a single MPI task + if (mpi_size > 1) { + // MPI all reduce nodal volume + this->template nodal_halo_exchange( + std::bind(&mpm::NodeBase::volume, std::placeholders::_1, + mpm::ParticlePhase::SinglePhase), + std::bind(&mpm::NodeBase::update_volume, std::placeholders::_1, + false, mpm::ParticlePhase::SinglePhase, + std::placeholders::_2)); + } +#endif + + // Compute nodal value of mass density + this->iterate_over_nodes_predicate( + std::bind(&mpm::NodeBase::compute_density, std::placeholders::_1), + std::bind(&mpm::NodeBase::status, std::placeholders::_1)); + + // Evaluate free surface particles + for (auto pitr = this->particles_.cbegin(); pitr != this->particles_.cend(); + ++pitr) { + bool status = (*pitr)->compute_free_surface_by_density(); + if (status) { + (*pitr)->assign_free_surface(status); + } + } + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + } + return status; +} + +//! Get free surface node set +template +std::set mpm::Mesh::free_surface_nodes() { + std::set id_set; + for (auto nitr = this->nodes_.cbegin(); nitr != this->nodes_.cend(); ++nitr) + if ((*nitr)->free_surface()) id_set.insert((*nitr)->id()); + return id_set; +} + +//! Get free surface cell set +template +std::set mpm::Mesh::free_surface_cells() { + std::set id_set; + for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); ++citr) + if ((*citr)->free_surface()) id_set.insert((*citr)->id()); + return id_set; +} + +//! Get free surface particle set +template +std::set mpm::Mesh::free_surface_particles() { + std::set id_set; + for (auto pitr = this->particles_.cbegin(); pitr != this->particles_.cend(); + ++pitr) + if ((*pitr)->free_surface()) id_set.insert((*pitr)->id()); + return id_set; +} + +//! Assign particle pore pressures +template +bool mpm::Mesh::assign_particles_pore_pressures( + const std::vector>& + particle_pore_pressures) { + bool status = true; + + try { + if (!particles_.size()) + throw std::runtime_error( + "No particles have been assigned in mesh, cannot assign pore " + "pressures"); + + for (const auto& particle_pore_pressure : particle_pore_pressures) { + // Particle id + mpm::Index pid = std::get<0>(particle_pore_pressure); + // Pore pressure + double pore_pressure = std::get<1>(particle_pore_pressure); + + if (map_particles_.find(pid) != map_particles_.end()) + map_particles_[pid]->assign_pressure(pore_pressure, + mpm::ParticlePhase::Liquid); + } + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} \ No newline at end of file diff --git a/include/node.h b/include/node.h index f3fee8da9..d5db61c62 100644 --- a/include/node.h +++ b/include/node.h @@ -294,8 +294,13 @@ class Node : public NodeBase { //! Compute multimaterial normal unit vector void compute_multimaterial_normal_unit_vector() override; - //! TwoPhase functions-------------------------------------------------------- + /** + * \defgroup TwoPhase Functions dealing with two-phase MPM + */ + /**@{*/ + //! Update internal force (body force / traction force) + //! \ingroup TwoPhase //! \param[in] update A boolean to update (true) or assign (false) //! \param[in] drag_force Drag force from the particles in a cell //! \retval status Update status @@ -303,21 +308,25 @@ class Node : public NodeBase { const VectorDim& drag_force) override; //! Compute acceleration and velocity for two phase + //! \ingroup TwoPhase //! \param[in] dt Timestep in analysis bool compute_acceleration_velocity_twophase_explicit( double dt) noexcept override; //! Compute acceleration and velocity for two phase with cundall damping + //! \ingroup TwoPhase //! \param[in] dt Timestep in analysis \param[in] damping_factor //! Damping factor bool compute_acceleration_velocity_twophase_explicit_cundall( double dt, double damping_factor) noexcept override; //! Return drag force at a given node + //! \ingroup TwoPhase VectorDim drag_force_coefficient() const override { return drag_force_coefficient_; } - //---------------------------------------------------------------------------- + + /**@}*/ private: //! Mutex diff --git a/include/node.tcc b/include/node.tcc index 5583ea056..2fd2bb162 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -183,17 +183,16 @@ bool mpm::Node::assign_pressure_constraint( bool status = true; try { // Constrain directions can take values between 0 and Tnphases - if (phase < Tnphases * 2) { - this->pressure_constraints_.insert(std::make_pair( - static_cast(phase), static_cast(pressure))); - // Assign pressure function - if (function != nullptr) - this->pressure_function_.insert( - std::make_pair>( - static_cast(phase), - static_cast>(function))); - } else - throw std::runtime_error("Pressure constraint phase is out of bounds"); + assert(phase < Tnphases * 2); + + this->pressure_constraints_.insert(std::make_pair( + static_cast(phase), static_cast(pressure))); + // Assign pressure function + if (function != nullptr) + this->pressure_function_.insert( + std::make_pair>( + static_cast(phase), + static_cast>(function))); } catch (std::exception& exception) { console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); diff --git a/include/node_base.h b/include/node_base.h index e5cf021c8..b8e206835 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -21,7 +21,7 @@ namespace mpm { enum NodePhase : unsigned int { NSolid = 0, NLiquid = 1, - nGas = 2, + NGas = 2, NMixture = 0 }; diff --git a/tests/mesh_test_2d.cc b/tests/mesh_test_2d.cc index c721b88d8..bdff4cad1 100644 --- a/tests/mesh_test_2d.cc +++ b/tests/mesh_test_2d.cc @@ -1236,13 +1236,6 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") { 0, pressure_constraints) == true); REQUIRE(constraints->assign_nodal_pressure_constraints( 1, pressure_constraints) == true); - // When constraints fail - REQUIRE(constraints->assign_nodal_pressure_constraints( - 4, pressure_constraints) == false); - - pressure_constraints.emplace_back(std::make_tuple(100, 0.0)); - REQUIRE(constraints->assign_nodal_pressure_constraints( - 0, pressure_constraints) == false); } // Test assign nodes concentrated_forces diff --git a/tests/mesh_test_3d.cc b/tests/mesh_test_3d.cc index 689f16ced..3b7a80b74 100644 --- a/tests/mesh_test_3d.cc +++ b/tests/mesh_test_3d.cc @@ -1347,13 +1347,6 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") { 0, pressure_constraints) == true); REQUIRE(constraints->assign_nodal_pressure_constraints( 1, pressure_constraints) == true); - // When constraints fail - REQUIRE(constraints->assign_nodal_pressure_constraints( - 4, pressure_constraints) == false); - - pressure_constraints.emplace_back(std::make_tuple(100, 0.0)); - REQUIRE(constraints->assign_nodal_pressure_constraints( - 0, pressure_constraints) == false); } // Test assign nodes concentrated_forces